Earthquake breakdown energy scaling despite constant fracture energy

In the quest to determine fault weakening processes that govern earthquake mechanics, it is common to infer the earthquake breakdown energy from seismological measurements. Breakdown energy is observed to scale with slip, which is often attributed to enhanced fault weakening with continued slip or at high slip rates, possibly caused by flash heating and thermal pressurization. However, seismologically inferred breakdown energy varies by more than six orders of magnitude and is frequently found to be negative-valued. This casts doubts about the common interpretation that breakdown energy is a proxy for the fracture energy, a material property which must be positive-valued and is generally observed to be relatively scale independent. Here, we present a dynamic model that demonstrates that breakdown energy scaling can occur despite constant fracture energy and does not require thermal pressurization or other enhanced weakening. Instead, earthquake breakdown energy scaling occurs simply due to scale-invariant stress drop overshoot, which may be affected more directly by the overall rupture mode – crack-like or pulse-like – rather than from a specific slip-weakening relationship.

I n an earthquake, strain energy stored elastically in the Earth's crust is quickly transformed into radiated seismic waves, new fractures and wear products, and heat. This process is controlled, to a large extent, by the fault constitutive law, which describes how a fault's strength evolves with slip or slip rate and how much energy is dissipated or available for continued rupture. Thus, the fault constitutive law and associated fault properties are key to better understand and possibly predict earthquake mechanics. However, these properties are difficult to measure in the Earth. A number of studies tried to constrain fault constitutive behavior using seismological observations of earthquakes, and, in particular, the way that earthquake parameters scale from small to large [1][2][3][4] . One significant seismological observation, known as the breakdown energy, is thought to be related to the slip weakening process, and is often assumed to be a proxy for fracture energy. However, the breakdown energy is frequently observed to be negative-valued 1,4,5 and, if positivevalued, to scale with slip. Hence, larger earthquakes with more total slip appear to dissipate more energy (per unit rupture area) than small earthquakes.
If breakdown energy was equivalent to fracture energy, as often assumed (e.g., ref. [6][7][8], it is expected to be a material or interfacial property that is positive-valued and, for most engineering materials, observed to be only mildly dependent on the rupture propagation 8 . Hence, it is not expected to scale over many orders of magnitude as seismologically observed. Various theories have been developed over the years to explain the discrepancy of breakdown energy scaling. For instance, it was suggested that frictional weakening distance could increase with increasing earthquake size 1 . Other studies suggested additional mechanisms that activate as the earthquake rupture grows larger and slip distances increase such as thermal pressurization (e.g., ref. 3,4,9,10 ) or off-fault energy sinks (e.g., ref. [11][12][13][14]. High-velocity friction experiments that show continued weakening with cumulative slip were also offered as support. However, the fracture energy measured from those experiments appears to plateau at around 2 m of slip 7 while seismologically estimated breakdown energy of natural earthquakes continues beyond this limit and scales across all sizes of earthquakes 1 . Furthermore, none of these theories and observations provide an explanation for negative-valued breakdown energy. Hence, a fully consistent theory explaining all observations of seismologically inferred breakdown energy remains missing.
A running earthquake arrests either because the rupture front enters a region of the fault that is stronger or has a more stable rheology compared to the nucleation region, i.e., a barrier 9,10,15,16 , or the front enters a region with low initial stress and hence subcritical driving force 17 . While the first scenario implies a larger fault fracture energy for a larger earthquake rupture, the second scenario does not.
Here, we will show that the observed breakdown energy scaling does not require any complex or scale dependent fault constitutive law. We will demonstrate a simple but plausible scenario where constant fault friction but non-uniform initial stress results in the same scaling of breakdown energy. The non-uniform initial stress is not required to produce scaling of seismologically estimated breakdown energy 1 G 0 ; a strong barrier model 18,19 that results in stress overshoot (described below) would produce the same result. However, the non-uniform initial stress enables our models to have constant fracture energy over the entire domain. We solve the three-dimensional model with fully dynamic numerical simulations that employ linear slip-weakening friction. From the models, we determine the breakdown energy and other seismic source parameters. The result is unambiguous: breakdown energy does not correspond to the locally imposed constant fracture energy; it results from a scale-invariant stress overshoot.
Stress overshoot occurs in crack-like ruptures when a fastpropagating rupture front arrests but parts of the fault continue to slip 19,20 . The implications of overshoot on breakdown energy were discussed in Abercrombie and Rice 1 and dismissed for crack-like ruptures by Viesca and Garagash 4 . However, we argue that overshoot plays a central role in breakdown energy scaling and that our model results are consistent with current estimates of apparent stress and stress drop for real earthquakes.

Results
Numerical model. We consider a two-dimensional planar fault (y = 0) embedded in a three-dimensional isotropic elastic medium with shear modulus μ = 12 GPa and shear wave speed c s = 2126 m/s (Fig. 1a). The fault is governed by a linear slipweakening friction law ( Fig. 1c) with peak strength τ p = 80 MPa, residual strength τ r = 60MPa and critical slip distance δ c = 10 −5 m. Hence, the fault is characterized by realistic stress levels for natural faults, and a well-defined and constant fracture energy G = 100 Jm −2 , which is larger than estimates on smooth, precut faults 21 , but smaller than intact rocks 22 . The initial stress has uniform amplitude α = 67.5 MPa within a circular region centered at the origin. Earthquake ruptures are nucleated at the origin by a slowly increasing area of reduced fault strength (see Supplementary Note S1). While the nucleation process is artificial, it is sufficiently small and slow to not affect the unstable propagation of the earthquake that occurs spontaneously. After nucleation, the earthquake rupture velocity quickly accelerates to the Rayleigh wave speed where it remains constant 23,24 . The earthquake is nearly axisymmetric but grows slightly faster in the direction of the applied shear load (i.e., x-direction).
We explore two different scenarios for earthquake scaling by imposing an initial stress distribution that depends on scaling factor χ. Both scenarios reproduce the scaling of breakdown energy despite constant fracture energy, and produce identical initial rupture growth within a circular region of radius a = 0.3χ m, but they differ in the way they arrest, and this provides insight into the role that arrest plays in calculated source parameters. In scaling case A, initial stress decreases outside the circular region (Fig. 1b) at spatial rate β = 75/χ MPa/m (cf. Fig. 1e and Fig. 1f and note scaling of x-axis). In scaling case B, we impose a constant stress falloff β = 300 MPa/m (cf. Fig. 1g and Fig. 1h). The χ = 2 −2 model is identical in both scaling cases. We highlight that the resolution of the numerical models, the nucleation procedure, and fault friction properties, including fracture energy G (Fig. 1c), are kept constant across all scales.
Scaling of earthquake source properties. The scaling of calculated earthquake source properties is shown in Fig. 2, and definitions for these parameters in the context of our models are described below. Earthquake rupture area A is the area of the slipped fault patch, i.e., A = ∫ Σ dS for Σ = {x, z ∈ U|δ(x, z) > 0} and U is the entire simulation fault surface domain with y = 0. The spatially averaged slip distance D is the final slip δ f averaged over the rupture area, i.e., MðtÞdt. Alternatively, the seismic moment could be computed as M 0 = μAD, which yields equivalent results.
The averaged stress drop Δτ is one of the most important inferred earthquake source properties. A common seismological approach to determine Δτ is to assume a uniform stress drop 25 and estimate it from seismic source parameters using . Despite the simplification of assumed uniform stress drop, which is clearly not satisfied in our model (see Δτ in Fig. 1e-h), the above formulation provides an accurate estimation of Δτ, and more sophisticated approaches 26,27 do not lead to significantly different results (see Supplementary Note S2). We observe that Δτ is nearly scale-invariant in our model ( Fig. 2e), where scaling case A and B bound closely a scaleinvariant behavior, consistent with seismologically inferred measurements [28][29][30][31] . With the self-similar initial stress distribution of scaling case A, the arrest zone width 17 w az = R − a is a nearly constant percentage of the average rupture radius R (Fig. 2d). For scaling case A, initial stress decreases more gradually (smaller β) for larger earthquakes than for smaller ones, and A scales with χ 2.1 rather than the expected χ 2.0 because large earthquakes rupture proportionally further into unfavorably stressed regions. For this reason M 0 scales with χ 3.1 instead of χ 3.0 . However, the average slip D scales as χ 1.0 , therefore stress drop ($ D= ffiffiffi ffi A p ) decreases slightly (Fig. 2e). In scaling case B, w az / 2R is smaller for larger earthquakes (Fig. 2d) and A scales as χ 1.9 (rather than χ 2.0 ) because larger and larger ruptures meet relatively steeper and steeper stress gradients and therefore rupture proportionally shorter and shorter distances into unfavorable stress. Similarly, M 0 scales with χ 2.9 and stress drop increases slightly with increasing earthquake size.
We also analyze the energy budget for our model earthquakes. Using the energy-considered averaged initial stress τ i and final stress τ f 26,27 , the averaged total released strain energy is which is a simplified formulation of the method proposed by Kostrov 32 (see Fig. 2g and Methods). We also compute the dissipated energy E D = E H + E G , which combines heat E H and total breakdown energy E G , by 2h). Finally, we estimate the radiated energy E R (see Methods and Supplementary Note S3), and we observe E R /A ∝ D 1 , consistent with a nearly constant apparent stress τ a = μE R /M 0 ≈ 2MPa (Fig. 2f, i; ref. 5,28,33,34 ).
From our scaling analysis, we conclude that our model reproduces the scaling behavior of all important seismic source properties. The expected scaling properties for self-similar rupture are either reproduced by the models or are bounded by scaling cases A and B (Fig. 2a-d). We specifically note that A ∝ D 2 , M 0 ∝ D 3 , and Δτ / D 0 , consistent with standard earthquake scaling 35 .
Scaling of breakdown energy. The dissipated energy E D is composed of the total breakdown energy E G and heat E H (Fig. 3a). The averaged breakdown energy G % E G =A is often interpreted as a proxy for the fracture energy G, which controls the dynamics of the fault rupture 23 , and hence is of great importance for seismology. Abercrombie and Rice 1 proposed a seismological approach to estimate G by which has been widely used (e.g., ref. 5,[36][37][38][39]. We compute G 0 for our simulations following Eq. (1) and observe that our scaling cases A and B follow G 0 / D 0:8 and G 0 / D 1:0 , respectively (Fig. 3c). Abercrombie and Rice 1 found G 0 / D 1:28 . Their exponent of 1.28 is greater than 1 as a result of magnitude-dependent   τ a and Δτ, which they argued was significant, but has not been supported by more recent studies 5,28,33,34 . While the scaling in our models is not an exact match, which could be attributed to additional aforementioned dissipative mechanisms, we argue that it is within the uncertainty of the observational data (Fig. S11), and emphasize that G 0 scales orders of magnitude in our simulations even though the actual fault property of fracture energy G is kept constant (Fig. 3). Using energy calculations proposed by previous studies 26, 27 , we rewrite G 0 as where Δτ OS ¼ τ r À τ f À Á is the energy-considered averaged stress overshoot (see Fig. 3a and Methods). Hence, the breakdown energy consists of the sum of the fracture energy and the overshoot energy, which we define as In our simulations, the averaged stress overshoot Δτ OS ¼ 1 À 2 MPa is positive and nearly scale-invariant (see Fig. S6 and Fig. S9). G 0 of large earthquakes with small G and scale invariant overshoot is dominated by the overshoot, hence: G 0 / D 1 (Eq. (2)). In Fig. 3c, we show the results of Eq. (2) for different levels of overshoot (Δτ OS ) compared to seismological estimates of G 0 for earthquakes of various sizes, and find Δτ OS % 1 MPa. Further, G serves as a lower bound of G 0 , as shown by the dotted curves at lower D in Fig. 3c. This implies the fault fracture energy is bounded by G ≤ 100 Jm −2 . Finally, we note that a model with a strong barrier 18,19 rather than an initial stress taper would present a similar scaling of G 0 due to a scale-invariant stress overshoot. Such a condition is most similar to scaling case B, where arrest occurs very abruptly, and hence it provides an upper bound for stress overshoot, which is Δτ OS ¼ 0:2Δτ (see Fig. S9). Abercrombie and Rice 1 considered whether G 0 scaling could be dominated by overshoot, as we suggest, but ultimately dismissed it because their observed τ a =Δτ ¼ 0:1 would then imply an overshoot of 0:4Δτ. Such a large overshoot would exceed the 0:15Δτ to 0:2Δτ of the Madariaga 19 model, which was considered an upper bound, due to its abrupt arrest. However, τ a =Δτ ¼ 0:1 is somewhat lower than found elsewhere 40 (τ a =Δτ ¼ 0:33). Our models produce τ a =Δτ ¼ 0:3 À 0:4, which is within the range of seismic estimates, and our observed overshoot (0:1Δτ to 0:2Δτ) is consistent with both the Madariaga 19 upper bound and the energy balance relation (ref. 20 Eq. 1c) used by Abercrombie and Rice 1 .
If our interpretation is correct, and G is indeed small (G ≤ 100 Jm −2 ), how can this be reconciled with estimates of G = 10 6 Jm −2 derived from finite-fault kinematic inversions (e.g., ref. 2 ), or with laboratory data 13 that indicates 1 m weakening distances? As for the kinematic inversions, we find that the minimum resolvable characteristic weakening lengthδ c could be limited by bandwidth 41 . For example, assuming a typical 0.5 s smoothing operator (e.g., ref. 2 ), minimum resolvable δ c is of order 500 mm. Assuming τ p − τ r = 10 MPa, this places minimum resolvable G = 2.5 MJ/m 2 . (See Supplementary Note S4, Fig. S4.) Considering the laboratory results, most of the observed weakening occurs while fault slip accelerates 13 , and 1 m weakening distances resulted from sluggish slip acceleration (6.5 m/s 2 ) compared with measurements 42 from dynamic rupture fronts (>20,000 m/s 2 ). Experiments 43 that imposed more abrupt loading (30 m/s 2 ) exhibited far smaller weakening distances (0.02 m). Thus, large inferred weakening distances and therefore large G may be an artifact of laboratory loading procedures that are slow compared to the fault acceleration imposed by a dynamic rupture front during an earthquake.
Finally, our interpretation helps reconcile observations of G 0 % 0 and negative G 0 , found for numerous earthquakes 1,4 , yet rarely discussed. Such cases would result from negligible overshoot or undershoot i.e., τ f > τ r , a condition typically associated with pulse-like earthquake ruptures 44 , where the fault is elastically reloaded after a short period of slip. Indeed, we find that a simple model with scale-invariant random overshoot that ranges between −1 and 2 MPa and negligibly small G produces G 0 / D 1:0 , fits the Abercrombie and Rice 1 data reasonably well (see Fig. S11 and dotted curves in Fig. 3c), and produces a catalog where one-third of all events exhibit G 0 ≤ 0, similar to the large earthquakes studied by Viesca and Garagash 4 . All in all, our simulations and proposed interpretation of G 0 reconciles the coexistence of breakdown energy scaling and negative G 0 .
Rupture growth and arrest from source time functions. We study the spontaneous growth and arrest of rupture in our dynamic model by comparing the earthquake source time function or moment rate function _ MðtÞ with natural observations. _ MðtÞ consists of three phases: (1) a self-similar growth phase (2) a divergence from self-similar growth near peak moment rate and (3) the post-peak decay. In the growth phase, we observe _ MðtÞ / t 2 , at all scales (see Fig. 4a), as expected for roughly circular earthquake ruptures propagating without bound 45,46 . Selfsimilar moment rate growth has also been observed for large earthquakesin some cases _ MðtÞ / t 2 ( Fig. 4c; ref. 47

) but in others _
MðtÞ / t 1 (Fig. 4b; ref. 48 ). Near peak _ MðtÞ, earthquake rupture begins to arrest and the boundaries to rupture growth are increasingly felt. Scaling case B produces earthquakes that arrest far too quickly compared to observations 47,48 . Scaling case A's gradual arrest is a reasonable fit to observations, however, all of our models exhibit rapid post-peak decay of _ MðtÞ and negative skew, while natural earthquakes show a more gradual post-peak behavior and slightly positive skew 47,48 . This suggests that our simulated earthquakes do not arrest slowly enough or that slip ceases too quickly after the rupture initially begins to arrest.
To better understand the above discrepancies, we compared the previously described simulations, where nucleation occurs in the center of the circular region of favorable initial stress, to a case where the earthquake nucleates close to the edge of the favorably stressed region, i.e., (x, y, z) = (0.7a, 0, 0.7a), denoted "edge nucleation" (see Fig. S10). The latter case quickly reaches unfavorable initial stress on NE side and then ruptures unilaterally to the SW. The difference in the source time function is striking (Fig. 4d). The edge nucleation quickly diverges from the _ MðtÞ / t 2 self-similar growth curve and instead _ MðtÞ grows nearly linearly. The growth phase is about 50 percent longer than in the symmetric case, though the decay is very similar.
While the edge nucleation model does not reconcile the abbreviated post-peak behavior (and actually makes the skewness worse), it offers a satisfying explanation for linear growth rate. The earliest part of the growth phase, when unbounded growth is expected, is difficult to resolve from kinematic models 48 . Proper resolution of early growth would require fault plane discretization size to depend on the distance from the hypocenter 49 . When normalized source time functions are plotted on a linear scale, as shown in Fig. 4b (or ref. 48 Fig. 2b-d), their form is dominated by the final increase in moment rate just before rupture arrest exceeds rupture growth. In this stage, rupture growth will, on average, be bounded on at least one side, and will produce the nearly linear increase in moment rate demonstrated by our edge nucleation model.

Discussion
Previous work assumed that breakdown energy G 0 was a proxy for fracture energy G, and was therefore dominated by the way strength evolved with slip. Our modeling offers an alternative interpretation. It demonstrates that scaling of G 0 can arise naturally from the dynamics of rupture growth and arrest, and it can occur even with constant G and does not require a scaledependent constitutive behavior such as enhanced weakening at large fault slip. As noted by Abercrombie and Rice 1 , G 0 depends strongly on overshoot and undershoot, which are affected by the rupture type (crack versus pulse) and the manner in which an earthquake arrests. Our observation of G 0 scaling is directly linked to the fact that our model exhibits scale-independent stress overshoot. This also suggests that negative G 0 is an indication of undershoot 44 , which is thought to occur for self-healing slip pulses, while the majority of earthquakes overshoot 20 , a property of crack-like ruptures. However, our models' overshoot-induced scaling (G 0 / D 0:8 and G 0 / D 1:0 ) is somewhat weaker than the G 0 / D 1:28 observed by Abercrombie and Rice 1 . Our model produces a self-similar source time function and offers an explanation for the seismologically observed linear growth phase 48 , but decays too quickly compared to natural observations.

Methods
Parametrization of initial stress distribution. The on-fault initial stress distribution τ i (x, y = 0, z) is applied in the x-direction and is parameterized as where r ¼ ffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi is the distance to the origin, H is the Heaviside step function, a is the radius of the stress plateau with amplitude α, and β is the spatial-rate of stress decay outside the stress plateau (Fig. 1b).
Numerical method. Numerical simulations are conducted with our implementation of the spectral boundary integral (SBI) method 50 to solve for a fully dynamic interface debonding problem. The method is based on previous methodological studies [51][52][53]  We verified that the simulated domain is large enough that the rupture and the reflected waves do not affect the ruptured region within the simulation duration. Ruptures were nucleated by a seed crack at the center of the fault, in which the peak strength is manually decreased to τ r and extended at 10% of the Rayleigh wave speed (see Supplementary Text S1).
Seismic source parameters. M 0 calculated through integration over _ MðtÞ ¼ μ R U _ δðx; z; tÞdS is theoretically identical to M 0 = μAD. We have verified that their results are the same as the amplitude at the low-frequency plateau of the moment rate spectrum where F denotes Fourier transformation.
Energy Calculations. The total strain energy released from an arrested earthquake rupture can be integrated solely on the ruptured domain Σ 32 : which can be rewritten into a simpler form with the energy-considered averaging method 26,27 : where τ E is the average stress weighted by the final slip Therefore, the energy-considered averaged stress drop is simply where τ E i ¼ τ E ðt ¼ 0Þ and τ E f ¼ τ E ðt ¼ t end Þ. We consider two approaches for radiated energy. First, we compute E N R from integrated on-fault stress and slip evolution. As shown by ref. 32 and Appendix C in ref. 55 , E R can be conveniently evaluated numerically from the result of numerical MðtÞ of numerical simulations and natural earthquakes. a All models follow an identical _ MðtÞ / t 2 growth independent of scaling parameter χ until reaching unfavorable stress conditions. The black dotted curve is a quadratic fit for growth phase of moment rate history. b Normalized source time functions, where Mdt and the area under _ MðtÞ is 1. Gray curves are natural earthquakes with M w ≥ 7 presented in ref. 48 . c Scaled source time functions with matching growth phase. Green and red curves are natural earthquakes at various source depths H adapted from ref. 47 . d Source time function of a rupture nucleated at the center of the stress plateau, and at the edge of the stress plateau. simulation through This approach is equivalent to an energy conservation equation (following ref. 1 ) if the τ i term in the second integration of Eq. (M7) is reorganized into Z t end 0 Z U τ i ðx; zÞ _ δðx; z; tÞdSdt ¼ and moved into the first term. Both approaches yield similar results with some differences for small χ, which we associate to small ruptures arresting before reaching Rayleigh wave speed and the effects of rupture initiation process in our numerical simulation. We use E N R to represent E R due to its superior numerical stability (see Supplementary Text S3).
Here we start from the energy conservation equation divided by A on both sides, where ΔW is the total strain energy released, A is rupture area, G 0 is the spatially averaged breakdown energy, E H is heat, and E R is radiated energy. Replacing With some rearrangements and replacing τ E i À τ E f ¼ Δτ E , G 0 can be expressed as which takes a very similar form as G 0 (Eq. (1)). The deduction above actually assumes that G 0 is Whereas the fracture energy should be the area above τ r , Similar to ref. 1 , when assuming G and τ r are constants, G 0 can be expressed as Note that the difference between G 0 and G involves the final-slip-weighted-average final stress τ E f , i.e., and the stress overshoot τ r À τ f ðx; zÞ À Á in our models appear to be larger at locations with larger slip, as shown in Supplementary Fig. S7. This highlights the spatially averaged stress overshoot cannot be used when evaluating the accuracy of G 0 , as the effect of stress overshoot is clearly amplified by its correlation with larger slip.

Data availability
The simulation data generated in this study have been deposited in the ETH Research Collection database under accession code ethz-b-000527677 [https://doi.org/10.3929/ ethz-b-000527677].