Nearly Magnitude-Invariant Stress Drops in Simulated Crack-Like Earthquake Sequences on Rate-and-State Faults with Thermal Pressurization of Pore Fluids

Stress drops, inferred to be magnitude-invariant, are a key characteristic used to describe natural earthquakes. Theoretical studies and laboratory experiments indicate that enhanced dynamic weakening, such as thermal pressurization of pore fluids, may be present on natural faults. At first glance, magnitude invariance of stress drops and enhanced dynamic weakening seem incompatible since larger events may experience greater weakening and should thus have lower final stresses and higher stress drops. We hypothesize that enhanced dynamic weakening can be reconciled with magnitude-invariant stress drops due to larger events having lower average prestress when compared to smaller events. We conduct numerical simulations of long-term earthquake sequences in fault models with rate-and-state friction and thermal pressurization, and in the parameter regime that results mostly in crack-like ruptures, we find that such models can explain both the observationally inferred stress drop invariance and increasing breakdown energy with event magnitude. Smaller events indeed have larger average initial stresses than medium-sized events, and we find nearly constant stress drops for events spanning up to two orders of magnitude in average slip, comparable to approximately six orders of magnitude in seismic moment. Segment-spanning events have more complex behavior, which depends on the properties of the arresting velocity-strengthening region at the edges of the faults.


Introduction
Stress drops and breakdown energy are important descriptors of natural earthquakes. Stress drops characterize the average change in stress state from before to after the dynamic event (Kanamori & Anderson, 1975;Knopoff, 1958;Kostrov, 1974). The stress drop distribution varies along the fault and can be averaged in several different ways in order to produce a single, representative value for an event (Section 3). There is a fair amount of scatter in the inferred average values of stress drops of natural earthquakes, from about 0.1 MPa up to values around 100 MPa (Baltay et al., 2011;Kanamori & Brodsky, 2004). However, the inferred values of stress drop are magnitude-invariant; most events have stress drops that fall between 1 MPa and 10 MPa, and this trend has been observed for events ranging nine orders of magnitude in seismic moment (Abercrombie & Rice, 2005;Allmann & Shearer, 2009;Cocco et al., 2016;Ide & Beroza, 2001). The generality of the inferred magnitude invariance of stress drops is still a topic of ongoing research, with some observations indicating that some individual earthquake sequences may exhibit mildly increasing trends in stress drop with increasing moment (e.g., Cocco et al., 2016;Viesca & Garagash, 2015). The interpretation and reliability of the stress drops estimates have been actively studied recently, with indications that the current standard methods of estimating stress drops can introduce some significant discrepancies between the actual and inferred stress drops (e.g., Kaneko & Shearer, 2014, 2015Lin & Lapusta, 2018;McGuire & Kaneko, 2018;Noda et al., 2013). However, there are no indications at present that the overall nearly magnitude-invariant trend should be questioned.
Breakdown energy, a quantity analogous to fracture energy from singular and cohesive zone models of fracture mechanics, is meant to capture the energy consumed near the rupture tip that controls the dynamics of the rupture front (Cocco et al., 2004;Palmer et al., 1973;Rice, 1980). Breakdown energy is a part of the overall energy budget of a seismic event, with the total strain energy released (ΔW) typically divided into the breakdown energy G, radiated energy E R , and other dissipation E D (Kanamori & Rivera, 2013). It is a more straightforward concept for shear stress versus slip behavior that follows slip weakening during dynamic Several numerical studies used these enhanced dynamic weakening effects to explain some observations for natural earthquakes. Thermal pressurization of pore fluids can explain the inferred increase in breakdown energy with the increasing event size (Rice, 2006;Viesca & Garagash, 2015); this has been shown using simplified theoretical arguments. Models with dynamic weakening have been successful in producing fault operation at low overall prestress and low heat production (Noda et al., 2009;Rice, 2006) as supported by several observations (Brune et al., 1969;Hickman & Zoback, 2004;Williams et al., 2004;Zoback et al., 1987).
However, it is not clear whether enhanced dynamic weakening is consistent with magnitude-invariant stress drops. In the following intuitive scenario, they are not. Let us assume that smaller and larger events nucleate at nearly the same level of average prestress. The smaller event has less slip and thus weakens a smaller amount. This results in a smaller breakdown energy (the dotted region) and a higher final stress. The larger event weakens more and has a larger breakdown energy and lower final stress. In this scenario, larger events would have systematically larger stress drops and larger breakdown energy ( Figure 1b). However, this intuitive scenario may be incorrect, due to the following hypothesis which is illustrated and supported by the simulations in this work. Both smaller and larger events would nucleate at locations with relatively high prestress, matching the quasi-static frictional strength. But we must consider the average initial stress of all points involved in the rupture, not just those involved in nucleation. Larger events would have larger slips and hence dynamically weaken more and may be able to propagate over areas of much less favorable (lower) prestress conditions. This means that the initial stress averaged over the entire rupture area may be lower for larger events than that for smaller events. Overall, larger events would dynamically weaken more and potentially arrest at a lower average final stress, but they would also have occurred with lower average initial stress. Thus, the average stress drop can be similar for smaller and larger events ( Figure 1c). However, the observed increase of the breakdown energy with event size is still preserved.
Here, we use fully dynamic simulations of earthquake sequences on rate-and-state faults to investigate this hypothesis and study if enhanced dynamic weakening can indeed be compatible with magnitude-invariant stress drops while also maintaining increasing breakdown energy with increasing event size. Different dynamic weakening mechanisms produce different weakening behaviors, but here we focus on thermal pressurization as a representative dynamic weakening mechanism that can lead to continuous fault weakening with earthquake-source slip. We consider the simplest scenario that allows us to explore this hypothesis, that of a seismogenic fault segment with uniform properties of quasi-static fault strength. For heterogeneous faults, the argument should still hold, since larger ruptures with larger slip and hence more pronounced weakening should be able to propagate over larger areas of locally unfavorable prestress, as compared to smaller ruptures, potentially still resulting in nearly magnitude-invariant stress drops, but with some scatter due to heterogeneity. Such scenarios will be investigated in future work.
We indeed find that the hypothesis of lower average initial stress before larger events holds for a wide range of events in our simulations that arrest within the seismogenic region, resulting in nearly magnitude-invariant 10.1029/2019JB018597 stress drops, at least for the range of parameters considered in this work that results in mostly crack-like ruptures. Our fully dynamic simulations also confirm the increase in breakdown energy with the event size consistent with observations. For seismogenic-region-spanning events, we find that the properties of the velocity-strengthening areas can have a profound impact on the stress drop. Models with large values of velocity strengthening do not allow ruptures to propagate much into the velocity-strengthening region, thus leading to higher stress drops, whereas models with smaller values of velocity strengthening allow farther propagation and thus lower stress drops.
For completeness, we start by investigating faults without enhanced dynamic weakening, with the Dieterich-Ruina rate-and-state friction only. Consistent with related findings of prior studies, we find that the stress drops are also magnitude-independent, but so is the breakdown energy. This is because Dieterich-Ruina rate-and-state friction resembles linear slip-weakening during dynamic rupture (Cocco & Bizzarri, 2002;Lapusta & Liu, 2009), which has prescribed and process-independent dynamic resistance and breakdown energy.
We also use our modeling to examine the accuracy of seismically estimated breakdown energies G SE , by comparing the values computed directly from the on-fault variables with inferred values G SE computed indirectly from seismically available observations.
Here we follow the assumption that most of the breakdown energy occurs on the shearing surface (e.g., Rice, 2006;Viesca & Garagash, 2015). While it is clear that some energy is dissipated in off-fault damage (Andrews, 2005;Poliakov et al., 2002), especially on rough, non-planar faults (Dieterich & Smith, 2009;Dunham et al., 2011), those amounts may be negligible compared to seismic estimates of breakdown energy, at least for relatively planar mature faults. The relative importance of the off-fault and on-fault dissipation during dynamic rupture is an important topic of ongoing studies which is beyond the scope of this work.

Fault Model Formulation
Our simulations are conducted following the methodological developments of Lapusta et al. (2000) and Noda and Lapusta (2010). In order to study long sequences of seismic events in simulations with enhanced dynamic weakening, we consider a mode III, two-dimensional (2-D) model with a one-dimensional (1-D) fault embedded into a 2-D uniform, isotropic, elastic medium ( Figure 2a). The earthquake sequences on the fault are simulated in their entirety: the nucleation process, the dynamic rupture propagation, postseismic slip that follows the event, and the interseismic period between events that can last up to tens or hundreds of years ( Figure 2b). In all models, the laboratory-derived rate-and-state friction (section 2.1) operates on the fault. Our 1-D fault ( Figure 2a) contains a velocity-weakening (VW) region surrounded by velocity-strengthening (VS) regions. The fault slip at the plate rate (V pl = 10 −9 m/s) is prescribed at the edges of the model. We begin with a standard rate-and-state model but then add thermal pressurization of pore fluids (section 2.2). Parameters for the specific models are listed in Tables 1 and 2. While many events arrest within the VW region, some span the entire VW region (Figure 2b). We refer to the events that span the entire VW region as "complete rupture" events and those that arrest within the VW region as "partial rupture" events.

Rate-and-State Friction
We use the laboratory-derived rate-and-state laws with the aging law proposed by Dieterich (1979) and Ruina (1983): where is the normal stress (constant in time), is the shear stress, is the friction coefficient, V is the slip velocity, p is the pore pressure, is the state variable, L is the characteristic slip for the evolution of the state variable, * is the reference friction coefficient corresponding to a reference slip rate V * , and a and b are constitutive parameters. At steady state (constant slip velocity), the values of and evolve to be their steady-state values ss and ss given by These steady-state relations show that the difference between the parameters a and b controls the fault behavior at steady state. If (a − b) > 0, then the fault has velocity-strengthening (VS) friction behavior in which increases in slip velocity result in increases in shear resistance. This leads to stable sliding on the fault  under steady loading. If (a − b) < 0, then the fault has velocity-weakening (VW) behavior. In this case, an increase in slip velocity leads to a decrease in shear resistance, making these regions of the fault potentially seismogenic (Rice & Andy, 1983;Rice et al., 2001;Rubin & Ampuero, 2005).
We first consider models with the standard rate-and-state formulation and no additional dynamic weakening, with parameters given in Tables 1 and 2.

Enhanced Dynamic Weakening due to Thermal Pressurization of Pore Fluids
Laboratory experiments have shown that the rate-and-state laws (Equations (1)-(4)) work well for relatively slow slip rates (10 −9 to 10 −3 m/s). However, at seismic rates of ∼1 m/s, additional dynamic weakening mechanisms, such as thermal pressurization, can be present. Thermal pressurization occurs when fluids within the fault heat up, expand, and pressurize during dynamic rupture, reducing the effective normal stress (Noda & Lapusta, 2010;Rice, 2006;Sibson, 1973). The thermal pressurization effect is governed in our model by the following coupled differential equations for pressure and temperature evolution (Noda & Lapusta, 2010): where T is the temperature of the pore fluid, h is the hydraulic diffusivity, th is the thermal diffusivity, V is the source of shear heating distributed over the shear zone of half-width w, c is the specific heat, is the distance normal to the fault plane, and Λ is the coupling coefficient that gives pore pressure change per unit temperature change under undrained conditions.
The efficiency of the thermal pressurization process depends on the interplay of several of these parameters. Shear heating, V, must be strong enough to raise the temperature, given both the specific heat of the rock, c, and the half-width of the shear zone, w. Furthermore, this heat generation must not be dissipated too quickly by the thermal diffusivity, th , of the system. If sufficient heat is generated, the temperature of the system increases, and this increase is coupled into an increase in pressure of the fluid. The fluid then pressurizes as long as the hydraulic diffusivity, h , is not too large. Several of these parameters are relatively well constrained from laboratory experiments: th = 10 −6 m/s, Λ = 0.1 MPa/K, and c = 2.7 MPa/K The stress versus slip evolution at three example points illustrate different behaviors along the fault (Row 4). Initial and final stresses during the event are marked (open circles) for each point, and some previous slip history is also shown preceding the initial stress marker. Point 3 shows the evolution in the nucleation zone, point 1 is in the region where the event arrests, and point 2 shows behavior of a well-established rupture (the majority of the ruptured points experience this behavior). Note that this event is crack-like, and the final stress is nearly equal to the dynamic level of stress for the three representative points. (Noda & Lapusta, 2010;Rempel & Rice, 2006;Wibberley & Shimamoto, 2005). Thus, the efficiency of the process is effectively controlled by the half-width w and hydraulic diffusivity h , which can vary significantly: w can vary from 10 −3 m to 10 −1 m and h can vary from 10 −2 m 2 /s to 10 −5 m 2 /s (e.g., Rice, 2006). Changing these two parameters within these ranges can make thermal pressurization either very efficient or completely negligible. The values we have chosen are motivated by prior studies (Rice, 2006;Noda & Lapusta, 2010) and are given in Tables 1 and 3.

Representative Simulated Events
Our simulations produce sequences of dynamic events together with interseismic creep, including aseismic nucleation processes ( Figure 2b). However, here, we focus on the properties of individual dynamic events. A sample dynamic event from our simulations is shown in Figure 3. In general, both slip throughout the event and final slip vary along the fault. The spatially varying initial and final shear stress distributions along the fault lead to a stress drop distribution that varies along the fault. Most of the ruptured area experiences a decrease in shear stress during the event, but both edges of the ruptured area in each of the events show an increase in stress (and hence negative stress drop). The shear stress versus slip evolution along the fault is illustrated for three representative locations. Locations near the nucleation region experience a small coseismic stress drop, with much of the stress change at these points achieved aseismically, during nucleation. Points near the middle of the ruptured area show the expected increase in stress to a peak value, followed by a drop, controlled by our rate-and-state laws, down to some near constant ( Figure 3) dynamic value. Where the event arrests, points only slip a small amount and do not completely weaken down to the expected dynamic level of shear stress.
Observations of natural earthquakes cannot resolve these fine variations in stress, slip, slip rate, etc. at all points along the fault as we are able to do in our simulations. Thus, natural events are often described by a single, average value for stress drop and average final slip. In the next sections, we discuss the condensing of heterogeneous slip and stress-drop distributions into average values for the seismic events.

Computation of Stress Drops and Breakdown Energy
We follow the averaging methodologies described in Noda et al. (2013), modified to fit our two-dimensional model, since our relevant variables are scalar fields rather than vector fields. The initial distribution of shear traction on the fault before an earthquake is denoted by i (x). An earthquake produces a slip distribution (x) and the traction along the fault changes to (x). The stress drop distribution is defined as

Averaging of Stress Drop Distribution Based on Seismic Moment
Seismologically estimated values of average stress drop are often based on the seismic moment M 0 of the event as well as the fault dimensions; the following formula is typically used (Kanamori & Anderson, 1975): where A is the ruptured area, = A 1∕2 is the characteristic spatial dimension, and C depends on the shape and aspect ratio of the ruptured domain: C = 2.44 for a circular ruptured area and increases for rectangular areas with larger aspect ratios (Noda et al., 2013).
If the actual stress drop is uniform over the ruptured domain Σ, then Δ M is exactly equal to that value. However, as evident for our example events (section 2.3), the stress drop across the fault is heterogeneous and given by the distribution Δ (x). In this case, Δ M is a weighted average of Δ (x). This average is weighted by the (elliptic) slip distribution E 12 that gives a uniform stress drop over the same ruptured domain (Madariaga, 1979):

Spatial Averaging of Stress Drop
The spatially averaged stress drop can be expressed as the integral of the stress drop of all ruptured points along the fault divided by the ruptured domain Σ: The stress change at every point has equal weighting of one in this averaging method, unlike Δ M where E 12 weights points differently along the fault. Similarly to Δ M , Δ A depends only on points in the ruptured domain. Considering the entire fault can result in severely underestimating the average stress drop of the event.
The ruptured domain Σ is defined as the region with non-zero slip (which is a line for our model, but a 2-D area in general):

10.1029/2019JB018597
However, it is difficult to precisely determine Σ for observed events due to non-uniqueness and smoothing when finding a solution to an inverse problem. Furthermore, in our models, the fault is prescribed to creep outside the locked, velocity-weakening region, and thus, there is small non-zero slip everywhere on the fault during every event. It is appropriate to only consider points where the inertial term becomes significant, but there is no exact quantitative criterion to define that, so we instead approximate this by defining the ruptured domain Σ to consist of locations that exceed a slip rate of 0.1 m/s during the event: Altering the seismic velocity threshold may change the effective rupture size. However, there is a sharp falloff in slip rate outside the ruptured area down to the creeping rate many orders of magnitude below the seismic slip rate. Thus, changing this threshold by even an order of magnitude does not change the rupture size appreciably.

Averaging of Stress Drop Distribution Based on Energy Considerations
The third method of averaging Δ (x) is consistent with energy partitioning (Noda & Lapusta, 2012;Noda et al., 2013). This stress drop is also part of the averaged shear stress versus slip evolution curve that conserves both the total strain energy released ΔW as well as the dissipated energy E D as discussed in section 3.4. Here, the final slip distribution (x) is used as the weighting function: In this method, the ruptured domain is implicitly defined by the slip distribution (x).
The three averaging methods ((9)-(10)) and (13) give similar but not identical results for the average stress drop for a given event. Noda et al. (2013) proved that Δ E ≥ Δ M and observed that Δ M ≥ Δ A in their simulations. Given that computing seismic moment on our 1-D faults requires additional assumptions of rupture aspect ratio and shape, we focus on computing the energy-based stress drop Δ E and the spatially averaged stress drop Δ A in this study, where the moment-based stress drop would be expected to lie in between these two values. For similar reasons, in this study, we present relationships between average stress drop and average slip, rather than moment. Examining these scaling relationships in 3-D calculations is a topic for future work.

Calculation of Energy Balance and Breakdown Energy G in Simulations
In our dynamic simulations, the slip and stress evolution is determined at every point along our fault at all times. As such, we are able to calculate the breakdown energy directly in our model. This can be done by integrating the breakdown energy along the fault for all ruptured points. Furthermore, we can construct a representative average curve for the event and use it to illustrate the breakdown energy.
In the earthquake energy budget per unit area, illustrated in Figure 4, the total strain energy released, ΔW∕A, is partitioned into dissipated energy per unit area, E D ∕A, which is the area underneath the stress-slip curve, and radiated energy E R ∕A: We write the balance per unit area because the breakdown work, G, is defined per unit area. The total strain energy released ΔW∕A is given by The average curve for the single event with enhanced dynamic weakening from Figure 9. The fault continues to weaken by more than an additional 5 MPa as it accumulates slip, leading to a larger breakdown energy. In both (a) and (b), the energy-based static stress drop Δ E is the difference between the average initial and final shear stresses. The total strain energy released ΔW∕A is outlined by the black dashed line; the associated trapezoid ends at the x axis (not shown). The dissipated energy E D ∕A is given by the total area underneath the stress versus slip curve (dotted + gray). Breakdown energy G is the subset of the dissipated energy labeled by the dotted area. Radiated energy can be calculated by subtracting total dissipated energy from the total strain energy released.
wherēis the average final slip for the event,̄i is the average initial shear stress weighted by the final slip, and̄is the average final shear stress. For our 1-D fault, let us define the edges of the ruptured domain Σ as L 1 and L 2 . Then the dissipated energy can be computed as The remainder of the total strain energy released is the radiated energy: The dissipated energy E D ∕A can further be partitioned into the breakdown energy G (Palmer et al., 1973;Rice, 1980) and frictionally dissipated energy E F ∕A which makes up the remainder (labeled as "other dissipation" in Figure 4): The breakdown energy G a is analogous to the fracture energy of fracture mechanics and can be calculated as where we use G a to indicate the "actual" or on-fault value of G.
One can illustrate the energy balance by a representative average shear stress versus slip curve ( Figure 4). We follow the averaging methodology of Noda and Lapusta (2012) to perform this calculation, which involves taking the stress versus slip evolution of every ruptured point and averaging them in slip rather than in time. Thus, this can only be done once the event is complete, and the stress versus slip evolution is known everywhere. The averaging method preserves total strain energy released ΔW∕A and total dissipated energy E D ∕A. Every ruptured fault location has, in general, a different amount of total slip (x), so the stress versus slip curves at each point are scaled in slip bȳ∕ (x) so that each point has the same average slip̄. Then the stress values are scaled by the factor of (x)∕̄, thus preserving the areas representing E D . Once all shear stress versus slip curves are scaled, the stress values at each value of slip are averaged among the curves. We can then calculate our energy quantities from this average curve. The strain energy released per unit area

10.1029/2019JB018597
ΔW∕A is given by the trapezoid indicated by the dashed line in Figure 4, and the dissipated energy per unit area E D ∕A is given by One can also compute the quantity motivated by the breakdown energy from the average curve, here titled as G curve : The average curve construction has been shown to preserve total strain energy released ΔW∕A and dissipated energy E D ∕A (Noda & Lapusta, 2012). However, it does not necessarily preserve the breakdown energy as the minimum shear stress of the average curve does not have a simple relation to the minima of the curves of each ruptured point. We later show that G a has a similar, but not identical, value to G curve for the crack-like ruptures considered in this study, and hence, G curve can be used to visualized G a .
Note that G a and G curve have units of energy per unit area, while ΔW, E D , and E R denote the energies per event and have units of energy. Representations of the type shown in Figure 4 show energies per unit area, and that is why we have been considering quantities ΔW∕A, E D ∕A, and E R ∕A. To compute the corresponding energies per event, one needs to multiply them by the total ruptured area.

Stress Drop and Breakdown Energy G from Observations
We seek to match the observed trends of magnitude-invariant stress drop and increasing breakdown energy G with increasing event size (Abercrombie & Rice, 2005;Ide & Beroza, 2001;Viesca & Garagash, 2015). However, as discussed earlier, these values cannot be directly measured in observed events and instead must be inferred from other observations. Stress drop is often calculated using the moment-based average (Equation (9)). For large events, the rupture shape and dimension is found from finite-fault inversions (Kanamori & Brodsky, 2004, and references therein). For small events, for which finite-fault inversions are not feasible, the spectral representation of the seismic waveforms is fitted by a model based on a circular crack with constant rupture speed to obtain the long-period displacement amplitude Ω 0 and corner frequency c measurements. These parameters are then used to calculate M 0 from Ω 0 (Brune, 1970) and the source radius r from c assuming a circular rupture and constant rupture velocity of 0.9c s (Madariaga, 1976).
The breakdown energy can be estimated from observations as follows (Abercrombie & Rice, 2005): where G ′ is the approximation for the breakdown energy G, Δ is the seismologically estimated (static) stress drop, is the shear modulus of the rock material,̄is the average slip of the event, M 0 is the seismic moment, and E R is the radiated energy. The relationship between of G ′ and the average breakdown energy assumes that (1) the initial stress is the peak stress and (2) that there is no stress overshoot or undershoot at the end of the event, making it potentially different from the actual G (see Figure 2 of Abercrombie & Rice, 2005). We refer to this G ′ as seismologically estimated breakdown energy G SE .

Theoretical Predictions for Breakdown Energy and Stress Drops on Rate-and-State Faults
Based on previous studies and theoretical considerations (Ampuero & Rubin, 2008;Cocco & Bizzarri, 2002;Lapusta & Liu, 2009;Rubin & Ampuero, 2005), we expect both the breakdown energy and the static stress drop to remain approximately the same for events of different sizes on a fault with uniform rate-and-state properties. This is because, at the rupture tip, the fault governed by the standard rate-and-state formulation behaves essentially as one governed by linear slip-weakening friction: 10.1029/2019JB018597 where LSW refers the linear slip weakening evolution of shear stress from the peak shear stress, p , at initial slip ini to the dynamic level of shear resistance d over the critical slip-weakening distance D c . The weakening rate W is defined as For the standard rate-and-state formulation, one can write the initial stress i from (1): As slip rate abruptly increases from near-zero V ini to dynamic V d n at the crack tip, stress will increase to some peak value p , which can be approximated by Assuming that the slip acceleration occurs at negligible slip and hence with no state evolution. As slip accumulates, the stress further evolves to a steady-state dynamic level given by This weakening effect occurs at weakening rate W: and hence, the evolution occurs over the effective critical slip-weakening distance D c given by If the final stress is approximately equal to the dynamic resistance, then we expect: These quantities depend on the dynamics of the process through V ini , V d n , and ini , but this is a weak dependence since they are contained within logarithms and changes of even an order of magnitude alter the final product by only a small amount. There is a much stronger dependence on the friction parameters a, b, and L, which are constant in a given model.

Dependence of G and
on Magnitude for Given a, b, and L Indeed, our simulations show that for uniform frictional parameters a, b, and L along the fault, both G a and Δ are nearly constant for events of different sizes. Both trends are evident in the accumulated slip profiles and average curves for three events of different sizes from the same simulation of earthquake sequences ( Figure 5). Larger events accumulate more slip and rupture longer fault stretches, but the breakdown energy (dotted area) and static stress drop are nearly equal for the three illustrated events.
There are some slight trends in G and Δ due to the dynamics of the process. Larger events tend to have lower average initial stresses, due to rupturing longer fault stretches, building more stress concentration, and entering slightly less favorably stressed regions. All events weaken down to approximately the same dynamic level, as expected. This leads to a mild decrease in the static stress drop from Δ E = 3.3 MPa for the smallest event down to Δ E = 1.8 MPa for the largest event. The peak stress p slightly increases with the event size, due to more stress concentration during the larger event and higher initial values for the state variable ini from longer recurrence times. The outcome is slightly higher breakdown energies as the The initial and final stresses are marked by circles; the breakdown energy is indicated by the dotted area. For a given value of L, the breakdown energy remains nearly constant. The stress drop slightly decreases, and the breakdown energy slightly increases with the event slip, as discussed in the text. These three events are marked with gray, downward-pointing triangles in Figure 6. event size increases. However, these two effects produce relatively small variations, within a factor of two, in both G and Δ .
We find that these trends extend for all events in our simulations (Figures 6 and 7a). For L = 250 μm (black circles), events differ by nearly an order of magnitude in slip, from 0.01 m to 0.1 m. The corresponding stress drops are nearly constant, around 2-3 MPa, with a slight decreasing trend with the increasing event size. The breakdown energies are also approximately constant, with a slight increasing but saturating trend (Figure 7a) for all events.  (a) Breakdown energies G a for events from several simulations with standard rate-and-state friction and L ranging from 0.125 mm to 4 mm (no enhanced dynamic weakening). Complete rupture events are marked with filled-in shapes. Increasing L leads to an increase in the breakdown energy. But breakdown energy only slightly increases and saturates for events with the same L. The two largest values of L lead to almost exclusively complete rupture events because the nucleation size is too large to produce small events given the size of the VW region. (b) Breakdown energies from simulated events overlaid on observational inferences for natural events from Rice (2006). The values are similar, though systematically lower, and the standard rate-and-state model produces breakdown energies that do not increase at the same rate as those inferred from natural events.
The example event discussed earlier (Figure 3) shows the expected behavior for the standard rate-and-state case. This event has the area-averaged stress drop of Δ A = 2.4 MPa, which matches well with the stress drop distribution seen in Figure 3. The entire ruptured domain is plotted in Figure 3, including penetration into the velocity-strengthening region. This is evident from the negative stress drops found at the edges of the event, greater than 3 km away from the center of the fault. Three representative points are chosen to show the variability of the stress versus slip evolution along the fault. The point at 2.4 km is in the nucleation zone and experiences mostly aseismic stress evolution (solid line preceding initial stress point) followed by little coseismic stress change with slip. The point in the arrest zone (−3.6 km) shows a very different behavior, with an increase to a peak level and a drop. However, the stress drop is negative (stress increase), owing to the velocity-strengthening properties of the fault at this point. The point at the center of the fault (0 km) is representative of the behavior of the majority of the fault. This point shows the typical rate-and-state behavior with an increase to a peak level of stress followed by a drop to a near-constant dynamic level of stress. This point experiences a stress drop similar to the average for the entire event. All of the points  Figure 6. (b) Three sample events with comparable stress drops, but varying final slips, from rate-and-state simulations with different L. (Row 1) Accumulated slip profiles and (Row 2) average stress versus slip curves. Increasing L increases both the slip weakening distance D c as well as the breakdown energy, but does not affect the average stress drop. These three events are marked with gray upward-pointing triangles in Figure 6. on the fault are averaged to create the illustrative average curve (Figure 4a). From the average curve, it is apparent that the majority of points follow the behavior qualitatively similar to the point at 0 km. Note that the energy-based stress drop from the average curve is Δ E = 4 MPa, which is higher than the Δ A = 2.4 MPa as expected (Noda & Lapusta, 2012).

Increasing G and Magnitude-Invariant
with Increasing Values of L Breakdown energy has a weak dependence on the dynamics in a standard rate-and-state fault model, but it has a stronger, quasi-linear, dependence on the characteristic slip distance L. One of the ways to reproduce an Figure 9. A representative event for the models with thermal pressurization. The plotting conventions are the same as on Figure 3. The three sample points exhibit decreasing dynamic stress with slip throughout the event, illustrating the effect of additional dynamic weakening due to thermal pressurization. increase in G a with average slip is to systematically increase L, which also systematically alters the effective critical slip-weakening distance D c (Figure 8). The peak stress of each event also increases, predominantly due to a longer recurrence time that results in fault strengthening. Increasing L increases the nucleation size of the event, and thus, a stress increase must penetrate further into the VW fault before an event nucleates, leading to a higher initial state variable ini , higher initial stress i , and higher peak stress p . This is even the case for events with the same amount of average slip (Figure 8a). However, the increase in the critical slip-weakening distance is clearly the main contributing factor to the increased G a . The dynamic levels of stress are nearly constant in all three cases as expected; this level does not directly depend on L. The stress drops increase with increasing L for these three events, due to the fact that we have chosen three events with very similar slips (Figure 6a, star symbols). Stress drops for the entire sequence of events do not change as Figure 10. Three sample partial rupture events from the same simulation with a 12-km-long velocity-weakening region and thermal pressurization. (Top) Accumulated slip profiles of the three events. (Bottom) Average shear stress versus slip curves. The initial and final stresses are marked by circles; the breakdown energy is indicated by the dotted area. As event size increases, both the average initial stress and average final stress decrease, so that the stress drops remains nearly constant at ∼7 MPa.
we increase L (Figure 6a). This is illustrated by selecting three other events that no longer have the same average final slip (Figure 8b), but do have comparable stress drops.
Varying L over an order of magnitude from 125 μm to 4 mm leads to a clear increase in breakdown energy (Figure 7a) that is much larger than the slight increasing trend we find for larger events of a given L. There are clear groups of events with similar breakdown energies, corresponding to simulations with each value of L. The values for the breakdown energies compare favorably to those from Rice (2006), though they are systematically lower, particularly at higher values of slip (Figure 7b). For a given L, the simulated breakdown energies level off and do not capture the observed trend. Even increasing L is not completely sufficient to match the observed trend.
Simulations with all values of L have comparable stress drops, determined by values of a and b. All of our calculated stress drops fall into the 1-3 MPa range which is consistent with inferred stress drops from natural events. We find two distinct trends when separating partial rupture from complete rupture events ( Figure 6). The first trend is that the partial rupture events show a slight decrease in stress drop with increasing slip. This is because all events arrest at similar levels of average final stress. However, as discussed in the previous section, larger events initiate with slightly lower average levels of prestress (unless they are complete ruptures, as discussed below) and thus have smaller stress drops. The second is for the complete rupture events; these events have the same ruptured domain, and the ones with larger slip correspond to larger stress drop, reflecting variability in the prestress level for complete rupture events. Note that, for each particular value of L, the decrease in stress drop with slip is within a factor of 3; however, the distribution of stress drops across the full set of events for all values of L is nearly magnitude-invariant, with a scatter well within that inferred for natural earthquakes.

Nearly Magnitude-Invariant Stress Drops and Increases in Breakdown Energy in Earthquake Sequence Simulations with Thermal Pressurization
We consider a 12-km-long VW segment surrounded by two 24-km-long VS sections and then increase our seismogenic zone from 12 km to 24 km in order to further expand the range of the simulated event sizes. Extending our fault to 24 km allows for a greater range of event sizes, with slips ranging from ∼0.07 m to ∼10 m.
Our simulations with a 12-km-long VW segment produce a range of events, with average slips of 0.1 m up to 5 m. One of the events is illustrated in more detail in Figure 9. It nucleates in an area of higher prestress and propagates along the fault until it reaches lower levels of prestress that are unfavorable enough to arrest the Figure 11. Stress drops Δ A (top) and Δ E (center) for events in the simulation with thermal pressurization and a 24-km-long velocity-weakening region. Events with complete ruptures are denoted by stars. (Bottom) Spatially averaged initial stress̄i and final stress̄in the simulation with thermal pressurization and a 24-km-long fault. Partial rupture events exhibit a decrease in both average initial and final stresses with increasing slip such that the change in average stress drop is relatively minor over a decade increase in average slip, resulting in nearly magnitude-invariant stress drops. Parameters from linear fits between the average stress drops and the logarithm of average slip are given in Table 4. event. The shear stress versus slip behavior is shown for three representative points. All three points show continuous weakening with slip, illustrating that thermal pressurization is acting effectively along the entire fault. The point in the nucleation zone (−5.25 km) again shows significant aseismic stress evolution (solid line preceding the initial stress point), followed by lesser coseismic stress change with slip. The other two points along the fault (−3.75 km and −2.4 km) show the expected behavior for most ruptured points with an initial increase and rapid decrease in stress (similar to the standard rate-and-state behavior) followed by a continuous decrease in stress with slip (due to dynamic weakening from thermal pressurization). The average curve for this event (Figure 4b) shows the behavior similar to the points outside the nucleation zone. Note that this event, as are others in our models, is largely crack-like, that is, have local durations of slip which are comparable to the overall rupture duration.
To illustrate how stress drop and breakdown energy vary with the event size, we consider three representative events with progressively larger average slip (Figure 10). The smallest event (Event 64) has the highest average prestress and also the highest average final stress. The intermediate-size event (Event 33) has a lower prestress, and it weakens more so it also has a lower final stress. The largest event (Event 20) has the lowest average initial stress, and it weakens the most, so it also has the lowest average final stress. As a result, all three events have approximately the same stress drop Δ E of 7 MPa. As the average slip of the events increases, so does the breakdown energy ( Figure 10). This increase in the breakdown energy is due to the additional dynamic weakening, as expected, based on considerations in Rice (2006).
Let us consider the stress drops for all events, using the 24-km fault models. Both energy-based stress drops Δ E and area-averaged stress drops Δ A are calculated (Figure 11). For the partial rupture events, the stress drops appear approximately constant, for the average slips ranging from 0.05 m to 2 m. The energy-based stress drops are higher than the area-averaged ones, consistent with Noda et al. (2013). We perform a linear fit between both the spatial and energy-based average stress drops and the logarithm of the average slip Table 4 Parameters From Linear Fit to Trends in Average Stress Drop and Log-10 Slip, as Shown in Figure 11 Model set Slope (MPa/log 10 (m)) Intercept (MPa) STD (MPa) for sets of partial and complete ruptures ( Table 4). The spatially averaged stress drops for partial ruptures in both models exhibit a mild trend with average slip, resulting in a 10% to 20% increase over a decade of average slip. We may illustrate how this would correspond to changes in seismic moment using the common approximation assuming that the stress drops are indeed magnitude-invariant and therefore that the average slip and rupture radius for a circular crack increase linearly with each other, resulting in a cubic relationship between moment and average slip. Thus, a decade of average slip corresponds approximately to three orders of magnitude in seismic moment or two units in moment magnitude. The energy-based average stress drop shows a stronger relationship with average slip, though the increase in stress drop with slip for partial ruptures in both models results in an increase of only around a factor of 1.5 over a decade of average slip. We consider these trends to exhibit near magnitude invariance, since the changes in average stress drop are relatively mild in comparison to the variation in average slip, with the resulting trend most likely not being discernable given the wide scatter and uncertainties in seismological inferences. Moreover, the overall weakening due to thermal pressurization increases far more substantially than the average stress drop with event size (Figure 10).
These findings confirm our hypothesis that larger events weaken more but also tend to occur at lower average initial stress, thus keeping stress drops relatively constant over a range of event sizes. In fact, for the entire sequence of partial rupture events, both average initial and average final stresses decrease with the increasing event size (or slip) ( Figure 11). The complete rupture events break the nearly magnitude-invariant trend, exhibiting average stress drops that increase more substantially with event size. However, the stress drops for these events are strongly affect by the properties of the VS region, as further discussed in section 6, and there is a range of VS properties for which these events also exhibit nearly magnitude-invariant stress drops.
Breakdown energy G a computed using the on-fault quantities from our simulations increases with increasing event size ( Figure 12) and matches estimates of breakdown energies for natural events, as expected from the simplified theoretical considerations in Rice (2006). We also compare the true breakdown energy G a and estimated value G SE for our simulated ruptures. The comparison ( Figure 13) shows that the actual and estimated values agree relatively well in the majority of cases, within a factor or two. This is because the ruptures are close to being crack-like, the case for which the estimate of G SE was developed. Moreover, despite the average initial stress not being the same as the peak stress in our simulated ruptures, the estimated value G SE still provides a reasonable representation of the actual average value G a . The strength excess increases the breakdown energy G a with expense to the radiated energy E R ∕A, so that the seismological estimate (24) still provides an adequate representation for the crack-like ruptures in our simulations. Our preliminary studies with stronger enhanced dynamic weakening that often leads to self-healing pulse-like ruptures (e.g., Noda et al., 2009) shows that G SE is a poor estimate in that case; an alternative estimate for the self-healing pulse-like case has been developed by Viesca and Garagash (2015). Next, comparing G a to breakdown energy calculated from the average curves G curve we see good, but not perfect agreement ( Figure 13). This is expected since the averaging process preserves the total strain energy release and the dissipated energy, but not the minimum dynamic level of stress. Therefore, the averaged curves provides a good illustration of G but not the exact value of it.  (Abercrombie & Rice, 2005;Rice, 2006). Our models are able to match the trend of the observed events quite well. That data set includes individual estimates for large earthquakes (triangles) from regional and teleseismic recordings, aftershocks from the 1994 Northridge earthquake (squares), averaged values from large earthquakes with several estimates of G ′ (ovals), and small earthquakes recorded at depth in the Cajon Pass borehole (circles with asterisks) and Long Valley borehole (diamonds).
The temperatures in our simulated shear zones should remain below melting, for self-consistency of the models, since melting and its consequences are not included in our constitutive relations. Figure 14 illustrates the evolution of the maximum temperature change measured within the 24-km velocity-weakening region. Both Models A and B assume relatively low effective normal stresses of 50 and 25 MPa, corresponding to substantial chronic fluid overpressurization. For Model A, with normal stress of 50 MPa and relatively mild thermal pressurization, the largest events increase the fault temperature by over 2000 K, well above the expected equlibrium melting temperature of 1000 • C for wet granitic compositions in the shallow crust (Rice, 2006). Note that the degree of shear heating during frictional sliding would be even more extreme for models incorporating only rate-and-state friction with comparable effective confining stress, as they would result in higher dynamic levels of shear resistance throughout slip. As we further reduce the effective normal stress and increase the efficiency of thermal pressurization, as in Model B, our models are able to accommodate more reasonable fault temperature fluctuations within 500 K, while maintaining the desired trends in magnitude-invariant stress drops and increasing breakdown energy with event size. Our future work will examine models with more efficient enhanced dynamic weakening with more localized shear, including conditions more consistent with slip on a plane as discussed by Rice (2006).

Complete Rupture Events and Effect on Stress Drop of Rupture-Arresting VS Regions
Complete rupture events that rupture all of the VW region tend to have different behavior from partial rupture events. These events do not encounter an area of unfavorable prestress within the VW region, but rather arrest due to the VS barriers. Their spatial extent is approximately equal to the length of the VW region, due to the relatively strong VS barriers adopted, but their slip varies. Hence, their stress drop, which approximately scales with slip divided by the rupture extent, scales with slip, as is evident in Figure 15.
This consideration implies that the stress drop of the complete rupture events can be altered if their extent can vary, due to different lengths of their penetration into the VS barriers. We explore how altering the properties of the velocity-strengthening barrier can affect the stress drops of the complete rupture events using six different models (VS1-VS6) with progressively less velocity-strengthening regions (Table 5). In other words, the VS regions surrounding the VW seismogenic zone become closer to velocity-neutral. We only alter the properties of the VS region; all other parameters match those from Tables 1 and 3. Each model is allowed to produce several complete rupture events and stress drops are plotted against rupture length (Figure 15 top) and average slip (Figure 15 bottom) for each event. We find that the stress drops of these complete rupture events indeed depend on the properties of the VS regions ( Figure 15). For models with moderate to relatively strong velocity strengthening regions (VS1-4), the stress drops for partial rupture events are magnitude-invariant over about one order of magnitude increase in slip. However, for models with stronger velocity-strengthening regions (VS1 and VS2), the largest complete rupture events continue to slip more but are unable to propagate appreciably further into the velocity-strengthening regions. As a result, for models with stronger VS regions, larger complete ruptures have increasingly larger stress drops with slip due to the larger degree of slip being confined in nearly the same spatial region.
As we decrease the degree of velocity strengthening in the VS regions, complete rupture events with larger slip propagate further into the VS region and their rupture length increases (Figure 15 top). Correspondingly, the stress drop of these largest complete rupture events decreases. In fact, for models with the least VS regions (VS4-6) the trend for the complete rupture events changes from that of stress drop increasing with their size to a decreasing trend. For less VS regions, the smaller partial rupture events are also able to Figure 14. Evolution of the maximum temperature change on the fault measured within the velocity-weakening domain in simulations for Model A (black) and Model B (red). Both models assume relatively low effective normal stresses (50 and 25 MPa, respectively) and hence substantial chronic fluid overpressurization, however to maintain reasonable fault temperatures to avoid melting our models also require efficient thermal pressurization, such as in Model B.
propagate further into the VS region, and thus, their average stress drops decrease as well. For the two models with the least VS regions (VS5 and VS6), we see stress drop slightly decrease with increasing event size for all events. The largest complete ruptures also have the lowest stress drops, close to ∼1 MPa. It is clear that the properties of the velocity-strengthening region can have a profound effect on the average stress drops. The exact nature of this effect is best studied in 3-D models with 2-D faults, where the relation of the VS boundary of events to their VW region can be different than in the 1-D faults considered in this work.

Conclusions
We have examined the variations of the average stress drop and breakdown energy with rupture size in fully dynamic simulations of earthquake sequences on rate-and-state faults with and without enhanced dynamic weakening due to the thermal pressurization of pore fluids.
Standard rate-and-state fault models are capable of reproducing realistic stress drops as well as the observationally inferred magnitude invariance in stress drops. However, the breakdown energies depend on the rate-and-state characteristic slip L and increase only slightly with increasing event size for models with a given value of L, before saturating. Simulations with larger L lead to larger values of breakdown energies. However, this alone is not sufficient to match the observed trend, because the nucleation size increases with large L and the models with large L are no longer able to produce small events. This problem can potentially be resolved by using a nonconstant value for L, perhaps one that evolves with slip or slip rate. One can physically motivate this by imagining that the characteristic slip distance evolves as the fault slips and undergoes physical changes including damage on the fault in the form of gouge and off the fault in the form of cracking. These processes may alter the "effective" characteristic slip distance on the fault during the dynamic event. Evolving L during the event may serve as a proxy for these additional phenomena.
Our simulations show that fault models with enhanced dynamic weakening due to thermal pressurization can explain both the increasing trend in breakdown energy with increasing event size as well as the near magnitude invariance of average stress drops. The simulated breakdown energies G a match well the inferred trend for natural events, and our stress drops are consistent with seismologically inferred values in the 1-10 MPa range for all of our event sizes, excluding the complete rupture events in some models. We find that, with enhanced dynamic weakening, larger partial events result in lower average levels of prestress, due to their penetration into lower-prestressed regions. These events also weaken the fault more than smaller events do and arrest at lower levels of final stress. Our simulations reproduce this effect for events ranging several orders of magnitude in size (two orders of slip and approximately four orders of magnitude in moment).
The thermal pressurization parameters assumed in this work, motivated by values from Noda and Lapusta (2010), result in moderate additional dynamic weakening and crack-like ruptures. Given the assumed frictional properties, in order to maintain reasonable fault temperatures that avoid wholesale melting of the shearing layer for self-consistency of the models, such models do require the assumption of relatively low effective confining stress and hence substantial chronic fluid overpressurization throughout seismogenic depths. Such fluid overpressure may be present on fault, for example, some subduction megathrusts. Models with more substantial dynamic weakening, examined for single rupture events (e.g., Noda et al., 2009), do show that sufficiently enhanced weakening can lead to reasonable fault temperatures even with hydrostatic values of pore pressure; examining such models in terms of earthquake sequences is the subject of ongoing work. Such models would result in relatively sharp self-healing pulses (e.g., Noda et al., 2009), which have been advocated as prevailing rupture modes in some observational studies (Heaton, 1990). Other observational studies inferred broader pulse-like ruptures (e.g., Ye et al., 2016), which could be an observational equivalent of crack-like ruptures with weak tails. Our future work will examine whether models with self-healing pulses also reproduce a range of available observations.
We also find that the properties of the arresting velocity-strengthening regions have an impact on the average stress drop of events that significantly propagate into these regions. This is most important for our complete rupture events. Partial ruptures encounter low levels of prestress which inhibit their propagation and lead to their arrest within the VW region. Complete ruptures do not encounter unfavorable prestresses which would inhibit their propagation more and instead are held to a limited rupture domain by the VS regions, no matter their slip. Arresting regions with higher values of VS inhibit rupture propagation and lead to increasing stress drops as larger events slip more but are unable to increase in their spatial extent. Lower values of VS allow for significant propagation into the arresting regions and can lead to decreasing stress drops as the rupture area increases.