# Effects of failure parameterization on pre- and co-seismic earthquake rupture

Ekaterina Bolotskaya, Bradford Hager, Jean-Paul Ampuero

► To cite this version:

Ekaterina Bolotskaya, Bradford Hager, Jean-Paul Ampuero. Effects of failure parameterization on pre-and co-seismic earthquake rupture. Geophysical Journal International, 2025, 240 (3), pp.1538-1550. 10.1093/gji/ggae462. hal-04885007

HAL Id: hal-04885007

https://hal.science/hal-04885007v1

Submitted on 4 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Effects of failure parametrization on pre- and coseismic earthquake rupture

Ekaterina Bolotskaya<fnref n="1" />, <fnref n="1" /> Bradford H. Hager<fnref n="1" /> and Jean-Paul Ampuero <fnref n="2" />

Accepted 2024 December 24. Received 2024 December 23; in original form 2024 August 16

## SUMMARY

This study investigates the impact of different slip-dependent fault failure parametrizations on the pre-seismic and coseismic phases of the earthquake cycle using 2-D finite element simulations. Various failure laws are considered: linear slip-weakening, double slip-weakening, parabolic cohesive zone and exponential cohesive zone laws. The pre-seismic phase is modelled using a quasi-static approach, while the coseismic phase is modelled dynamically. Results demonstrate the importance of accounting for the different shapes of failure laws during the pre-seismic phase, as they lead to variations in pre-seismic phase duration, associated pre-seismic slip and lateral extent of slip front. Failure laws with gentler initial slopes require more time to initiate dynamic rupture compared to those with steeper initial weakening slopes. The presence of an initial strengthening segment in the failure law affects the amount and lateral extent of pre-seismic slip. During dynamic rupture propagation, the specific details of the failure law become less significant due to the Lorentz contraction of the process zone. However, variations in the duration of the pre-seismic phases result in different dynamic stress drops and, consequently, variations in rupture acceleration times and earthquake magnitudes for different failure laws.

**Key words:** Numerical modelling; Earthquake dynamics; Earthquake parametrization; Rheology and friction of fault zones.

## 1 INTRODUCTION

Different physical processes associated with material failure affect earthquake rupture nucleation and propagation. These processes include frictional sliding along faults (Brace & Byerlee 1966; Scholz 1998), fracturing of intact rock sections connecting pre-existing planes of weakness (fault step-overs) or cemented fault segments (Aki 1979, 1984; Ohnaka 2003), fault branching and “wing-crack” formation (Archuleta 1984; Aochi & Fukuyama 2002) and inelastic deformation in the fault damage zone (Rudnicki 1980; Wilson et al. 2003), among others.

In earthquake modelling, it is necessary to combine all these processes into simplified mathematical representations known as failure laws. Failure laws typically relate the evolution of shear strength of the fault to slip, slip rate and/or fault state variables. The fracture energy or breakdown work, represented by the area under the shear strength versus slip curve and above the residual shear strength, is a crucial characteristic of the failure law and significantly influences earthquake rupture. The most commonly used failure laws for earthquake modelling are linear slip-weakening friction (hereafter referred to as linear slip-weakening or slip-weakening failure law) (Barenblatt 1959; Ida 1972; Palmer & Rice 1973), where shear

resistance linearly decreases with slip, and rate-and-state friction (Dieterich 1979; Ruina 1983), which is based on laboratory observations of rock friction at relatively slow slip velocities and can exhibit stick-slip behaviour, fault healing and memory effects (Rice & Tse 1986; Marone 1998).

While rate-and-state friction laws are appropriate for modelling mature faults characterized by frictional sliding, these laws are not expected to provide an adequate description of the onset of slip across consolidated interfaces. Slip-dependent laws are thought to be better suited for describing the initial loss of cohesion on faults locked for long periods and healed through hydrochemical processes like cementation (Brantut & Viesca 2015). Although there is a scarcity of laboratory experiments focused on cohesion loss, this process has been previously modelled as slip-dependent rather than rate-dependent (Brantut & Viesca 2015). Thus, slip-dependent laws are appropriate to represent the mechanical behaviour of healed faults during earthquake nucleation.

In this study, we consider several different slip-dependent failure parametrizations (different classes of slip-dependent strength curve shapes) and compare the results of earthquake rupture simulations obtained using those. We consider linear slip-weakening failure law as a reference due to its widespread use and simplicity.

<footnote n="1">¹Massachusetts Institute of Technology, Postal code: 02139, Cambridge, USA. E-mail: bolee@mit.edu</footnote>

<footnote n="2">²Université Côte d'Azur, Observatoire de la Côte d'Azur, IRD, CNRS, Geoazur. Postal code: 06560, Valbonne, France.</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Model geometry (fault length not to scale). Dashed line (blue) is the fault and solid line (red) is the loading patch. $\sigma_n$ is the normal stress on the fault and $\tau_0$ is the initial shear stress on the fault.</figcaption></figure>

Additionally, we investigate double-slip-weakening failure law, inspired by recent experimental results on frictional sliding in polymethyl methacrylate (PMMA) (Paglialunga et al. 2022), and two failure laws inspired by the fracturing of intact rocks: parabolic and exponential cohesive zone laws (Ohnaka et al. 1997; Alfano & Crisfield 2001). While we acknowledge that the various physical processes associated with fault failure may correspond to different magnitudes of fracture energy, it is challenging to accurately constrain this parameter. Therefore, in this study, we set the fracture energy to be the same for all failure laws. By doing so, we aim to explore the variations arising from the different shapes of the failure laws during various stages of earthquake rupture propagation.

In summary, we utilize the 2-D finite element method (FEM) (Section 2) to conduct earthquake simulations in a linear elastic material with different failure laws prescribed along the fault: slip-weakening failure law, double-slip-weakening failure law, parabolic cohesive zone law, and exponential cohesive zone law (Section 3). Quasi-static simulations are used for the pre-seismic phase of the earthquake cycle, while dynamic simulations are used for earthquake rupture propagation. The slip and stress values obtained from the quasi-static portion of the simulation serve as inputs for the dynamic phase. From the quasi-static models, we find that different failure laws play a significant role in the duration of the pre-seismic phase of the earthquake cycle, the amount of accumulated pre-seismic slip and the lateral extent of slip propagation (Section 4.1). The varying durations of the pre-seismic phases for the different

failure laws, coupled with a constant background stressing rate, lead to different dynamic stress drops during the dynamic phases of the earthquakes. This discrepancy results in variations in rupture acceleration times and earthquake magnitudes for simulations with different failure laws (Section 4.2). We compare our 2-D results with those obtained using 1-D spring-slider models (Section 5). We discuss the implications of our findings for earthquake rupture modelling in Section 6.1 and offer a potential mechanism for slow slip preceding earthquakes in Section 6.2. Finally, we present our conclusions in Section 7.

## 2 MODEL

We conduct numerical simulations using the finite element method implemented in the PyLith software (Aagaard et al. 2013). Our study focuses on two fundamental aspects of earthquake behaviour: the slow propagation of slip into the seismogenic portion of the fault (known as the pre-seismic phase of the earthquake cycle) and the rapid rupture process during the earthquake (referred to as the coseismic phase). The former is modelled using a quasi-static approach without radiation damping, neglecting the inertial term in the momentum conservation equation, considering background stressing rate over a timescale of months. The latter is modelled fully dynamically, including the inertial term, considering timescales on the order of milliseconds. The quasi-static to dynamic transition criterion has an objective, quantitative basis and a physical meaning: under a slowly increasing tectonic load, the quasi-static simulation ends when a static equilibrium solution no longer exists. When this solution ceases to exist, the quasi-static solver fails to converge. The slip and stress values obtained from the last converged time step of the quasi-static solver are used as inputs for the dynamic solver.

Our model consists of a horizontal fault embedded within a 2-D rectangular domain (Fig. 1). The fault does not span the full domain width, preventing wave reflections from the boundaries from affecting the dynamic rupture. The bulk rheology of the material is assumed to be homogeneous, isotropic and linear elastic under 2-D plane strain conditions. Material parameters and model geometry are presented in Table 1 and Fig. 1, respectively.

We impose uniform initial stresses throughout the domain (both $\sigma_{xx}$ and $\sigma_{yy}$ are present), resulting in uniform initial fault shear stress $\tau_0$ and fault normal stress $\sigma_n$. In the quasi-static portion of the simulation, the bottom boundary is fixed, a uniform and constant shear stressing rate is applied along the top boundary and

<table><caption>**Table 1.** Model parameters.</caption><thead><tr><th></th><th>Geometry</th><th></th></tr></thead><tbody><tr><td rowspan="3">L</td><td>Total fault length</td><td>12 m</td></tr><tr><td>Loading patch length</td><td>1 m</td></tr><tr><td>Domain size</td><td>20 m × 30 m</td></tr><tr><td></td><td><strong>Bulk properties</strong></td><td></td></tr><tr><td>v<sub>S</sub></td><td>S-wave speed</td><td>3000 m s<sup>-1</sup></td></tr><tr><td>v<sub>P</sub></td><td>P-wave speed</td><td>4500 m s<sup>-1</sup></td></tr><tr><td>ρ</td><td>Density</td><td>2400 kg m<sup>-3</sup></td></tr><tr><td>G</td><td>Shear modulus</td><td>21.6 GPa</td></tr><tr><td>E</td><td>Young's modulus</td><td>47.52 GPa</td></tr><tr><td>ν</td><td>Poisson's ratio</td><td>0.1</td></tr><tr><td></td><td><strong>Stress state</strong></td><td></td></tr><tr><td>τ<sub>0</sub></td><td>Fault initial shear stress</td><td>44.95 MPa</td></tr><tr><td>σ<sub>n</sub></td><td>Fault normal stress</td><td>75 MPa</td></tr><tr><td>σ&#773;<sub>xy</sub></td><td>Shear stressing rate</td><td>0.5 MPa per month</td></tr></tbody></table>
<table><caption>**Table 2.** Failure law parameters.</caption><thead><tr><th></th><th>SW</th><th>DSW</th><th>PCZ</th><th>ECZ</th></tr></thead><tbody><tr><td>μ<sub>s</sub> | τ<sub>s</sub></td><td>0.650 | 48.75 MPa</td><td>0.650 | 48.75 MPa</td><td>0.650 | 48.75 MPa</td><td>0.650 | 48.75 MPa</td></tr><tr><td>μ<sub>t</sub> | τ<sub>t</sub></td><td></td><td>0.630 | 47.25 MPa</td><td></td><td></td></tr><tr><td>μ<sub>d</sub> | τ<sub>d</sub></td><td>0.600 | 45.00 MPa</td><td>0.600 | 45.00 MPa</td><td>0.600 | 45.00 MPa</td><td>0.600 | 45.00 MPa</td></tr><tr><td>μ<sub>i</sub> | τ<sub>i</sub></td><td></td><td></td><td>0.641 | 48.07 MPa</td><td>0.641 | 45.36 MPa</td></tr><tr><td>D<sub>c</sub></td><td>120.0 μm</td><td></td><td></td><td></td></tr><tr><td>D<sub>w1</sub></td><td></td><td>24.0 μm</td><td></td><td></td></tr><tr><td>D<sub>w2</sub></td><td></td><td>136.0 μm</td><td></td><td></td></tr><tr><td>D<sub>1</sub></td><td></td><td></td><td></td><td>0.8 μm</td></tr><tr><td>D<sub>2</sub></td><td></td><td></td><td></td><td>22.1 μm</td></tr><tr><td>D<sub>s</sub></td><td></td><td></td><td>24.0 μm</td><td>21.3 μm</td></tr><tr><td>D<sub>w</sub></td><td></td><td></td><td>56.2 μm</td><td>∞ μm</td></tr><tr><td>G<sub>c</sub></td><td>225 Pa m</td><td>225 Pa m</td><td>225 Pa m</td><td>225 Pa m</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>**Figure 2.** Shear strength as a function of slip for the four failure laws.</figcaption></figure>

a uniform and constant normal stress is applied along the side boundaries to preserve the initially applied $\sigma_{xx}$. In the dynamic simulation, all boundaries are absorbing, and the internal stress state is carried over from the final time step of the quasi-static simulation.

A loading patch with constant friction coefficient, $\mu_d$ (Table 2), is prescribed on the left end of the fault; it serves a purpose similar to a velocity-strengthening region in other earthquake rupture models (e.g. Liu & Rice 2007; Rubin 2008): it slides stably and its induced stress loads the locked portion of the fault, until an earthquake is triggered.

The mesh is quadrilateral, with element size increasing from the fault toward the domain boundaries. The element size along the fault, $\Delta x$, is 1 cm for quasi-static simulations and 1 mm for dynamic simulations. The static process zone size for mode II rupture with slip-weakening (SW) failure law is $\Lambda_0 = \frac{9\pi G D_c}{32 (1-\nu)(\tau_s - \tau_d)}$ (Day et al. 2005). $\Lambda_0 = 0.68$ m in our case, with 680 elements within the process zone at the onset of the dynamic simulation, ensuring sufficient resolution. The time step is 0.003 months for quasi-static simulations and 0.1 µs for dynamic simulations. The dynamic time step satisfies the Courant–Friedrichs–Lewy stability condition: CFL = $\frac{\nu_s \Delta t}{\Delta x} = 0.3$.

## 3 FAILURE LAWS

We consider the following four failure parametrizations, each providing a different representation of the slip-dependent evolution of fault strength during an earthquake:

(i) SW: linear slip-weakening failure law, which is the most common simplification of other more complex failure laws (Barenblatt

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Shear strength versus slip schematics for the different failure parametrizations. $D_s$ and $D_w$ for ECZ and PCZ denote strengthening and weakening distances, respectively. For ECZ these are a function of $D_1$ and $D_2$ given in Table 2. More details are provided in Supplementary Section S1 (Supporting Information).</figcaption></figure>

1959; Ida 1972; Palmer & Rice 1973; Andrews 1976b, a) and is frequently used in earthquake modelling.

(ii) DSW: double slip-weakening model, featuring two linear weakening portions with different slip-weakening rate (slope). Bilinear and polylinear failure laws have been used in the literature to model fracturing of cohesive surfaces (Hillerborg et al. 1976; Mi et al. 1998; Alfano & Crisfield 2001; Alfano 2006).

(iii) PCZ: parabolic cohesive zone, which (along with linear-parabolic variations of it) is frequently used to model surface debonding and fracturing processes (Allix et al. 1995; Allix & Corigliano 1996; Alfano 2006).

(iv) ECZ: exponential cohesive zone, which is one of many descriptions of fracture process, inspired by experiments on fracturing of cohesive surfaces (Needleman 1990; Ohnaka 1996; Chandra et al. 2002; Alfano 2006; Park & Paulino 2011).

The governing equations of these failure parametrizations are provided in Supplementary Section S1 (Supporting Information) and their shear strength versus slip curves are shown in Fig. 2.
<figure><img src="image_5.png" /><figcaption>**Figure 4.** Cumulative fault slip in mm (a) and shear stress in MPa (b) in a simulation with SW failure law. Each curve corresponds to a different time, indicated by the parula colour scale during the quasi-static phase (in months) and by the plasma colour scale during the dynamic phase (in ms). Vertical dotted line shows the boundary of the loading patch on the left. Red box in panel (a) is the zoomed in area shown in Fig. 5.</figcaption></figure>

The area under the shear strength versus slip curve and above the residual strength represents the fracture energy, $G_c$, of the failure law (Fig. 3). In this study, we are interested in the effects of the different failure law shapes on earthquake properties, thus, we maintain the same fracture energy and strength drop ($\tau_s - \tau_d$) across all four failure parametrizations. We use the concept of fracture energy only in the context of dynamic ruptures. Since pre-seismic slip is minimal or non-existent for most points along the fault, both strengthening and weakening occur dynamically, making it more appropriate to compute fracture energy cumulatively under the entire strength curve (Fig. 3). The parameter values for each law are provided in Table 2. The selected values align with laboratory data from shear experiments on cemented fault surfaces and intact rocks (Karner et al. 1997; Ohnaka et al. 1997; Ohnaka 2013).

As pointed out by Xu et al. (2019), some features of the failure laws we use—such as the initial strengthening and the long weakening tail—could reflect an artifact due to off-fault rather than on-fault measurements. We acknowledge this limitation but maintain that comparing these four failure laws is useful, as all are applicable to earthquake modelling and understanding their differences is important.

# 4 RESULTS

Fig. 4 displays the simulation results for the slip-weakening (SW) failure law on a 12-m long fault. Panel (a) presents the cumulative fault slip at different times, while panel (b) illustrates the distribution of shear stress at different times, obtained from both the quasi-static simulation (with time represented by parula colours) and the dynamic simulation (with time represented by plasma colours).

During the pre-seismic phase, the shear stress far from the loading patch increases linearly with time, while the slip front slowly propagates from the edge of the loading patch into the seismogenic portion of the fault. The duration of the pre-seismic phase is approximately

six months, and the contour spacing is 0.39 months. In the dynamic phase, we observe unilateral crack-like rupture propagation along the fault, that is, with non-zero slip rates behind the rupture front, reaching the end of the fault. The duration of the dynamic phase is 9 ms, and the contour spacing is 0.15 ms. In the subsequent sections, we provide a detailed analysis of both portions of the simulation, discussing the results for the different failure parametrizations under consideration.

## 4.1 Pre-seismic phase

In the quasi-static portion of the simulation, we initiate the loading with a shear stress below the threshold required for sliding even along the loading patch ($\tau_0 < \tau_d$). The block is then continuously loaded with a constant shear stressing rate until the loading patch starts sliding. The loading patch undergoes stable sliding, accumulating slip and generating a stress concentration of increasing magnitude near its tip (Fig. 4b—parula contours). Once the magnitude of the stress concentration surpasses $\tau_s$, slip begins to propagate into the originally locked section of the fault.

Fig. 5(a) provides a zoomed-in view of Fig. 4(a), highlighting the pre-seismic slip evolution. Correspondingly, Figs 5(b)–(d) offer zoomed-in views for the other failure laws: DSW, PCZ and ECZ. There are differences in the slip contours depending on the failure law used. PCZ (panel c) and ECZ (panel d) exhibit greater lateral extent of the pre-seismic slip front compared to SW (panel a) and DSW (panel b), as seen by comparing the dark blue contours that correspond to the last converged time step of the quasi-static simulations. This difference is attributed to the presence of initial strengthening in the PCZ and ECZ failure laws (Fig. 2). ECZ demonstrates the most pronounced effect of initial strengthening: pre-seismic sliding extends over the entire fault because the background shear stress eventually exceeds the initial strength of the failure law ($\tau_i$ in Fig. 3).
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Zoom in of Fig. 4 for fault slip (in mm) focusing on the quasi-static portion of the simulation for the different failure laws (a–d). Panel (d) shows both the area near the loading patch and the rightmost end of the fault to illustrate pre-seismic slip front propagation to the end of the fault.</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Figure 6.** Mean slip rate in mm per month (a), mean slip in mm (b), and slip front location in m (c) as a function of time in months along the initially locked section of the fault. Vertical dash-dotted lines indicate the end of the pre-seismic phase for the different failure laws.</figcaption></figure>

Fig. 6 displays the mean slip rate (a), mean slip (b), and slip front location (c) along the initially locked section of the fault during the pre-seismic phase. Notably, the ECZ failure law (purple lines) exhibits the highest slip rate and slip magnitudes, with an average pre-seismic slip of approximately 1 per cent of the eventual average coseismic slip. Panel (c) shows the lateral extent of

the quasi-static slip front and highlights that, for ECZ, the pre-seismic slip front propagates across the entire fault. The abrupt step in mean slip velocity observed for ECZ at approximately 1.7 months (the purple line in panel a) corresponds to the time when the pre-seismic slip front propagation along the whole fault occurs.
<figure><img src="image_8.png" /><figcaption>**Figure 7.** Coseismic fault slip rate in m s⁻¹ (colour) versus time and position along the fault for the four failure laws (a) to (d). Solid red lines show the Rayleigh wave speed, v<sub>R</sub>. Dashed red lines show the S-wave speed, v<sub>S</sub>.</figcaption></figure>

The maximal extent of quasi-static slip propagation into the initially locked portion of the fault can be considered the “nucleation length”. Uenishi & Rice (2003) showed that the nucleation length in the case of a linear slip-weakening failure law in an unbounded fault is

$$h_n \approx 1.158 \frac{\mu^* D_c}{\tau_s - \tau_d}, \qquad (1)$$

where $\mu^* = \frac{G}{1-\nu}$ in our case of mode II rupture. With our model and failure law parameters, the theoretical nucleation length is $h_n = 0.89$ m. The SW nucleation length obtained in our simulations is shorter, 0.47 m. However, this difference is expected due to the different boundary conditions: the presence of a loading patch with accumulated slip next to the nucleation area in our simulations should reduce the nucleation length. This expectation is consistent with our modelling results and is quantitatively validated in Appendix A by comparison to the results of a theoretical analysis analogous to the eigenvalue analysis of Uenishi & Rice (2003).

Fig. 6 also demonstrates the differences in the duration of the pre-seismic phase (vertical dash-dotted lines in all panels), with the DSW failure law triggering an earthquake in the shortest time and the PCZ failure law taking the longest. The pre-seismic phase durations between PCZ and DSW are different by almost one month. This behaviour can be attributed to the initial weakening slope of the two failure laws: the steeper slope of DSW causes it to become unstable faster, while the gentle initial weakening slope of PCZ, coupled with the presence of strengthening (Fig. 2), stabilizes it and delays the occurrence of the earthquake. These findings align with pre-seismic behaviours observed in spring-slider models (Section 5).

An important implication of the different durations of the pre-seismic phases is the variation in the background stress reached at

the onset of the coseismic phase. Since we are applying a constant stress rate in the domain, the failure laws that remain stable for a longer duration lead to a larger coseismic stress drop, such as approximately 77 per cent of $\tau_s - \tau_d$ for PCZ and 65 per cent of $\tau_s - \tau_d$ for DSW.

## 4.2 Coseismic phase

As previously mentioned, the stress and slip results from the last converged time step of the quasi-static simulations serve as the initial conditions for the dynamic simulations. Thus, the coseismic stress drop varies among the different failure laws due to the differing durations of the pre-seismic phases.

Fig. 7 presents the fault slip rate as a function of time and position along the fault for the dynamic portion of the simulation for the four failure laws. Unilateral crack-like ruptures are observed, accelerating and propagating to the end of the fault. Distinct differences are observed among the four failure laws during the onset of the rupture, including differences in the rupture acceleration time (the time taken for the rupture to accelerate to the Rayleigh wave speed), shape and width of the rupture front. Similar differences during coseismic rupture nucleation have been reported previously by Rice & Uenishi (2010), where the initial weakening rate, governed by the exponent of power-law slip-weakening friction, affects rupture acceleration rate. However, a direct comparison to this study is not possible: in our case, acceleration rates are mainly influenced by coseismic stress drop differences. Specifically, PCZ, which has the lowest initial weakening rate, nucleates the fastest due to its high background stress, while DSW, with the highest rate, nucleates the slowest. Supplementary Section S4 (Supporting Information) presents results with the same coseismic stress drop for the four failure laws, where acceleration rates
align with initial weakening rates and are consistent with previous findings.

The observed differences at rupture onset become insignificant as the ruptures grow. This is expected when a rupture is much larger than its process zone, the region near the rupture front where weakening occurs. In this case, the small-scale yielding assumption of linear elastic fracture mechanics theory is satisfied. According to that theory, fracture energy (rather than the detailed shape of the failure law) controls subshear rupture dynamics (Freund 1990). This effect is enhanced as the ruptures accelerate toward the Rayleigh wave speed, the limiting speed for subshear mode II earthquake ruptures in elastic media, because this leads to shrinkage of the process zone (Supplementary Section S3, Supporting Information), a phenomenon known as Lorentz contraction (Ida 1973; Freund 1990).

An interesting observation pertains to the slip rate evolution of the DSW law (panel b). Initially, the slip rate front consists of two distinct peaks: the first one is associated with the change in slope of the failure law, and the second one with reaching the residual strength (Fig. 2 and Supplementary Section S2, Supporting Information). Paglialunga et al. (2022) and Ferry & Molinari (2025) reported similar results, showing two successive increases in slip rate in their numerical models with the same failure law. As the rupture propagates, these peaks gradually converge and eventually merge into a single peak (Supplementary Section S2, Supporting Information), as also observed in their studies. This merging is attributed to differences in the energy release rate during the two stages of linear weakening, as explained in section 5 of the supplementary material of Paglialunga et al. (2022).

While near-field earthquake rupture parameters can be obtained from dynamic rupture laboratory experiments with appropriate instrumentation, they are seldom recoverable from field data due to the scarcity of near-fault recordings. Fig. 8 illustrates some proxies of the far-field observables derived from dynamic simulations using different failure laws, which can be retrieved from seismic data and are commonly used for earthquake analysis. These proxies include source-time functions, rupture velocity evolution and spectra.

The recorded seismogram is the convolution of the source-time function, the elastic Green's function connecting the source and receiver, including the attenuation operator acting along the ray path, and the receiver response (Shearer 2019). As this study primarily focuses on modelling the source and does not address wave propagation in the medium, we examine directly the source-time functions (moment rate) obtained from our simulations. The seismic moment rate is (Kanamori & Brodsky 2004):

$$ \dot{M}_0(t) = G \bar{V}(t) S, \qquad (2) $$

where G is the shear modulus, S is the fault area and $\bar{V}$ is the mean fault slip rate. Moment rate cannot be meaningfully discussed for a 1-D fault, thus, we approximately convert our results to 3-D assuming radial rupture symmetry using (e.g. Lapusta & Rice 2003)

$$ \dot{M}_0(t) = 2\pi G \int_0^L V(r, t) r \, dr, \qquad (3) $$

where r is the distance relative to the left end of the fault and $V(r, t)$ is the fault slip rate. Fig. 8(a) shows the source-time functions for the different failure laws. They exhibit similar shapes, nearly triangular but with smooth onset and tail, although there are variations in the

moment acceleration time (time to reach the peak moment rate) and final moment (area under the moment rate curve).

Fig. 8(b) depicts earthquake rupture velocities. For all failure laws, the ruptures accelerate to the Rayleigh wave speed and propagate at velocities close to the Rayleigh wave speed until reaching the fault boundary. As the rupture velocity approaches the Rayleigh wave speed, the process zone shrinks due to Lorentz contraction (Supplementary Section S3, Supporting Information), and the specifics of the failure law used in these simulations become less significant. However, different rupture acceleration times are evident in the rupture velocity evolution among the different failure laws.

Another far-field observable that we consider is the spectrum of the source-time function, which is proportional to the far-field displacement spectra commonly used to estimate source dimensions and stress drop (Brune 1970). We show the moment rate spectrum in Fig. 8(c) for the four failure laws. The spectra exhibit similar high-frequency decay rates that are between $f^{-1}$ and $f^{-2}$, but differ in the plateau amplitude at low frequencies, which corresponds to the different seismic moment obtained for the different failure laws. The seismic moments, $M_0$, for the four cases are shown in Fig. 8(a), with the largest $M_0 = 1.48 \times 10^9$ N·m corresponding to PCZ and the smallest $M_0 = 1.19 \times 10^9$ N·m corresponding to DSW.

The discrepancies observed in the far-field observables discussed earlier are attributed to the differing coseismic stress drops. Different failure laws result in different rupture acceleration times (Fig. 8a or b), which are inversely proportional to the stress drop. PCZ exhibits the highest stress drop, while DSW has the lowest. Both the source-time functions (area under the curves in panel a) and the spectra (plateau levels in panel c) manifest differences in the mean slip (hence seismic moment—see Fig. 8a) obtained from the various failure laws. To demonstrate that most of these differences are solely a consequence of the different stress drops, we provide a set of simulations in the Supplementary Section S4 (Supporting Information) with the same four failure laws and the same coseismic stress drop, resulting in nearly indistinguishable far-field observables.

## 5 COMPARISON WITH 1-D SPRING-SLIDER RESULTS

A 1-D spring-slider model has been widely used in earthquake research as a simplified representation of strain accumulation and release throughout an earthquake cycle. In this section, we compare the results obtained from our 2-D simulations with those obtained using spring-slider models. Model parameters for the spring-slider model are given in the Supplementary Section S5 (Supporting Information). We could not match all the model parameters perfectly between 1-D and 2-D, thus, we compare the results qualitatively rather than quantitatively. Failure law parameters are perfectly matched.

The 1-D spring-slider model consists of a block pulled along a surface by a spring, with friction resisting the sliding of the block (Burridge & Knopoff 1967) (representing the different failure laws of interest in our case). The stretching of the spring corresponds to the elastic strain in the medium surrounding a fault segment. Under certain conditions, this system exhibits stick-slip behaviour: the block remains stationary as long as the frictional resistance exceeds the force from the stretched spring.
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Moment rate along the seismogenic portion of the fault (source-time function) (a), rupture velocity (b) and moment rate spectra (c). Horizontal black line in panel (b) shows the Rayleigh wave speed.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 9.** Logarithm of slip rate versus logarithm of reversed time for the four failure laws. The coseismic phase is on the left, and the interseismic phase is on the right of the plot. $V_0$ is the loading velocity, $V_{ei}$ is the velocity of early- to late-interseismic phase transition and $V_{p-c}$ is the velocity of pre-seismic to coseismic phase transition.</figcaption></figure>

Once the static frictional resistance is surpassed, it decreases to its dynamic value, and the block slips unstably, simulating an earthquake.

Fig. 9 displays the logarithm of slip rate plotted against the logarithm of reversed time for a single seismic cycle modelled with the four failure laws. The slip rate evolution curves are aligned to intersect at the loading velocity, $V_0$. Block velocity reaching $V_0$ almost coincides with the loading force reaching the maximum strength of the interface, $\tau_s$. Thus, we consider this moment as

the onset of the pre-seismic phase. More detailed information regarding the model, equations used, and a comprehensive comparison of the 1-D spring-slider earthquake cycle with different failure laws can be found in Bolotskaya & Hager (2022) and Bolotskaya (2023).

The main distinction between the 2-D FEM and 1-D spring-slider models is that, in the case of the 1-D spring-slider model, each phase of the earthquake cycle is governed by different stages of the failure law. This is due to a single node representing the entire
fault, whereas in 2-D models, different fault nodes are located at different points along the failure law curve during each stage of the earthquake cycle.

Nevertheless, we do observe similarities in the behaviour of the two systems. For instance, in the 2-D models, the inclusion of a strengthening segment in the failure law leads to an increased amount of pre-seismic slip (area under the slip rate curve before the crossing—compare with Fig. 6b) for the ECZ and PCZ failure laws.

In our 2-D FEM models, failure laws with gentler initial slopes, such as the PCZ, take longer to nucleate compared to failure laws with steeper initial weakening slopes. These findings align with the durations of the pre-seismic phase obtained from 1-D spring-slider models (compare the durations of acceleration from $V_0$ to $V_{p-c}$ in Fig. 9 with Fig. 6a). Different pre-seismic phase durations also result in varying full cycle durations (recurrence intervals) in the 1-D spring-slider models, with longer intervals observed for failure laws that include strengthening.

We observe nearly identical results for the coseismic phase in both the 2-D and 1-D spring-slider models. However, in the case of the 1-D spring-slider model, the similarity arises due to the residual strength segment of the failure laws (at $\tau_d$) being their common feature. In the 2-D models, this similarity is a consequence of the Lorentz contraction of the process zone. The disparity in earthquake magnitude observed in the 2-D models cannot be observed with 1-D spring-slider models since, in the 1-D spring-slider model, the stress drop is equal to the fault strength drop, $\tau_s - \tau_d$, which is the same for all four failure laws.

# 6 DISCUSSION

## 6.1 Rupture speed and implications for earthquake modelling

In dynamic rupture modelling, it is commonly assumed that once the earthquake rupture accelerates, the size of the process zone becomes small enough to apply a small-scale yielding approximation. In this scenario, the specific details of the failure law can be disregarded, and rupture propagation is solely governed by the balance between the energy release rate flux towards the rupture tip and fracture energy (Ampuero et al. 2006; Ripperger et al. 2007).

Our results support the validity of this assumption for dynamic ruptures that are sufficiently accelerated. As the process zone contracts, the specific details of shear strength and slip evolution within the process zone, governed by the failure law, become less significant. Instead, the dynamic ruptures are primarily controlled by the magnitude of fracture energy, $G_c$, and stress drop (Section 4.2 and Supplementary Section S4, Supporting Information). Thus, it appears appropriate to neglect the details of the failure law shape for ruptures approaching the Rayleigh wave speed, provided that they are modelled with appropriate stress drops, which differ for the different failure laws due to variations in the duration of their respective pre-seismic phases.

It is important to point out, however, that our findings are limited to subshear ruptures; Dunham (2007) previously demonstrated that variations in the initial weakening rates of failure laws influence the transition from subshear to supershear rupture speeds. Furthermore, few earthquakes have been observed to propagate at velocities close to the Rayleigh wave speed (Weng & Ampuero 2020). In the absence of acceleration up to the Rayleigh wave speed, the contraction of the process zone becomes less pronounced, and the details

of the failure law may regain significance during dynamic rupture propagation. Various mechanisms can limit the speed of earthquake rupture, including off-fault damage, fault branching (Washabaugh & Knauss 1994), heterogeneous surface roughness affecting frictional properties, increasing size of the process zone during rupture propagation (Fineberg et al. 1991, 1992; Sharon & Fineberg 1996; Yang & Ravi-Chandar 1996), heterogeneous initial stress state, or velocity-dependent fracture energy. Conducting dynamic rupture simulations incorporating different failure laws and any of these speed limiting mechanisms would be valuable to investigate the hypothesis that different failure law shapes can yield different results for dynamic rupture propagation with less pronounced Lorentz contraction.

Similar considerations apply to modelling the pre-seismic phase of earthquake rupture and earthquake nucleation. As demonstrated in Section 4.1, it is crucial to account for the various shapes of the failure laws during the pre-seismic phase of an earthquake, as they lead to significant differences in pre-seismic phase duration, associated slip, and slip rate depending on the failure law used.

## 6.2 Pre-seismic slip and earthquake observations

The question of earthquake nucleation and the transition from slow to fast slip in earthquakes is a widely debated topic in the literature. Two models have been proposed to explain the initiation of seismic rupture (Dodge et al. 1996). One model assumes that the accelerated moment release observed before large earthquakes (Bowman & King 2001) is triggered by transient aseismic slip on the fault interface (Bouchon et al. 2013; Ruiz et al. 2014). In the other model, a cascade of failures eventually triggers the main shock (Dodge et al. 1996).

There are laboratory observations reporting the occurrence of a slow sliding phase preceding the main shock (Wu & McLaskey 2019; Gvirtzman & Fineberg 2021). Additionally, slow slip events have been detected using geodetic data prior to large earthquakes (Socquet et al. 2017), as well as in the context of icequakes (Barcheck et al. 2021). Although the accumulation of observations of slow slip events is increasing, there is still ongoing debate regarding a plausible mechanism that can explain this behaviour. Several mechanisms proposed in the literature include a change in frictional behaviour from velocity weakening to velocity strengthening (Shibazaki & Iio 2003; Shibazaki & Shimamoto 2007), rate-and-state friction near the stability boundary (Liu & Rice 2005, 2007) with specific constraints on fault length, pore pressure and other parameters, as well as dilatancy-related strengthening and the transition from drained to undrained conditions that quench the instability (Segall et al. 2010).

We present an alternative mechanism that can give rise to slow slip preceding an earthquake, without requiring specific conditions related to pore pressure, fault stiffness or slip rate dependence of the failure law. This mechanism involves the presence of initial strengthening in a slip-dependent failure law, associated with a significant magnitude of slip and stress increase. Similar models were previously considered by Shibazaki & Matsu'ura (1992) and Matsu'ura et al. (1992). Such a failure law is particularly appropriate for modelling earthquakes on cemented faults that have regained cohesion due to microcrack healing and mineralization of pore space and cracks (Smith et al. 2013). Fig. 5(d) illustrates the prominent effect of initial strengthening observed in the ECZ failure law, where the background shear stress eventually surpasses the initial strength
of the failure law. As a result, the entire fault undergoes a slow quasi-static sliding process months before the occurrence of the earthquake.

# 7 CONCLUSIONS

Here, we address the impact of different slip-dependent failure parametrizations prescribed along the fault on the pre-seismic and coseismic phases of the earthquake cycle through numerical simulations using the finite element method. The pre-seismic phase is modelled using a quasi-static approach, considering background stressing rate over a timescale of months, while the coseismic phase is modelled dynamically, considering timescales on the order of milliseconds. The slip and stress values obtained from the last converged time step of the quasi-static simulation are used as inputs for the dynamic simulation.

Both model portions employ the same bulk rheology, model geometry and set of failure laws prescribed along the fault. The failure laws considered include linear SW, DSW, parabolic cohesive zone (PCZ), and ECZ laws. The strength drop and fracture energy parameters are set to be the same for all failure laws, but the shapes and effective slip-weakening distances vary.

Our results demonstrate the importance of accounting for the different shapes of the failure laws during the pre-seismic phase. Failure law shape variations lead to significant differences in pre-seismic phase duration, associated slip and lateral extent of slip front. Failure laws with gentler initial weakening slopes, such as PCZ, require significantly more time to initiate dynamic rupture compared to failure laws with steeper initial weakening slopes, such as DSW. The presence of a strengthening segment in the failure law increases the amount of slip obtained before the rupture becomes unstable and influences the lateral extent of slip propagation. In cases, where the initial strength of the failure law is below the background stress reached during fault loading, such as ECZ, we observe quasi-static propagation of the slip front to the end of the fault.

During dynamic rupture propagation, due to the elastic nature of our bulk rheology, we observe Lorentz contraction of the process zone as the rupture speed approaches the Rayleigh wave speed. As the process zone shrinks, the specific details of shear strength and slip evolution within the process zone become insignificant. Consequently, we obtain dynamic ruptures that are almost indistinguishable from each other in terms of source-time functions and spectra, despite being governed by different failure laws. This holds true as long as we model them with the same coseismic stress drops (Supplementary Section S4, Supporting Information).

However, an important implication of the different durations of the pre-seismic phases for the different failure laws is the variation in the background stress reached prior to the onset of the coseismic phase. The different resulting coseismic stress drops lead to ruptures with varying acceleration times and earthquake magnitudes. For example, PCZ inherits the highest stress drop from the quasi-static simulation, resulting in the shortest rupture acceleration time and largest earthquake magnitude, while DSW exhibits the lowest stress drop, leading to the slowest rupture onset and smallest earthquake magnitude.

Our results support the validity of the small-scale yielding approximation for sufficiently accelerated dynamic ruptures. In our case, dynamic rupture propagation is primarily controlled by the

magnitude of fracture energy, $G_c$, and dynamic stress drop. Thus, it is appropriate to overlook the specific details of the failure law shape for ruptures approaching the Rayleigh wave speed, as long as they are modelled with the appropriate dynamic stress drops, which differ for the different failure laws due to variations in their respective pre-seismic phases.

# ACKNOWLEDGMENTS

The authors thank the editorial team and the three anonymous reviewers for their constructive comments, which helped improve the article. This study was supported by ExxonMobil Upstream Research Company (Grant Number EM10654).

# SUPPORTING INFORMATION

Supplementary data are available at *GJIRAS* online.

**Section S1:** Equations for the different failure laws used.

**Figure S1.1.** Shear strength versus slip curves for the different failure parametrizations.

**Section S2:** Double-peak rupture in the case of the DSW failure law.

**Figure S2.1.** Dynamic simulation results for the DSW failure law showing two peak-slip-velocity fronts: spatial distribution of (A) slip rate and (B) shear stress, at two times: $t = 3.2$ ms and $t = 5.2$ ms. Vertical dashed lines indicate slip rate peaks.

**Section S3:** Lorentz contraction of the process zone illustrated with the dynamic rupture simulation results.

**Figure S3.1.** Shear stress versus fault coordinate for SW failure law (A) from the dynamic simulation. Colour-time. Single shear stress peak during the onset of the rupture for the different failure laws (B). Dashed arrows mark the process zone length. Single shear stress peak for an accelerated rupture for the different failure laws (C).

**Figure S3.2.** Dynamic process zone size versus time from the simulation data and the analytical prediction for the SW failure law.

**Section S4:** Dynamic rupture simulations with different failure laws and the same stress drop.

**Figure S4.1.** Mean slip rate along the seismogenic portion of the fault (source-time function) (A), rupture velocity (B) and moment rate spectra (C).

**Section S5:** 1D spring-slider model parameters.

**Table S5.1.** 1-D spring-slider model parameters.

# DATA AVAILABILITY

No data were used in this article.

The codes for the custom failure laws implemented for PyLith v2 are available at: https://doi.org/10.5281/zenodo.14047393

Fully dynamic equations for 1-D spring-slider models for different failure laws and the codes used to produce Fig. 9 are available at: https://doi.org/10.5281/zenodo.7035922

# REFERENCES

* Aagaard, B.T., Knepley, M.G. & Williams, C.A., 2013. A domain decomposition approach to implementing fault slip in finite-element models of quasi-static and dynamic crustal deformation, *J. Geophys. Res. Solid Earth*, **118**(6), 3059–3079.

* Aki, K., 1979. Characterization of barriers on an earthquake fault, *J. Geophys. Res. Solid Earth*, **84**(B11), 6140–6148.
Aki, K., 1984. Asperities, barriers, characteristic earthquakes and strong motion prediction, *J. Geophys. Res. Solid Earth*, **89**(B7), 5867–5872.

Alfano, G., 2006. On the influence of the shape of the interface law on the application of cohesive-zone models, *Compos. Sci. Technol.*, **66**(6), 723–730.

Alfano, G. & Crisfield, M., 2001. Finite element interface models for the delamination analysis of laminated composites: mechanical and computational issues, *Int. J. Numer. Methods Eng.*, **50**(7), 1701–1736.

Allix, O. & Corigliano, A., 1996. Modeling and simulation of crack propagation in mixed-modes interlaminar fracture specimens, *Int. J. Fract.*, **77**(2), 111–140.

Allix, O., Ladeveze, P. & Corigliano, A., 1995. Damage analysis of interlaminar fracture specimens, *Compos. Struct.*, **31**(1), 61–74.

Ampuero, J.-P., Ripperger, J. & Mai, P.M., 2006. *Properties of Dynamic Earthquake Ruptures with Heterogeneous Stress Drop*, pp. 255–261, American Geophysical Union (AGU).

Andrews, D., 1976a. Rupture propagation with finite stress in antiplane strain, *J. geophys. Res.*, **81**(20), 3575–3582.

Andrews, D., 1976b. Rupture velocity of plane strain shear cracks, *J. geophys. Res.*, **81**(32), 5679–5687.

Aochi, H. & Fukuyama, E., 2002. Three-dimensional nonplanar simulation of the 1992 landers earthquake, *J. Geophys. Res. Solid Earth*, **107**(B2), ESE-4.

Archuleta, R.J., 1984. A faulting model for the 1979 imperial valley earthquake, *J. Geophys. Res. Solid Earth*, **89**(B6), 4559–4585.

Barcheck, G., Brodsky, E., Fulton, P., King, M., Siegfried, M. & Tulaczyk, S., 2021. Migratory earthquake precursors are dominant on an ice stream fault, *Sci. Adv.*, **7**(6), eabd0105, doi: 10.1126/sciadv.abd0105.

Barenblatt, G.I., 1959. The formation of equilibrium cracks during brittle fracture. general ideas and hypotheses. axially-symmetric cracks, *J. Appl. Math. Mech.*, **23**(3), 622–636.

Bolotskaya, E., 2023. *Effects of fault failure parameterization and bulk rheology on earthquake rupture, Doctoral dissertation*, Massachusetts Institute of Technology.

Bolotskaya, E. & Hager, B.H., 2022. A 1D spring-slider model with a simple poly-linear failure law produces rich variations in slip behaviour, *Bull. seism. Soc. Am.*, **112**(6), doi: 10.1785/0120220052.

Bouchon, M., Durand, V., Marsan, D., Karabulut, H. & Schmittbuhl, J., 2013. The long precursory phase of most large interplate earthquakes, *Nat. Geosci.*, **6**(4), 299–302.

Bowman, D.D. & King, G.C., 2001. Accelerating seismicity and stress accumulation before large earthquakes, *Geophys. Res. Lett.*, **28**(21), 4039–4042.

Brace, W. & Byerlee, J., 1966. Stick-slip as a mechanism for earthquakes, *Science*, **153**(3739), 990–992.

Brantut, N. & Viesca, R.C., 2015. Earthquake nucleation in intact or healed rocks, *J. Geophys. Res. Solid Earth*, **120**(1), 191–209.

Brune, J.N., 1970. Tectonic stress and the spectra of seismic shear waves from earthquakes, *J. geophys. Res.*, **75**(26), 4997–5009.

Burridge, R. & Knopoff, L., 1967. Model and theoretical seismicity, *Bull. seism. Soc. Am.*, **57**(3), 341–371.

Chandra, N., Li, H., Shet, C. & Ghonem, H., 2002. Some issues in the application of cohesive zone models for metal–ceramic interfaces, *Int. J. Solids Struct.*, **39**(10), 2827–2855.

Day, S.M., Dalguer, L.A., Lapusta, N. & Liu, Y., 2005. Comparison of finite difference and boundary integral solutions to three-dimensional spontaneous rupture, *J. Geophys. Res. Solid Earth*, **110**(B12), doi: 10.1029/2005JB003813.

Dieterich, J.H., 1979. Modeling of rock friction: experimental results and constitutive equations, *J. Geophys. Res. Solid Earth*, **84**(B5), 2161–2168.

Dieterich, J.H., 1992. Earthquake nucleation on faults with rate- and state-dependent strength, *Tectonophysics*, **211**(1-4), 115–134.

Dodge, D.A., Beroza, G.C. & Ellsworth, W., 1996. Detailed observations of california foreshock sequences: implications for the earthquake initiation process, *J. Geophys. Res. Solid Earth*, **101**(B10), 22 371–22 392.

Dunham, E.M., 2007. Conditions governing the occurrence of supershear ruptures under slip-weakening friction, *J. Geophys. Res. Solid Earth*, **112**(B7), doi: 10.1029/2006JB004717.

Ferry, R. & Molinari, J.-F., 2025. Do slip-weakening laws shapes influence rupture dynamics?, *Extreme Mech. Lett.*, **74**, 102 272, doi: 10.1016/j.eml.2024.102272.

Fineberg, J., Gross, S.P., Marder, M. & Swinney, H.L., 1991. Instability in dynamic fracture, *Phys. Rev. Lett.*, **67**(4), 457, doi: 10.1103/Phys-RevLett.67.457.

Fineberg, J., Gross, S.P., Marder, M. & Swinney, H.L., 1992. Instability in the propagation of fast cracks, *Phys. Rev. B*, **45**(10), 5146, doi: 10.1103/Phys-RevB.45.5146.

Freund, L.B., 1990. *Dynamic Fracture Mechanics*, Cambridge University Press, Cambridge.

Gvirtzman, S. & Fineberg, J., 2021. Nucleation fronts ignite the interface rupture that initiates frictional motion, *Nat. Phys.*, **17**(9), 1037–1042.

Hillerborg, A., Modéer, M. & Petersson, P.-E., 1976. Analysis of crack formation and crack growth in concrete by means of fracture mechanics and finite elements, *Cem. Concr. Res.*, **6**(6), 773–781.

Ida, Y., 1972. Cohesive force across the tip of a longitudinal-shear crack and griffith's specific surface energy, *J. geophys. Res.*, **77**(20), 3796–3805.

Ida, Y., 1973. Stress concentration and unsteady propagation of longitudinal shear cracks, *J. geophys. Res.*, **78**(17), 3418–3429.

Kanamori, H. & Brodsky, E.E., 2004. The physics of earthquakes, *Rep. Prog. Phys.*, **67**(8), 1429, doi: 10.1088/0034-4885/67/8/R03.

Karner, S.L., Marone, C. & Evans, B., 1997. Laboratory study of fault healing and lithification in simulated fault gouge under hydrothermal conditions, *Tectonophysics*, **277**(1-3), 41–55.

Lapusta, N. & Rice, J.R., 2003. Nucleation and early seismic propagation of small and large events in a crustal earthquake model, *J. Geophys. Res. Solid Earth*, **108**(B4), doi: 10.1029/2001JB000793.

Liu, Y. & Rice, J.R., 2005. Aseismic slip transients emerge spontaneously in three-dimensional rate and state modeling of subduction earthquake sequences, *J. Geophys. Res. Solid Earth*, **110**(B8), doi: 10.1029/2004JB003424.

Liu, Y. & Rice, J.R., 2007. Spontaneous and triggered aseismic deformation transients in a subduction fault model, *J. Geophys. Res. Solid Earth*, **112**(B9), doi: 10.1029/2007JB004930.

Marone, C., 1998. Laboratory-derived friction laws and their application to seismic faulting, *Annu. Rev. Earth Planet. Sci.*, **26**(1), 643–696.

Matsuura, M., Kataoka, H. & Shibazaki, B., 1992. Slip-dependent friction law and nucleation processes in earthquake rupture, *Tectonophysics*, **211**(1-4), 135–148.

Mi, Y., Crisfield, M., Davies, G. & Hellweg, H., 1998. Progressive delamination using interface elements, *J. Compos. Mater.*, **32**(14), 1246–1272.

Needleman, A., 1990. An analysis of tensile decohesion along an interface, *J. Mech. Phys. Solids*, **38**(3), 289–324.

Ohnaka, M., 1996. Nonuniformity of the constitutive law parameters for shear rupture and quasistatic nucleation to dynamic rupture: a physical model of earthquake generation processes., *Proc. Natl. Acad. Sci. USA*, **93**(9), 3795–3802.

Ohnaka, M., 2003. A constitutive scaling law and a unified comprehension for frictional slip failure, shear fracture of intact rock, and earthquake rupture, *J. Geophys. Res.: Solid Earth*, **108**(B2), doi: 10.1029/2000JB000123.

Ohnaka, M., 2013. *The Physics of Rock Failure and Earthquakes*, Cambridge University Press, Cambridge.

Ohnaka, M., Akatsu, M., Mochizuki, H., Odedra, A., Tagashira, F. & Yamamoto, Y., 1997. A constitutive law for the shear failure of rock under lithospheric conditions, *Tectonophysics*, **277**(1-3), 1–27.

Paglialunga, F., Passelègue, F.X., Brantut, N., Barras, F., Lebihain, M. & Violay, M., 2022. On the scale dependence in the dynamics of frictional rupture: constant fracture energy versus size-dependent breakdown work, *Earth planet. Sci. Lett.*, **584**, 117442, doi: 10.1016/j.epsl.2022.117442.

Palmer, A.C. & Rice, J.R., 1973. The growth of slip surfaces in the progressive failure of over-consolidated clay, *Proc. R. Soc. Lond. A. Math. Phys. Sci.*, **332**(1591), 527–548.
Park, K. & Paulino, G.H., 2011. Cohesive zone models: a critical review of traction-separation relationships across fracture surfaces, *Appl. Mech. Rev.*, **64**(6), doi: 10.1115/1.4023110.

Rice, J.R. & Tse, S.T., 1986. Dynamic motion of a single degree of freedom system following a rate and state dependent friction law, *J. Geophys. Res. Solid Earth*, **91**(B1), 521–530.

Rice, J.R. & Uenishi, K., 2010. Rupture nucleation on an interface with a power-law relation between stress and displacement discontinuity, *Int. J. Fract.*, **163**, 1–13.

Ripperger, J., Ampuero, J.-P., Mai, P. & Giardini, D., 2007. Earthquake source characteristics from dynamic rupture with constrained stochastic fault stress, *J. Geophys. Res. Solid Earth*, **112**(B4), doi: 10.1029/2006JB004515.

Rubin, A.M., 2008. Episodic slow slip events and rate-and-state friction, *J. Geophys. Res. Solid Earth*, **113**(B11), doi: 10.1029/2008JB0056 42.

Rudnicki, J.W., 1980. Fracture mechanics applied to the earth's crust, *Annu. Rev. Earth Planet. Sci.*, **8**, 489, doi: 10.1146/an-nurev.ea.08.050180.002421.

Ruina, A., 1983. Slip instability and state variable friction laws, *J. Geophys. Res. Solid Earth*, **88**(B12), 10 359–10 370.

Ruiz, S. et al., 2014. Intense foreshocks and a slow slip event preceded the 2014 iquique m w 8.1 earthquake, *Science*, **345**(6201), 1165–1169.

Scholz, C.H., 1998. Earthquakes and friction laws, *Nature*, **391**(6662), 37– 42.

Segall, P., Rubin, A.M., Bradley, A.M. & Rice, J.R., 2010. Dilatant strengthening as a mechanism for slow slip events, *J. Geophys. Res. Solid Earth*, **115**(B12), doi: 10.1029/2010JB007449.

Sharon, E. & Fineberg, J., 1996. Microbranching instability and the dynamic fracture of brittle materials, *Phys. Rev. B*, **54**(10), 7128, doi: 10.1103/Phys-RevB.54.7128.

Shearer, P.M., 2019. *Introduction to Seismology*, Cambridge University Press, Cambridge.

Shibazaki, B. & Iio, Y., 2003. On the physical mechanism of silent slip events along the deeper part of the seismogenic zone, *Geophys. Res. Lett.*, **30**(9), doi: 10.1029/2003GL017047.

Shibazaki, B. & Matsu'ura, M., 1992. Spontaneous processes for nucleation, dynamic propagation, and stop of earthquake rupture, *Geophys. Res. Lett.*, **19**(12), 1189–1192.

Shibazaki, B. & Shimamoto, T., 2007. Modelling of short-interval silent slip events in deeper subduction interfaces considering the frictional properties at the unstable—stable transition regime, *Geophys. J. Int.*, **171**(1), 191–205.

Smith, S.A., Bistacchi, A., Mitchell, T.M., Mittempergher, S. & Di Toro, G., 2013. The structure of an exhumed intraplate seismogenic fault in crystalline basement, *Tectonophysics*, **599**, 29–44.

Socquet, A. et al., 2017. An 8 month slow slip event triggers progressive nucleation of the 2014 chile megathrust, *Geophys. Res. Lett.*, **44**(9), 4046– 4053.

Uenishi, K. & Rice, J.R., 2003. Universal nucleation length for slip-weakening rupture instability under nonuniform fault loading, *J. Geophys. Res. Solid Earth*, **108**(B1), doi: 10.1029/2001JB001681.

Washabaugh, P.D. & Knauss, W., 1994. A reconciliation of dynamic crack velocity and rayleigh wave speed in isotropic brittle solids, *Int. J. Fract.*, **65**(2), 97–114.

Weng, H. & Ampuero, J.-P., 2020. Continuum of earthquake rupture speeds enabled by oblique slip, *Nat. Geosci.*, **13**(12), 817–821.

Wilson, J., Chester, J. & Chester, F., 2003. Microfracture analysis of fault growth and wear processes, punchbowl fault, san andreas system, california, *J. Struct. Geol.*, **25**(11), 1855–1873.

Wu, B.S. & McLaskey, G.C., 2019. Contained laboratory earthquakes ranging from slow to fast, *J. Geophys. Res. Solid Earth*, **124**(10), 10 270– 10 291.

Xu, S., Fukuyama, E. & Yamashita, F., 2019. Robust estimation of rupture properties at propagating front of laboratory earthquakes, *J. Geophys. Res. Solid Earth*, **124**(1), 766–787.

Yang, B. & Ravi-Chandar, K., 1996. On the role of the process zone in dynamic fracture, *J. Mech. Phys. Solids*, **44**(12), 1955–1976.

## APPENDIX: A

Here, we numerically compute the nucleation length for linear slip-weakening failure law as the solution of an eigenvalue problem, following Uenishi & Rice (2003) while including a loading patch with constant friction coefficient, $\mu_d$ (Fig. 1). We only outline the main steps of the method. The shear stress along the fault in static equilibrium is given by the following boundary integral

$$ \tau(x) = \tau_{bg} - \frac{\mu^*}{2\pi} \int_{a_{-}}^{a_{+}} \frac{\partial D}{\partial \xi} d\xi, \quad (A1) $$

where $\tau_{bg}$ is the uniform background load, $\mu^* = \frac{G}{1-\nu}$ in our case of mode II rupture, $a_{-}$ and $a_{+}$ are the tips of the slip zone and $D$ is slip. At equilibrium, the shear stress must satisfy the failure law in the slip-weakening zone and the constant strength condition on the loading patch. Assuming $D(x) < D_c$ during nucleation, as verified in our simulations, we write these two conditions as

$$ \tau(x) = \tau_s - W(x)D(x) \quad (A2) $$

within the slip zone ($a_{-} < x < a_{+}$), where $W(x) = \frac{\tau_s - \tau_d}{D_c}$ inside the slip-weakening zone and $W(x) = 0$ inside the loading patch. Finite stress conditions are also assumed near the edges of the slip zone. Taking the time derivative of the equilibrium equation and applying a change of variable that maps the slip zone ($x \in [a_{-}, a_{+}]$) onto a fixed reference segment ($X \in [-1, +1]$), we obtain

$$ -W(X)V(X) = \dot{\tau}_{bg} - \frac{\mu^*}{2\pi a} \int_{-1}^{+1} \frac{\partial V}{\partial s} ds, \quad (A3) $$

where $V$ is slip velocity and $a = \frac{a_{+}-a_{-}}{2}$ is the half-size of the slip zone. At the approach of instability, the slip changes due to an increment of external load become very large, and the equilibrium equation becomes dominated by its slip-dependent terms. Keeping only those dominant terms leads to

$$ \lambda \Lambda(X) V(X) = \frac{1}{2\pi a} \int_{-1}^{+1} \frac{\partial V}{\partial s} ds, \quad (A4) $$

where $\lambda = \frac{a(\tau_s - \tau_d)}{\mu^* D_c}$ is the normalized slip zone size, and $\Lambda(X) = 1$ within the slip-weakening zone and $\Lambda(X) = 0$ within the loading patch. For a given spatial distribution of $\Lambda(X)$, the equation above is an eigenvalue problem: it has solutions $V(X)$ (eigenfunctions) only for a discrete set of values of $\lambda$ (eigenvalues). In particular, the smallest eigenvalue $\lambda$ corresponds to non-negative slip rate and yields the critical slip zone size at which instability starts, which is associated with the nucleation size. For instance, in the case where $\Lambda(X) = 1$ along the whole slip zone, the problem considered by Uenishi & Rice (2003), the smallest eigenvalue is $\lambda = 0.579$, which leads to the nucleation length given by eq. (1) (UR length). In practice, we numerically solve a discretized version of the eigenvalue problem (eq. A4): we split the reference segment $[-1, +1]$ into $N = 4000$ elements of equal size, and assume constant slip within each element; the kernel of the integral is given by eqs (32) and (33) of Dieterich (1992).

We consider two geometries: asymmetric (a) and symmetric (b) cases. In the asymmetric case, the loading patch is located at the left boundary of the fault, and slip propagates in only one direction (the same geometry as used in the main article and shown in Fig. 1). In the symmetric case, the loading patch is located at the centre of the fault, allowing slip to propagate in both directions from its boundaries into the slip-weakening sections of the fault. Fig. A1 illustrates
<figure><img src="image_11.png" /><figcaption>**Figure A1.** Slip contours coloured by time for (a)—asymmetric loading case and (b)—symmetric loading case. Loading patch length = 0.4 m.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure A2.** Normalized nucleation length versus normalized loading patch half length. Lines—analytical solution, markers—simulation data.</figcaption></figure>

both geometries used and shows how the loading patch length and nucleation length are measured. Fig. A2 shows the numerical solutions of eq. (A4) for both geometries (solid lines). As expected, as the loading patch length approaches zero, our solution converges to the nucleation length reported by Uenishi & Rice (2003): the UR length for the asymmetric geometry and half UR length for the symmetric geometry. In the symmetric case, we measure only one side of the nucleation (Fig. A1), hence the half. For large loading patch lengths, both solutions become independent of the loading patch length and asymptotically approach approximately 0.63 UR length.

To verify the analytical solution for earthquake nucleation length on faults with a pre-existing loading patch, we use a quasi-static modelling approach (Section 2) to model earthquake nucleation

on faults with varying loading patch lengths and slip-weakening failure law across the remainder of the fault. As discussed in the main manuscript, the quasi-static simulation terminates when the inertial term becomes too large for the solver to converge. We use the last converged time step to measure the nucleation length (Fig. A1). The results of 14 simulations (7 for each geometry) with loading patch half-length varying from 0.01 to 0.8 m are shown in Fig. A2 (markers). We observe good correlation between the simulation results and the numerical solution of a discretized version of the eigenvalue problem (eq. A4) (solid lines). The discrepancies are likely due to the constant time step used in the finite element simulations, which affects the time of numerical non-convergence, and mesh resolution, which affects the tolerance in measuring the nucleation length.