AbstractStress transfer between earthquakes is recognized as a fundamental mechanism governing aftershock sequences. A common approach to relate stress changes to seismicity rate changes is the rate‐and‐state constitutive law developed by Dieterich: these elements are the foundation of Coulomb‐rate‐and‐state (CRS) models. Despite the successes of Coulomb hypothesis and of the rate‐and‐state formulation, such models perform worse than statistical models in an operational forecasting context: one reason is that Coulomb stress is subject to large uncertainties and intrinsic spatial heterogeneity. In this study, we characterize the uncertainties in Coulomb stress inherited from different physical quantities and assess their effect on CRS models. We use a Monte Carlo method and focus on the following aspects: the existence of multiple receiver faults; the stress heterogeneity within grid cells, due to their finite size; and errors inherited from the coseismic slip model. We study two well‐recorded sequences from different tectonic settings: the Mw = 6.0 Parkfield and the Mw= 9.0 Tohoku earthquakes. We find that the existence of multiple receiver faults is the most important source of intrinsic stress heterogeneity, and CRS models perform significantly better when this variability is taken into account. The choice of slip model also generates large uncertainties. We construct an ensemble model based on published slip models and find that it outperforms individual models. Our findings highlight the importance of identifying sources of errors and quantifying confidence boundaries in the forecasts; moreover, we demonstrate that consideration of stress heterogeneity and epistemic uncertainty has the potential to improve the performance of operational forecasting models.