Resolving Simulated Sequences of Earthquakes and Fault Interactions: Implications for Physics‐Based Seismic Hazard Assessment

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 the outcomes of numerical simulations of sequences of earthquakes and aseismic slip are 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 fault models with different properties and modeling assumptions can produce similar frequency‐magnitude statistics and static stress drops but have different rates of two‐segment ruptures. 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 similar among adequately resolved simulations, whereas others, such as the rate of two‐segment ruptures, can be highly sensitive to numerical procedures and modeling assumptions. While it is possible that the response of models with additional ingredients ‐Realistic fault geometry, fluid effects, etc. ‐Would be less sensitive to numerical procedures, 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.


of 33
The 2016 Mw 7.8 Kaikoura earthquake ruptured at least 21 segments of the Marlborough fault system (Ulrich et al., 2019). Increasingly, seismological observations show that it is not uncommon to see ruptures navigating and triggering subsequent ruptures within fault networks, including the recent 2019 Mw 6.4 and 7.1 Ridgecrest earthquakes (Ross et al., 2019), and the 2012 Mw 8.6 and 8.2 Indian Ocean earthquakes (Wei et al., 2013), the largest and second-largest recorded strike-slip earthquakes to date. Yet, in any given seismogenic region, the record of past large events is not long enough to forecast the behavior of ruptures with respect to the existing fault segments, specifically how likely would the rupture be to jump between nearby segments, prompting the discussion on whether and how physics-based models may inform this and other questions important for seismic hazard assessment (Field, 2019).
Determining what conditions allow a dynamic rupture to propagate or arrest are key to understanding the maximum potential magnitude of an earthquake. Previous modeling of single fully dynamic ruptures have shown great success in investigating earthquake propagation in nonplanar and multi-segment fault models, including step-overs and branched geometries (Ando & Kaneko, 2018;Douilly et al., 2015;Duan & Oglesby, 2006;Dunham et al., 2011b;Galvez et al., 2014;Harris & Day, 1993, 1999Harris et al., 1991;Hu et al., 2016;Kame et al., 2003;Lozos et al., 2015;Ulrich et al., 2019;Withers et al., 2018;Wollherr et al., 2019). In particular, such modeling has shown that the ability of a rupture to propagate across segments depends on the stresses before the rupture and shear resistance assumptions, as well as the geometry of the fault system. However, single-rupture simulations need to select initial conditions and need additional assumptions to incorporate the effect of previous seismic and aseismic slip.
Fault processes involve both sequences of dynamic events and complex patterns of quasi-static slip. Simulating this behavior in its entirety is a fascinating scientific problem. However, even for the more pragmatic goal of physics-based predictive modeling of destructive large dynamic events, it is still important to consider sequences of earthquakes and aseismic slip (SEAS), since prior slip events, including aseismic slip, may determine where earthquakes would nucleate as well as modify stress and other initial conditions before dynamic rupture. Furthermore, such simulations provide a framework for determining physical properties consistent with a range of observations including geodetically recorded surface motions, microseismicity, past (including paleoseismic) events, and thermal constraints, and hence may inform us about the current state of a fault segment or system and potential future rupture scenarios (e.g., Allison & Dunham, 2018;Barbot et al., 2012;Ben-Zion & Rice, 1997;Cattania, 2019;Chen & Lapusta, 2009;Erickson & Day, 2016;Erickson & Dunham, 2014;Jiang & Lapusta, 2016;Kaneko et al., 2010;Lambert & Barbot, 2016;Lambert, Lapusta, Perry, 2021;Lapusta & Rice, 2003;Lapusta et al., 2000;Lin & Lapusta, 2018;Liu & Rice, 2005;Perry et al., 2020;Segall et al., 2010). However, simulating long-term slip histories is quite challenging because of the variety of temporal and spatial scales involved.
Recently, several earthquake simulators have been developed with the goal of simulating millions of earthquake ruptures over regional fault networks for tens of thousands of years Shaw et al., 2018;Tullis et al., 2012). The term "simulators" typically refers to approaches that employ significant simplifications, compared to most SEAS simulations, in solution procedures and physical processes, in order to simulate earthquake sequences on complex, regional scale 3-D fault networks for long periods of time. For example, earthquake simulators typically account only for the quasi-static stress transfer due to earthquake events, ignoring wave-mediated stress changes, aseismic slip/deformation, and fluid effects; employ approximate rule-based update schemes for earthquake progression instead of solutions of the governing continuum mechanics equations; and use oversized numerical cells. Such simplifications are currently necessary to permit simulations of hundreds of thousands of events over hundreds of fault segments that comprise the regional networks (Shaw et al., 2018). Earthquake simulators have matched a number of regional-scale statistical relations, including the Gutenberg-Richter frequency-magnitude scaling (Shaw et al., 2018), and highlighted the importance of large-scale fault and rupture interactions.
There is growing interest in using earthquake simulators to directly determine quantities of interest for seismic hazard, such as the probability of an earthquake rupture jumping from one fault segment to another (Field, 2019;Shaw et al., 2018). However, assessing the predictive power of numerical models remains a topic of active research. Determining how sensitive simulated outcomes may be to modeling choices and how reliably they can be determined from a given numerical model are topics of great importance for physics-based hazard assessment. Here, we examine the sensitivity of the long-term interaction of fault segments 3 of 33 to choices in numerical discretization and representations of inertial effects in simulated sequences of earthquakes and aseismic slip, using a relatively simple 2-D model of two co-planar strike-slip fault segments separated by a velocity-strengthening (VS) barrier. We explore how considerations for adequate numerical resolution and convergence depend on the physical assumptions and complexity of earthquake sequences as well as on the modeling outcome of interest. We especially focus on the rate of earthquake ruptures jumping across the VS barrier and examine whether reproducing comparable earthquake frequency-magnitude statistics and static stress drops provides sufficient predictive power for the jump rate, a quantity of interest to seismic hazard studies (Field, 2019).

Model Setup and Numerical Resolution
Our simulations are conducted following the methodological developments of Lapusta et al. (2000), Noda and Lapusta (2010), and Lambert, Lapusta, Perry (2021). We consider a one-dimensional (1-D) fault embedded into a 2-D uniform, isotropic, elastic medium ( Figure 1). The 2-D model approximates a faulted crustal plate coupled to a moving substrate using the idea of a constrained continuum (Johnson, 1992;Lehner et al., 1981). Fault slip may vary spatially along-strike but it is depth-averaged through a prescribed seismogenic thickness = 15 km, beneath which the elastic domain is coupled to a substrate moving at the prescribed loading rate ( pl = 10 −9 m/s). The elastodynamic equation for the depth-averaged displacement along-strike ( ) is given by (Kaneko & Lapusta, 2008;Lehner et al., 1981): where eff = ( ∕4) S and = 1∕(1 − ) , with being the Poisson's ratio. The effective wave speed alongstrike for the crustal plane model is = , where is the shear wave speed. The along-strike slip is then given by ( ) =̄( = 0 + ) −̄( − ) .
Our simulations resolve sequences of earthquakes and aseismic slip (SEAS) in their entirety, including the gradual development of frictional instability and spontaneous nucleation, dynamic rupture propagation, post-seismic slip that follows the event, and the interseismic period between events ( Figure 2). In all models, frictional resistance along the fault interface is governed by the standard laboratory-derived rate-andstate friction law with the state evolution described by the aging law (Dieterich, 1979;Ruina, 1983): where = ( − ) is the effective normal stress, is the normal stress, is the pore pressure, is the shear stress, is the friction coefficient, is the slip velocity, is the state variable, RS is the characteristic slip for the evolution of the state variable, * is the reference steady-state friction coefficient corresponding to a reference slip rate * , and and are the direct and evolution effect constitutive parameters, respectively. At steady-state (constant slip velocity), the shear stress and state variable evolve to their steady-state values and given by: In our simulations, we use a 2D approximation of the problem with a 1D along-strike depth-averaged fault, in which the fault is assumed to be creeping at the loading plate rate pl = 10 −9 m/s below the depth of = 15 km.

of 33
The combination of frictional properties such that ( − ) > 0 results in steady-state velocity-strengthening (VS) behavior, where the shear resistance increases with an increase in slip velocity and where stable slip is expected. If ( − ) < 0 then the fault exhibits velocity-weakening (VW) behavior, in which case an increase in slip velocity leads to a decrease in shear resistance, making these regions of the fault potentially seismogenic if their size exceeds a critical nucleation size. Interaction of two co-planar fault segments in well-resolved simulations of model M1 demonstrating convergence of simulated earthquake sequences. (a and b) History of cumulative slip over 4,000 years in well-resolved fully dynamic simulations of fault model M1 with initial conditions S1 using (a) 12.5-m and (b) 25-m cell size. Contours for seismic slip are plotted every 0.5 s, with ruptures that jump across the VS barrier colored blue. The simulated fault behavior is virtually indistinguishable between the two resolutions. (c) Frequency-magnitude histograms of events, on top of each other for the two resolutions. The well-resolved simulations produce the same relatively simple and quasi-periodic behavior. (d and e) The evolution of local shear stress and slip velocity at a point ( = −20.5 km, shown by star in (a and b), practically indistinguishable even after over 3,800 years of simulated time. (f-h) Spatial distribution of shear stress at the rupture front for three locations ( = −20 km, 5 and 20 km) throughout the first rupture in (a and b). While the quasi-static estimate of the cohesive zone Λ0 is about 1.1 km, the actual size of the cohesive zone varies with the local rupture speed throughout the rupture. In these wellresolved simulations, the cohesive zone is always resolved by at least 10 cells.

of 33
Two theoretical estimates of the nucleation size in mode II are (Rice & Ruina, 1983;Rubin & Ampuero, 2005): where is the shear modulus. The estimate * was derived from the linear stability analysis of steady frictional sliding by Rice and Ruina (1983). It also represents the critical cell size for steady-state quasi-static sliding such that larger cells can become unstable on their own. Thus * represents a key length scale to resolve for slow interseismic processes and earthquake nucleation (Lapusta et al., 2000;Rice & Ruina, 1983). The estimate * was determined in the parameter regime ∕ 0.5 using the energy balance of a quasi-statically expanding crack (Rubin & Ampuero, 2005), and provides an estimate of the minimum size for a slipping region that releases enough stored energy to result in the radiation of waves.
We aim to explore the impact of numerical resolution on the long-term simulated slip behavior of sequences of earthquakes and aseismic slip. The nucleation size, * , estimated by either * or * from Equation 6, is one length-scale that clearly needs to be well resolved. Early resolution studies for sequences of events showed that resolution of the nucleation scale * by 20-40 cells is required for stable numerical results (Lapusta et al., 2000). Later, the need to resolve the nucleation size by at least 20 cells was shown to be due to the more stringent criterion of resolving the region where shear resistance breaks down at the rupture front, often referred to as the cohesive zone. The cohesive zone can be an order of magnitude smaller than the nucleation size, depending on the constitutive description (Day et al., 2005;Lapusta & Liu, 2009). The size of the cohesive zone depends on the weakening rate of shear stress with slip associated with the constitutive law. The quasi-static estimate Λ0 of the cohesive zone size at near-zero rupture speed and constant is given by: where 1 is a constant, ′ = for mode III, and ′ = ∕(1 − ) for mode II (Rice, 1980). For standard rate-and-state friction with the aging form of the state variable evolution, the weakening rate is given by = RS∕( ) (Lapusta & Liu, 2009) and: If one assumes that the traction distribution within the cohesive zone is linear, then the constant 1 can be approximated as 1 = 9 ∕32 (Rice, 1980). For fully dynamic rupture simulations, continuously resolving the breakdown process at the rupture front becomes even more challenging as the cohesive zone size Λ exhibits a contraction with increasing rupture speed (e.g., Rice, 1980): . Note that −1 (0 + ) = 1 , giving the quasi-static cohesive zone estimate Λ0 when = 0 + . As the rupture speed approaches the limiting wave speed, → (Rayleigh wave speed) for mode II and → (shear wave speed) for mode III, one has −1 ( ) → 0 and the width of the breakdown region approaches zero. Hence, it becomes increasingly more challenging to resolve the rupture front during fully dynamic simulations if the rupture accelerates toward the limiting speeds. Such acceleration typically occurs during long enough propagation of dynamic rupture over favorable prestress, unless impeded by additional factors such as unfavorable prestress or situations with increasing effective breakdown energy, for example, due to off-fault inelasticity or navigating fault roughness (Andrews, 2005;Dunham et al., 2011b;Okubo et al., 2019;Perry et al., 2020;Poliakov et al., 2002;Rice, 2006). Simulations of faults with rate-and-state friction and the aging form of the state variable evolution embedded in elastic bulk result in ruptures with near-constant breakdown energy (Perry et al., 2020) and this holds for most cases considered in this study. In Section 7, we show that adding an approximation of off-fault inelasticity to our simulations that reduces the rupture speeds does not alter our conclusions.
In our model, the fault contains a frictional domain consisting of two VW regions of length = 32 km that are separated by a 2-km-long VS region that impedes rupture propagation. We select large enough 6 of 33 values of the velocity strengthening in the central VS region so that the region acts like a barrier, requiring ruptures to jump/renucleate on the other side of the barrier to propagate over the second segment. This region is a proxy for what would be a gap in the fault connectivity, at least at the surface, requiring the ruptures to jump across. The remainder of the frictional region surrounding the VW segments has more mild VS properties (Figure 1). At the edges of the model, outside of the frictional domain, fault slip is prescribed at the loading plate rate. Values for the model parameters used in our simulations are provided in Tables 1 and 2. We first examine models with lower instability ratio ∕ℎ * that result in quasi-periodic sequences of events, and then consider models with higher instability ratios that result in more complex earthquake sequences and qualitatively different convergence behavior.

Resolving Quasi-Periodic Fully Dynamic Sequences of Earthquakes and Aseismic Slip (SEAS)
Let us consider the simulated slip behavior of fault model M1 with instability ratio ∕ℎ * = 21 (Table 2). Its quasi-static cohesive zone ( Λ0 = 1.1 km) should be well-resolved by cell sizes of 12.5 and 25 m, with 88 and 44 cells over Λ0 , respectively; the nucleation size is even larger and hence also well-resolved. Consistently with these considerations, these two well-discretized simulations produce the same relatively simple quasi-periodic sequences of earthquake events that periodically jump across the VS barrier (Figures 2a and 2b). We clearly see that the results are the same for the two simulations with different resolutions, including the local evolution of slip rate and shear stress during ruptures late in the earthquake sequence (Figures 2d and 2e). Note that the cohesive zone evolves throughout the rupture process, shrinking with the increasing rupture speed by 3-4 times in these simulations (Figures 2f-2h) and the spatial discretization is fine enough to adequately characterize the rupture front throughout the entire dynamic process. The jump rate for both simulations is 0.54; we define this rate of ruptures jumping across the VS barrier within a given time period as the total number of ruptures that propagate toward the barrier and result in seismic slip on both fault segments divided by the total number of ruptures that propagate toward the barrier and span at least the segment on which they nucleate.
The variability between different ruptures in fault model M1 is generally mild, as shown by their frequency-magnitude histograms ( Figure 2c). To create the frequency-magnitude histograms, we compute the moment for each simulated event in our 2-D models as =̄ where the rupture area is defined with  The quasi-periodic nature of events observed over the first 4,000 years in well-resolved simulations of fault model M1 persists in longer-duration simulations over 20,000 years, resulting in similar long-term jump rates of 0.48-0.54 depending on the time interval considered (Figure 3). We also examine simulations of fault model M1 with different initial shear stress conditions and find that the long-term sequences of events converge to the same quasi-periodic behavior upon adequate discretization, despite the initial few events being different (Figure 3a vs. 3b; details of initial shear stress distributions S1 and S2 are provided in the Supporting Information S1). Simulations of fault model M1 thus exhibit long-term numerical convergence upon adequate discretization, producing virtually indistinguishable long-term slip behavior and a consistent rate of two-segment ruptures among simulations with differing initial conditions, after a sufficiently large initial sequence of events. Convergence of well-resolved simulated earthquake sequences in model M1 for longer-term simulations and different initial conditions. (a and b) Cumulative slip over 0-4,000 years and 16,000-20,000 years in two well-resolved fully dynamic simulations of fault model M1 with two different initial conditions, S1 and S2. Contours of seismic slip are plotted every 0.5 s with ruptures that jump across the VS barrier colored blue. The quasi-periodic behavior seen in the first 4,000 years in well-resolved simulations, including the rate of ruptures jumping across the VS barrier, remains generally consistent throughout longer-term simulations over 20,000 years (right). Simulations using different initial shear stress conditions produce different initial sequences of events; however, the simulated sequences converge to the same slip behavior and have the same long-term rates of two-segment ruptures (0.50 over 2,000-20,000 years). (c and d) Normalized frequency-magnitude histograms for events from (a) to (b), respectively, over 4,000 and 20,0000 years, illustrating that the population statistics in this relatively simple system is the same, apart from the initial start-up period.

of 33
Let us now consider simulations that use larger computational cells. The cell sizes of 250 and 125 m resolve the quasi-static cohesive zone Λ0 with 4.5-9 cells ( Figure 4). While this resolution seems adequate, one can anticipate that the dynamic shrinking of the cohesive zone size by 3-4 times would result in a more marginal resolution of 1-3 cells. Indeed, we see that the simulated long-term sequences of events and jump rates differ substantially from those of the well-resolved simulations (Figures 2a and 2b vs. 4a and 4b).
Considering even larger cell sizes of 500 and 1,000 m brings further differences in the event sequences and jump rates ( Figure 5), with the earthquake sequences that look plausible and not obviously numerically compromised even for the largest cell sizes ( Figure S1 in Supporting Information S1). Note that the jump rate in simulations with marginal and oversized cells is neither systematically larger nor smaller than the range 0.48-0.54 from the well-resolved cases, but varies from 0.25 to 0.95 depending on the choice of numerical discretization.
Increasingly poor resolution of the dynamic cohesive zone at the rupture front and, for the largest cell sizes, of the nucleation zone results in an increasing abundance of small events ( Figure 5), as had been shown in previous studies (Lapusta & Liu, 2009;Rice, 1993;Rice & Ben-Zion, 1996). Inadequate resolution of the dynamic rupture front prevents simulating the actual stress concentration and promotes event arrest. Inadequate resolution of the nucleation size enables individual cells or small number of cells to fail 9 of 33 independently due to the inadequate resolution of the stress interactions (Lapusta & Liu, 2009;Rice, 1993;Rice & Ben-Zion, 1996). Using sufficiently oversized cells can result in power-law statistics in terms of the frequency-magnitude distribution of simulated earthquake ruptures (Figures 5e-5j, Rice, 1993;Rice & Ben-Zion, 1996).
Note that the suggested minimum average resolution of three cells of the (variable) cohesive zone from the dynamic rupture study by Day et al. (2005) is not adequate for convergent results in these earthquake sequence simulations. That criterion would be achieved in this model for a cell size between the 250 and 125 m. Yet, the simulated long-term behavior for those cell sizes is clearly different from the better-resolved 10 of 33 and convergent results with the cell sizes of 25 and 12.5 m. At the same time, the criterion by Day et al. (2005) works well for a single dynamic rupture as intended, since the first dynamic events in simulations with cell sizes 12.5, 25, 125, and 250 m are quite similar to each other ( Figure S2 in Supporting Information S1). The events are not identical, however; for example, the average slip with the resolution of 12.5 and 125 m differs by 0.7%. Clearly, these differences, acceptable for a single event, accumulate in these highly nonlinear solutions, resulting in different event statistics and jump rate ( Figure 5).
We find that our fully dynamic 2-D simulations of fault model M1, which include uniform VW properties with relatively mild weakening due to standard rate-and-state friction, converge when the quasi-static cohesive zone estimate Λ0 is discretized by at least 22 cells, which translates to the average resolution of the dynamically variable cohesive zone size of 10-15 cells. Fault models with additional or different ingredients, such as fault heterogeneity/roughness, more efficient weakening, 3D elastodynamics with 3D faults, or different instability ratio, would require further considerations for resolution requirements that result in convergent simulations. For example, as we discuss in Section 6, the convergence and resolution properties of models with higher instability ratios, which result in more complex earthquake sequences, are qualitatively different.
In the more complicated earthquake sequences observed in under-resolved simulations of fault model M1, some statistics, such as the rate of two-segment ruptures, depends on the specific period that one considers throughout the simulation. To explore the variability in the event statistics and jump rate across the VS barrier in models with different numerical resolution, we examine the jump rate over different 2000-year periods throughout longer term simulations of 20,000 years, using a sliding window of 1,000 years starting at the beginning of the simulation (19 periods total; Figure 5). The choice of a 2000-year period allows us to have a sufficient number ( ∼ 20) of large earthquakes within a period to estimate jump rates. We also consider the outcomes for two different initial conditions S1 and S2, as before. For the well-resolved simulations exhibiting long-term convergence, the frequency-magnitude and 2000-year jump rate statistics for simulations with different initiation conditions are comparable, with the jump rate for all 2000-year periods being consistent with the overall 20,000 years jump rate (Figures 5a and 5b). As the numerical resolution worsens, the sequences of events become more complex with greater variability in rupture sizes and increased production of smaller events (Figures 5c-5j). The jump rate during any 2000-year period also becomes more variable in marginally resolved simulations and can considerably differ from both the 20,000-year jump rate of the same simulation as well as from the true jump rate in the well-resolved simulations. Note that despite being clearly affected by numerical resolution, the frequency-magnitude and jump-rate distributions of inadequately resolved simulations can appear generally consistent among simulations with similar cell sizes and different initial conditions ( Figure 5 left vs. right columns). In other words, even if simulations using marginal or oversized cells produce comparable statistical properties for different initial conditions, these characteristics do not necessarily represent robust features of the physical system but rather may still be numerical artifacts.

Interaction of Fault Segments in Simulations With Quasi-Dynamic Approximation for Inertial Effects
Many numerical studies of long-term fault behavior utilize quasi-dynamic solutions to the equations of motion, in which the wave-mediated stress transfers during the coseismic phase are replaced with a radiation damping approximation (Rice, 1993). The quasi-dynamic approximation substantially reduces the computational expense of the simulation, as the consideration of stress redistribution by waves requires substantial additional storage and computational expense. Considerable insight into fault mechanics has been derived from studies using quasi-dynamic formulations, particularly when such approximations are used to incorporate new physical effects that may otherwise result in prohibitive computational expense, as well as in scenarios where it may be argued that inertial effects are relatively mild, such as during earthquake nucleation or during aseismic slip transients (Allison & Dunham, 2018;Erickson et al., 2017;Lambert & Barbot, 2016;Liu, 2014;Liu & Rice, 2005, 2007Rice, 1993;Rubin & Ampuero, 2005;Segall & Rice, 1995;Segall et al., 2010). However, as with all approximations, it is important to be aware of how such simplifications modify the outcome of study (Thomas et al., 2014).

of 33
Let us review the quasi-dynamic approximation for inertial effects during sliding and study their implications for the long-term interaction of two fault segments. In the 2D boundary integral formulation, the elastodynamic shear stress along a 1D fault plane, can be expressed as (Cochard & Madariaga, 1994;Perrin et al., 1995): where 0 ( ) are the "loading" tractions (i.e., the stress induced on the fault plane if it were constrained against any slip), static( ) and dynamic( ) represent the static and dynamic contributions to the stress transfer along the fault, respectively, and the last term represents radiation damping ( = ∕(2 ) for mode III).
The static solution for the equations of motion would only contain static , which depends only on the current values of slip along the fault. However, the static solution does not exist during dynamic rupture when inertial effects becomes important.
dynamic and both arise due to the inertial effects. dynamic represents the wave-mediated stress interactions along the interface and this term is challenging to compute as it requires calculating convolutions on time and storing the history of deformation. Radiation damping is much easier to incorporate as it depends on the current slip rate, and represents part of the radiated energy (Rice, 1993). The quasi-dynamic approximation, in which dynamic is ignored and only is included, allows the solution to exist during inertially controlled dynamic rupture. However, the solution is altered from the true elastodynamic representation.
Let us consider the long-term behavior of fault model M1, as examined in Section 3, but now using the quasi-dynamic approximation. For well-resolved quasi-dynamic simulations of fault model M1, we find that the long-term slip behavior of the two-fault segment system is even simpler than for the fully dynamic case, with ruptures being exclusively isolated to individual segments and the jump rate being zero (Figure 6a). For simulations with the increasing cell size, and thus decreasing spatial resolution, we see increased variability in the size of the individual ruptures, to the point where some marginally resolved simulations produce ruptures that jump across the VS barrier, whereas well-resolved simulations of the same fault model never do (Figures 6b and 6c). The increasing cell size also leads to increased production of smaller events and more complicated fault behavior, similarly to the fully dynamic simulations (Figures 6d-6f).
In addition to substantially reducing the computational expense associated with calculating the wave-mediated stress transfers, quasi-dynamic simulations place milder constraints on the spatial resolution since the cohesive zone always remains near the quasi-static estimate, Λ ≈ Λ0 = Λ( = 0 + ) (Figures 6g and 6h). This is because the stress transfer calculated for the ruptures is always quasi-static, and the much stronger stress transfer due to waves is ignored ( Figure 7e). As a result, the quasi-dynamic simulations produce significantly smaller slip velocities and rupture speeds than the fully dynamic ones 7a-7c, consistent with previous studies (Lapusta & Liu, 2009;Thomas et al., 2014).
One can attempt to enhance the slip rates and rupture speeds in the quasi-dynamic simulations by reducing the radiation damping term ; this can be interpreted as increasing the effective shear wave speed in the radiation damping term enh. = , thus allowing for higher slip rates (Lapusta & Liu, 2009). We compare the enhanced quasi-dynamic simulations ( = 3 ) with the standard quasi-dynamic ( = 1 ) and fully dynamic simulations of fault model M1 ( Figure S3 in Supporting Information S1). Decreasing the radiation damping increases the effective rupture speed and slip rate (Figure 7a-7c) in comparison to the standard quasi-dynamic simulation, however, for the parameters considered, it does not substantially alter the longterm interactions of the two fault segments, nor match the rate of ruptures jumping across the VS barrier in the fully dynamic case ( Figure S3 in Supporting Information S1).
In comparing the three simulations with different treatment of the inertial effects, it is clear that the fully dynamic ruptures result in higher slip rates and narrowing of the cohesive zone (Figure 7). For simulations with standard rate-and-state friction, the peak shear stresses vary mildly from fully dynamic versus quasi-dynamic representations, as they are limited by the shear resistance of the fault, which has a relatively mild logarithmic dependence on slip rate. However, the stress transfer along the fault substantially differs for fully dynamic versus quasi-dynamic representations (Figure 7e). The difference between the stress transfer term and the shear stress is accommodated by the radiation damping , which results in higher slip rates to balance the larger dynamic stress transfers (Figures 7c-7e). Hence while the resolved peak 12 of 33 shear stresses along the fault may be comparable due to the specific choice of the constitutive relationship, the rupture dynamics and kinematics, as seen through the stress transfer, slip rate, and rupture speed along the fault, differ considerably with and without the inclusion of full inertial effects.
These larger dynamic stress transfers facilitate the triggering and continued propagation of slip on the neighboring fault segment, rather than leaving the rupture to always be arrested by the creeping barrier, as in the well-resolved quasi-dynamic simulations ( Figure S3 in Supporting Information S1). Decreasing 13 of 33 the radiation damping term allows for somewhat higher slip rates and arbitrarily higher rupture speeds, but it does not mimic the full effects in the dynamic stress transfer, particularly at the rupture front. As the result, the fully dynamic simulations have higher jump rates. The differences between fully dynamic and quasi-dynamic approximations can be even more substantial for models with enhanced weakening at seismic slip rates from the flash heating of contact asperities or the thermal pressurization of pore fluids (Thomas et al., 2014). The fully dynamic rupture accelerates to a rupture speed close to the limiting wave speed of ≈ 4.4 km/s throughout the rupture, whereas the quasi-dynamic ruptures maintain lower effective rupture speeds. Decreasing the radiation damping term for quasidynamic ruptures increases the slip rate and rupture speed, but does not truly mimic the acceleration of the fully dynamic rupture. (c and d) A closer look at the spatial distribution of (c) slip velocity and (d) shear stress at a given time highlights how full consideration of inertial effects leads to much higher slip velocities and a more localized stress concentration at the rupture front, which facilitates rupture propagation. Enhancing the quasi-dynamic ruptures with lower radiation damping increases the slip rate but maintains the same quasi-static spatial pattern of stress at the rupture front. (e) The corresponding values of the stress transfer functional near the rupture front. The radiation damping approximation of the inertial effects results in dramatically reduced stress transfer along the fault. The larger total stress transfer in the fully dynamic simulations is balanced by higher slip rates, as shown in (c).

Rupture Jump Rates in Simulations With Similar Power-Law Frequency-Magnitude Statistics and Earthquake Stress Drops
Two common observations about natural earthquakes and regional seismicity are the average static stress drops between 1 and 10 MPa independently of the event magnitude (e.g., Allmann & Shearer, 2009;Ye et al., 2016) as well as the frequency-magnitude statistics of earthquakes within a region, which commonly follow the Gutenberg-Richter power law relation (Field et al., 2013). Some earthquake simulators are capable of matching these observations (Shaw et al., 2018). An important question is whether matching these constraints allows simulators to be predictive for other quantities of interest to seismic hazard assessment, such as the probability of multiple fault-segment ruptures, despite using approximations for inertial effects and oversized computational cells.
Let us consider this question using simulations of earthquake sequences in five fault models with the same (simple) fault geometry but different friction properties and different assumptions about inertial effects, and one additional model in which the effective seismogenic depth is slightly reduced from 15 to 14 km ( Figure 8, Table 2). All six models have comparable nucleation and quasi-static cohesive zone sizes (Table 2) and use oversized cells of Δ = 1,000 m. (An example of well-resolved simulations with similar conclusions is given in Section 7.) The six simulations produce comparable frequency-magnitude distributions when interpreted as power-law, characterized by b-values of 0.3-0.4 for 4,000 years of the simulated time. All six simulations also produce ruptures with comparable average static stress drops ( Figure S4 in Supporting Information S1), with values typically between 1 and 10 MPa, as commonly inferred for natural earthquakes (Allmann & Shearer, 2009;Ye et al., 2016).
However, the probability of a rupture jumping across the VS barrier varies among the six simulations, ranging from 0 to near 1. This substantial variability in jump rate for simulations with comparable power-law frequency-magnitude statistics persists in longer-duration simulations over 20,000 years, where both the 20,000-year jump rate and distributions of jump rates within individual 2000-year periods can substantially differ (Figure 9). In particular, fault model M1 results in a jump rate of 0 for the quasi-dynamic simulation and near 1 for the fully dynamic simulation (Figures 8 and 9a vs. 9d). This case illustrates how using approximations for inertial effects may considerably bias estimates of the actual rate of multisegment ruptures, even if the power-law frequency-magnitude statistics and static stress drops are comparable. In addition, the suite of simulations suggest that the probability of ruptures jumping across the VS barrier is sensitive to variations in the frictional parameters, effective normal stress, as well as minor changes in the seismogenic depth.
Note that one-segment and two-segment events would not appear that different on the logarithmic moment-magnitude scale. In our 2D crustal fault model, slip is limited by the seismogenic thickness, resulting in (geometric) pulse-like rupture propagation for large events. If one assumes that large ruptures have comparable average slip given the specific geometry of our simple model, then ruptures that span a single segment or two segments differ by about a factor of two based on the difference in rupture area. A factor of 2 in seismic moment changes the moment magnitude only modestly. For example, a single-segment rupture with a seismic moment of 0.5 × 10 20 Nm corresponds to 7.06 whereas a two-segment rupture with a seismic moment of 10 20 Nm corresponds to 7.26; such moment magnitudes would be binned together for the choice of 0.5 bin width in Figure 8. Furthermore, the actual amount of slip in individual ruptures can vary; in our simulations, some single-segment ruptures exhibit greater average slip than two-segment ruptures (Figure 8), resulting in more variability in seismic moment among large ruptures and greater overlap in the moment of single and two-segment ruptures. The distinction in the seismic moment of large single-segment and multisegment ruptures may be even less obvious in more realistic models of fault networks with fault segments of varying surface area.
The results from our simple 2-D modeling suggest that reproducing static stress drops and the b-value of the frequency-magnitude event distribution interpreted as a power-law does not clearly provide sufficient predictive power for rupture jump rates, although the results could be different for models with additional physics, as further discussed in Section 8. This finding is perhaps not surprising given that many combinations of rate-and-state properties and effective normal stress may produce ruptures with comparable static stress changes ( Figure S4 in Supporting Information S1), but different overall levels of shear resistance.

of 33
Similarly, a number of studies have demonstrated that power-law frequency-magnitude statistics can be reproduced in many models, including discrete fault models (Bak & Tang, 1989;Burridge & Knopoff, 1967;Olami et al., 1992), continuum fault models that are inadequately resolved and therefore numerically discrete (Ben-Zion & , continuum models with spatially heterogeneous frictional properties (Hillers et al., 2006(Hillers et al., , 2007, and continuum models of single or multiple fault segments with larger instability ratio Figure 8. Models with comparable power-law frequency-magnitude statistics and static stress drops but very different rate of two-segment ruptures. (a-f) Cumulative frequency-magnitude histograms (top) and history of cumulative slip (bottom) over 4,000 years in (a-c) fully dynamic and (d-f) quasi-dynamic SEAS simulations. Contours for seismic slip are plotted every 0.5 s, with ruptures that jump across the VS barrier colored blue. The simulations assume different physical conditions described in the text. All six simulations produce comparable average static stress drops ( Figure S4 in Supporting Information S1) and comparable population statistics with a b-value around 0.33. However, the rate of two-segment ruptures varies from 0 to 1.  Figure 8. The six simulations have comparable frequency-magnitude statistics but the 20,000-year rate of two-segment ruptures varies from 0 to 0.91. The distribution of 2000-year jump rates is also highly variable among the six simulations. (Barbot, 2021;Cattania, 2019;Wu & Chen, 2014). Therefore, Gutenberg-Richter statistics may be compatible with a range of fault properties, and may not pose a considerable physical constraint on its own. In particular, the presence of smaller-event complexity can be due to factors unrelated to the potential for ruptures jumping between segments, such as small-scale heterogeneity or heterogeneous loading conditions (Hillers et al., 2006(Hillers et al., , 2007Rice, 1993).
Note that, while the simulations of Figures 8 and 9 have similar power-law frequency-magnitude statistics with comparable b-values, the underlying frequency-magnitude distributions are different, especially for the larger events ( Figure S5 in Supporting Information S1). Some of these distributions are more consistent with the characteristic-earthquake statistics suggested by some studies for individual fault segments (Bakun & McEvilly, 1984;Parsons, 2008;Schwartz & Coppersmith, 1984). Further study is needed to determine whether reproducing these more detailed frequency-magnitude distributions, especially for the largest events and perhaps with finer magnitude binning ( Figure S6 in Supporting Information S1), would allow simulations to capture the rates of single and two-segment ruptures, and how close the agreement would have to be. However, the practical implications of such a study could be limited due to limited statistical information about large earthquakes for a given fault network of interest, due to their infrequent occurrence and relatively short periods of reliable earthquake records.

Resolution and Convergence of SEAS Simulations of Faults With Higher Instability Ratios
As discussed in Section 3, we find that the discretization required to achieve long-term numerical convergence in simulations of fault model M1, with instability ratio of ∕ℎ * ≈ 21 , is more stringent than the current standards based on simulations of single dynamic ruptures and shorter SEAS simulations with lower instability ratios (Day et al., 2005;Lapusta & Liu, 2009). It has been demonstrated that fault models with relatively low instability ratios can result in quasi-periodic behavior, as seen in fault model M1 (Figure 2), whereas increasing the instability ratio can lead to more variable sequences of events with partial-segment ruptures of different rupture size, potentially consistent with Gutenberg-Richter scaling (e.g., Cattania, 2019; Lapusta & Rice, 2003;Lapusta et al., 2000;Michel et al., 2017;Wu & Chen, 2014). As simulations with higher instability ratios can produce ruptures with a wider variety of rupture sizes, with the rupture size depending on the prestress conditions before rupture nucleation, one could hypothesize that simulations of fault models with higher instability ratios may be more sensitive to how the evolution of shear stress is resolved over long-term fault behavior.
To test that, let us consider sequences of events in fault model M5 (Table 2), which has smaller characteristic slip distance, hence smaller nucleation size ( * ≈ 603 m), and larger instability ratio ( ∕ℎ * = 53 vs. 21 in M1). Interestingly, we find that the long-term sequence of simulated events in this model is not the same for finely discretized simulations with cell sizes of 25, 12.5, and 6.25 m (Figure 10), in which the quasi-static cohesive zone Λ0 is resolved by 18, 36, and 72 cells, respectively. The simulations produce nearly identical fault behavior for the first several hundred years of simulated time, but then eventually begin to differ (Figures 10a-10c).
Let us consider the first event in the three simulations of model M5 with fine discretization (Figures 10a-10c), which all have the same initial conditions. If we examine the local evolution of shear stress versus slip at two spatial points in the simulations, the results are virtually identical (Figures S7a and S7b in Supporting Information S1), suggesting that a single dynamic rupture in these finely discretized simulations is adequately resolved. The evolution of shear stress and slip rate at the rupture front with time is also well resolved for each individual simulation. While the different spatial resolutions result in small variations in the timing and magnitude of the resolved properties at specified locations ( Figure S7c-S7f in Supporting Information S1), these differences are well within of what is considered well resolved and convergent in prior studies (e.g., Day et al., 2005). Early in the rupture, shortly after nucleation (near = 30 km), the rupture front is almost identical in the three simulations. (Figures S7c and S7e in Supporting Information S1). As the rupture continues, small numerical differences for different resolutions result in minor differences in the rupture, such as less than 0.08% difference in the rupture arrival time and 2% difference in the peak slip 19 of 33 rate between the two best-resolved simulations at the location close to the end of the rupture (Figures S7d and S7f in Supporting Information S1). Such minor differences arise even for fine resolutions due to cumulative effects of slightly different representations of the solution by the discrete cells; for example, the fixed computational cells sample slightly different portions of the passing rupture front, leading to small accumulating differences in the magnitude of the shear stress and slip rate.
These small differences, that do not substantially alter the resulting rupture characteristics of individual events, do eventually alter the resulting earthquake sequences. For several ruptures early on in finely discretized simulations, the slip and shear stress distributions before and after individual events are virtually indistinguishable (Figures 11a and 11b). However, eventually the small variations accumulate, resulting in enough differences in prestress conditions to cause more substantial differences in rupture lengths and amounts of slip within individual events, as well as changes in timing and location of earthquake nucleations (Figures 11c-11e). As a result, the long-term history of sequences of slip events is altered (Figures 11f), including the rate of ruptures that jump across the VS barrier. We hypothesize that this alteration occurs for higher but not lower instability ratios due to more complex earthquake sequences in the latter case, although this issue requires further study.
Despite the specific sequences of events being different in the finely discretized simulations shown in Figures 10a-10c, we do find that certain outcomes are quite similar between these simulations, such as relationships between average static stress drop and seismic moment, average slip and rupture length, and breakdown energy and average slip, as well as general characteristics of the evolution of average shear stress and shear heating with time ( Figure 12). Other parameters, such as the rate of ruptures jumping from one fault segment to another, are sensitive to numerical resolution even in these finely resolved simulations, although they have relatively similar values (from 0.64 to 0.78). This highlights how the criteria for adequate discretization in numerical simulations can depend on both the physical problem being considered and the outcome of interest. Note that while it is plausible that further discretization of fault model M5 would result in eventual convergence, and thus potentially a true rate of two-segment ruptures, the spatial discretization considered in this study is already much finer than those considered in most numerical SEAS studies, especially in more realistic models of 2D faults in 3D media which are often challenged to resolve Λ0 by even three cells.
While the specific rate of ruptures jumping across the VS barrier varies among these finely discretized simulations of fault model M5, it is possible that some broader statistical features of the jump rate are more robust. We examine the frequency-magnitude and 2000-year jump rate statistics for the long-term sequences of events in simulations of model M5 with different discretization. While the distributions mildly vary among finely discretized simulations with differing cell sizes (12.5 and 25 m), they are comparable (Figure 13 and Figure S8 in Supporting Information S1). Thus, one can ascertain information about the probability distribution for the rate of multi-segment ruptures, even if specific results vary due to numerical discretization. Such small numerical perturbations could potentially be considered representative of various sources of physical perturbations on natural faults, and the statistical consistency of the distributions could be explored by producing ensembles of simulations with varying initial conditions. However, our results suggest that it is still important to sufficiently resolve the rupture process as the statistical distributions for rupture properties in simulations using oversized cells can be more substantially impacted by numerical artifacts and considerably vary from simulations with finer discretization (Figures 13 and 14). Figure 10. Sequences of earthquakes and rates of two-segment ruptures over 4,000 years in fully dynamic simulations with different resolution of fault model M5 with higher instability ratio. Seismic slip is contoured every 0.5 s with ruptures jumping across the VS barrier colored blue. (a-c) Slip history for increasingly better-discretized simulations. While the initial 1,000 years of simulated behavior appear well resolved and comparable, longer-term simulations begin to diverge due to the compounded effects of small numerical differences, leading to similar but inconsistent jump rates across the barrier. (d and e) The spatial distribution of shear stress at the rupture front. For well-resolved simulations (d), the cohesive zone is resolved by several cells, but is resolved by less than even one cell for poorly resolved simulations (e). (f and g) Simulations with decreasing numerical resolution can exhibit additional artificial complexity and substantially different long-term fault behavior, including different rates of two-segment ruptures.

Resolution and Convergence in SEAS Simulations With Moderate Rupture Speeds Due To an Approximation for Off-Fault Plasticity
While the 2-D fault models discussed in this study can be considered relatively simple, in some ways they can be particularly challenging to resolve. In fault models with purely elastic bulk, dynamic ruptures are able to accelerate to rupture speeds close to the limiting values (e.g., Figure 7 for fault model M1), making (c-e) Eventually, small differences in shear stress before events build up due to different numerical approximations, resulting in small differences in slip and rupture length for individual events, as well as the location and timing for the nucleation of smaller events. (f & g) The differences in shear stress accumulate over sequences of events, resulting in noticeable variations of slip in larger events after 800 years of simulated time and, eventually, different histories of large segmentspanning events between the two well-resolved simulations, as shown in Figure 10. 21 of 33 it difficult to resolve the significantly shrinking cohesive zone Λ . For example, during fully dynamic ruptures in simulations of fault model M5, the rupture speed approaches 0.99 and the cohesive zone shrinks more than seven times to about 63.5 m. In real rocks, high slip rates and hence high strain rates associated with dynamic rupture would be mitigated by off-fault inelastic behavior around the rupture front, which would contribute to limiting the rupture speed (Andrews, 2004;Dunham et al., 2011a).
In order to examine how conditions for resolution and convergence may differ in long-term SEAS simulations with more moderate rupture speeds, we approximate the effects of off-fault yielding by employing a limit on the slip velocity, as suggested by Andrews (2004) and discussed in detail in Lambert, Lapusta, Perry (2021). Such slip velocity limit is meant to capture the effect on rupture propagation of limiting off-fault elastic strains around the rupture front, as would be expected with plastic yielding. We consider long-term fully dynamic simulations of fault model M5 with the slip velocity limited to 2 m/s in order to maintain rupture speeds around 0.8 , consistent with the cohesive zone shrinking by about a factor of 2 from the quasi-static estimate.
Surprisingly, the finely discretized simulations of fault model M5 with limited rupture speed still produce differing sequences of events, despite the rupture front and local behavior being well resolved and nearly identical for cell sizes of 6.25-25 m (Figure 15 and Figure S9 in Supporting Information S1). As with the standard fully dynamic simulations without the plasticity approximation, well-resolved simulations of fault model M5 with the velocity limit are nearly identical for the initial few sequences (Figure S10a-S10b in Supporting Information S1). However, the sequences of events begin to differ due to slight differences in how the evolution of shear stress is resolved during a slow-slip transient within the nucleation region of an impending rupture, resulting in a 3-year delay between the nucleation of the subsequent rupture in each simulation ( Figure S10c and S10d in Supporting Information S1). As discussed earlier for the standard fully dynamic simulations, the small differences in prestress lead to mild differences in slip and rupture size in subsequent events, which eventually compound to produce more substantial variations in the long-term sequences of events (Figure S10e-S10h in Supporting Information S1). These results once again illustrate the extreme sensitivity of the long-term sequences of events, and rates of two-segment ruptures, in this (c-f) Simulations with marginal or inadequate resolution have enhanced production of smaller events, as small groups of cells nucleate into ruptures but fail to propagate substantially due to poorly resolved stress concentration at the rupture front. The 20,000-year jump rates and 2000-year jump rate distributions substantially vary for simulations using oversized cells compared to the well-resolved simulations. 23 of 33 highly nonlinear problem, as well as the significance of resolving how aseismic processes load, relax, and redistribute stress along faults.
It is possible that sharp edges in the slip-rate function due to the velocity-limit approximation for off-fault plasticity may introduce higher-frequency contributions to the elastic stress interactions that are challenging to resolve. Numerical studies have also shown that the explicit consideration of off-fault plasticity in numerical earthquake models requires special consideration; for example, some common plasticity models, such as rate-independent Drucker-Prager plasticity, can be mathematically ill-posed for some stress states and require a regularization to determine unique solutions with sufficient numerical refinement (Andrews, 2005;Dunham et al., 2011a;Erickson et al., 2017;Rudnicki & Rice, 1975). Thus, while incorporating additional sources of dissipation like off-fault plasticity may assist in limiting the rupture speed and making the dynamic cohesive zone easier to resolve, such additional physics may require its own special consideration to achieve numerical convergence. Further investigation of appropriate approximations for off-fault plasticity and other off-fault behavior using detailed dynamic rupture simulations would facilitate the efficient incorporation of such physics into long-term SEAS simulations.
Interestingly, we see similar lack of convergence in quasi-dynamic simulations of fault model M5, where long-term sequences, including the rate of two-segment ruptures, differ in seemingly well-resolved simulations due to the compounded effects of small numerical differences (Figures S11 and S12 in Supporting Information S1). Moreover, despite the rupture front being better resolved in the quasi-dynamic simulations and in fully dynamic simulations with the plasticity approximation than in the standard fully dynamic simulations, the sequences of events begin to diverge earlier. Specifically, while the standard fully dynamic 25 of 33 simulations of fault model M5 with cell sizes of Δ = 6.25 and Δ = 12.5 m have the same event sequences through ∼600-700 years of simulated time, fully dynamic simulations with the plasticity approximation begin to substantially differ between 200 and 300 years, and quasi-dynamic simulations begin to noticeably differ between 100 and 200 years.
A potential explanation for this finding is that both the quasi-dynamic approximation and strong limitation on slip rate for fully dynamic simulations also limit the magnitude of the stress transfer along the fault ( Figure S13 in Supporting Information S1), making the simulations more sensitive to small numerical differences. Thus, while the lower stress concentrations in both cases facilitate maintaining slower ruptures and resolving the breakdown of shear resistance at the rupture front, the smaller magnitudes for the stress transfer along the fault makes rupture propagation more sensitive to variations in the preexisting shear stress ahead of the rupture front. Note that while the approximation for off-fault plasticity substantially limits the peak slip rate and magnitude of the stress transfer along the fault, the overall stress transfer for the fully dynamic rupture including the plasticity approximation is still more pronounced than that of the quasi-dynamic ruptures, and remains more pronounced well behind the rupture front due to the continued arrival of waves from ongoing slip in already-ruptured regions. Both the quasi-dynamic simulations and the fully dynamic simulations with the plasticity approximation produce comparable static stress drops and frequency-magnitude statistics to the standard fully dynamic simulations (Figures S8, S11, and S14 in Supporting Information S1). However, the rupture speeds and rates of two-segment ruptures are consistently higher for the fully dynamic simulations due to the substantially larger stress transfer. These results emphasize the significance of inertial effects when considering how ruptures navigate various forms of fault heterogeneity. The simulations of model M5, without and with the plasticity approximation, provide another example of how earthquake sequences with similar frequency-magnitude statistics can result in different jump rates across the velocity-strengthening barrier. While the simulations with cell sizes of 6.25, 12.5, and 25 m have well-resolved cohesive zones (Figures 10 and 15) and similar event statistics (Figures S8 and S14 in Supporting Information S1), they have jump rates ranging from 0.7 to 0.8 without the plasticity approximation to 0.3-0.5 with the plasticity approximation (Figures 10 and 15).

Conclusions and Discussion
We have investigated the sensitivity of numerical simulations of long-term sequences of earthquakes and aseismic slip (SEAS) to numerical discretization and treatment of inertial effects, using a simplified 2-D model of a 1D fault with two co-planar seismogenic, VW segments separated by a VS barrier. Our focus is, in part, on the resulting rate of rupture jumps across the barrier.
We find that the convergence of long-term simulated earthquake sequences with increasing numerical resolution may not always be achievable, at least practically given current computational capabilities. Even if simulations are sufficiently discretized to produce consistent modeling results for individual ruptures or short sequences of events, they can still produce different long-term sequences due to compounded effects of small numerical differences over many events. We have achieved the convergence for fault models with lower instability ratios ∕ℎ * , that is, lower fault lengths in comparison to the nucleation size ( Figure 3). In contrast, models with higher instability ratios exhibit different long-term behavior even in simulations that are well discretized by standard metrics (Day et al., 2005;Lapusta & Liu, 2009), including different specific sequences of earthquakes and different probability of ruptures jumping across the VS barrier. In the cases with convergent long-term behavior, the criteria for numerical resolution that leads to the same evolution of slip are more stringent than those for individual dynamic ruptures, that is, the dynamic cohesive zone size needs to be discretized by more cells. Figure 15. Sequences of earthquakes and rate of two-segment ruptures over 4,000 years in fully dynamic simulations with different resolution of fault model M5 and an approximation of off-fault plasticity. The rupture speed reduces to 0.8 due to the approximation using a velocity limit of lim = 2 m/s. Seismic slip is contoured every 0.5 s with ruptures jumping across the VS barrier colored blue. (a-c) Slip history for increasingly well-resolved simulations. The initial few sequences of events appear comparable among well-resolved simulations, however the sequences begin to differ due to the compounded effects of small numerical differences. (d and e) The cohesive zone shrinks by only about a factor of two for rupture speeds below 0.8 , so the rupture front is very wellresolved. (f and g) Simulations with decreasing numerical resolution exhibit additional artificial complexity and substantially different long-term fault behavior, including rates of two-segment ruptures.

of 33
Our results show that numerical convergence in SEAS simulations depends not only on how well important length-scales are discretized but also on the sensitivity of the specific physical problem to small numerical perturbations. In particular, our results suggest that faults with higher instability ratios are more sensitive to accumulating numerical perturbations (Figure 16), although that conclusion is reached here for simulations governed by only standard rate-and-state friction and requires further study in models with different/ additional fault physics. While quasi-dynamic simulations are easier to resolve and thus should result in smaller numerical discrepancies for sufficiently small cell sizes, the milder stress transfer compared to fully dynamic ruptures can make long-term quasi-dynamic simulations more sensitive to small perturbations in shear stress, as occurs in fault model M5. In contrast, while earthquake sequence simulations including enhanced dynamic weakening can require finer discretization to resolve the rapid weakening rates at the rupture front compared to standard rate-and-state friction, the larger dynamic stress changes in fault models with more efficient weakening can make rupture propagation less sensitive to different levels of fault heterogeneity, promoting the propagation of predominantly larger ruptures and hence resulting in relatively simple earthquake sequences even for models with high instability ratios (Lambert, Lapusta, Perry, 2021;Lambert, Lapusta, & Faulkner, 2021). Hence empirical discretization criteria, such as those of (Day et al., 2005), should be seen as guidelines that may not be universally applicable to all physical models and outcomes of interest. Moreover, for some models, there is the possibility that numerical convergence of Figure 16. Conceptual diagram illustrating potentially convergent versus divergent numerical behavior depending on resolution and model complexity, parameterized by the instability ratio as an example. Well-resolved fault models with low enough instability ratio may potentially be numerically deterministic where adequate discretization results in virtually indistinguishable numerical outcomes. Fault models with higher instability ratio may either have more stringent requirements for numerical discretization in order to achieve long-term convergence, or such convergence may be impossible; either way, achieving numerical convergence in simulations of sufficiently complex fault models, such as with higher instability ratios, would be impractical. In such cases, it may still be possible to achieve statistical consistency among some outcomes within well-resolved simulations, though other properties of the system may be highly sensitive to numerical precision and considerably vary depending on the numerical procedures.

10.1029/2021JB022193
27 of 33 long-term slip may not be achievable, though statistical consistency may hold for some modeling results but not others (Figure 16).
For the fault models considered, we find that the specific sequences of slip events and rate of earthquake ruptures jumping across a VS barrier are sensitive to the numerical resolution, representation of inertial effects, as well as minor changes in physical properties, such frictional parameters, confining stress, seismogenic depth, and barrier size. This suggests that, at least in this relatively simple model, the rate of ruptures jumping across a VS barrier is not a stable outcome that can always be reliably estimated from numerical simulations, unless the barrier is so large or small that the rate is reliably zero or 1 (Figure 17). The sensitivity of the specific timing and sequences of slip events in models with higher instability ratios is suggestive of deterministically chaotic long-term fault behavior, as has been suggested by theoretical studies of coupled spring-sliders and continuum models with increasing instability ratios (Barbot, 2019;Huang & Turcotte, 1990;Kato, 2014). The sensitivity of rupture jump rates to small changes in our models suggests that the jump rates across barriers that serve as earthquake gates may also be highly sensitive to small physical perturbations on natural faults, and thus may be impractical to estimate in a reliable manner.
However, even for the models that do not achieve deterministic convergence with finer resolution, we find that some characteristics of well-resolved simulations are preserved, qualitatively and quantitatively. The characteristics include ranges of average source properties such as the average static stress drop, quantities related to energy partitioning such as the average breakdown energy, as well as general features of the average shear stress and shear heating evolution throughout time (Figure 12). These results suggest that some 28 of 33 aspects of physical systems may be reliably determined from a given physics-based model, while others perhaps cannot, in the sense that they are very sensitive to numerical procedures and initial conditions, and even well-resolved models produce different outcomes with respect to those quantities. Our findings also suggest that it may be possible to discern some statistical aspects of the probability distribution for multisegment ruptures from well-formulated numerical models, even if they do not exhibit convergence of long-term behavior with numerical resolution. However, as the jump rate appears to be sensitive to small perturbations in numerical and physical properties, it would be prudent to examine the statistical consistency of the jump rate distribution through large ensembles of models, particularly given the uncertainty in values for physical model parameters. Another route for examining plausible rupture scenarios for large earthquakes navigating key sections of fault networks would be to conduct detailed dynamic rupture simulations that can handle more realistic fault geometries with full treatment of inertial effects (e.g., Ulrich et al., 2019;Wollherr et al., 2019), and produce large ensembles of dynamic rupture scenarios with variations in initial conditions inspired by SEAS simulations.
Our results confirm that quasi-dynamic simulations that ignore wave-mediated stress transfer during dynamic rupture can lead to qualitative differences in the resolved rupture behavior and long-term sequences of slip events. The wave-mediated stress redistribution not only facilitates long-range interactions among portions of a fault and neighboring segments, but also alters the state of stress at the rupture front, promoting higher slip rates and more focused stress concentrations. In particular, the relatively small static stress transfer in quasi-dynamic simulations makes the rupture front more susceptible to unfavorable conditions, such as those one may expect from frictional heterogeneity, fault roughness, and regions of unfavorably low prestress. In contrast, the larger wave-mediated dynamic stresses in fully dynamic ruptures may assist rupture propagation in navigating unfavorable fault conditions and geometric irregularities (Ando & Kaneko, 2018;Douilly et al., 2015;Duan & Oglesby, 2006;Dunham et al., 2011b;Galvez et al., 2014;Harris & Day, 1993, 1999Harris et al., 1991;Kame et al., 2003;Lozos et al., 2015;Thomas et al., 2014;Ulrich et al., 2019;Withers et al., 2018;Wollherr et al., 2019). Moreover, the spatial pattern for dynamic stresses, which affects the preferential direction for ruptures to branch or jump to neighboring faults, rotates as a function of the rupture speed, and hence can be considerably different from a quasi-dynamic rupture (Kame et al., 2003). Thus, considering full inertial effects during individual dynamic ruptures and longterm sequences of slip events is particularly important when considering the interaction of multiple fault segments and the likelihood of ruptures propagating through fault portions with potentially unfavorable conditions.
Our results also confirm that using increasingly oversized cells, with or without wave-mediated stress transfers, results in a progressively more complex slip response, with broader distributions of event sizes, consistent with conclusions from prior studies (Ben-Zion & Rice, 1995). Using oversized cells and/or ignoring wave-mediated stress transfer significantly modifies the probability of two-segment ruptures, as well as the resulting earthquakes sequences.
Finally, we find that simulations with similar b-values of the cumulative frequency-magnitude statistics and similar average static stress drops can have substantially different rates of ruptures that jump between two segments. We caution that models with more realistic geometry and/or additional physics may behave differently; for example, more realistic non-planar fault geometry or enhanced dynamic weakening may dominate the model outcomes (e.g., the rupture jump rates) to the point of making them insensitive to certain variations in model assumptions or parameters. However, our findings suggest that such insensitivity of model outcomes cannot be assumed and needs to be established through suites of simulations with different plausible properties and assumptions.
In order to constrain the most relevant fault physics and property ranges, and hence to maximize the predictive power of numerical earthquake models for future fault behavior, numerical simulations of fault models can explore how choices in physical mechanisms and conditions are reflected in a range of observations taken together, including, but not limited to, geodetic observations of interseismic and postseismic motions (e.g., Barbot et al., 2012), geological records of past earthquakes (e.g., Sieh, 1978;Biasi & Wesnousky, 2017), heat flow constraints (e.g., Henyey & Wasserburg, 1971;Lachenbruch & Sass, 1980;Fulton et al., 2013), and seismologically determined properties of earthquakes such as average static stress drops, radiated energy, and the inferred increase in average breakdown energy with rupture size (Abercrombie & Rice, 2005;Lambert, Lapusta, Perry, 2021;Perry et al., 2020;Rice, 2006;Viesca & Garagash, 2015). As modeling capabilities continue to develop and capture realistic fault geometries tailored for specific fault segments, such models can also aim to reproduce details of more local seismicity patterns, including the interplay between earthquakes and slow-slip transients (e.g., Ito et al., 2013;Kato et al., 2012;Ruiz et al., 2014), repeating earthquakes (e.g., Chen & Lapusta, 2009;Mele Veedu & Barbot, 2016), seismic swarm activity (e.g., Ross et al., 2020), and spatiotemporal variations in frequency-magnitude statistics (Ishibe & Shimazaki, 2012;Schorlemmer & Wiemer, 2005). While the frequency-magnitude distribution of seismicity over broad regions like California are generally consistent with Gutenberg-Richter scaling with typical b-values near unity (Field et al., 2013), it remains a topic of active debate as to whether such scaling applies to individual fault segments and their immediate surroundings (Field et al., 2017;Ishibe & Shimazaki, 2012;Kagan et al., 2012;Page & Felzer, 2015;Page & van der Elst, 2018;Schwartz & Coppersmith, 1984;Stirling et al., 1996;Wesnousky, 1994). In particular, some segments of major mature faults, such as the Cholame and Carrizo segments of the San Andreas Fault, exhibit notable deviations from typical Gutenberg-Richter scaling and are nearly absent of microseismicity (Bouchon & Karabulut, 2008;Hauksson et al., 2012;Jiang & Lapusta, 2016;Michailos et al., 2019;Sieh, 1978;Wesnousky, 1994). Determining plausible physical conditions compatible with seismicity patterns on such large mature fault segments, including observed seismic quiescence, is of particular importance for seismic hazard assessment, since the segments have historically hosted some of the largest recorded continental earthquakes.
Such physics-based modeling can then inform laboratory and observational scientists about the most strategic measurements that can improve constraints on plausible models and help assess the resolution and sources of potential bias in geophysical inferences (e.g., Abercrombie, 2021;Lin & Lapusta, 2018;. Note that a number of physical properties not included in our simplified 2-D models may qualitatively alter the behavior and hence interaction of neighboring fault segments, such as the explicit consideration of depth variations in slip and the depth extent to which ruptures propagate (e.g., Jiang & Lapusta, 2016;Ulrich et al., 2019;Wollherr et al., 2019), time-dependent variations in loading from distributed deformation at depth (Allison & Dunham, 2018;Lambert & Barbot, 2016), enhanced dynamic weakening at seismic slip rates (Di Toro et al., 2011;Dunham et al., 2011b;Lambert, Lapusta, Perry, 2021;Perry et al., 2020;Tullis, 2007), evolving damage and fault geometry (Lyakhovsky & Ben-Zion, 2009Okubo et al., 2019), especially near fault junctions, and hydromechanical coupling due to poroelasticity, fault-zone dilatancy and other fluid effects (Dunham & Rice, 2008;Segall & Rice, 1995;Zhu et al., 2020). These are just a few physical ingredients that merit detailed study in the long-term interaction of fault segments.
Overall, our findings highlight the importance for numerical studies to examine the sensitivity of their outcomes of interest to variations in their modeling parameters, particularly when assessing their predictive value for seismic hazard assessment. Community initiatives, such as the Southern California Earthquake Center (SCEC) code comparisons for dynamic rupture simulations and simulations of sequences of seismic and aseismic slip (Barall & Harris, 2014;Erickson et al., 2020;Harris et al., 2009Harris et al., , 2018, can provide further insight into how numerically derived results for different physical quantities may depend on numerical methodologies and computational practices. The significant sensitivity of the rate of multi-segment ruptures to small changes in our numerical models implies that this hazard parameter may also be sensitive to physical perturbations on natural faults. This consideration motivates further evaluation of metrics for describing long-term fault behavior and assessing seismic hazard, tasks for which physics-based modeling is well-suited.

Data Availability Statement
Details about the fault models and numerical parameters for calculations are provided in the main and supplementary text, and Tables 1-2. Data were not used, nor created for this research. This study was supported by the National Science Foundation (grants EAR 1142183 and 1520907) and the Southern California Earthquake Center (SCEC), contribution No. 10089. SCEC is funded by NSF Cooperative Agreement EAR-1033462 and USGS Cooperative Agreement G12AC20038. Numerical simulations for this study were carried out on the High Performance Computing Center cluster of the California Institute of Technology. This study was motivated by insightful discussions within the SCEC community.