Scale Dependence of Earthquake Rupture Prestress in Models With Enhanced Weakening: Implications for Event Statistics and Inferences of Fault Stress

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 (a) illustrate that large earthquakes can occur on faults that appear not to be critically stressed compared to stresses required for slip nucleation; (b) highlight the importance of finite‐fault modeling in relating the local friction behavior determined in the lab to the field scale; and (c) 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.


Introduction
Determining the absolute level and controlling factors of the stress state on faults has profound implications for earthquake physics, seismic hazard assessment, and the role of faulting in plate tectonics and geodynamics. Numerous lines of field evidence suggest that the average shear stress acting on mature faults must be low, 20 MPa or less, in comparison to the expected shear resistance of 100-200 MPa averaged over the seismogenic depth, given rock overburden and hydrostatic pore fluid pressure, along with typical quasi-static friction coefficients of 0.6-0.85 (aka "Byerlee friction") measured in laboratory experiments (Brune Abstract 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 (a) illustrate that large earthquakes can occur on faults that appear not to be critically stressed compared to stresses required for slip nucleation; (b) highlight the importance of finite-fault modeling in relating the local friction behavior determined in the lab to the field scale; and (c) 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.

Plain Language Summary
Understanding the evolution of stress on faults over periods of slow and fast motion is crucial for assessing how earthquakes start, grow, and ultimately stop. Here we use computer models to explore the stress conditions required for simulated earthquake ruptures to occur. We find that the critical stress conditions for rupture propagation depend on the size of the rupture and how efficiently the fault shear resistance weakens during fast slip. In particular, larger earthquakes on faults that experience increasingly more efficient weakening during rupture can propagate under systematically lower stress conditions than those required for rupture nucleation. As a result, we find that faults that exhibit efficient weakening can host predominantly large earthquakes at the expense of smaller earthquakes. Our findings illustrate how large earthquakes can occur on faults that appear to be understressed compared to expected conditions for rupture nucleation. Moreover, our findings support a body of work suggesting that the scarcity of small earthquakes on some major mature fault segments, like the central section of the San Andreas Fault, may indicate that they experience substantial weakening during large earthquakes, a consideration that may be particularly useful for earthquake early warning systems.
LAMBERT ET AL.
A relatively straightforward explanation for the low-stress operation of mature faults is that they may be persistently weak (Figure 1), due to the presence of anomalously low quasi-static friction coefficients and/ or low effective normal stress from pervasive fluid overpressure (Bangs et al., 2009;Brown et al., 2003;Collettini et al., 2009;Faulkner et al., 2006;Lockner et al., 2011). Mature fault zones typically have a layer of fine gouge material at their core (e.g., Faulkner et al., 2006). However, most materials with low quasi-static friction coefficients (less than 0.5) under laboratory conditions tend to exhibit velocity-strengthening behavior (Ikari et al., 2011), which would preclude spontaneous nucleation of dynamic ruptures. Moreover, while evidence of substantial fluid overpressure has been documented for many subduction zones (Bangs et al., 2009;Brown et al., 2003), there remains much debate over the ubiquity of chronic near-lithostatic fluid overpressurization along faults in other tectonic settings, such as continental faults, with some borehole Figure 1. Field observations suggest that the average effective friction on mature faults must be low ( E  0.1). One explanation for this inferred low effective friction would be that mature faults are persistently weak, such as from the presence of fault materials with persistently low friction coefficients / ( ) E p An alternative hypothesis for explaining such low-stress, low-heat operation is that mature faults are indeed strong at slow, quasi-static sliding rates but undergo considerable enhanced dynamic weakening at seismic slip rates, which has been widely hypothesized in theoretical studies and documented in laboratory experiments ( Figure 1, dashed black line; Acosta et al., 2018;Di Toro et al., 2011;Noda et al., 2009;Rice, 2006;Sibson, 1973;Tsutsumi & Shimamoto, 1997;Wibberley et al., 2008). The presence of enhanced dynamic weakening on natural faults can be questioned by the expectation that enhanced dynamic weakening would result in magnitude-dependent static stress drops (Beeler et al., 2012;Gao & Wang, 2014), with larger ruptures resulting in higher stress drops than typical values of 1-10 MPa inferred from earthquakes on natural faults (Allmann & Shearer, 2009;Ye et al., 2016b). The expectation is based on a common assumption that the shear prestress over the entire rupture area should be near the static strength of the fault while the final shear stress should be near the dynamic resistance of the fault, resulting in a large static stress change for more efficient weakening. However, a number of numerical and laboratory studies have demonstrated that, once nucleated, dynamic ruptures can propagate under regions with prestress conditions that are well below the expected static strength, based on prescribed or measured quasi-static friction coefficients and confining conditions (Dunham et al., 2011;Fineberg & Bouchbinder, 2015;Gabriel et al., 2012;Lu et al., 2010;Noda et al., 2009;Schmitt et al., 2015;Zheng & Rice, 1998) while the final shear stress could be higher than dynamic shear stress for pulse-like ruptures, with both inferences promoting reasonable stress drops. Such studies have often considered a single dynamic rupture nucleated artificially and propagating over uniform prestress conditions.
Recent numerical studies of earthquake sequences have shown that fault models with a combination of both hypotheses for low-stress operation, including some chronic fluid overpressure as well as mild-to-moderate enhanced dynamic weakening due to the thermal pressurization of pore fluids, work well for reproducing a range of observations (Lambert et al., 2021;Perry et al., 2020). These include reasonable static stress drops between 1 and 10 MPa nearly independent of earthquake magnitude, the seismologically inferred increase in average breakdown energy with rupture size, the radiation ratios between 0.1 and 1 inferred for natural events, and the heat flow constraints. The simulations produce mainly crack-like or mild pulse-like ruptures, with no significant undershoot. The near magnitude-invariance of average static stress drop arises in these fault models because enhanced dynamic weakening results in both lower average prestress and lower average final shear stress for larger ruptures with larger slip, with the average static stress drops being nearly magnitude-independent. These studies suggest that distinguishing between the conditions required for rupture nucleation and propagation is important for assessing the relationship between laboratory friction measurements, seismological observations and the absolute stress conditions on faults.
Here, we use and expand upon the set of numerical models from Perry et al. (2020) and Lambert et al. (2021) to document the variability of prestress on a fault that arises from the history of previous ruptures, and to study the relation between the size of dynamic rupture events and the average shear prestress over the rupture area. We also examine how the complexity of earthquake sequences, in terms of the variability of rupture size, differs with the efficiency of dynamic weakening. We study these behaviors in the context of simulations of sequences of earthquakes and slow slip, which allow the prestress conditions before earthquakes to be set by the loading conditions, evolving fault shear resistance (including weakening and healing), and stress redistribution by prior slip, as would occur on natural faults. Moreover, our simulations resolve the spontaneous nucleation process with the natural acceleration of slow unsteady slip prior to dynamic rupture. The constitutive relations for the evolving fault resistance and healing adopted in our models have been formulated as a result of a large body of laboratory, field and theoretical work (e.g., Dieterich, 1979;Di Toro et al., 2011;Rice, 2006;Ruina, 1983;Sibson, 1973;Wibberley et al., 2008). Indeed, laboratory experiments of fault shear resistance at both slow and fast slip rates have been indispensible for our understanding of fault behavior and for formulating fault models such as the ones used in this study. The modeling allows us to examine the implications of the laboratory-derived constitutive behaviors for the larger-scale behavior of faults, and we compare our inferences of average shear prestress from relatively large-scale finite-fault modeling to field measurements of crustal stresses acting on mature faults and small-scale laboratory measurements of the shear resistance of typical fault materials.
LAMBERT ET AL.

Building on Laboratory Constraints to Model Larger-Scale Fault Behavior
Laboratory experiments have been instrumental for exploring aspects of fault resistance during both slow and fast sliding ( 9 10 E  -1 m/s, Figure 1). Experiments with slow sliding velocities ( 3 10 E   m/s) are critical for formulating fault constitutive laws that form the basis for understanding the nucleation of earthquake ruptures. High-velocity laboratory friction experiments have demonstrated enhanced dynamic weakening of faults and elucidated a range of mechanisms by which this dynamic weakening can occur (e.g., Acosta et al., 2018;De Paola et al., 2015;Di Toro et al., 2011;Faulkner et al., 2011;Goldsby & Tullis, 2011;Han et al., 2007;Wibberley et al., 2008). Most slow-and high-velocity experiments measure or infer the relevant quantities-slip, slip rate, shear stress etc.-averaged over the sample and examine the evolution of shear resistance corresponding to a particular history of loading, such as imposed variations in the displacement rate of the loading piston, and the particular fault conditions (normal stress, temperature, pore fluid pressure, etc.). Some experimental studies imposed the expected sliding motion during earthquakes in order to directly relate laboratory stress measurements to seismological quantities, such as static stress drop and breakdown energy (e.g., Fukuyama & Mizoguchi, 2010;Nielsen et al., 2016;Sone & Shimamoto, 2009).
To understand the full implications of the evolution of shear resistance measured in small-scale experiments for slip at larger scales along natural faults, they are synthesized into mathematical formulations and used in numerical modeling, for the following reasons. During slipping events on a finite fault over scales of tens of meters to kilometres-much larger than the experimental scale-the fault does not slip uniformly with a predetermined slip-rate history. Rather, the slip event initiates on a portion of the fault and then spreads along the fault, with varying slip-rate histories and final slips at different points along the fault. This is captured in inversions of large earthquakes (e.g., Heaton, 1990;Simons et al., 2011;Tinti et al., 2016;Ye et al., 2016a) and, to a degree, in larger-scale experiments, sometimes involving analog materials (Lu et al., 2010;McLaskey et al., 2014;Rubino et al., 2017;Svetlizky & Fineberg, 2014;Yamashita et al., 2015). In the process, the slip (a) transfers stress to the more locked portions of the fault and (b) enters portions of the fault with different conditions-such as levels of shear pre-stress, pore fluid pressure, etc.-and potentially different friction and hydraulic properties. Hence the resulting coupled evolution of shear resistance ( , ) ln ln , 10.1029/2021JB021886 5 of 29 models are governed by a form of the laboratory-derived Dieterich-Ruina rate-and-state friction law regularized for zero and negative slip rates (Lapusta et al., 2000;Noda & Lapusta, 2010). The evolution of the state variable can be described by various evolution laws; we employ the aging law (Ruina, 1983): which describes evolution during sliding as well as time-dependent healing in near-stationary contact. In our models, the shear resistance and shear stress also change due to the evolution of pore fluid pressure E p .
We conduct numerical simulations following the methodological developments of Lapusta et al. (2000), Noda and Lapusta (2010), and Lambert et al. (2021) in order to solve the elastodynamic equations of motion with the fault boundary conditions, including the evolution of pore fluid pressure and temperature on the fault coupled with off-fault diffusion. The simulations solve for mode III slip on a 1-D fault embedded into a 2-D uniform, isotropic, elastic medium ( Figure 2). The potential types of slip on the fault include sequences of earthquakes and aseismic slip (SEAS) and they are simulated in their entirety, including the nucleation process, dynamic rupture propagation, postseismic slip that follows the event, and interseismic period between seismic events that can last up to tens or hundreds of years and host steady and transient slow slip (Figure 2).
The simulated fault in our models contains a 24-km-long segment with velocity-weakening (VW) frictional properties where earthquake ruptures may nucleate and propagate, surrounded by velocity-strengthening (VS) segments that inhibit rupture nucleation and propagation. Our simulations include enhanced dynamic weakening due to the thermal pressurization of pore fluids, which occurs when pore fluids within the fault shearing layer heat up and pressurize during dynamic rupture, reducing the effective normal stress and shear resistance (Noda & Lapusta, 2010;Rice, 2006;Sibson, 1973). Thermal pressurization is one potential mechanism for enhanced weakening; qualitatively similar results should hold for models with other types of enhanced dynamic weakening. We follow the thermal pressurization formulation of Noda and Lapusta (2010) (Supporting Information S1). We approximate the effects of off-fault yielding by employing a limit on the slip velocity max 15 E V  m/s, which is motivated by detailed dynamic rupture simulations with offfault yielding (Andrews, 2004) and discussed in detail in Lambert et al. (2021).
For the purpose of comparing local frictional behavior with the average prestress for dynamic ruptures of varying sizes, we focus this study on simulated ruptures that arrest within the VW region, where the friction properties are uniform with a quasi-static reference friction of 0.6, consistent with many materials exhibiting VW behavior in laboratory experiments (Ikari et al., 2011). We examine the evolution of the apparent friction coefficient, or the ratio of the current shear stress E  to the interseismic drained effective normal stress where int E p is the interseismic drained value of the pore pressure. The term "drained" refers to the effective stress with ambient pore pressure unaffected by slip processes such as dilatancy, compaction, or thermal pressurization. The values used for the direct and evolution effect parameters, E a and E b respectively, within the VW region are consistent with typical laboratory values (Tables 1-3; Blanpied et al., 1991Blanpied et al., , 1995Marone, 1998). Within the VS regions. We use higher values of the direct effect E a in order to more efficiently stop VW-region-spanning ruptures, which helps to maintain a reasonable domain size for computational expense. The properties of the VS regions should not appreciably alter the conclusions of this work, as we focus on the average stress measures for "partial" ruptures that are arrested predominantly in the VW region. The properties of larger ruptures that span the entire VW region and that arrest in the fault regions with VS properties are more sensitive to the VS properties, as discussed in Perry et al. (2020).
We examine fault models with varying levels of ambient fluid overpressure in terms of the effective normal stress, as well as varying degrees of efficiency in enhanced weakening due to thermal pressurization.
The parameter values we have chosen (Tables 1-3) are motivated by prior studies that have reproduced a range of seismological observations as well as low-stress, low-heat operation of mature faults (Lambert et al., 2021;Perry et al., 2020). The parameter values also facilitate our goal of examining ruptures in fault models with a range of efficiency in enhanced dynamic weakening. We refer to the weakening behavior of fault models as being mild, moderate or efficient in comparison to the underlying weakening of the standard rate-and-state friction. This can be approximately quantified as the difference between the steadystate frictional resistance at the plate rate and seismic slip rates, equal to ( , which corresponds to about 10% E times the interseismic effective normal stress for seis 1 E V  m/s. Mild weakening 6 of 29 corresponds to additional dynamic weakening comparable to that of standard rate-and-state friction (e.g., Fault model TP 1 in Table 2) and efficient corresponds to weakening by 70 100% E  of the prestress, all the way to near-zero shear resistance during seismic slip. Moderate weakening refers to the regime in between (e.g., Fault models TP 2-4 in Table 2).  (Table 2). Seismic events are illustrated by red lines plotted every 0.5 s while aseismic slip is shown by black lines plotted every 10 yr. (d-g) Evolution of local slip rate with time and slip at points representative of nucleation and typical rupture propagation behavior within a crack-like rupture (colored blue in c). Points throughout rupture propagation (e, g) are initially locked and are driven to rupture by the concentration of dynamic stresses at the rupture front, thus experiencing more rapid acceleration of slip compared to points within the nucleation region (d and f). (h-i) The difference in local slip rate history contributes to a difference in the evolution of shear stress with slip. (h) Evolution of the apparent coefficient of friction   / int ( )  p with slip in the nucleation region is consistent with the laboratory notion of quasi-statically strong, dynamically weak behavior, with the apparent friction coefficient initially close to the reference value of 0.6 and dropping to a low dynamic resistance below 0.2 with slip. (i) Evolution of the apparent friction coefficient at points throughout rupture propagation is more complicated as the scaled prestress can be much lower than the reference friction before the arrival of the dynamic stress concentration.

Evolution of Local Slip and Shear Resistance and Notions of Failure
Our simulations capture the evolution of motion and shear stress across the fault over sequences of earthquakes spanning several thousands of years ( Figure 2c). The initial distributions of shear stress and other quantities such as the slip rate are assumed to be uniform along most of the VW region of the fault at the start of our simulations, other than a small region of initially high prestress near the VW-VS boundary to nucleate the first rupture in the earthquake sequence. The distributions of shear stress and slip along the fault evolve to become highly variable throughout periods of fast earthquake-producing slip as well as slow aseismic slip and fault locking. Below we review how the rate-and-state friction framework allows the model to represent creeping, locked, and seismically slipping fault areas as well as transitions between these different styles of slip, a key feature of SEAS simulations (e.g., Cattania, 2019;T. Chen & Lapusta, 2009;Erickson et al., 2020;Lapusta & Rice, 2003;Michel et al., 2017;Rice, 1993;Wu & Chen, 2014).
During dynamic rupture, the evolution of slip rate and shear stress can be particularly complex and variable along the fault. At points where individual ruptures nucleate, the slip rate gradually accelerates toward seismic slip rates and shear stress at the beginning of rupture, ini E t , is relatively high, with the apparent friction coefficient   / int ( )  p close to the quasi-static reference friction of 0.6. As seismic slip rates are reached,   / int ( )  p drops substantially due to thermal pressurization of pore fluids in a manner qualitatively consistent with the enhanced dynamic weakening observed in high-velocity laboratory friction experiments ( Figure 2h). The evolution of slip rate and shear stress outside of the nucleation region is even more complicated: The shear stress at ini E t , prior to the arrival of the rupture front, can be much lower than the shear stress levels where the rupture nucleates, then increases to a higher peak shear stress (set by the balance between the dynamic stress concentration at the rupture front due to inertial effects and peak shear resistance due to interseismic fault healing and the rate-and-state direct effect), and then decreases due to weakening with seismic slip (Figure 2h vs. 2i). Consistently, the slip rate rapidly increases to seismic values at the beginning of slip and then decreases, as in a  typical Yoffe-like behavior for dynamic ruptures (Figure 2g; e.g., Tinti et al., 2005). Thus, even with the uniform normal stress and uniform parameters of the assumed friction and pore pressure equations within the seismogenic VW region, the prestress conditions throughout the rupture area can be highly variable and, in part, substantially different between regions of rupture nucleation and rupture propagation (Barbot, 2019;Cattania, 2019;Galis et al., 2017;Lapusta & Rice, 2003;Rice, 1993;Schmitt et al., 2015;Wu & Chen, 2014).
Note that the peak shear stress during dynamic rupture of fault locations outside the nucleation zone can correspond to much higher apparent friction coefficient (e.g., 0.95 in Figure 2i) than the reference friction coefficient ( * 0.6 E f  in this study). This is due to both the direct effect at the rupture tip and the high, interseismically "healed" value of the state variable E  , as discussed in Lambert and Lapusta (2020) and Equation S3. As follows from the first line of Equation S3, the difference between the peak friction coefficient and * E f due to the direct effect of a V V ln( ) * peak / would be 0.14-0.16 for peak 1 E V  to 10 m/s and other parameters of our model, with the rest due to the much larger value of the "healed" state variable than that for sliding at the reference sliding rate.
The local evolution of shear stress throughout the VW seismogenic zone differs among points based on the long-term history of motion, including both local slip as well as slip across the entire fault. For example, a point at the center of the VW region (z = 0 km) of one of our simulations (fault model TP 3 in Table 2, as shown in Figure 2c) experiences substantial slip only during the largest earthquake ruptures that span the entire VW domain, resulting in a relatively simple and quasi-repetitive pattern of stress accumulation and weakening over sequences of earthquakes (Figures 3a and 3c). In contrast, another point in the VW region closer to the VS boundary (z = −9.6 km) experiences different amounts of slip during dynamic ruptures of varying size, resulting in a more complicated evolution of shear stress with accumulating slip (Figures 3b and 3d).
In between individual earthquakes, the VS regions of the fault creep (i.e., slowly slip) with the slip rate close to the prescribed tectonic plate rate, due to that rate being imposed on the fault areas nearby, with occasional quasi-static accelerations due to post-seismic slip ( Figure 4, left column). The creep penetrates into the VW regions nearby, creating fault areas prone to earthquake nucleation (Jiang & Lapusta, 2016;Lapusta & Rice, 2003;Michel et al., 2017) (Figure 4, right column). These points of the VW region close to the VS region (within one or so nucleation length) are reloaded due to creep and post-seismic slip from previous rupture within the VS regions. The loading rate at these points near the VS-VW boundary varies over time depending on the rate of motion in the VS region, which in turn depends on the previous history of co-seismic slip during dynamic ruptures in the VW region.
The slip rate and apparent friction at points close to the VW-VS boundary are typically brought to near steady conditions around the loading plate rate, however both exhibit small oscillations as these points continue to be loaded by creep in the VS region, resulting in further acceleration, slip and weakening, and thus the transmission of stress further into the VW region until a sufficiently large area is loaded to sustain rupture nucleation and acceleration to seismic slip rates (Figures 4e-4g). This oscillatory behavior is consistent with predictions from the stability analysis of a single degree-of-freedom spring-slider undergoing frictional slip, where the amplitude of the oscillations is expected to grow as the spring stiffness decreases below a critical stiffness value (Gu et al., 1984). The effective stiffness of the slipping fault zone in a continuum model is inversely proportional to the slipping zone size (Rice & Ruina, 1983), decreasing with the  increasing slipping region. Note that this rate-and-state nucleation process has been used to explain the period-dependent response of microseismicity to periodic stress perturbations in Nepal, where seismicity shows significant variations in response to annual monsoon-induced stress variations but not to semidiurnal tidal stresses of the same magnitude (Ader et al., 2014).
In contrast, much of the VW region further away from the VS regions is essentially locked, which is expressed in the rate-and-state formulation as sliding at very low, but still non-zero, slip rates that are many orders of magnitude smaller than the loading rate (Figures 5a and 5b). This differential motion between the VS and VW regions loads points in the VW region (Figures 5c and 5d), gradually increasing shear stress there (e.g., between 700 and 800 yr in Figure 5c). Note that the interseismic stressing rate is higher at locations closer to the creeping regions than further away from it ( Figure 5c vs. 5d vs. 4f), as one would expect (Cattania, 2019;Herrendörfer et al., 2015;Lapusta & Rice, 2003;Michel et al., 2017;Rice, 1993;Wu & Chen, 2014). At the same time, the essentially locked points within the VW region experience time-dependent healing of the local shear resistance encapsulated in the increase of the state variable E  (Figures 5e and 5f). One of the manifestations of this healing is that larger interseismic increases in the state variable generally lead to higher peak shear stress during dynamic rupture propagation . Despite the increase in the state variable, its value is far below the steady-state one for the very low interseismic slip rates, consistent with continuing healing prior to dynamic rupture (Figures 5g and 5h). Depending on whether the local shear stressing rate (which increases the shear stress E  on the left of Equation 1) is larger or smaller than the rate of healing (expressed by the last, E  term on the right hand side of Equation 2), the local slip rate (that enters the second term of Equation 2) increases (as between 700 and 800 yr in Figure 5a) or decreases, that is, the fault is accelerating toward failure or becomes even more locked. However, most of  Table 2). The stars denote instances in the earthquake sequence in Figure 2c, with pink stars marking the initiation of the three large modelspanning ruptures, the blue and red stars denoting the beginning and end of the moderate-sized rupture illustrated by blue contours, respectively. The yellow stars denote small to moderate-sized ruptures occurring along the VW-VS boundary at 12 E z   km. (a, c) The point in the center of the VW region (z = 0 km) ruptures and experiences substantial slip only in large ruptures. The point exhibits an increase in shear stress over time due to the stress transfer from smaller ruptures that do not penetrate into the center of the VW region (such as the rupture colored blue in Figure 2c). (b and d) Points closer to the boundary between the VW and VS regions can rupture during both smaller and large ruptures depending on the prestress conditions when ruptures arrive, resulting in a more complicated evolution of shear stress with accumulating slip. For both points in the VW region, the shear stress is brought to the peak stress and failure during ruptures by the dynamic stresses at the rupture front.  Figure 2c, with pink stars marking the initiation of the first two large model-spanning ruptures, the blue star denoting the beginning of the moderate-sized rupture illustrated by blue contours and the yellow stars denoting smaller ruptures. (a) Points within the VS region typically slip near the loading plate rate but can experience transient accelerated slip during and following ruptures occurring within the VW region. (b-d) The apparent friction coefficient and state variable in the VS region is typically near steady state, except during accelerated slip. (e and f) Slow slip penetrates into the VW region, driving points near the VW-VS boundary close to the loading slip rate, with the apparent friction coefficient being close to the corresponding steady-state value ss pl ( ) E f V . The slip rate and apparent friction exhibit small oscillations as the points near the VW-VS boundary continue to be loaded by slow slip in the VS region, accelerate, and weaken, thus transmitting stress further into the VW region until a sufficiently large region is loaded to sustain rupture nucleation and acceleration to seismic slip rates. The loading rate of the VW region also depends on the amount of accelerated slip in the VS due to previous ruptures (e.g., a and e around 650 vs. 750 yr). (g and h) Following dynamic rupture, the state variable heals close to the steady-state value around the prescribed loading rate LAMBERT ET AL.
10.1029/2021JB021886 11 of 29 more involved formulations of the state evolution law exist, for example, the slip law and various composite laws (Kato & Tullis, 2001;Perrin et al., 1995;Ruina, 1983), some of which provide a better (although still imperfect) match to laboratory experiments (Bhattacharya et al., 2015). Note that aging and slip laws resulted in qualitatively similar results for simulations of repeating earthquake sequences (T. Chen & Lapusta, 2019). Exploring the behavior of our models with different state evolution laws and incorporating more , corresponding to the current local slip rate E V , is much smaller than 1 during the interseismic periods, indicating the continued healing of shear resistance prior to rupture. As the slip rate rapidly accelerates during dynamic rupture, the state variable temporarily exceeds the new much lower steady-state values corresponding to the dynamic slip rate , then evolves to this lower steady-state value, and then falls to values below steady-state during the interseismic periods, indicating fault healing.
12 of 29 realistic healing into shear resistance formulations and numerical modeling is an important goal for future work. This can be done by modifying the evolution of the state variable E  or adding other state variables that would encode healing. Yet, qualitatively, additional healing mechanisms would have similar effects on the simulations as the current rate-and-state healing, in that the healing would modify the peak shear resistance and the subsequent evolution of the resistance based on the interseismic fault state, potentially further amplifying differences in shear resistance evolution for different points along the fault (e.g., nucleation points vs. locked points) that our simulations already highlight.
The presence of time-dependent healing as well as persistent, potentially unperceivable, slow (quasi-static) motion and its acceleration under variable levels of shear stress illustrate how the concepts of failure, and hence strength, are not easily defined for frictional sliding. For realistic frictional interfaces, the precise value of a static friction coefficient is ill-defined, since no interface loaded in shear is perfectly static; rather creep processes occur at slow, unperceivable slip rates at any level of shear loading (Bhattacharya et al., 2017;Dieterich & Kilgore, 1994) and/or over parts of the contacting interfaces (Ben-David et al., 2010;Rubinstein et al., 2004Rubinstein et al., , 2006. Hence the transition from locked interfaces to detectable slip is always a gradual process (although it may be occurring faster than the time scales of interest/observation in many applications). This reality is reflected in lab-derived fault constitutive relations such as rate-and-state friction.
Since failure typically refers to the presence of irreversible or inelastic deformation, frictional interfaces may be considered failing under any style or rate of motion, be it during slow steady sliding, transient slow slip, or dynamic rupture. Therefore, any meaningful notion of strength first requires definition of the failure of interest, for example, reaching seismic slip rates of the order of 1 m/s. Without such explicit definition, failure is then implicitly defined as transition from locked to slipping and corresponds to sliding with a detectable velocity; for laboratory experiments or observational studies, this would imply that whether the interface is locked or slipping depends on the instrumental precision for detectable motion.
In this study, we would like to compare the shear stress values required for aseismic slip nucleation and for dynamic rupture propagation. During spontaneous aseismic slip nucleation, the slip rates evolve from very low to seismic, passing in the process through the slip rate equal to the tectonic loading rate pl E V . In the standard rate-and-state friction, at each fixed sliding rate E V , the friction coefficient eventually evolves to a steady-state value ( ) ; for very small slip rates, the regularized formulation of Equation S5 needs to be considered). Under slow loading, aseismic earthquake nucleation on a finite fault is typically a gradual process, with many points within the nucleation zone being close to the steady state (Figure 4; Kaneko & Lapusta, 2008;Rubin & Ampuero, 2005). While the steady-state values of friction depend on the sliding rate, the dependence is relatively minor at the low, quasi-static slip rates between the plate rate of approximately 9 10 E  m/s and sub-seismic slip rates of 3 10 E   m/s ( Figure 1) which are relevant for fault creep and earthquake nucleation, and for which the standard rate-and-state formulation is (approximately) valid. The product of this collection of steady-state quasi-static friction coefficients and the interseismic drained effective stress gives the shear resistance of faults at sustained slow sliding rates, which we call the steadystate quasi-static fault shear resistance (referred to in short as local SSQS shear resistance). As the representative value of such local SSQS shear resistance, we choose the shear resistance of the fault steadily creeping at the prescribed long-term tectonic plate rate pl E V (which the fault would have long-term if it were slipping stably), with the interseismic drained value of the pore pressure int E p : In the following section, we compare this representative value of local SSQS shear resistance to the spatial distribution of shear stress prior to dynamic ruptures in our simulations. Note that the local SSQS shear resistance is similar to what is typically viewed as "frictional fault strength" in the sense of Byerlee (1978), that is, this is the resistance that needs to be met for noticeable quasi-static slip with the loading rate or another reference rate.
LAMBERT ET AL.

Larger Ruptures Associated With Lower Shear Prestress Over the Rupture Scale but Higher Prestress Over Smaller Scales Near Nucleation
The interseismic periods in between individual earthquake ruptures in our simulations vary from months to decades, depending on the size of the rupture and the stress state resulting from the history of prior slip along the fault. Our earthquake sequence simulations produce a wide variety of rupture sizes due to heterogeneous prestress conditions along the fault that spontaneously arise in our models.
Let us consider the evolution of slip and shear stress in representative simulated spontaneous ruptures of increasing sizes within the same simulation ( Figure 6). Over sequences of rupture events, the shear stress conditions prior to and after individual dynamic ruptures become spatially heterogeneous. This stress heterogeneity is due in part to the history of spatially variable slip and local static stress drop produced in previous ruptures, as well as stress relaxation and redistribution due to aseismic slip. In addition, while our simulated fault models are loaded by a constant long-term loading rate of pl E V , the effective loading conditions along the fault interface vary in space and time due to differences in slip rate along the fault. Ruptures nucleate preferentially in regions with the highest shear prestress, which in our models occur near the creeping regions as discussed in Section 3 ( Figure 6; Barbot, 2019;Lapusta et al., 2000;Michel et al., 2017). The ruptures then propagate into the less stressed areas of the fault. Put another way, the average prestress over the nucleation region is higher than the average prestress over the entire ruptured region (Figure 7a vs. 7b), as we quantify in the following.
We compute the average shear prestress right before a dynamic rupture event over the entire future rupture area (which we do as post-processing of data in our simulation). We also compute the average shear prestress over the slow-slip nucleation zone, which we call the nucleation stress. We compare these average shear stress measures with the local steady-state quasi-static (SSQS) fault shear resistance pl V ss E  , which is related to the local fault constitutive properties during slow slip and given by Equation 4. Averaging of spatially variable stress fields can be done in several different ways (Noda & Lapusta, 2012;Noda et al., 2013). The simplest definition of the average shear prestress over the rupture region E  is the spatially averaged prestress ini A E  acting in the overall slip direction at the beginning of the rupture ini E t , given by: We can similarly define the spatially averaged nucleation stress nucl A E  within the nucleation region. We define the nucleation region to be the fault segment between the expanding stress fronts at the initiation of dynamic rupture; the size of the nucleation regions in our simulations is comparable to the theoretical nucleation size estimate * RA E h of Rubin and Ampuero (2005) (Equation S6; Figure S1).
Not surprisingly and consistent with prior studies, we find that the spatially averaged nucleation stress . As a consequence, it does not significantly depend on the ultimate rupture size or slip. Since the nucleation stress here is computed at the beginning of dynamic rupture, it is then the shear stress within the nucleation zone at the end of the nucleation, when parts of the zone slip with near-dynamic slip rates approaching 2 10 E  m/s. That is why the nucleation stress is systematically slightly lower than the local SSQS shear resistance defined as the steady-state shear resistance to slip with the (lower) plate rate. The difference between the nucleation stress and local SSQS shear resistance could be more substantial if dynamic weakening were efficient enough to affect some portion of the earthquake nucleation region (Segall & Rice, 2006).
In contrast, the spatially averaged prestress over the entire ruptured area ini A E  tends to decrease with the rupture size and increasingly deviate from the local SSQS shear resistance and nucleation stress for increasingly efficient dynamic weakening (Figures 6 and 7b). Such behavior is also true for another average prestress measure, the energy-based average prestress   being slightly larger ( Figure S2).
For the range of simulated rupture sizes considered, we find that the decrease in average shear prestress with increasing rupture size in our 2-D models with 1-D faults is consistent with a power-law relationship  Table 2). Slip contours are plotted every 0.25 s. The purple and gray shading illustrates the extent of the nucleation and ruptured regions, respectively, over which the prestress is averaged. While the ruptures nucleate in regions with stress levels near the local steady-state quasi-static shear resistance (dashed orange line), larger ruptures propagate over lower prestressed areas, resulting in lower average prestress and lower average coefficients of . The shear stress distribution for a typical moment during rupture propagation is shown in black, demonstrating the stress concentration at the rupture front that brings the fault stress to values comparable to the SSQS shear resistance. The peak stress is even higher since the fault is initially dynamically stronger due to the rate-and-state direct effect. (d and e) Significant differences in local evolution of slip and stress at the same fault location ( 9.6 E z  km) for different ruptures that depend on the prestress conditions due to previous slip events and the dynamic stress interactions during the individual ruptures.
LAMBERT ET AL.
where 1 E c and 2 E c are fit parameters (Table 4). We find that 2 E c is comparable to unity for all our simulations (Table 4), while 1 E c depends on the efficiency of weakening in each simulation, increasing by about a factor of 5 from 0.07 in our fault models with no enhanced weakening (RS1 and 2) to around 0.3 in models with moderate enhanced weakening (TP3 and 4). Thus, for a decade increase in rupture length rupt E  with respect to nucleation size * RA E h , the decrease in the spatially averaged prestress with increasing rupture size differs by a factor of 5 between our models with decrease with increasing rupture size; the effect is more pronounced with increasing efficiency of weakening. The three ruptures shown in Figure 6 are denoted by red stars. (c) The decrease in average prestress with rupture size for our simulated ruptures is generally consistent with a power-law relationship with the ratio of the rupture length and nucleation size (Equation 7;   Figure 7c) LAMBERT ET AL.

10.1029/2021JB021886
16 of 29 no to moderately efficient enhanced dynamic weakening. We find comparable results for the energy-based prestress ini E E  (Equations 9 and 10; Figure S3 and Table S1).
The finding that larger ruptures are associated with smaller average shear prestress over the ruptured area may appear counterintuitive. Why do smaller ruptures not become larger if they are more favorably prestressed? To understand this behavior, let us consider the prestress averaged over several fixed scales around the nucleation region for ruptures of different sizes. We locate the VW-VS boundary next to which each of our simulated ruptures nucleate and average the prestress along the VW region over fixed distances (1, 2, 4, 8, 12, and 16 km) from the corresponding VW-VS boundary (Figure 8; shown for fault model TP4 from Table 2). While the spatially-averaged prestress over the entire rupture length decreases with increasing rupture size, we see that the prestress spatially-averaged over smaller fixed scales is generally higher for larger ruptures than for smaller ruptures (Figure 8 warmer vs. cooler colored triangles). For smaller ruptures, the average shear stress over scales just larger than their total rupture length is lower than the average prestress of larger ruptures with comparable length to the fixed averaging scales (Figure 8, triangles below the circles). This confirms that the smaller ruptures arrest because the prestress conditions ahead of the rupture are too low to sustain further rupture propagation. For larger ruptures, the average prestress levels at scales smaller than their total rupture length are generally higher or comparable to the average prestress over smaller ruptures with the length comparable to the fixed averaging scales (Figure 8, triangles above the circles). This finding suggests that larger ruptures have higher, more favorable average prestress conditions at smaller scales compared to smaller ruptures, which facilitates continued rupture propagation. Hence we find that the shear prestress prior to our simulated ruptures of varying sizes self-organizes into a spatial distribution of scale-dependent average shear stress that governs the rupture occurrence.  Figure 7, the spatially-averaged prestress over the total rupture area ini A E  (circles) decreases considerably with rupture size in fault model TP4 from Table 2 with moderate enhanced dynamic weakening. However, larger ruptures have generally higher average shear stresses over smaller fixed scales around the nucleation region compared to smaller ruptures (red vs. blue triangles). The spatially-averaged shear stress over 1 km from the VW-VS boundary near the nucleation region of ruptures (triangles on the far-left) is relatively high (comparable to the local SSQS resistance) for both small and large ruptures, indicating that ruptures nucleate in regions of relatively high prestress compared to the average prestress over the entire rupture area (circles). For smaller ruptures, the average prestress at the fixed scales just larger than their total rupture length is lower than the average prestress of ruptures with comparable length to the fixed scale, suggesting that the prestress levels were too low to sustain further rupture propagation.
LAMBERT ET AL.

10.1029/2021JB021886
17 of 29 As discussed in Perry et al. (2020), the decrease in average prestress with increasing rupture size can result in nearly magnitude-invariant static stress drops. Note that, as the static stress drop for points throughout typical rupture propagation for our simulated ruptures with efficient enhanced dynamic weakening is lower than the static stress drop within the nucleation region (Figure 2h vs. 2i), the average static stress drop can decrease as the rupture length increases beyond the nucleation length ( Figure S3; Ampuero et al., 2006;Perry et al., 2020). The effect is more noticeable for the energy-based stress drop which weights the local static stress drop by the local slip, but becomes less pronounced for our simulated ruptures that are considerably larger than the nucleation length. We refer to these ruptures as having nearly magnitude-invariant static stress drops, since the variation of the stress drop with rupture size is relatively mild and would likely fall within the scatter and uncertainty of seismological observations (Allmann & Shearer, 2009;Ye et al., 2016b).

Role of Dynamic Stress Transfers and Motion-Dependent Local Shear Resistance
Such scale-and motion-dependent average fault shear prestress before ruptures results from two related and interacting factors. First, as dynamic rupture propagates, some of the released energy is carried by waves along the fault, creating a substantial stress concentration near the rupture tip that is a well-known feature of dynamic rupture (e.g., Freund, 1990). The stress concentration enables rupture propagation over regions where the prestress is lower than the local SSQS shear resistance, drawing the local shear stress up to the peak stress before the subsequent stress drop due to local weakening (black lines in Figure 6). The dynamic stress concentration increases with the rupture dimension and/or slip and thus allows larger ruptures to continue propagating over regions with lower, and hence less favorable, prestress conditions ( Figure 6). This is illustrated in this work for largely crack-like ruptures that occur in the presented models with mild to moderate enhanced dynamic weakening (Lambert et al., 2021), but similar conclusions would be reached for pulse-like ruptures provided that they satisfy the observational constraint of magnitude-independent stress drops, which implies that ruptures with larger magnitudes would have larger average slip and hence larger stress concentrations. Note that a pulse-like rupture with the same or similar spatial distribution of the slip rate (and hence the same local slip) propagating along the fault would result in a similar stress concentration at the rupture tip regardless of the rupture length; however, in that scenario, pulses with larger rupture propagation lengths would have systematically lower static stress drops, as the stress drops would be proportional to the (uniform) pulse slip divided by ever increasing propagation lengths.
Second, the evolving local shear resistance substantially depends on both the prior history of slip events on the fault through fault prestress and on the motion during the current rupture event through dynamic stress transfers that add substantial time-dependent loading. This pronounced dependence is due to strong coupling between the evolving motion and shear resistance, which is coupled to the resulting shear heating in the case of thermally-activated weakening mechanisms like thermal pressurization. As a result, the evolution of local slip rate and local shear resistance (a) significantly differs at different fault locations of each rupture (despite uniform constitutive properties) and (b) significantly differs at the same fault location for different ruptures (Figures 2d-2i and 6d-6e).
These two factors create a substantial positive feedback, in which larger ruptures with more slip generate larger stress concentrations, leading to faster and larger slip, which dynamically causes more fault weakening, which in turn promotes more/faster slip, more energy release, larger stress concentrations, and increasing rupture sizes.
The result that larger ruptures are associated with lower average prestress indicates the need for increasingly less favorable stress conditions to arrest growing ruptures. For a given rupture size, if the prestress ahead of the rupture is favorable, then the rupture would continue to grow until it experiences sufficiently unfavorable prestress conditions, thus lowering the overall average prestress. Alternatively, the rupture may be forcibly arrested by other means such as strong geometric or rheological barriers. For example, ruptures propagating over higher prestress conditions within the VW region can be arrested by fault regions with VS properties; in those cases, the overall average prestress conditions would depend on the properties of the VS regions (Perry et al., 2020). Detailed study of the implications of fault geometry and heterogeneity for rupture arrest and the average stress conditions prior to rupture is an important topic for future work.
LAMBERT ET AL.

Comparison of Finite-Fault Modeling to Single-Degree-of-Freedom Representations
As captured in field observations of natural earthquakes and reflected in our simulations, sufficiently large earthquake ruptures nucleate on a subsection of the fault and then propagate through other sections of the fault. Capturing such space-dependent behavior is typically called "finite-fault" modeling, in contrast to the point source that considers a spatially averaged representation of an event, as if it occurs at one "point". A typical numerical model of a point source is the single-degree-of-freedom system (SDOF) of a slider with friction pulled by a spring (e.g., Dieterich, 1979;Rice & Ruina, 1983;Ruina, 1983). Small-scale laboratory experiments often measure properties averaged over a sample and are typically modeled as a SDOF spring-slider systems.
The significant role of spatially varying prestress conditions and dynamic stress transfers during rupture propagation in determining the rupture behavior implies that capturing the finite-fault nature of the process is essential for determining the stress evolution characteristic of dynamic rupture. For example, several laboratory studies applied variable slip rates histories inferred from natural earthquakes to rock samples, measured the resulting shear resistance, and then related laboratory stress measurements to seismological source properties such as breakdown energy and stress drops (e.g., Fukuyama & Mizoguchi, 2010;Nielsen et al., 2016;Sone & Shimamoto, 2009). Such experiments have provided invaluable data about the local shear resistance of faults, specifically enhanced dynamic weakening, that have informed theoretical and numerical modeling of finite faults (e.g., Dunham et al., 2011;Gabriel et al., 2012;Lambert et al., 2021;Noda & Lapusta, 2010;Noda et al., 2009;Perry et al., 2020;Rice, 2006;Zheng & Rice, 1998), including the current study. However, the interpretation of such experiments needs to take into account their SDOF nature. For example, to improve alignment etc, the experiments often impose pre-sliding at slow slip rates (of the order of micron/s) prior to imitating seismic motion. That procedure results in the shear prestress before seismic slip comparable to the local SSQS shear resistance (Equation 4) and near steady-state values of the state variable, as appropriate for a location within a nucleation zone. In contrast, our simulations show that most points on a fault through which the rupture propagates have much lower shear prestress and much larger values of the state variable corresponding to well-healed fault (Figures 6 and 9b). Furthermore, the experiments often apply smoothened slip-rate histories obtained from finite-fault inversions, while the stress concentration at the tip of dynamic rupture makes the slip rate variation much more dramatic.
To illustrate the differences for the shear resistance evolution obtained with such experimental procedures versus the one from our simulated finite-fault models, let us compare the local fault behavior during one of our dynamic ruptures with a SDOF calculation. In the SDOF calculation, we use the same fault properties (Equations 3, S4, S7, and S8) and same parameter values as in the finite-fault VW regions but apply quasi-static presliding and modified, smoothened slip rates motivated by the laboratory procedures of Fukuyama and Mizoguchi (2010) (further details in Supporting Information S1). We conduct the comparison for two fault locations, one in the nucleation region and one within dynamic rupture propagation region (Figure 9). These SDOF calculations are successful at reproducing the presence of the enhanced dynamic weakening with slip as occurs during dynamic ruptures and generally capture the more moderate slip evolution and behavior of points within the nucleation region of our simulated dynamic ruptures. At the same time, the overall shear stress evolution during typical propagation of the dynamic rupture substantially differs from that of the SDOF calculation, with notable features including the low initial stress (which depends on prior slip history) relative to the SSQS shear resistance, the much more dramatic increase in shear stress associated with the dynamic rupture front (which arises due to the more healed fault coupled with the dynamic stress concentration), and the shear stress evolution at the end of slip (which depends on the final slip distribution over the entire finite fault) (Figure 9).

Implications for Earthquake Statistics
A notable feature of the scale dependence of average prestress before dynamic rupture is that, as an earthquake grows larger, the prestress needed for further propagation decreases (Figure 7b). In addition, the higher the weakening rate, the easier it should be for a rupture to have favorable prestress conditions to continue growing, rather than arresting as a smaller earthquake. Hence one could hypothesize that the LAMBERT ET AL.

10.1029/2021JB021886
19 of 29 more efficient the enhanced dynamic weakening, the smaller the complexity of the resulting earthquake sequences, with increasing representation of larger events at the expense of smaller events.
This is exactly what our modeling shows (Figure 10). The fault models with increasingly more efficient weakening produce earthquake sequences with increasingly fewer small events and decreasing b-values of the cumulative size distribution (Figure 10; details for calculating seismic moment in Text S1). Fault models with even more efficient dynamic weakening than considered in this study, such as those that produce sharp self-healing pulses, result in relatively simple earthquake sequences consisting of only large events (Lambert et al., 2021). The fault models governed by relatively mild to more moderate weakening as considered in this work develop a wider range of earthquake sizes, due to a feedback loop of more likely rupture arrest due to milder weakening creating stress heterogeneity that in turn makes rupture arrest more likely. This result is consistent with those of previous quasi-dynamic earthquake sequence simulations demonstrating Figure 9. Comparison of the results of our dynamic modeling with a single-degree of freedom (SDOF) model that might more closely represent what would be obtained in laboratory experiments given the same constitutive properties and typical lab procedures. (a) Comparison of the local slip rate during nucleation ( 11.5 E z  km, yellow) and typical propagation ( 9.6 E z  km, black) of the simulated dynamic rupture of Figure 6b with the slip rate evolution that could be imposed in lab experiments represented by two regularized Yoffe functions (Tinti et al., 2005) with peak slip rate of 2 m/s and comparable slip to the point at 9.6 E z  km. The imposed regularized Yoffe functions are generally comparable to the evolution of slip within the nucleation region ( 11.5 E z  km), however they do not capture the rapid acceleration of slip associated with the arrival of the rupture front at points of typical propagation, as observed at 9.6 E z  km. (b) Comparison of the state variable evolution from our finite-fault dynamic simulation and the SDOF simulation of the lab experiment. The simulated lab experiment starts with the steady-state conditions for 0.1 mm/s based on the experiments of Fukuyama and Mizoguchi (2010), which results in a much lower initial state value compared to the point 9.6 E z  km in our simulations which, prior to dynamic rupture, had negligible motion over a 20-yr interseismic period. (c and d) Evolution of the local apparent coefficient of friction with time and slip for the point in our simulated finite-fault dynamic rupture and simulated SDOF lab experiments. The dynamic weakening is generally comparable between the points in the finite rupture and the simulated SDOF experiments, however the evolution of shear stress substantially differ with regards to the much lower prestress at 9.6 E z  km before the finite dynamic rupture and the abrupt increase and then decrease in stress due to the arrival of the dynamic rupture front and the associated rapid weakening.
LAMBERT ET AL.

10.1029/2021JB021886
20 of 29 complex earthquake sequences with higher b-values around 0.75 on faults governed by standard rate-andstate friction only and milder quasi-dynamic stress transfer, as well as relatively large fault lengths compared to the nucleation size,  VW RA /h *  100 (Cattania, 2019;Wu & Chen, 2014). Our results show b-values around 0.5 for fully dynamic simulations with instability ratios of  VW RA /h *  44 and without enhanced dynamic weakening, which further decrease to 0.25 or so for the most efficient weakening considered in this study. Note that simulations with higher instability ratios and even more efficient dynamic weakening than considered in this work, such as that consistent with the propagation of sharp self-healing pulses, can result in only large earthquakes (Lambert et al., 2021), consistent with our observation of decreasing variability of rupture sizes for fault models with more efficient dynamic weakening, even in cases with relatively high instability ratios. While the frequency-magnitude distribution of seismicity over relatively large regions, such as Northern or Southern California, is generally well-described by Gutenberg-Richter scaling with typical b-values near unity (Field et al., 2013), whether such scaling applies to individual fault segments and/or their immediate surroundings is a topic of active research (Field et al., 2017;Ishibe & Shimazaki, 2012;Kagan et al., 2012;Page & Felzer, 2015;Page & van der Elst, 2018;Wesnousky, 1994). Estimates of b-values associated with individual fault segments can exhibit considerable variability (e.g., between 0.5 and 1.5 along faults in California; Schorlemmer & Wiemer, 2005;Tormann et al., 2014), and are sensitive to a number of factors, including the magnitude of completeness of the relevant earthquake catalog and the choice of observation region and time window (Ishibe & Shimazaki, 2012;Page & Felzer, 2015;Page & van der Elst, 2018;Tormann et al., 2014). A number of studies suggest that the rate of large earthquakes on major faults, such as the San Andreas Fault, is elevated above what would be expected given typical Gutenberg-Richter scaling from smaller magnitude events (Field et al., 2017;Schwartz & Coppersmith, 1984). In particular, some mature fault segments that have historically hosted large earthquakes, such as the Cholame and Carrizo segments of the San Andreas Fault, exhibit substantial deviations from typical Gutenberg-Richter scaling, being nearly absent of small earthquakes (Bouchon & Karabulut, 2008;Hauksson et al., 2012;Jiang & Lapusta, 2016;Michailos et al., 2019;Sieh, 1978;Wesnousky, 1994). Our findings suggest that the paucity of microseismicity on such mature fault segments may indicate that they undergo substantial dynamic weakening during earthquakes ruptures.

Discussion
Our simulations reemphasize that the average shear prestress required for rupture propagation can be considerably lower than the average shear stress required for the rupture nucleation (Barbot, 2019;Cattania, 2019;Galis et al., 2017;Lapusta & Rice, 2003;Rice, 1993;Schmitt et al., 2015;Wu & Chen, 2014). This is because the quasi-static nucleation process is governed by relatively small stress changes and hence requires favorable prestress conditions-close to the local steady-state quasi-static shear resistance-to proceed. In contrast, during dynamic rupture, the rupture front is driven by larger wave-mediated dynamic stress concentrations, which are more substantial for larger ruptures and facilitate rupture propagation over less favorably stressed regions, resulting in the spatially-averaged prestress over the ruptured area being much lower than the average local SSQS shear resistance. More efficient weakening facilitates larger dynamic stress changes at the rupture front, allowing propagation over even less favorable prestress conditions. Our results highlight the significance of heterogeneity in prestress, or shear resistance, for the nucleation and ultimate arrest of finite ruptures, even in fault models that have otherwise uniform material and confining properties.
The decrease in averaged prestress with rupture length can be interpreted as a decrease in the average quasi-static friction coefficient  with rupture size (Figure 7). The average quasi-static friction coefficients for ruptures on the scale of the nucleation size are consistent with the prescribed quasi-static reference friction coefficient near typical Byerlee values. However, as we average the prestress over larger rupture lengths, the average quasi-static friction coefficient can considerably decrease depending on the efficiency in weakening.
The presence of enhanced dynamic weakening draws the average shear stress along larger regions of the fault below the local SSQS consistent with earthquake nucleation, resulting in lower average shear stress conditions in terms of both the average prestress for larger ruptures and the average dynamic resistance associated with shear heating during ruptures ( Figure 11). The models presented in this study with mild-to-moderate enhanced weakening include considerable persistent fluid overpressurization to maintain low-heat, low-stress conditions with average dynamic shear resistance during seismic slip rates below 10 MPa; however the degree of fluid overpressure required to maintain low-heat conditions is less than that with comparable rate-and-state properties but no enhanced weakening. The presence of some enhanced dynamic weakening is also needed for persistently weak fault models due to chronic fluid overpressure in order to ensure that static stress drops are not too small, as they would otherwise be with low effective stress and small changes in the friction coefficient due to standard rate-and-state laws (Figures 11 and S3;Lambert et al., 2021). Fault models with more efficient dynamic weakening have been shown to be able to reproduce low-stress operation and reasonable static stress drops with quasi-static friction coefficients 22 of 29 around Byerlee values and higher effective normal stress (e.g., 100 E  MPa; Dunham et al., 2011;Lambert et al., 2021;Noda et al., 2009). Earthquake sequence simulations of such fault models typically consist of only large ruptures (Lambert et al., 2021), consistent with the notion that large fault areas governed by efficient weakening maintain substantially lower average shear stresses than that required for nucleation. These findings further strengthen the conclusion of prior studies that enhanced dynamic weakening can help explain the discrepancy between laboratory values of (quasi-static) friction coefficients around 0.6 and geophysical inferences of low effective coefficients of friction ( E  0.2), along with mild average static stress drops of 1-10 MPa, over fault areas that host large earthquakes (e.g., Allmann & Shearer, 2009;Dunham et al., 2011;Gao & Wang, 2014;Ikari et al., 2011;Lambert et al., 2021;Marone, 1998;Noda et al., 2009;Perry et al., 2020;Suppe, 2007;Ye et al., 2016b).  (black line) over earthquakes sequences. (a, b) Standard rate-and-state friction results in modest changes in shear resistance from the average local steady-state quasi-static (SSQS) shear resistance (orange line). Ruptures on persistently strong faults produce realistic static stress drops (a); however, the fault temperature would increase by more than 3,000°C during a dynamic event for a shearzone half-width of 10 mm. (b) Persistently weak fault models due to low effective normal stress but with no enhanced weakening (RS two of The scale dependence of average prestress before ruptures can also be interpreted as a scale dependence of average fault strength, since the average prestress represents a measure of how much shear stress that fault region can hold before failing in a rupture. Given this interpretation, our simulations suggest that faults maintain lower average shear stresses, and hence appear weaker (or understressed with respect to quasi-static failure), at larger scales than at smaller scales. This interpretation is conceptually consistent with laboratory measurements of scale-dependent yield stress for rocks and a number of engineering materials, which demonstrate decreasing material strength with increasing scale (Bandis et al., 1981;Greer et al., 2005;Jaeger & Cook, 1976;Pharr et al., 2010;Thom et al., 2017;Uchic et al., 2004;Yamashita et al., 2015). Note that our larger simulated ruptures, even with more efficient weakening, still require higher average shear stresses over smaller scales in order to nucleate and grow. Thus the lower average prestress levels that allow continued failure in dynamic ruptures at larger scales only become relevant once the rupture event has already nucleated and sufficiently grown over smaller scales. Hence, while the faults appear to be weaker on larger-scales, they are clearly not truly so, in the sense that the "stronger" small scales have to fail quasi-statically before the larger scales can fail dynamically, and that small-scale quasi-static strength can be high everywhere on the fault. This consideration suggests that the critical stress conditions for rupture occurrence are governed not by a single stress quantity but by a distribution of scale-dependent stress criteria for rupture nucleation and continued propagation. An important implication of our findings is that the critical stress for earthquake occurrence may not be governed by a simple condition such as a certain level of Coulomb stress. Given our findings, in order to reason about the stress conditions critical for a rupture to occur, it is important to consider both the size of the rupture and the weakening behavior, and hence the style of motion, that may occur throughout rupture propagation.
The scale dependence of fault material strength has also been hypothesized to explain the measured scaling of roughness on natural fault surfaces (Brodsky et al., 2016). Dynamic rupture simulations on geometrically irregular faults motivated by such roughness measurements have indicated an additional contribution to fault shear resistance arising from roughness drag during rupture propagation (Fang & Dunham, 2013). Further examination of the scale dependence of average shear resistance across faults including realistic fault geometry is an important topic for future work.
A common assumption when considering the average motion along a fault in relation to a SDOF system is that the shear prestress over the entire ruptured area must be near the average local static (or quasi-static) strength, comparable to the SSQS shear resistance discussed in this study. Our results reemphasize that the assumption is not necessarily valid for finite ruptures; moreover, faults with enhanced dynamic weakening and history of large earthquake ruptures would, in fact, be expected to have low average shear stress over large enough scales. At the same time, the state of stress needs to be heterogeneous, with the average stresses over small scales (comparable to earthquake nucleation) being close to the (much higher) local SSQS shear resistance in some places. Thus, while individual measurements of low resolved shear stress onto a fault may suggest that those locations appear to not be critically stressed for quasi-static failure, those regions, and much of the fault, may be sufficiently stressed to sustain dynamic rupture propagation and hence large earthquake ruptures.
In addition, our findings suggest that inferences of stress levels on faults may differ if they are obtained over different scales or influenced by different rupture processes. For example, low-stress conditions on mature faults from observations of low heat flow may not only represent average shear stress conditions over large fault segments as a whole but also be dominated by low dynamic resistance during fast slip, whereas averages over smaller scales would be expected to reflect the heterogeneity of the underlying prestress distribution, as perhaps reflected in varying stress rotations inferred over scales of tens of kilometers (Hardebeck, 2015;Hardebeck & Hauksson, 1999. This assertion fits with the notion of "strong" intraplate faults where crustal stress measurements suggest stress levels commensurate with incipient failure on faults with Byerlee friction (Townend & Zoback, 2000). In these regions, the smaller intraplate faults require higher prestress levels across the entirety of the fault surface not only for nucleation, but also propagation.
Our modeling shows that increasingly efficient dynamic weakening leads to different earthquake statistics, with fewer small events and increasing number of large events. Another factor that can significantly affect the ability of earthquake ruptures to propagate is fault heterogeneity, including variations in the rate-andstate frictional properties, effective normal stress, and fault roughness (e.g., Ampuero et al., 2006;Cattania 24 of 29 & Segall, 2021;Heimisson, 2020;Hillers et al., 2006Hillers et al., , 2007Schaal & Lapusta, 2019). Some dynamic heterogeneity in shear stress spontaneously develops in our simulations, leading to a broad distribution of event sizes for cases with mild to moderate enhanced dynamic weakening. Our findings suggest that the effects of pre-existing types of fault heterogeneity need to be considered with respect to the size of the rupture and weakening behavior on the fault. For example, faults that experience more substantial weakening would require the presence of larger amplitudes of small-wavelength heterogeneity in shear stress or resistance to produce small events. Examining the relationship between earthquake sequence complexity and varying levels of fault heterogeneity and enhanced dynamic weakening is an important topic for future work.
The models presented in this work consider changes in the pore fluid pressure due to the thermal expansion of pore fluids and off-fault diffusion of heat and fluids, but do not account for additional pore-scale fluid effects such as poroelastic stress changes or dilatancy. Such effects should not substantially alter the conclusions of this work provided that the enhanced dynamic weakening-either due to thermal pressurization of pore fluids as considered in this work, or due to other mechanisms-survives their addition. Dilatancy of the pore space during the initiation of unstable slip can stabilize slip by decreasing the pore fluid pressure, and hence increasing shear resistance. However, thermal pressurization is expected to overwhelm dilatancy at seismic slip rates during dynamic rupture (e.g., Segall et al., 2010), so the dilatancy should not substantially alter the results for our models, unless the dilatancy is much greater than accounted for by the commonly used formulation for dilatancy (Marone et al., 1990;Segall & Rice, 1995). The effect of poroelasticity on rupture propagation is less known; however, it can create an effective bimaterial response promoting slip in certain directions and discouraging it in others (Dunham & Rice, 2008;Heimisson et al., 2019). The poroelastic effect in studies so far is smaller than that of moderate or efficient thermal pressurization; for example, the poroelastic changes in normal stress are 10% of friction change or less in Dunham and Rice (2008). Both dilatancy and poroelastic effects warrant further study in future work, including laboratory studies of dilatancy in well-slipped gouge layers typical of mature natural faults.

Conclusions
Our modeling of faults with rate-and-state friction and enhanced dynamic weakening indicates that average shear prestress before dynamic rupture-which can serve as a measure of average fault strength-can be scale-dependent and decrease with the increasing rupture size. Such decrease is more prominent for faults with more efficient dynamic weakening. The finding holds for faults with the standard rate-andstate friction only, without any additional dynamic weakening, although the dependence is relatively unremarkable in that case (Figures 7 and S4). However, the scale-dependent decrease in average prestress is quite pronounced even for fault models with mild to moderate enhanced dynamic weakening that satisfy a number of other field inferences, including nearly magnitude-invariant static stress drops of 1-10 MPa, increasing average breakdown energy with rupture size, radiation ratios between 0.1 and 1.0, and low-heat fault operation (Lambert et al., 2021;Perry et al., 2020).
Our simulations illustrate that both critical fault stress required for rupture propagation and static stress drops are products of complex finite-fault interactions, including wave-mediated stress concentrations at the rupture front and redistribution of stress post-rupture by dynamic waves. Hence it is important to keep in mind the finite-fault effects-and their consequences in terms of the spatially variable fault prestress, slip rate, and shear stress evolution-when interpreting single-degree-of-freedom representations, such as spring-slider models and small-scale laboratory measurements. This consideration highlights the need to continue developing a better physical understanding of faulting at various scales through a combination and interaction of small-scale and intermediate-scale lab and field experiments, constitutive relations formulated based on such experiments, and finite-fault numerical modeling constrained by inferences from large-scale field observations. Our comparison of local fault behavior in SDOF and dynamic rupture simulations also demonstrate how small-scale experiments can be used in conjunction with finite-fault modeling to improve our understanding of the earthquake source: the finite-fault modeling can suggest the initial conditions and slip-rate histories for the small-scale experiments to impose, and then the shear stress evolution from the small-scale experiments can be compared to the numerically obtained ones, which would allow to further improve the constitutive laws used in finite-fault modeling.
We find that increasingly efficient dynamic weakening leads to different earthquake statistics, with fewer small events and increasingly more large events. This finding adds to the body of work suggesting that enhanced dynamic weakening may be responsible for deviations-inferred for large, mature fault segments-of earthquake statistics from the Gutenberg-Richter scaling (Bouchon & Karabulut, 2008;Hauksson et al., 2012;Jiang & Lapusta, 2016;Michailos et al., 2019;Sieh, 1978). For example, fault models with efficient dynamic weakening are consistent with mature faults that have historically hosted large earthquakes but otherwise appear seismically quiescent, such as the Cholame and Carrizo segments of the San Andreas Fault, which hosted the 1857 Fort Tejon earthquake (Jiang & Lapusta, 2016). In contrast, the presence of a wider range of rupture sizes and styles of slip transients on other faults, such as the Japan trench (e.g., Ito et al., 2013), may suggest that they undergo more mild to moderate enhanced weakening during dynamic ruptures, and/or exhibit more pronounced fault heterogeneity.
Such considerations may be useful for earthquake early warning systems, which currently do not take into account the potential physics-based differences in the event-size distribution. Under the assumption of Gutenberg-Richter statistics, the probability that a smaller, Mw 5 or 6 event becomes a much larger earthquake is not great; however, that probability may be substantially larger on mature faults if they are indeed governed by enhanced dynamic weakening.
Our results indicate that critical stress conditions for earthquake occurrence cannot be described by a single number but rather present a complex spatial distribution with scale-dependent averages. When considering the critical stress conditions, it is essential to take into account both the size of the rupture and the weakening behavior, and hence the style of motion, that may occur throughout rupture propagation. These results warrant further investigation, specifically how the weakening behavior during dynamic rupture would interact with different degrees of fault heterogeneity as well as implications for earthquake early warning.

Data Availability Statement
The data supporting the analysis and conclusions is given in figures and tables, in the main text and Supporting Information S1. Data is accessible through the CaltechDATA repository (https://data.caltech.edu/ records/1612). We thank Tom Heaton, Hiroo Kanamori, Emily Brodsky, Mark Simons, Jean-Philippe Avouac and Zhongwen Zhan for helpful discussions and the Associate Editor and two anonymous reviewers for thoughtful suggestions that helped improve the manuscript.