Physics-based numerical modeling of earthquake source processes strives to predict quantities of interest for seismic hazard, such as the probability of an earthquake rupture jumping between fault segments. How to assess the predictive power of numerical models remains a topic of ongoing debate. Here, we investigate how sensitive are the outcomes of numerical simulations of sequences of earthquakes and aseismic slip to choices in numerical discretization and treatment of inertial effects, using a simplified 2-D crustal fault model with two co-planar segments separated by a creeping barrier. Our simulations demonstrate that simplifying inertial effects and using oversized cells significantly affects the resulting earthquake sequences, including the rate of two-segment ruptures. We find that a number of fault models with different properties and modeling assumptions can produce comparable frequency-magnitude statistics and static stress drops but have rates of two-segment ruptures ranging from 0 (single-segment ruptures only) to 1 (two-segment ruptures only). For sufficiently long faults, we find that long-term sequences of events can substantially differ even among simulations that are well-resolved by standard considerations. In such simulations, some outcomes, such as static stress drops, are stable among adequately-resolved simulations, whereas others, such as the rate of two-segment ruptures, can be highly sensitive to numerical procedures and physical assumptions, and hence cannot be reliably inferred. Our results emphasize the need to examine the potential dependence of simulation outcomes on the modeling procedures and resolution, particularly when assessing their predictive value for seismic hazard assessment.

Determining conditions for earthquake slip on faults is a key goal of fault mechanics highly relevant to seismic hazard. Previous studies have demonstrated that enhanced dynamic weakening (EDW) can lead to dynamic rupture of faults with much lower shear stress than required for rupture nucleation. We study the stress conditions before earthquake ruptures of different sizes that spontaneously evolve in numerical simulations of earthquake sequences on rate-and-state faults with EDW due to thermal pressurization of pore fluids. We find that average shear stress right before dynamic rupture (aka shear prestress) systematically varies with the rupture size. The smallest ruptures have prestress comparable to the local shear stress required for nucleation. Larger ruptures weaken the fault more, propagate over increasingly under-stressed areas due to dynamic stress concentration, and result in progressively lower average prestress over the entire rupture. The effect is more significant in fault models with more efficient EDW. We find that, as a result, fault models with more efficient weakening produce fewer small events and result in systematically lower b-values of the frequency-magnitude event distributions. The findings 1) illustrate that large earthquakes can occur on faults that appear not to be critically stressed compared to stresses required for slip nucleation; 2) highlight the importance of finite-fault modeling in relating the local friction behavior determined in the lab to the field scale; and 3) suggest that paucity of small events or seismic quiescence may be the observational indication of mature faults that operate under low shear stress due to EDW.