# A proxy implementation of thermal pressurization for earthquake cycle modelling on rate-and-state faults

Marco T Herrera, Jean P Ampuero, Jorge G F Crempien

► To cite this version:

Marco T Herrera, Jean P Ampuero, Jorge G F Crempien. A proxy implementation of thermal pressurization for earthquake cycle modelling on rate-and-state faults. Geophysical Journal International, 2024, 237 (3), pp.1432-1441. 10.1093/gji/ggae113. hal-04593792

HAL Id: hal-04593792

https://hal.science/hal-04593792

Submitted on 30 May 2024

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
# A proxy implementation of thermal pressurization for earthquake cycle modelling on rate-and-state faults

Marco T Herrera, Jean P Ampuero, Jorge G F Crempien

► To cite this version:

Marco T Herrera, Jean P Ampuero, Jorge G F Crempien. A proxy implementation of thermal pressurization for earthquake cycle modelling on rate-and-state faults. Geophysical Journal International, 2024, 237 (3), pp.1432-1441. 10.1093/gji/ggae113. hal-04593792

HAL Id: hal-04593792

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

Submitted on 30 May 2024

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
# A proxy implementation of thermal pressurization for earthquake cycle modelling on rate-and-state faults

Marco T. Herrera<fnref n="1" />,<fnref n="2" />, Jean P. Ampuero<fnref n="1" />,<fnref n="3" /> and Jorge G. F. Crempien<fnref n="1" />,<fnref n="2" />

Accepted 2024 March 12. Received 2024 January 30; in original form 2023 August 21

## SUMMARY

The reduction of effective normal stress during earthquake slip due to thermal pressurization of fault zone pore fluids is a significant fault weakening mechanism. Explicit incorporation of this process into frictional fault models involves solving the diffusion equations for fluid pressure and temperature outside the fault at each time step, which significantly increases the computational complexity. Here, we propose a proxy for thermal pressurization implemented through a modification of the rate-and-state friction law. This approach is designed to emulate the fault weakening and the relationship between breakdown energy and slip resulting from thermal pressurization and is appropriate for fully dynamic simulations of multiple earthquake cycles. It preserves the computational efficiency of conventional rate-and-state friction models, which in turn can enable systematic studies to advance our understanding of the effects of fault weakening on earthquake mechanics. In 2.5-D simulations of pulse-like ruptures on faults with finite seismogenic width, based on our thermal pressurization proxy, we find that the spatial distribution of slip velocity near the rupture front is consistent with the conventional square-root singularity, despite continued slip-weakening within the pulse, once the rupture has propagated a distance larger than the rupture width. An unconventional singularity appears only at shorter rupture distances. We further derive and verify numerically a theoretical estimate of the breakdown energy dissipated by our implementation of thermal pressurization. These results support the use of fracture mechanics theory to understand the propagation and arrest of very large earthquakes.

**Key words:** Friction; Numerical modelling; Earthquake dynamics; Subduction zone processes.

## 1 INTRODUCTION

Understanding how faults weaken during rapid earthquake slip is important for constraining the minimum level of fault stress required to generate large ruptures, which can help us understand the state of stress in Earth's crust and improve earthquake hazard assessment (Abercrombie & Rice 2005; Tinti et al. 2005; Noda & Lapusta 2013; Viesca & Garagash 2015; Perry et al. 2020). Earthquake models with enhanced dynamic weakening have been successful in reproducing the low apparent strength of faults interseismically and during coseismic slip (i.e. low average ratio of shear stress to normal stress) and their low heat production (Thomas et al. 2014; Viesca & Garagash 2015; Perry et al. 2020; Lambert & Lapusta 2023). Two important mechanisms for dramatic fault weakening incorporated in such models are flash heating (FH, Mase & Smith 1987; Goldsby & Tullis 2011; Noda et al. 2011) and thermal pressurization (TP,

Rempel & Rice 2006; Garagash 2012; Noda & Lapusta 2013; Perry et al. 2020). Both arise from frictional heating, the conversion of frictional work into heat. FH is a reduction of the friction coefficient resulting from temperature rise at the microscopic contacts of rough fault surfaces or gouge, often modelled as a friction decrease inversely proportional to slip rate at high slip rates above ~0.1 m s⁻<fnref n="1" /> (Rice 2006; Noda et al. 2011; Thomas et al. 2014; Lambert et al. 2021). TP is a reduction of the effective normal stress due to thermal expansion of fault zone pore fluids and their consequent pressure increase (Rempel & Rice 2006; Garagash 2012; Noda & Lapusta 2013; Viesca & Garagash 2015; Perry et al. 2020).

Incorporating a realistic scaling of breakdown energy in earthquake models is essential to capture first-order aspects of rupture propagation and arrest (Weng & Ampuero 2022). As shown by Viesca & Garagash (2015), TP can account for important aspects of

<footnote n="1">¹Department of Structural and Geotechnical Engineering, Pontificia Universidad Católica de Chile. Postal code: 7820436, Santiago, Chile. E-mail: mtherrera2@uc.cl.</footnote>

<footnote n="2">²National Research Center for Integrated Natural Disaster Management (CIGIDEN). Postal code: 7820436, Santiago, Chile</footnote>

<footnote n="3">³Université Côte d'Azur, Observatoire de la Côte d'Azur, IRD, CNRS, Geoazur. Postal code: 06560, Valbonne, France</footnote>
# A proxy implementation of thermal pressurization for earthquake cycle modelling on rate-and-state faults

Marco T. Herrera<fnref n="1" />,<fnref n="1" />,<fnref n="2" /> Jean P. Ampuero<fnref n="1" />,<fnref n="3" /> and Jorge G. F. Crempien<fnref n="1" />,<fnref n="1" />,<fnref n="2" />

Accepted 2024 March 12. Received 2024 January 30; in original form 2023 August 21

## SUMMARY

The reduction of effective normal stress during earthquake slip due to thermal pressurization of fault zone pore fluids is a significant fault weakening mechanism. Explicit incorporation of this process into frictional fault models involves solving the diffusion equations for fluid pressure and temperature outside the fault at each time step, which significantly increases the computational complexity. Here, we propose a proxy for thermal pressurization implemented through a modification of the rate-and-state friction law. This approach is designed to emulate the fault weakening and the relationship between breakdown energy and slip resulting from thermal pressurization and is appropriate for fully dynamic simulations of multiple earthquake cycles. It preserves the computational efficiency of conventional rate-and-state friction models, which in turn can enable systematic studies to advance our understanding of the effects of fault weakening on earthquake mechanics. In 2.5-D simulations of pulse-like ruptures on faults with finite seismogenic width, based on our thermal pressurization proxy, we find that the spatial distribution of slip velocity near the rupture front is consistent with the conventional square-root singularity, despite continued slip-weakening within the pulse, once the rupture has propagated a distance larger than the rupture width. An unconventional singularity appears only at shorter rupture distances. We further derive and verify numerically a theoretical estimate of the breakdown energy dissipated by our implementation of thermal pressurization. These results support the use of fracture mechanics theory to understand the propagation and arrest of very large earthquakes.

**Key words:** Friction; Numerical modelling; Earthquake dynamics; Subduction zone processes.

## 1 INTRODUCTION

Understanding how faults weaken during rapid earthquake slip is important for constraining the minimum level of fault stress required to generate large ruptures, which can help us understand the state of stress in Earth's crust and improve earthquake hazard assessment (Abercrombie & Rice 2005; Tinti et al. 2005; Noda & Lapusta 2013; Viesca & Garagash 2015; Perry et al. 2020). Earthquake models with enhanced dynamic weakening have been successful in reproducing the low apparent strength of faults interseismically and during coseismic slip (i.e. low average ratio of shear stress to normal stress) and their low heat production (Thomas et al. 2014; Viesca & Garagash 2015; Perry et al. 2020; Lambert & Lapusta 2023). Two important mechanisms for dramatic fault weakening incorporated in such models are flash heating (FH, Mase & Smith 1987; Goldsby & Tullis 2011; Noda et al. 2011) and thermal pressurization (TP,

Rempel & Rice 2006; Garagash 2012; Noda & Lapusta 2013; Perry et al. 2020). Both arise from frictional heating, the conversion of frictional work into heat. FH is a reduction of the friction coefficient resulting from temperature rise at the microscopic contacts of rough fault surfaces or gouge, often modelled as a friction decrease inversely proportional to slip rate at high slip rates above ~0.1 m s⁻<fnref n="1" /> (Rice 2006; Noda et al. 2011; Thomas et al. 2014; Lambert et al. 2021). TP is a reduction of the effective normal stress due to thermal expansion of fault zone pore fluids and their consequent pressure increase (Rempel & Rice 2006; Garagash 2012; Noda & Lapusta 2013; Viesca & Garagash 2015; Perry et al. 2020).

Incorporating a realistic scaling of breakdown energy in earthquake models is essential to capture first-order aspects of rupture propagation and arrest (Weng & Ampuero 2022). As shown by Viesca & Garagash (2015), TP can account for important aspects of

<footnote n="1">¹Department of Structural and Geotechnical Engineering, Pontificia Universidad Católica de Chile. Postal code: 7820436, Santiago, Chile. E-mail: mtherrera2@uc.cl.</footnote>

<footnote n="2">²National Research Center for Integrated Natural Disaster Management (CIGIDEN). Postal code: 7820436, Santiago, Chile</footnote>

<footnote n="3">³Université Côte d'Azur, Observatoire de la Côte d'Azur, IRD, CNRS, Geoazur. Postal code: 06560, Valbonne, France</footnote>
the scaling of breakdown energy with slip inferred from seismo-logical observations over seven orders of slip magnitude, spanning small to large earthquakes. This suggests that TP is a widespread and prominent process for fault weakening.

Models with enhanced dynamic weakening due to TP can explain both the increasing trend in breakdown energy with increasing event size and the near magnitude invariance of stress drops (Viesca & Garagash 2015; Perry et al. 2020). Breakdown energy $G_c$ can be constrained through scaling relations, depending on final slip $D$, which have been compiled over a wide range of earthquake sizes by combining different approaches: dynamic earthquake modelling (Weng & Yang 2018; Gallovic et al. 2019), laboratory experiments (Nielsen et al. 2016) and seismological methods such as kinematic source inversion (Tinti et al. 2009; Viesca & Garagash 2015). As illustrated in Fig. 1, the widely used rate-and-state (R-S, Dieterich 1978, 1979; Ruina 1983) and linear slip-weakening (Ida 1972; Andrews 1976) friction models have a distinct residual strength and produce a breakdown energy $G_c$ that does not depend on slip $D$. The linear slip-weakening friction model prescribes a linear decrease of the friction coefficient $f$ as a function of slip from a static value $f_s$ to a dynamic value $f_d$ reached at a characteristic slip distance $D_c$ (Ida 1972; Andrews 1976). It features a constant $G_c = \sigma(f_s - f_d)D_c/2$, where $\sigma$ is the effective normal stress. The R-S friction model prescribes a dependence of friction on slip rate $V$ and a state variable that represents memory effects (Dieterich 1978, 1979; Ruina 1983), and is adequate at relatively low slip rates ($10^{-9}-10^{-3}$ m s$^{-1}$, Rubin & Ampuero 2005; Lapusta & Liu 2009; Allison & Dunham 2018; Cattania 2019). It yields a $G_c$ independent on slip and weakly (log-arithmically) dependent on slip rate (Rubin & Ampuero 2005). In contrast, TP models feature continued slip-weakening, leading to $G_c \propto D^{2/3}$ in the diffusion-dominated regime at large slip and $G_c \propto D^2$ in the undrained-adiabatic regime at small slip $\ll 0.1$ m (Viesca & Garagash 2015).

Current seismic cycle models combine R-S friction and FH by modifying the friction law (Noda et al. 2011; Noda & Lapusta 2013; Thomas et al. 2014; Perry et al. 2020; Lambert et al. 2021). In contrast, incorporating TP requires solving the diffusion equations for fluid pressure and temperature within the fault zone, coupled to fault slip (Noda & Lapusta 2010, 2013; Mavrommatis et al. 2017; Perry et al. 2020). This implies a higher computational cost and complexity, which limits the capacity to conduct large sets of simulations, especially with realistically large ratios of fault size to the length scales arising from friction and TP. Therefore, in Section 2, we propose a TP proxy through a modification of the R-S friction law, which is convenient for implementation in earthquake simulators, has a low computational cost, and is designed to mimic the slip-weakening and scaling between breakdown energy and slip that emerge from TP. We incorporate the thermal weakening effect via the friction coefficient and keep the effective normal stress constant in time. The proposed proxy implementation has similar complexity and computational cost as conventional R-S friction.

In Section 3, we present results of earthquake cycle simulations with the TP proxy and FH, implemented in a fully dynamic 2.5-D earthquake cycle simulator that accounts for the effect of a finite seismogenic width $W$. We show that the asymptotic behaviour of slip rate, as a function of distance behind the rupture front, supports the applicability of fracture mechanics theory for very large earthquakes with rupture lengths $> W$ despite continued weakening by TP. To enable such fracture mechanics analyses, in Section 4 we verify that our TP proxy produces the scaling between breakdown energy and slip intended in its design.

# 2 THERMAL PRESSURIZATION PROXY

## 2.1 Rate-and-state friction law

We take as starting point the formulation of R-S friction (Dieterich 1978, 1979; Rice & Ruina 1983; Ruina 1983) with FH used by (Harris et al. 2018; Noda et al. 2011; Thomas et al. 2014). This formulation is summarized in this section.

The shear stress $\tau$ on the fault satisfies

$$ \tau = \sigma f(V, \Psi), \qquad (1) $$

where $\sigma$ is the effective normal stress on the fault and $f$ the friction coefficient, which depends on slip velocity, $V$, and a fault state variable, $\Psi$, as

$$ f(V, \Psi) = a \operatorname{arcsinh} \left( \frac{V}{2V_0} \exp \frac{\Psi}{a} \right), \qquad (2) $$

where $a$ is an R-S parameter and $V_0$ a reference velocity. This is a regularized version of the classical R-S friction law, designed to prevent a singularity at zero velocity (Ben-Zion & Rice 1997; Lapusta et al. 2000). We consider the slip law for the state evolution:

$$ \dot{\Psi} = -\frac{V}{L}(\Psi - \Psi_{ss}(V)), \qquad (3) $$

where

$$ \Psi_{ss}(V) = f_0 + b \log \frac{V_0}{V}, \qquad (4) $$

$L$ is a characteristic slip distance, $f_0$ a reference friction coefficient and $b$ an R-S parameter. Eqs (1)–(4) completely define the frictional fault strength.

The state variable, $\Psi$, is a re-formulation of the classical R-S state variable $\theta$:

$$ \Psi = f_0 + b \log \frac{V_0 \theta}{L}. \qquad (5) $$

This re-formulation is computationally advantageous, since it improves the stability of the numerical friction solver.

At steady state, when $\dot{\Psi} = 0$, the state variable $\Psi$ is equal to $\Psi_{ss}(V)$ and the friction coefficient is

$$ f_{ss}(V) = a \operatorname{arcsinh} \left( \frac{V}{2V_0} \exp \frac{\Psi_{ss}(V)}{a} \right). \qquad (6) $$

The corresponding classical state variable at steady state is $\theta_{ss} = L/V$.

## 2.2 Rate-and-state friction law with flash heating

The abrupt velocity-dependent reduction of the friction coefficient due to FH is introduced by re-defining the steady-state friction as

$$ f_{ss}^{FH}(V) = \frac{f_{ss}(V) - f_w}{1 + \frac{V}{V_w}} + f_w, \qquad (7) $$

where $f_w$ is the residual friction coefficient and $V_w$ the threshold velocity for the activation of FH. The steady-state values of friction and state variable including FH are related by an equation analogous to eq. (6):

$$ f_{ss}^{FH}(V) = a \operatorname{arcsinh} \left( \frac{V}{2V_0} \exp \frac{\Psi_{ss}^{FH}(V)}{a} \right). \qquad (8) $$

Combining eqs (7) and (8), we obtain an expression for $\Psi_{ss}^{FH}(V)$, which we then use in the state evolution law as

$$ \dot{\Psi} = -\frac{V}{L}(\Psi - \Psi_{ss}^{FH}(V)). \qquad (9) $$
the scaling of breakdown energy with slip inferred from seismo-logical observations over seven orders of slip magnitude, spanning small to large earthquakes. This suggests that TP is a widespread and prominent process for fault weakening.

Models with enhanced dynamic weakening due to TP can explain both the increasing trend in breakdown energy with increasing event size and the near magnitude invariance of stress drops (Viesca & Garagash 2015; Perry et al. 2020). Breakdown energy $G_c$ can be constrained through scaling relations, depending on final slip $D$, which have been compiled over a wide range of earthquake sizes by combining different approaches: dynamic earthquake modelling (Weng & Yang 2018; Gallovic et al. 2019), laboratory experiments (Nielsen et al. 2016) and seismological methods such as kinematic source inversion (Tinti et al. 2009; Viesca & Garagash 2015). As illustrated in Fig. 1, the widely used rate-and-state (R-S, Dieterich 1978, 1979; Ruina 1983) and linear slip-weakening (Ida 1972; Andrews 1976) friction models have a distinct residual strength and produce a breakdown energy $G_c$ that does not depend on slip $D$. The linear slip-weakening friction model prescribes a linear decrease of the friction coefficient $f$ as a function of slip from a static value $f_s$ to a dynamic value $f_d$ reached at a characteristic slip distance $D_c$ (Ida 1972; Andrews 1976). It features a constant $G_c = \sigma(f_s - f_d)D_c/2$, where $\sigma$ is the effective normal stress. The R-S friction model prescribes a dependence of friction on slip rate $V$ and a state variable that represents memory effects (Dieterich 1978, 1979; Ruina 1983), and is adequate at relatively low slip rates ($10^{-9}-10^{-3}$ m s$^{-1}$, Rubin & Ampuero 2005; Lapusta & Liu 2009; Allison & Dunham 2018; Cattania 2019). It yields a $G_c$ independent on slip and weakly (logarithmically) dependent on slip rate (Rubin & Ampuero 2005). In contrast, TP models feature continued slip-weakening, leading to $G_c \propto D^{2/3}$ in the diffusion-dominated regime at large slip and $G_c \propto D^2$ in the undrained-adiabatic regime at small slip $\ll 0.1$ m (Viesca & Garagash 2015).

Current seismic cycle models combine R-S friction and FH by modifying the friction law (Noda et al. 2011; Noda & Lapusta 2013; Thomas et al. 2014; Perry et al. 2020; Lambert et al. 2021). In contrast, incorporating TP requires solving the diffusion equations for fluid pressure and temperature within the fault zone, coupled to fault slip (Noda & Lapusta 2010, 2013; Mavrommatis et al. 2017; Perry et al. 2020). This implies a higher computational cost and complexity, which limits the capacity to conduct large sets of simulations, especially with realistically large ratios of fault size to the length scales arising from friction and TP. Therefore, in Section 2, we propose a TP proxy through a modification of the R-S friction law, which is convenient for implementation in earthquake simulators, has a low computational cost, and is designed to mimic the slip-weakening and scaling between breakdown energy and slip that emerge from TP. We incorporate the thermal weakening effect via the friction coefficient and keep the effective normal stress constant in time. The proposed proxy implementation has similar complexity and computational cost as conventional R-S friction.

In Section 3, we present results of earthquake cycle simulations with the TP proxy and FH, implemented in a fully dynamic 2.5-D earthquake cycle simulator that accounts for the effect of a finite seismogenic width $W$. We show that the asymptotic behaviour of slip rate, as a function of distance behind the rupture front, supports the applicability of fracture mechanics theory for very large earthquakes with rupture lengths $> W$ despite continued weakening by TP. To enable such fracture mechanics analyses, in Section 4 we verify that our TP proxy produces the scaling between breakdown energy and slip intended in its design.

# 2 THERMAL PRESSURIZATION PROXY

## 2.1 Rate-and-state friction law

We take as starting point the formulation of R-S friction (Dieterich 1978, 1979; Rice & Ruina 1983; Ruina 1983) with FH used by (Harris et al. 2018; Noda et al. 2011; Thomas et al. 2014). This formulation is summarized in this section.

The shear stress $\tau$ on the fault satisfies

$$ \tau = \sigma f(V, \Psi), \qquad (1) $$

where $\sigma$ is the effective normal stress on the fault and $f$ the friction coefficient, which depends on slip velocity, $V$, and a fault state variable, $\Psi$, as

$$ f(V, \Psi) = a \operatorname{arcsinh} \left( \frac{V}{2V_0} \exp \frac{\Psi}{a} \right), \qquad (2) $$

where $a$ is an R-S parameter and $V_0$ a reference velocity. This is a regularized version of the classical R-S friction law, designed to prevent a singularity at zero velocity (Ben-Zion & Rice 1997; Lapusta et al. 2000). We consider the slip law for the state evolution:

$$ \dot{\Psi} = -\frac{V}{L}(\Psi - \Psi_{ss}(V)), \qquad (3) $$

where

$$ \Psi_{ss}(V) = f_0 + b \log \frac{V_0}{V}, \qquad (4) $$

$L$ is a characteristic slip distance, $f_0$ a reference friction coefficient and $b$ an R-S parameter. Eqs (1)–(4) completely define the frictional fault strength.

The state variable, $\Psi$, is a re-formulation of the classical R-S state variable $\theta$:

$$ \Psi = f_0 + b \log \frac{V_0 \theta}{L}. \qquad (5) $$

This re-formulation is computationally advantageous, since it improves the stability of the numerical friction solver.

At steady state, when $\dot{\Psi} = 0$, the state variable $\Psi$ is equal to $\Psi_{ss}(V)$ and the friction coefficient is

$$ f_{ss}(V) = a \operatorname{arcsinh} \left( \frac{V}{2V_0} \exp \frac{\Psi_{ss}(V)}{a} \right). \qquad (6) $$

The corresponding classical state variable at steady state is $\theta_{ss} = L/V$.

## 2.2 Rate-and-state friction law with flash heating

The abrupt velocity-dependent reduction of the friction coefficient due to FH is introduced by re-defining the steady-state friction as

$$ f_{ss}^{FH}(V) = \frac{f_{ss}(V) - f_w}{1 + \frac{V}{V_w}} + f_w, \qquad (7) $$

where $f_w$ is the residual friction coefficient and $V_w$ the threshold velocity for the activation of FH. The steady-state values of friction and state variable including FH are related by an equation analogous to eq. (6):

$$ f_{ss}^{FH}(V) = a \operatorname{arcsinh} \left( \frac{V}{2V_0} \exp \frac{\Psi_{ss}^{FH}(V)}{a} \right). \qquad (8) $$

Combining eqs (7) and (8), we obtain an expression for $\Psi_{ss}^{FH}(V)$, which we then use in the state evolution law as

$$ \dot{\Psi} = -\frac{V}{L}(\Psi - \Psi_{ss}^{FH}(V)). \qquad (9) $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Slip-weakening curves and breakdown energy $G_c$ for (a) classical R-S friction, (b) linear slip weakening friction and (c) TP weakening in the diffusion-dominated regime. The curves are the shear stress versus slip response at a given point along a rupture, in examples with little or no stress overshot/undershot. Stars indicate the initial and final states of an earthquake. Final slip $D$ and stress drop $\Delta\tau$ are indicated.</figcaption></figure>

The frictional fault strength including FH is completely defined by eqs (1), (2) and (7)–(9).

## 2.3 Rate-and-state friction law with TP proxy

To implement a TP proxy, we introduce a new friction coefficient under steady-state conditions, denoted as $f_{ss}^{TP}$. To maintain simplicity, we implement the concept of strong weakening by modifying the friction coefficient rather than by modifying the effective normal stress $\sigma$. This approach facilitates the derivation of a scaling relationship between $G_c$ and slip, following the approach outlined by Viesca & Garagash (2015). Our objective is not to fully replicate the intricate physics of TP, but rather to satisfy its resulting scaling relationship between breakdown energy and slip.

We define a new state variable $\phi$, that approximates the slip only in the coseismic phase. Following the procedure shown in Beeler et al. (2008), we design the evolution of $\phi$ to meet two physical constraints. During earthquakes, $\phi$ should approximate the coseismic slip. Thus, its derivative should approximate the slip rate:

$$ \dot{\phi} \approx V. \qquad (10) $$

During interseismic phases, $\phi$ should reset to zero, thus

$$ \dot{\phi} \approx -\phi/T^*, \qquad (11) $$

where $T^*$ is a characteristic time scale for the state variable to return to zero after the coseismic phase. We set that value at 0.1 s, shorter than the typical rise time of large earthquakes. We define the $\phi$ evolution equation as a weighted sum of both expressions:

$$ \dot{\phi} = \Gamma V - (1 - \Gamma) \frac{\phi}{T^*}. \qquad (12) $$

As $\phi$ behaves differently, depending on the velocity (coseismic or interseismic), we include a weight $\Gamma(V)$ to separate both regimes. The weight is a sigmoid function, with velocity threshold $V_{th}$, whose smoothness favours numerical stability:

$$ \Gamma(V) = \frac{1}{1 + \exp\left(-40\left(\frac{V}{V_{th}} - 1\right)\right)}. \qquad (13) $$

Our objective is to model large earthquakes, with slip values >0.1 m. Following Viesca & Garagash (2015), the relevant TP regime for large slips is the diffusion-dominated regime, in which the breakdown energy,

$$ G_c = \int_0^D [\tau(\delta) - \tau_r] \, d\delta, \qquad (14) $$

scales with coseismic slip $D$ as $G_c \propto D^{2/3}$. To satisfy this scaling, $\tau(\delta) - \tau(0)$ must scale with $\delta^{-1/3}$, as shown in Section 4. In our proxy TP implementation, $\phi$ tends to $\delta$ during earthquakes, and we introduce the $\delta^{-1/3}$ scaling by modifying the steady-state friction as:

$$ f_{ss}^{TP}(V, \phi) = \frac{f_{ss}(V)}{\left(1 + \frac{\phi}{L^*}\right)^{1/3}}, \qquad (15) $$

where $L^*$ is a characteristic slip distance for TP (Rempel & Rice 2006; Rice 2006):

$$ L^* = \left( \frac{2\rho c}{f\Lambda} \right)^2 \frac{\sqrt{\alpha_{hy}} + \sqrt{\alpha_{th}}}{V_*}^2, \qquad (16) $$

where $\rho c$ is the specific heat, $\alpha_{hy}$ is the hydraulic diffusivity, $\alpha_{th}$ is the thermal diffusivity, $f$ is the friction coefficient prior the TP, $V_* = f\sigma c_s/\mu$ is a characteristic elastodynamic slip rate, with $c_s$ the shear wave speed; and $\Lambda$ is the TP coefficient relating increments of pore fluid pressure to increments in temperature. Note that in eq. (15), the state variable, $\phi$, is not necessarily at steady state. Nonetheless, for historical reasons and simplicity, we keep the terminology as $f_{ss}^{TP}$. The steady-state friction and state variables, including TP, are related by an equation analogous to eq. (6):

$$ f_{ss}^{TP}(V, \phi) = a \operatorname{arcsinh} \left( \frac{V}{2V_0} \exp \frac{\Psi_{ss}^{TP}(V, \phi)}{a} \right). \qquad (17) $$

Combining eqs (15) and (17), we obtain an expression for $\Psi_{ss}^{TP}(V, \phi)$, which we use in the state evolution law as

$$ \dot{\Psi} = -\frac{V}{L}(\Psi - \Psi_{ss}^{TP}(V, \phi)). \qquad (18) $$

To include both weakening mechanisms, FH and TP, we adopt the reasonable assumption that FH operates at the beginning of slip, while the normal stress changes induced by TP in the diffusion-dominated regime occur at much larger slips (Rice 2006; Viesca & Garagash 2015; Lambert et al. 2021). Thus, we replace in eq. (15) $f_{ss}$ by $f_{ss}^{FH}$ as:

$$ f_{ss}^{FH+TP} = \frac{f_{ss}^{FH}(V)}{\left(1 + \frac{\phi}{L^*}\right)^{1/3}}. \qquad (19) $$

We equate this expression with

$$ f_{ss}^{FH+TP} = a \operatorname{arcsinh} \left( \frac{V}{2V_0} \exp \frac{\Psi_{ss}^{FH+TP}}{a} \right), \qquad (20) $$
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Slip-weakening curves and breakdown energy $G_c$ for (a) classical R-S friction, (b) linear slip weakening friction and (c) TP weakening in the diffusion-dominated regime. The curves are the shear stress versus slip response at a given point along a rupture, in examples with little or no stress overshot/undershot. Stars indicate the initial and final states of an earthquake. Final slip $D$ and stress drop $\Delta\tau$ are indicated.</figcaption></figure>

The frictional fault strength including FH is completely defined by eqs (1), (2) and (7)–(9).

## 2.3 Rate-and-state friction law with TP proxy

To implement a TP proxy, we introduce a new friction coefficient under steady-state conditions, denoted as $f_{ss}^{TP}$. To maintain simplicity, we implement the concept of strong weakening by modifying the friction coefficient rather than by modifying the effective normal stress $\sigma$. This approach facilitates the derivation of a scaling relationship between $G_c$ and slip, following the approach outlined by Viesca & Garagash (2015). Our objective is not to fully replicate the intricate physics of TP, but rather to satisfy its resulting scaling relationship between breakdown energy and slip.

We define a new state variable $\phi$, that approximates the slip only in the coseismic phase. Following the procedure shown in Beeler et al. (2008), we design the evolution of $\phi$ to meet two physical constraints. During earthquakes, $\phi$ should approximate the coseismic slip. Thus, its derivative should approximate the slip rate:

$$ \dot{\phi} \approx V. \qquad (10) $$

During interseismic phases, $\phi$ should reset to zero, thus

$$ \dot{\phi} \approx -\phi/T^*, \qquad (11) $$

where $T^*$ is a characteristic time scale for the state variable to return to zero after the coseismic phase. We set that value at 0.1 s, shorter than the typical rise time of large earthquakes. We define the $\phi$ evolution equation as a weighted sum of both expressions:

$$ \dot{\phi} = \Gamma V - (1 - \Gamma) \frac{\phi}{T^*}. \qquad (12) $$

As $\phi$ behaves differently, depending on the velocity (coseismic or interseismic), we include a weight $\Gamma(V)$ to separate both regimes. The weight is a sigmoid function, with velocity threshold $V_{th}$, whose smoothness favours numerical stability:

$$ \Gamma(V) = \frac{1}{1 + \exp\left(-40\left(\frac{V}{V_{th}} - 1\right)\right)}. \qquad (13) $$

Our objective is to model large earthquakes, with slip values >0.1 m. Following Viesca & Garagash (2015), the relevant TP regime for large slips is the diffusion-dominated regime, in which the breakdown energy,

$$ G_c = \int_0^D [\tau(\delta) - \tau_r] \, d\delta, \qquad (14) $$

scales with coseismic slip $D$ as $G_c \propto D^{2/3}$. To satisfy this scaling, $\tau(\delta) - \tau(0)$ must scale with $\delta^{-1/3}$, as shown in Section 4. In our proxy TP implementation, $\phi$ tends to $\delta$ during earthquakes, and we introduce the $\delta^{-1/3}$ scaling by modifying the steady-state friction as:

$$ f_{ss}^{TP}(V, \phi) = \frac{f_{ss}(V)}{\left(1 + \frac{\phi}{L^*}\right)^{1/3}}, \qquad (15) $$

where $L^*$ is a characteristic slip distance for TP (Rempel & Rice 2006; Rice 2006):

$$ L^* = \left( \frac{2\rho c}{f\Lambda} \right)^2 \frac{\sqrt{\alpha_{hy}} + \sqrt{\alpha_{th}}}{V_*}, \qquad (16) $$

where $\rho c$ is the specific heat, $\alpha_{hy}$ is the hydraulic diffusivity, $\alpha_{th}$ is the thermal diffusivity, $f$ is the friction coefficient prior the TP, $V_* = f\sigma c_s/\mu$ is a characteristic elastodynamic slip rate, with $c_s$ the shear wave speed; and $\Lambda$ is the TP coefficient relating increments of pore fluid pressure to increments in temperature. Note that in eq. (15), the state variable, $\phi$, is not necessarily at steady state. Nonetheless, for historical reasons and simplicity, we keep the terminology as $f_{ss}^{TP}$. The steady-state friction and state variables, including TP, are related by an equation analogous to eq. (6):

$$ f_{ss}^{TP}(V, \phi) = a \operatorname{arcsinh} \left( \frac{V}{2V_0} \exp \frac{\Psi_{ss}^{TP}(V, \phi)}{a} \right). \qquad (17) $$

Combining eqs (15) and (17), we obtain an expression for $\Psi_{ss}^{TP}(V, \phi)$, which we use in the state evolution law as

$$ \dot{\Psi} = -\frac{V}{L}(\Psi - \Psi_{ss}^{TP}(V, \phi)). \qquad (18) $$

To include both weakening mechanisms, FH and TP, we adopt the reasonable assumption that FH operates at the beginning of slip, while the normal stress changes induced by TP in the diffusion-dominated regime occur at much larger slips (Rice 2006; Viesca & Garagash 2015; Lambert et al. 2021). Thus, we replace in eq. (15) $f_{ss}$ by $f_{ss}^{FH}$ as:

$$ f_{ss}^{FH+TP} = \frac{f_{ss}^{FH}(V)}{\left(1 + \frac{\phi}{L^*}\right)^{1/3}}. \qquad (19) $$

We equate this expression with

$$ f_{ss}^{FH+TP} = a \operatorname{arcsinh} \left( \frac{V}{2V_0} \exp \frac{\Psi_{ss}^{FH+TP}}{a} \right), \qquad (20) $$
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Results of a 2.5-D earthquake cycle simulation on a homogeneous dip-slip shallow fault. Friction is applied over a fault length of 100 km, including a central VW segment of length 80 km surrounded by VS segments, and steady slip at velocity $V_{PL} = 10^{-9}$ m s$^{-1}$ is applied elsewhere along the fault. Model parameters are $W = 20$ km, $\sigma = 50$ MPa, $L = 0.0075$ m, $L^* = 0.1$ m, $a_{VW} = 0.01$, $a_{VS} = 0.015$, $b_{VW} = 0.015$, $b_{VS} = 0.005$, $V_{th} = 10^{-3}$ m s$^{-1}$, $f_0 = 0.6$, $V_0 = 10^{-6}$ m s$^{-1}$, $\mu = 32$ GPa and $c_s = 3464$ m s$^{-1}$. (a) Slip evolution along the whole fault plotted every 0.5 s during earthquakes (blue curves) and every 5 yr in between earthquakes (red curves). (b) Shear stress as a function of slip, (c) friction coefficient as a function of time and (d) slip and TP state variable $\phi$ as a function of time at a representative point near the middle of the fault.</figcaption></figure>

to obtain $\Psi_{ss}^{FH+TP}$, which we then use in the state evolution law:

$$ \dot{\Psi} = -\frac{V}{L}(\Psi - \Psi_{ss}^{FH+TP}). \qquad (21) $$

The computational cost of the combined FH and TP model is comparable to that of a classical R-S fault model.

## 3 SAMPLE SIMULATION RESULTS

We implemented the TP proxy and FH in the fully dynamic 2.5-D earthquake cycle simulation software SEM2DPACK, based on the spectral element method (Ampuero 2002, 2012; Kaneko et al. 2008, 2011; Liang et al. 2022). The 2.5-D formulation is an approximation of the 3-D problem, accounting for the finite seismogenic width

$W$ (Weng & Ampuero 2019, 2020; Liang et al. 2022). The slip on a deeply buried fault with relatively uniform stress drop has a semi-elliptical depth profile, which can be approximated as a half cosine of wavelength $2W$. The 2.5-D model is thus derived by assuming the displacement field has a sinusoidal depth dependence with wavelength $2W$ (Weng & Ampuero 2019). In Appendix A, we calibrate the value of a geometric coefficient involved in the 2.5-D model so that it produces results consistent with 3-D dip-slip earthquake cycle simulations.

An example of the fault response in a 2.5-D seismic cycle simulation with FH and TP is shown in Fig. 2. The fault comprises a central segment with velocity-weakening (VW) R-S properties, surrounded by segments with velocity-strengthening (VS) properties. Results are plotted during five cycles, starting at the second
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Results of a 2.5-D earthquake cycle simulation on a homogeneous dip-slip shallow fault. Friction is applied over a fault length of 100 km, including a central VW segment of length 80 km surrounded by VS segments, and steady slip at velocity $V_{PL} = 10^{-9}$ m s$^{-1}$ is applied elsewhere along the fault. Model parameters are $W = 20$ km, $\sigma = 50$ MPa, $L = 0.0075$ m, $L^* = 0.1$ m, $a_{VW} = 0.01$, $a_{VS} = 0.015$, $b_{VW} = 0.015$, $b_{VS} = 0.005$, $V_{th} = 10^{-3}$ m s$^{-1}$, $f_0 = 0.6$, $V_0 = 10^{-6}$ m s$^{-1}$, $\mu = 32$ GPa and $c_s = 3464$ m s$^{-1}$. (a) Slip evolution along the whole fault plotted every 0.5 s during earthquakes (blue curves) and every 5 yr in between earthquakes (red curves). (b) Shear stress as a function of slip, (c) friction coefficient as a function of time and (d) slip and TP state variable $\phi$ as a function of time at a representative point near the middle of the fault.</figcaption></figure>

to obtain $\Psi_{ss}^{FH+TP}$, which we then use in the state evolution law:

$$ \dot{\Psi} = -\frac{V}{L}(\Psi - \Psi_{ss}^{FH+TP}). \qquad (21) $$

The computational cost of the combined FH and TP model is comparable to that of a classical R-S fault model.

## 3 SAMPLE SIMULATION RESULTS

We implemented the TP proxy and FH in the fully dynamic 2.5-D earthquake cycle simulation software SEM2DPACK, based on the spectral element method (Ampuero 2002, 2012; Kaneko et al. 2008, 2011; Liang et al. 2022). The 2.5-D formulation is an approximation of the 3-D problem, accounting for the finite seismogenic width

$W$ (Weng & Ampuero 2019, 2020; Liang et al. 2022). The slip on a deeply buried fault with relatively uniform stress drop has a semi-elliptical depth profile, which can be approximated as a half cosine of wavelength $2W$. The 2.5-D model is thus derived by assuming the displacement field has a sinusoidal depth dependence with wavelength $2W$ (Weng & Ampuero 2019). In Appendix A, we calibrate the value of a geometric coefficient involved in the 2.5-D model so that it produces results consistent with 3-D dip-slip earthquake cycle simulations.

An example of the fault response in a 2.5-D seismic cycle simulation with FH and TP is shown in Fig. 2. The fault comprises a central segment with velocity-weakening (VW) R-S properties, surrounded by segments with velocity-strengthening (VS) properties. Results are plotted during five cycles, starting at the second
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Dynamic response of a 2.5-D homogeneous dip-slip buried fault for cases with only R-S friction, only TP, only FH and including both (FH+TP). Model parameters are $W = 10$ km, $σ = 50$ MPa, $f_0 = 0.6$, $a = 0.0125$, $b = 0.0172$, $L = 0.015$ m, $L^* = 1$ m, $f_w = 0.2$, $V_{th} = 10^{-3}$ m s$^{-1}$, $V_0 = 10^{-6}$ m s$^{-1}$, $V_w = 0.1$ m s$^{-1}$, $μ = 32$ GPa and $c_s = 3464$ m s$^{-1}$. (a) Initial shear stress, including a nucleation patch with higher stress at the centre of the fault. (b) Shear stress as a function of slip at a point located at 10 km from the nucleation patch.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Spatial distribution of slip rate at two different times, when the rupture front has propagated up to positions $X = W$ and $4W$. Model parameters are the same as in the TP model shown in Fig. 3. (a) Slip rate as a function of distance $x$ relative to the location of the peak slip rate. The slip rate has a pulse-like shape. (b) Spatial decay of slip rate, with a focus on short distances from the rupture front. This log-log scale plot shows that the asymptotic behaviour is similar to a conventional singularity, $v ∝ x^{-0.5}$, at large rupture propagation distance ($X = 4W$, red curve), and to an unconventional singularity, $v ∝ x^{-0.25}$, at small rupture propagation distance ($X = W$, blue curve). (c) Same as (b), but with semi-logarithmic scale, showing an exponential decay of slip rate at larger distances behind the rupture front, $log(v) ∝ -x$.</figcaption></figure>

cycle to avoid effects of initial conditions. The simulation produces ‘characteristic earthquake’ behaviour: events spanning the whole fault occur with regular recurrence time and similar spatial and temporal distributions of slip (Fig. 2a). The dependence of shear stress on slip follows the desired scaling $τ ∝ δ^{-1/3}$ at large slip for each event (Fig. 2b). The shear stress prior to each event ($τ/σ ∼ 0.3$ near the centre of the fault, Fig. 2c) is well below the peak friction coefficient (~0.6). As designed, the state variable $φ$ mimics well the coseismic slip of each event and resets to zero rapidly after each event (Fig. 2d).

We now compare the slip-weakening response of models with and without FH (Fig. 3). We set equivalent initial conditions (Fig. 3a) to compare the shear stress evolution as a function of slip (Fig. 3b). For reference, we also show a model with only FH (no TP) and a model with only R-S friction. For the latter, we had to set up a higher initial shear stress $τ_0$ to achieve a runaway rupture. Although FH alone reduces the shear strength substantially (red curve in Fig. 3b), it does not produce continued slip-weakening. When combining both mechanisms (magenta curve in Fig. 3b), weakening at short slip values is dominated by FH, while at large slip values it is dominated by TP. Only the models that include TP result in a $δ^{-1/3}$ decay of
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Dynamic response of a 2.5-D homogeneous dip-slip buried fault for cases with only R-S friction, only TP, only FH and including both (FH+TP). Model parameters are $W = 10$ km, $σ = 50$ MPa, $f_0 = 0.6$, $a = 0.0125$, $b = 0.0172$, $L = 0.015$ m, $L^* = 1$ m, $f_w = 0.2$, $V_{th} = 10^{-3}$ m s$^{-1}$, $V_0 = 10^{-6}$ m s$^{-1}$, $V_w = 0.1$ m s$^{-1}$, $μ = 32$ GPa and $c_s = 3464$ m s$^{-1}$. (a) Initial shear stress, including a nucleation patch with higher stress at the centre of the fault. (b) Shear stress as a function of slip at a point located at 10 km from the nucleation patch.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Spatial distribution of slip rate at two different times, when the rupture front has propagated up to positions $X = W$ and $4W$. Model parameters are the same as in the TP model shown in Fig. 3. (a) Slip rate as a function of distance $x$ relative to the location of the peak slip rate. The slip rate has a pulse-like shape. (b) Spatial decay of slip rate, with a focus on short distances from the rupture front. This log-log scale plot shows that the asymptotic behaviour is similar to a conventional singularity, $v ∝ x^{-0.5}$, at large rupture propagation distance ($X = 4W$, red curve), and to an unconventional singularity, $v ∝ x^{-0.25}$, at small rupture propagation distance ($X = W$, blue curve). (c) Same as (b), but with semi-logarithmic scale, showing an exponential decay of slip rate at larger distances behind the rupture front, $log(v) ∝ -x$.</figcaption></figure>

cycle to avoid effects of initial conditions. The simulation produces ‘characteristic earthquake’ behaviour: events spanning the whole fault occur with regular recurrence time and similar spatial and temporal distributions of slip (Fig. 2a). The dependence of shear stress on slip follows the desired scaling $τ ∝ δ^{-1/3}$ at large slip for each event (Fig. 2b). The shear stress prior to each event ($τ/σ ∼ 0.3$ near the centre of the fault, Fig. 2c) is well below the peak friction coefficient (~0.6). As designed, the state variable $φ$ mimics well the coseismic slip of each event and resets to zero rapidly after each event (Fig. 2d).

We now compare the slip-weakening response of models with and without FH (Fig. 3). We set equivalent initial conditions (Fig. 3a) to compare the shear stress evolution as a function of slip (Fig. 3b). For reference, we also show a model with only FH (no TP) and a model with only R-S friction. For the latter, we had to set up a higher initial shear stress $τ_0$ to achieve a runaway rupture. Although FH alone reduces the shear strength substantially (red curve in Fig. 3b), it does not produce continued slip-weakening. When combining both mechanisms (magenta curve in Fig. 3b), weakening at short slip values is dominated by FH, while at large slip values it is dominated by TP. Only the models that include TP result in a $δ^{-1/3}$ decay of
<figure><img src="image_6.png" /><figcaption>**Figure 5.** Slip-weakening response and breakdown energy $G_c$ on an R-S fault with TP proxy, from numerical simulation and theoretical predictions.</figcaption></figure>

shear stress with slip. We observe stress undershoot in models that include FH, due to their fast healing process, which results in slip pulses with shorter rise time than in models without FH.

We further verify that our large earthquake models with TP are compatible with fracture mechanics theory. Large earthquakes with rupture lengths >W are pulse-like, due to the 2.5-D effect of the finite seismogenic width W (Fig. 4a). When the rupture has propagated over distances significantly larger than W, the slip rate decays with distance x behind the rupture front as $v \propto x^{-0.5}$ (red curve in Fig. 4b). This asymptotic behaviour emerges behind the peak slip velocity location, at distances sufficiently larger than the slip acceleration zone size (analogous to a cohesive zone or process zone size) and smaller than the slip pulse size. This square root decay is the same as in singular pulse models with uniform stress drop (Freund 1979; Nielsen & Madariaga 2003). This conventional singularity asymptotics supports the use of fracture mechanics theory to understand the propagation and arrest of very large earthquakes (larger than wide), despite the continued weakening induced by TP. Only when the rupture has propagated over distances shorter than or comparable to W, the slip rate decays as $v \propto x^{-0.25}$ (blue curve in Fig. 4b), which is the expected ‘unconventional singularity’ in 2-D rupture models with continued slip-weakening due to TP (Viesca & Garagash 2015; Brener & Bouchbinder 2021; Paglialunga et al. 2024). At further distance into the tail of the pulse, the slip rate decays exponentially (Fig. 4c), which is the expected behaviour induced by the finite rupture width W in 2.5-D models (fig. S6e and f of Weng & Ampuero 2022).

## 4 BREAKDOWN ENERGY OF THE TP PROXY

Our implementation is designed to satisfy the relation between $G_c$ and slip $\delta$ derived for TP in the diffusion-dominated regime at slips exceeding ~0.1 m by Viesca & Garagash (2015), which they validated by comparison to seismological observations. Here, we verify numerically that our implementation indeed produces the intended $G_c(\delta)$ relation. A different decay consistent with the undrained-adiabatic regime of TP at short slip, $G_c \propto \delta^2$ (Viesca & Garagash

2015), could be similarly introduced in the model by modifying eq. (15).

To estimate $G_c$, we consider an earthquake with initial shear stress $\tau_0$ right before nucleation and residual stress $\tau_r \approx \tau(D)$ when slip stops (negligible undershoot or overshoot). During seismic slip, while TP dominates, friction is primarily a function of slip $\delta$ (eqs 15 and 19 and Fig. 5), thus

$$ \tau(\delta) \approx \frac{f_1 \sigma}{(1 + \frac{\delta}{L^*})^{1/3}}, \qquad (22) $$

where $f_1$ is a reduced reference friction coefficient that accounts for the early dependence of $\tau$ on velocity due to R-S friction and FH. Numerical results without FH show that for $f_0 = 0.6$ and typical coseismic velocities in the range of 1–20 m s$^{-1}$, $f_1 \approx 0.51–0.53$. The logarithmic dependence on velocity in R-S friction does not produce a significant reduction of the steady-state friction. When including FH, $f_1 \approx f_w$. Fig. 5 shows that eq. (22) approximates well the slip-weakening curve resulting from dynamic rupture simulations with the TP proxy. Then, for final slip D:

$$ \tau_r = \tau(D) \approx \frac{f_1 \sigma}{(1 + \frac{D}{L^*})^{1/3}}. \qquad (23) $$

Using eqs (22) and (23), we solve the integral defining the breakdown energy in eq. (14) and obtain:

$$ G_c \approx \left\{ \frac{3}{2} L^* \left[ \left( 1 + \frac{D}{L^*} \right)^{2/3} - 1 \right] - \frac{D}{\left( 1 + \frac{D}{L^*} \right)^{1/3}} \right\} f_1 \sigma. \qquad (24) $$

For large slip, $D \gg L^*$, we obtain

$$ G_c \approx \frac{1}{2} L^{*1/3} D^{2/3} f_1 \sigma, \qquad (25) $$

which is the intended scaling $G_c \propto D^{2/3}$. Fig. 6 compares values of $G_c$ calculated from the stress and slip outputs of 2.5-D dynamic simulations to the theoretical estimates given by eqs (24) and (25). The simulations reach various values of final slip D and have two different values of $L^*$. The comparison verifies that eq. (24) is a very accurate approximation, implying that $G_c$ is dominated by the contribution from TP. For slip larger than 1 m, eq. (25) also predicts accurately $G_c$. Comparing Figs 6(a) and (b), we verify that this approximation improves at shorter values of $L^*$ for given D; it generally improves at larger values of $D/L^*$.

In some studies, the average stress-slip relation is inferred from seismological observations (Abercrombie & Rice 2005). This average relation can be different from the local relation evaluated at a particular point on the fault (Tinti et al. 2009). While Fig. 5 shows the stress-slip curve at a selected point on the fault, in our simulations the scaling between average shear stress and average slip is similar, because our 2.5-D models with rupture lengths much larger than the seismogenic width W produce almost-steady pulse-like ruptures, with almost uniform slip and shear stress along the fault.

## 5 DISCUSSION AND CONCLUSIONS

We have introduced a proxy model for TP on faults during earthquakes, implemented through a modified R-S friction law. This approach is designed to emulate the fault slip-dependent weakening and the relationship between breakdown energy and slip resulting from TP in the diffusion-dominated regime relevant at large slip, in a way that is amenable for earthquake cycle simulations. Compared to the classical R-S model, an additional state variable is introduced
<figure><img src="image_6.png" /><figcaption>Figure 5. Slip-weakening response and breakdown energy $G_c$ on an R-S fault with TP proxy, from numerical simulation and theoretical predictions.</figcaption></figure>

shear stress with slip. We observe stress undershoot in models that include FH, due to their fast healing process, which results in slip pulses with shorter rise time than in models without FH.

We further verify that our large earthquake models with TP are compatible with fracture mechanics theory. Large earthquakes with rupture lengths >W are pulse-like, due to the 2.5-D effect of the finite seismogenic width W (Fig. 4a). When the rupture has propagated over distances significantly larger than W, the slip rate decays with distance x behind the rupture front as $v \propto x^{-0.5}$ (red curve in Fig. 4b). This asymptotic behaviour emerges behind the peak slip velocity location, at distances sufficiently larger than the slip acceleration zone size (analogous to a cohesive zone or process zone size) and smaller than the slip pulse size. This square root decay is the same as in singular pulse models with uniform stress drop (Freund 1979; Nielsen & Madariaga 2003). This conventional singularity asymptotics supports the use of fracture mechanics theory to understand the propagation and arrest of very large earthquakes (larger than wide), despite the continued weakening induced by TP. Only when the rupture has propagated over distances shorter than or comparable to W, the slip rate decays as $v \propto x^{-0.25}$ (blue curve in Fig. 4b), which is the expected ‘unconventional singularity’ in 2-D rupture models with continued slip-weakening due to TP (Viesca & Garagash 2015; Brener & Bouchbinder 2021; Paglialunga et al. 2024). At further distance into the tail of the pulse, the slip rate decays exponentially (Fig. 4c), which is the expected behaviour induced by the finite rupture width W in 2.5-D models (fig. S6e and f of Weng & Ampuero 2022).

## 4 BREAKDOWN ENERGY OF THE TP PROXY

Our implementation is designed to satisfy the relation between $G_c$ and slip $\delta$ derived for TP in the diffusion-dominated regime at slips exceeding ~0.1 m by Viesca & Garagash (2015), which they validated by comparison to seismological observations. Here, we verify numerically that our implementation indeed produces the intended $G_c(\delta)$ relation. A different decay consistent with the undrained-adiabatic regime of TP at short slip, $G_c \propto \delta^2$ (Viesca & Garagash

2015), could be similarly introduced in the model by modifying eq. (15).

To estimate $G_c$, we consider an earthquake with initial shear stress $\tau_0$ right before nucleation and residual stress $\tau_r \approx \tau(D)$ when slip stops (negligible undershoot or overshoot). During seismic slip, while TP dominates, friction is primarily a function of slip $\delta$ (eqs 15 and 19 and Fig. 5), thus

$$ \tau(\delta) \approx \frac{f_1 \sigma}{(1 + \frac{\delta}{L^*})^{1/3}}, \qquad (22) $$

where $f_1$ is a reduced reference friction coefficient that accounts for the early dependence of $\tau$ on velocity due to R-S friction and FH. Numerical results without FH show that for $f_0 = 0.6$ and typical coseismic velocities in the range of 1–20 m s$^{-1}$, $f_1 \approx 0.51–0.53$. The logarithmic dependence on velocity in R-S friction does not produce a significant reduction of the steady-state friction. When including FH, $f_1 \approx f_w$. Fig. 5 shows that eq. (22) approximates well the slip-weakening curve resulting from dynamic rupture simulations with the TP proxy. Then, for final slip D:

$$ \tau_r = \tau(D) \approx \frac{f_1 \sigma}{(1 + \frac{D}{L^*})^{1/3}}. \qquad (23) $$

Using eqs (22) and (23), we solve the integral defining the breakdown energy in eq. (14) and obtain:

$$ G_c \approx \left\{ \frac{3}{2} L^* \left[ \left( 1 + \frac{D}{L^*} \right)^{2/3} - 1 \right] - \frac{D}{\left( 1 + \frac{D}{L^*} \right)^{1/3}} \right\} f_1 \sigma. \qquad (24) $$

For large slip, $D \gg L^*$, we obtain

$$ G_c \approx \frac{1}{2} L^{*1/3} D^{2/3} f_1 \sigma, \qquad (25) $$

which is the intended scaling $G_c \propto D^{2/3}$. Fig. 6 compares values of $G_c$ calculated from the stress and slip outputs of 2.5-D dynamic simulations to the theoretical estimates given by eqs (24) and (25). The simulations reach various values of final slip D and have two different values of $L^*$. The comparison verifies that eq. (24) is a very accurate approximation, implying that $G_c$ is dominated by the contribution from TP. For slip larger than 1 m, eq. (25) also predicts accurately $G_c$. Comparing Figs 6(a) and (b), we verify that this approximation improves at shorter values of $L^*$ for given D; it generally improves at larger values of $D/L^*$.

In some studies, the average stress-slip relation is inferred from seismological observations (Abercrombie & Rice 2005). This average relation can be different from the local relation evaluated at a particular point on the fault (Tinti et al. 2009). While Fig. 5 shows the stress-slip curve at a selected point on the fault, in our simulations the scaling between average shear stress and average slip is similar, because our 2.5-D models with rupture lengths much larger than the seismogenic width W produce almost-steady pulse-like ruptures, with almost uniform slip and shear stress along the fault.

## 5 DISCUSSION AND CONCLUSIONS

We have introduced a proxy model for TP on faults during earthquakes, implemented through a modified R-S friction law. This approach is designed to emulate the fault slip-dependent weakening and the relationship between breakdown energy and slip resulting from TP in the diffusion-dominated regime relevant at large slip, in a way that is amenable for earthquake cycle simulations. Compared to the classical R-S model, an additional state variable is introduced
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Breakdown energy $G_c$ on 2.5-D fully dynamic simulations, compared to theoretical estimates. Simulations span a range of values of averaged final slip $D$, obtained by setting different $W$ values from 1 to 200 km and $L$ values from 0.001 to 0.1 m. The $G_c$ values of simulations are computed from the shear stress and slip outputs along the fault, then averaged spatially. Two sets of simulations are shown, with (a) $L^* = 0.1$ m and (b) $L^* = 0.01$ m. Other model parameters are $\sigma = 25$ MPa, $f_0 = 0.6$, $V_0 = 10^{-6}$ m s$^{-1}$, $V_{th} = 10^{-3}$ m s$^{-1}$, $a = 0.01$, $b = 0.015$, $\mu = 32$ GPa and $c_s = 3464$ m s$^{-1}$.</figcaption></figure>

and one additional state evolution equation is solved, whose calculation time is negligible. Therefore, the additional complexity and computational cost are negligible, unlike complete implementations of TP, which require solving fluid pressure and thermal diffusion equations. Note that we do not aim to reproduce the full physics of TP with our proposed proxy model. By bypassing the diffusion equations for temperature and pressure, we chose not to model explicitly the temperature and pore pressure evolution during the seismic cycle, nor the normal stress changes induced by TP. Developing a different proxy model by implementing TP via changes of normal stress and comparing to the proxy introduced here could be a topic of future research. Our approach focuses on first order effects of weakening and energy dissipation associated with TP.

Although we demonstrate the concept in a 2.5-D earthquake cycle simulator based on the spectral element method, the TP proxy only involves modifications of the friction solver, thus it can be

readily implemented in 1-D, 2-D and 3-D, and in simulators based on other numerical methods, such as the finite-element method or the boundary element method. Moreover, with a similar implementation by modifying eq. (15), it is possible to reproduce different forms of slip-weakening, such as the undrained-adiabatic regime of TP that could be dominant at shorter slip values (Viesca & Garagash 2015).

Our proxy encapsulates the effects of TP in a single parameter, the length scale $L^*$ defined in eq. (16). It depends on the fault zone physical parameters involved in the thermal and hydraulic diffusion equations. Varying $L^*$, earthquake models can be tuned between strong and weak TP effects, or account for both cases on the same fault through spatial variations of $L^*$ (Noda & Lapusta 2013).

Carrying 2.5-D simulations of large earthquakes, we find that, despite the continued slip-weakening behind the rupture front produced by TP, the asymptotic behaviour of slip rate is consistent with the conventional singularity of fracture mechanics theory, as soon as
<figure><img src="image_7.png" /><figcaption>**Figure 6.** Breakdown energy $G_c$ on 2.5-D fully dynamic simulations, compared to theoretical estimates. Simulations span a range of values of averaged final slip $D$, obtained by setting different $W$ values from 1 to 200 km and $L$ values from 0.001 to 0.1 m. The $G_c$ values of simulations are computed from the shear stress and slip outputs along the fault, then averaged spatially. Two sets of simulations are shown, with (a) $L^* = 0.1$ m and (b) $L^* = 0.01$ m. Other model parameters are $\sigma = 25$ MPa, $f_0 = 0.6$, $V_0 = 10^{-6}$ m s$^{-1}$, $V_{th} = 10^{-3}$ m s$^{-1}$, $a = 0.01$, $b = 0.015$, $\mu = 32$ GPa and $c_s = 3464$ m s$^{-1}$.</figcaption></figure>

and one additional state evolution equation is solved, whose calculation time is negligible. Therefore, the additional complexity and computational cost are negligible, unlike complete implementations of TP, which require solving fluid pressure and thermal diffusion equations. Note that we do not aim to reproduce the full physics of TP with our proposed proxy model. By bypassing the diffusion equations for temperature and pressure, we chose not to model explicitly the temperature and pore pressure evolution during the seismic cycle, nor the normal stress changes induced by TP. Developing a different proxy model by implementing TP via changes of normal stress and comparing to the proxy introduced here could be a topic of future research. Our approach focuses on first order effects of weakening and energy dissipation associated with TP.

Although we demonstrate the concept in a 2.5-D earthquake cycle simulator based on the spectral element method, the TP proxy only involves modifications of the friction solver, thus it can be

readily implemented in 1-D, 2-D and 3-D, and in simulators based on other numerical methods, such as the finite-element method or the boundary element method. Moreover, with a similar implementation by modifying eq. (15), it is possible to reproduce different forms of slip-weakening, such as the undrained-adiabatic regime of TP that could be dominant at shorter slip values (Viesca & Garagash 2015).

Our proxy encapsulates the effects of TP in a single parameter, the length scale $L^*$ defined in eq. (16). It depends on the fault zone physical parameters involved in the thermal and hydraulic diffusion equations. Varying $L^*$, earthquake models can be tuned between strong and weak TP effects, or account for both cases on the same fault through spatial variations of $L^*$ (Noda & Lapusta 2013).

Carrying 2.5-D simulations of large earthquakes, we find that, despite the continued slip-weakening behind the rupture front produced by TP, the asymptotic behaviour of slip rate is consistent with the conventional singularity of fracture mechanics theory, as soon as
the rupture has propagated a distance larger than the rupture width. This result supports the applicability of fracture mechanics theory to understand the propagation and arrest of large earthquakes. Consequently, we posit that the breakdown energy derived from our model is equivalent to the fracture energy described in cohesive-zone models of fracture mechanics, which are a powerful framework to model the dynamics of very large earthquakes (Weng & Ampuero 2019, 2020).

We derive and verify an expression for the breakdown energy $G_c$ (eq. 24) in earthquake cycle models with TP, which agrees asymptotically with a previous theoretical prediction (Viesca & Garagash 2015) by design, and which can be used to evaluate the rupture potential of fault segments (Weng & Ampuero 2019, 2020). Thus, these results can have important implications for understanding earthquake mechanics.

## DATA AVAILABILITY

The code employed in this research is SEM2DPACK, an open-source spectral element method code, available to download in https://github.com/jpampuero/sem2dpack/tree/cycle_thermal pressurization. The parameters used to reproduce the results are described in the main text and in Supporting Information. For the calibration of the 2.5-D model in APPENDIX A, we use QDYN, a boundary element software, available to download in https://github.com/ydluo/qdyn.

## ACKNOWLEDGMENTS

MTH acknowledges the support by ANID 'Agencia Nacional de Investigación y Desarrollo Subdirección de Capital Humano/Doctorado Nacional/2020' National PhD scholarship 21201352. JPA was supported by the French government through the UCAJEDI Investments in the Future project (ANR-15-IDEX-01) managed by the National Research Agency (ANR). The authors are grateful to the OPAL infrastructure and the Université Côte d'Azur's Center for High-Performance Computing for providing resources and support. MTH and JGFC acknowledge the funding support by ANID/FONDAP/15110017 (Centro de Investigación para la Gestión Integrada del Riesgo de Desastres [CIGIDEN]). JGFC acknowledges support from ANID/FONDECYT/11201180 (Dynamic Rupture on Faults with Heterogeneous Frictional Properties).

## SUPPORTING INFORMATION

Supplementary data are available at GJI online. In the supporting information we included a comparison of the proxy implementation with the full implementation of thermal pressurization and the parameters used to reproduce the results shown in the main text.

### suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## REFERENCES

* Abercrombie, R.E. & Rice, J.R., 2005. Can observations of earthquake scaling constrain slip weakening?, *Geophys. J. Int.*, **162**(2), 406–424.

* Aki, K., 1966. Generation and propagation of G-waves from the Niigata earthquake of June 16, 1964, 2. Estimation of earthquake moment, released energy, and stress-strain drop from G-wave spectrum, *Bull. Earthq. Res. Inst. Tokyo Univ.*, **44**, 73–88.

* Allison, K.L. & Dunham, E.M., 2018. Earthquake cycle simulations with rate-and-state friction and power-law viscoelasticity, *Tectonophysics*, **733**, 232–256.

* Ampuero, J.-P., 2002. *Étude physique et numérique de la nucléation des séismes*, PhD thesis, University of Paris VII.

* Ampuero, J.-P., 2012. SEM2DPACK, a spectral element software for 2D seismic wave propagation and earthquake source dynamics, v2.3.8, Zenodo, doi:10.5281/zenodo.230363.

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

* Beeler, N.M., Tullis, T.E. & Goldsby, D.L., 2008. Constitutive relationships and physical basis of fault strength due to flash heating, *J. Geophys. Res.*, **113**(B1), B01401.

* Ben-Zion, Y. & Rice, J.R., 1997. Dynamic simulations of slip on a smooth fault in an elastic solid, *J. Geophys. Res.*, **102**(B8), 17771–17784.

* Brener, E. & Bouchbinder, E., 2021. Unconventional singularities and energy balance in frictional rupture, *Nat. Commun.*, **12**, 2585.

* Cattania, C., 2019. Complex earthquake sequences on simple faults, *Geophys. Res. Lett.*, **46**(17-18), 10384–10393.

* Dieterich, J., 1978. Time-dependent friction and the mechanics of stick-slip, *Pure Appl. Geophys.*, **116**(4), 790–806.

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

* Freund, L.B., 1979. The mechanics of dynamic shear crack propagation, *J. Geophys. Res.*, **84**(B5), 2199–2209.

* Gallovic, F., Valentova, L., Ampuero, J.-P. & Gabriel, A.-A., 2019. Bayesian dynamic finite-fault inversion: 2. Application to the 2016 mw 6.2 ama-trice, italy, earthquake, *J. Geophys. Res.: Solid Earth*, **124**(7), 6970–6988.

* Garagash, D.I., 2012. Seismic and aseismic slip pulses driven by thermal pressurization of pore fluid, *J. Geophys. Res.*, **117**(B4), B04314.

* Goldsby, D.L. & Tullis, T.E., 2011. Flash heating leads to low frictional strength of crustal rocks at earthquake slip rates, *Science*, **334**(6053), 216–218.

* Harris, R.A. et al., 2018. A suite of exercises for verifying dynamic earthquake rupture codes, *Seismol. Res. Lett.*, **89**(3), 1146–1162.

* 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.

* Kanamori, H. & Anderson, D., 1975. Theoretical basis of some empirical relations in seismology, *Bull. Seism. Soc. Am.*, **65**(5), 1073–1095.

* Kaneko, Y., Ampuero, J.-P. & Lapusta, N., 2011. Spectral-element simulations of long-term fault slip: effect of low-rigidity layers on earthquake-cycle dynamics, *J. Geophys. Res.*, **116**(B10), B10313.

* Kaneko, Y., Lapusta, N. & Ampuero, J.-P., 2008. Spectral element modeling of spontaneous earthquake rupture on rate and state faults: Effect of velocity-strengthening friction at shallow depths, *J. Geophys. Res.*, **113**(B9), B09317.

* Lambert, V. & Lapusta, N., 2023. Absolute stress levels in models of low-heat faults: Links to geophysical observables and differences for crack-like ruptures and self-healing pulses, *Earth Planet. Sci. Lett.*, **618**, 118277.

* Lambert, V., Lapusta, N. & Perry, S., 2021. Propagation of large earthquakes as self-healing pulses or mild cracks, *Nature*, **591**(7849), 252–258.

* Lapusta, N. & Liu, Y., 2009. Three-dimensional boundary integral modeling of spontaneous earthquake sequences and aseismic slip, *J. Geophys. Res.*, **114**(B9), B09303.

* Lapusta, N., Rice, J., Ben-Zion, Y. & Zheng, G., 2000. Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate-and state-dependent friction, *J. Geophys. Res.*, **105**(B10), 23765–23789.

* Liang, C., Ampuero, J.-P. & Muñoz, D., 2022. The paucity of supershear earthquakes on large faults governed by rate and state friction, *Geophys. Res. Lett.*, **49**(22), e2022GL099749.
the rupture has propagated a distance larger than the rupture width. This result supports the applicability of fracture mechanics theory to understand the propagation and arrest of large earthquakes. Consequently, we posit that the breakdown energy derived from our model is equivalent to the fracture energy described in cohesive-zone models of fracture mechanics, which are a powerful framework to model the dynamics of very large earthquakes (Weng & Ampuero 2019, 2020).

We derive and verify an expression for the breakdown energy $G_c$ (eq. 24) in earthquake cycle models with TP, which agrees asymptotically with a previous theoretical prediction (Viesca & Garagash 2015) by design, and which can be used to evaluate the rupture potential of fault segments (Weng & Ampuero 2019, 2020). Thus, these results can have important implications for understanding earthquake mechanics.

## DATA AVAILABILITY

The code employed in this research is SEM2DPACK, an open-source spectral element method code, available to download in https://github.com/jpampuero/sem2dpack/tree/cycle_thermal pressurization. The parameters used to reproduce the results are described in the main text and in Supporting Information. For the calibration of the 2.5-D model in APPENDIX A, we use QDYN, a boundary element software, available to download in https://github.com/ydluo/qdyn.

## ACKNOWLEDGMENTS

MTH acknowledges the support by ANID 'Agencia Nacional de Investigación y Desarrollo Subdirección de Capital Humano/Doctorado Nacional/2020' National PhD scholarship 21201352. JPA was supported by the French government through the UCAJEDI Investments in the Future project (ANR-15-IDEX-01) managed by the National Research Agency (ANR). The authors are grateful to the OPAL infrastructure and the Université Côte d'Azur's Center for High-Performance Computing for providing resources and support. MTH and JGFC acknowledge the funding support by ANID/FONDAP/15110017 (Centro de Investigación para la Gestión Integrada del Riesgo de Desastres [CIGIDEN]). JGFC acknowledges support from ANID/FONDECYT/11201180 (Dynamic Rupture on Faults with Heterogeneous Frictional Properties).

## SUPPORTING INFORMATION

Supplementary data are available at GJI online. In the supporting information we included a comparison of the proxy implementation with the full implementation of thermal pressurization and the parameters used to reproduce the results shown in the main text.

### suppl_data

Please note: Oxford University Press is not responsible for the content or functionality of any supporting materials supplied by the authors. Any queries (other than missing material) should be directed to the corresponding author for the paper.

## REFERENCES

* Abercrombie, R.E. & Rice, J.R., 2005. Can observations of earthquake scaling constrain slip weakening?, *Geophys. J. Int.*, **162**(2), 406–424.

* Aki, K., 1966. Generation and propagation of G-waves from the Niigata earthquake of June 16, 1964, 2. Estimation of earthquake moment, released energy, and stress-strain drop from G-wave spectrum, *Bull. Earthq. Res. Inst. Tokyo Univ.*, **44**, 73–88.

* Allison, K.L. & Dunham, E.M., 2018. Earthquake cycle simulations with rate-and-state friction and power-law viscoelasticity, *Tectonophysics*, **733**, 232–256.

* Ampuero, J.-P., 2002. *Étude physique et numérique de la nucléation des séismes*, PhD thesis, University of Paris VII.

* Ampuero, J.-P., 2012. SEM2DPACK, a spectral element software for 2D seismic wave propagation and earthquake source dynamics, v2.3.8, Zenodo, doi:10.5281/zenodo.230363.

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

* Beeler, N.M., Tullis, T.E. & Goldsby, D.L., 2008. Constitutive relationships and physical basis of fault strength due to flash heating, *J. Geophys. Res.*, **113**(B1), B01401.

* Ben-Zion, Y. & Rice, J.R., 1997. Dynamic simulations of slip on a smooth fault in an elastic solid, *J. Geophys. Res.*, **102**(B8), 17771–17784.

* Brener, E. & Bouchbinder, E., 2021. Unconventional singularities and energy balance in frictional rupture, *Nat. Commun.*, **12**, 2585.

* Cattania, C., 2019. Complex earthquake sequences on simple faults, *Geophys. Res. Lett.*, **46**(17-18), 10384–10393.

* Dieterich, J., 1978. Time-dependent friction and the mechanics of stick-slip, *Pure Appl. Geophys.*, **116**(4), 790–806.

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

* Freund, L.B., 1979. The mechanics of dynamic shear crack propagation, *J. Geophys. Res.*, **84**(B5), 2199–2209.

* Gallovic, F., Valentova, L., Ampuero, J.-P. & Gabriel, A.-A., 2019. Bayesian dynamic finite-fault inversion: 2. Application to the 2016 mw 6.2 ama-trice, italy, earthquake, *J. Geophys. Res.: Solid Earth*, **124**(7), 6970–6988.

* Garagash, D.I., 2012. Seismic and aseismic slip pulses driven by thermal pressurization of pore fluid, *J. Geophys. Res.*, **117**(B4), B04314.

* Goldsby, D.L. & Tullis, T.E., 2011. Flash heating leads to low frictional strength of crustal rocks at earthquake slip rates, *Science*, **334**(6053), 216–218.

* Harris, R.A. et al., 2018. A suite of exercises for verifying dynamic earthquake rupture codes, *Seismol. Res. Lett.*, **89**(3), 1146–1162.

* 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.

* Kanamori, H. & Anderson, D., 1975. Theoretical basis of some empirical relations in seismology, *Bull. Seism. Soc. Am.*, **65**(5), 1073–1095.

* Kaneko, Y., Ampuero, J.-P. & Lapusta, N., 2011. Spectral-element simulations of long-term fault slip: effect of low-rigidity layers on earthquake-cycle dynamics, *J. Geophys. Res.*, **116**(B10), B10313.

* Kaneko, Y., Lapusta, N. & Ampuero, J.-P., 2008. Spectral element modeling of spontaneous earthquake rupture on rate and state faults: Effect of velocity-strengthening friction at shallow depths, *J. Geophys. Res.*, **113**(B9), B09317.

* Lambert, V. & Lapusta, N., 2023. Absolute stress levels in models of low-heat faults: Links to geophysical observables and differences for crack-like ruptures and self-healing pulses, *Earth Planet. Sci. Lett.*, **618**, 118277.

* Lambert, V., Lapusta, N. & Perry, S., 2021. Propagation of large earthquakes as self-healing pulses or mild cracks, *Nature*, **591**(7849), 252–258.

* Lapusta, N. & Liu, Y., 2009. Three-dimensional boundary integral modeling of spontaneous earthquake sequences and aseismic slip, *J. Geophys. Res.*, **114**(B9), B09303.

* Lapusta, N., Rice, J., Ben-Zion, Y. & Zheng, G., 2000. Elastodynamic analysis for slow tectonic loading with spontaneous rupture episodes on faults with rate-and state-dependent friction, *J. Geophys. Res.*, **105**(B10), 23765–23789.

* Liang, C., Ampuero, J.-P. & Muñoz, D., 2022. The paucity of supershear earthquakes on large faults governed by rate and state friction, *Geophys. Res. Lett.*, **49**(22), e2022GL099749.
Luo, Y., Ampuero, J.-P., Galvez, P., van den Ende, M. & Idini, B., 2017. QDYN: a Quasi-DYNamic earthquake simulator (v1.1), Zenodo, doi:10.5281/zenodo.322459.

Mase, C.W. & Smith, L., 1987. Effects of frictional heating on the thermal, hydrologic, and mechanical response of a fault, *J. Geophys. Res.*, **92**(B7), 6249–6272.

Mavrommatis, A.P., Segall, P. & Johnson, K., 2017. A physical model for interseismic erosion of locked fault asperities, *J. Geophys. Res.: Solid Earth*, **122**(10),8326–8346.

Nielsen, S. & Madariaga, R., 2003. On the self-healing fracture mode, *Bull. Seism. Soc. Am.*, **93**(6), 2375–2388.

Nielsen, S., Spagnuolo, E., Smith, S. A.F., Violay, M., Di Toro, G. & Bis-tacchi, A., 2016. Scaling in natural and laboratory earthquakes, *Geophys. Res. Lett.*, **43**(4), 1504–1510.

Noda, H. & Lapusta, N., 2010. Three-dimensional earthquake sequence simulations with evolving temperature and pore pressure due to shear heating: Effect of heterogeneous hydraulic diffusivity, *J. geophys. Res.*, **115**(B12), B12314.

Noda, H. & Lapusta, N., 2013. Stable creeping fault segments can become destructive as a result of dynamic weakening, *Nature*, **493**(7433), 518–521.

Noda, H., Lapusta, N. & Rice, J., 2011. Earthquake sequence calculations with dynamic weakening mechanisms, in Borja, R.I.(ed.), *Multiscale and Multiphysics Processes in Geomechanics, Springer Series in Geomechanics and Geoengineering*, Springer, Berlin, Heidelberg, 149–152.

Paglialunga, F., Passelègue, F., Lebihain, M. & Violay, M., 2024. Frictional weakening leads to unconventional singularities during dynamic rupture propagation, *Earth Planet. Sci. Lett.*, **626**, 118550.

Perry, S., Lambert, V. & Lapusta, N., 2020. Nearly magnitude-invariant stress drops in simulated crack-like earthquake sequences on rate-and-state faults with thermal pressurization of pore fluids, *J. Geophys. Res.: Solid Earth*, **125**(3), e2019JB018597.

Rempel, A.W. & Rice, J.R., 2006. Thermal pressurization and onset of melting in fault zones, *J. Geophys. Res.*, **111**(B9), B09314.

Rice, J. & Ruina, A., 1983. Stability of Steady Frictional Slipping, *J. Appl. Mech.*, **50**(2), 343–349.

Rice, J.R., 2006. Heating and weakening of faults during earthquake slip, *J. Geophys. Res.*, **111**(B05), B05311.

Rubin, A. & Ampuero, J.-P., 2005. Earthquake nucleation on (aging) rate and state faults, *J. Geophys. Res.*, **110**(B11), B11312.

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

Starr, A., 1928. Slip in a crystal and rupture in a solid due to shear, *Math. Proc. Cambridge Phil. Soc.*, **24**(4), 489–500.

Thomas, M.Y., Lapusta, N. & Avouac, J.-P., 2014. Quasi-dynamic versus fully dynamic simulations of earthquakes and aseismic slip with and without enhanced coseismic weakening, *J. Geophys. Res.: Solid Earth*, **119**(3), 1986–2004.

Tinti, E., Cocco, M., Fukuyama, E. & Piatanesi, A., 2009. Dependence of slip weakening distance (Dc) on final slip during dynamic rupture of earthquakes, *Geophys. J. Int.*, **177**(3), 1205–1220.

Tinti, E., Spudich, P. & Cocco, M., 2005. Earthquake fracture energy inferred from kinematic rupture models on extended faults, *J. Geophys. Res.*, **110**(B11), B12303.

Viesca, R. & Garagash, D., 2015. Ubiquitous weakening of faults due to thermal pressurization, *Nat. Geosci.*, **8**(11), 875–879.

Weng, H. & Ampuero, J.-P., 2019. The dynamics of elongated earthquake ruptures, *J. Geophys. Res.: Solid Earth*, **124**(8), 8584–8610.

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

Weng, H. & Ampuero, J.-P., 2022. Integrated rupture mechanics for slow slip events and earthquakes., *Nat. Commun.*, **13**, 7327.

Weng, H. & Yang, H., 2018. Constraining frictional properties on fault by dynamic rupture simulations and near-field observations, *J. Geophys. Res.: Solid Earth*, **123**(8), 6658–6670.

## APPENDIX A: CALIBRATION OF THE 2.5-D MODEL

The slip in 2.5-D models corresponds to the peak slip across the seismogenic depth of a 3-D model. The slip on a deeply buried fault with uniform stress drop can be crudely approximated as one-half of a cosine of wavelength $2W$. Similarly, in a shallow fault on a half-space, the slip is maximal at the surface and zero at the bottom of the rupture, and the depth profile can be approximated as one quarter of a cosine of wavelength $4W$ (Weng & Ampuero 2019). However, as these approximations are crude, the $W_{2.5D}$ used in the 2.5-D code might differ from an equivalent $W_{3D}$ in 3-D. Here, we determine the relation between $W_{2.5D}$ and $W_{3D}$ in such a way that 2.5-D and 3-D models with same stress drop produce the same slip.

For the 3-D model, a theoretical relation between stress drop $\Delta\tau$ and slip is (Kanamori & Anderson 1975):

$$ \Delta\tau = C\mu\frac{\bar{\delta}}{L}. \qquad (A.1) $$

For a shallow infinitely long dip-slip fault, $\bar{L} = W$ (width) and:

$$ C = \frac{4(\lambda + \mu)}{\pi(\lambda + 2\mu)}, \qquad (A.2) $$

where $\lambda$ is the Lamé constant (Starr 1928; Aki 1966). Then, we obtain:

$$ \frac{\Delta\tau}{\bar{\delta}} = \frac{\frac{4(\lambda+\mu)}{\pi(\lambda+2\mu)}\mu}{W_{3D}}. \qquad (A.3) $$

From fig. 1 of Kanamori & Anderson (1975), noting that average slip $\bar{\delta}$ and peak slip $D_{max}$ in a crack are related by $\bar{\delta} = \frac{\pi}{4} D_{max}$:

$$ \frac{\Delta\tau}{D_{max}} = \frac{\frac{(\lambda+\mu)}{(\lambda+2\mu)}\mu}{W_{3D}}. \qquad (A.4) $$

Replacing $\lambda = \frac{2\mu\nu}{1-2\nu}$, we get:

$$ D_{max} = \frac{2(1 - \nu)\Delta\tau W_{3D}}{\mu}. \qquad (A.5) $$

For the 2.5-D model without free-surface effect:

$$ D = \frac{\Delta\tau W_{2.5D}}{\pi \mu}. \qquad (A.6) $$

As $D$ in the 2.5-D model corresponds to $D_{max}$ in a 3-D equivalent model with similar $\Delta\tau$, equating eqs (A.5) and (A.6) gives:

$$ W_{2.5D} = 2\pi(1 - \nu)W_{3D}. \qquad (A.7) $$

For $\nu = 1/4$, we must set $W_{2.5D} = 4.71 W_{3D}$ to obtain an equivalent model for a shallow fault. For a buried fault, following a similar procedure we obtain $W_{2.5D} = 2.36 W_{3D}$. Both values were validated numerically (see Fig. A1) by comparing 3-D simulations with 2.5-D simulations using QDYN, a Quasi-DYNamic boundary element method simulator (Luo et al. 2017).
Luo, Y., Ampuero, J.-P., Galvez, P., van den Ende, M. & Idini, B., 2017. QDYN: a Quasi-DYNamic earthquake simulator (v1.1), Zenodo, doi:10.5281/zenodo.322459.

Mase, C.W. & Smith, L., 1987. Effects of frictional heating on the thermal, hydrologic, and mechanical response of a fault, *J. Geophys. Res.*, **92**(B7), 6249–6272.

Mavrommatis, A.P., Segall, P. & Johnson, K., 2017. A physical model for interseismic erosion of locked fault asperities, *J. Geophys. Res.: Solid Earth*, **122**(10),8326–8346.

Nielsen, S. & Madariaga, R., 2003. On the self-healing fracture mode, *Bull. Seism. Soc. Am.*, **93**(6), 2375–2388.

Nielsen, S., Spagnuolo, E., Smith, S. A.F., Violay, M., Di Toro, G. & Bis-tacchi, A., 2016. Scaling in natural and laboratory earthquakes, *Geophys. Res. Lett.*, **43**(4), 1504–1510.

Noda, H. & Lapusta, N., 2010. Three-dimensional earthquake sequence simulations with evolving temperature and pore pressure due to shear heating: Effect of heterogeneous hydraulic diffusivity, *J. geophys. Res.*, **115**(B12), B12314.

Noda, H. & Lapusta, N., 2013. Stable creeping fault segments can become destructive as a result of dynamic weakening, *Nature*, **493**(7433), 518–521.

Noda, H., Lapusta, N. & Rice, J., 2011. Earthquake sequence calculations with dynamic weakening mechanisms, in Borja, R.I.(ed.), *Multiscale and Multiphysics Processes in Geomechanics, Springer Series in Geomechanics and Geoengineering*, Springer, Berlin, Heidelberg, 149–152.

Paglialunga, F., Passelègue, F., Lebihain, M. & Violay, M., 2024. Frictional weakening leads to unconventional singularities during dynamic rupture propagation, *Earth Planet. Sci. Lett.*, **626**, 118550.

Perry, S., Lambert, V. & Lapusta, N., 2020. Nearly magnitude-invariant stress drops in simulated crack-like earthquake sequences on rate-and-state faults with thermal pressurization of pore fluids, *J. Geophys. Res.: Solid Earth*, **125**(3), e2019JB018597.

Rempel, A.W. & Rice, J.R., 2006. Thermal pressurization and onset of melting in fault zones, *J. Geophys. Res.*, **111**(B9), B09314.

Rice, J. & Ruina, A., 1983. Stability of Steady Frictional Slipping, *J. Appl. Mech.*, **50**(2), 343–349.

Rice, J.R., 2006. Heating and weakening of faults during earthquake slip, *J. Geophys. Res.*, **111**(B05), B05311.

Rubin, A. & Ampuero, J.-P., 2005. Earthquake nucleation on (aging) rate and state faults, *J. Geophys. Res.*, **110**(B11), B11312.

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

Starr, A., 1928. Slip in a crystal and rupture in a solid due to shear, *Math. Proc. Cambridge Phil. Soc.*, **24**(4), 489–500.

Thomas, M.Y., Lapusta, N. & Avouac, J.-P., 2014. Quasi-dynamic versus fully dynamic simulations of earthquakes and aseismic slip with and without enhanced coseismic weakening, *J. Geophys. Res.: Solid Earth*, **119**(3), 1986–2004.

Tinti, E., Cocco, M., Fukuyama, E. & Piatanesi, A., 2009. Dependence of slip weakening distance (Dc) on final slip during dynamic rupture of earthquakes, *Geophys. J. Int.*, **177**(3), 1205–1220.

Tinti, E., Spudich, P. & Cocco, M., 2005. Earthquake fracture energy inferred from kinematic rupture models on extended faults, *J. Geophys. Res.*, **110**(B11), B12303.

Viesca, R. & Garagash, D., 2015. Ubiquitous weakening of faults due to thermal pressurization, *Nat. Geosci.*, **8**(11), 875–879.

Weng, H. & Ampuero, J.-P., 2019. The dynamics of elongated earthquake ruptures, *J. Geophys. Res.: Solid Earth*, **124**(8), 8584–8610.

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

Weng, H. & Ampuero, J.-P., 2022. Integrated rupture mechanics for slow slip events and earthquakes., *Nat. Commun.*, **13**, 7327.

Weng, H. & Yang, H., 2018. Constraining frictional properties on fault by dynamic rupture simulations and near-field observations, *J. Geophys. Res.: Solid Earth*, **123**(8), 6658–6670.

## APPENDIX A: CALIBRATION OF THE 2.5-D MODEL

The slip in 2.5-D models corresponds to the peak slip across the seismogenic depth of a 3-D model. The slip on a deeply buried fault with uniform stress drop can be crudely approximated as one-half of a cosine of wavelength $2W$. Similarly, in a shallow fault on a half-space, the slip is maximal at the surface and zero at the bottom of the rupture, and the depth profile can be approximated as one quarter of a cosine of wavelength $4W$ (Weng & Ampuero 2019). However, as these approximations are crude, the $W_{2.5D}$ used in the 2.5-D code might differ from an equivalent $W_{3D}$ in 3-D. Here, we determine the relation between $W_{2.5D}$ and $W_{3D}$ in such a way that 2.5-D and 3-D models with same stress drop produce the same slip.

For the 3-D model, a theoretical relation between stress drop $\Delta\tau$ and slip is (Kanamori & Anderson 1975):

$$ \Delta\tau = C\mu\frac{\bar{\delta}}{L}. \qquad (A.1) $$

For a shallow infinitely long dip-slip fault, $\bar{L} = W$ (width) and:

$$ C = \frac{4(\lambda + \mu)}{\pi(\lambda + 2\mu)}, \qquad (A.2) $$

where $\lambda$ is the Lamé constant (Starr 1928; Aki 1966). Then, we obtain:

$$ \frac{\Delta\tau}{\bar{\delta}} = \frac{\frac{4(\lambda+\mu)}{\pi(\lambda+2\mu)}\mu}{W_{3D}}. \qquad (A.3) $$

From fig. 1 of Kanamori & Anderson (1975), noting that average slip $\bar{\delta}$ and peak slip $D_{max}$ in a crack are related by $\bar{\delta} = \frac{\pi}{4} D_{max}$:

$$ \frac{\Delta\tau}{D_{max}} = \frac{\frac{(\lambda+\mu)}{(\lambda+2\mu)}\mu}{W_{3D}}. \qquad (A.4) $$

Replacing $\lambda = \frac{2\mu\nu}{1-2\nu}$, we get:

$$ D_{max} = \frac{2(1 - \nu)\Delta\tau W_{3D}}{\mu}. \qquad (A.5) $$

For the 2.5-D model without free-surface effect:

$$ D = \frac{\Delta\tau W_{2.5D}}{\pi \mu}. \qquad (A.6) $$

As $D$ in the 2.5-D model corresponds to $D_{max}$ in a 3-D equivalent model with similar $\Delta\tau$, equating eqs (A.5) and (A.6) gives:

$$ W_{2.5D} = 2\pi(1 - \nu)W_{3D}. \qquad (A.7) $$

For $\nu = 1/4$, we must set $W_{2.5D} = 4.71 W_{3D}$ to obtain an equivalent model for a shallow fault. For a buried fault, following a similar procedure we obtain $W_{2.5D} = 2.36 W_{3D}$. Both values were validated numerically (see Fig. A1) by comparing 3-D simulations with 2.5-D simulations using QDYN, a Quasi-DYNamic boundary element method simulator (Luo et al. 2017).
<figure><img src="image_8.png" /><figcaption>**Figure A1.** Calibration of *W* for 2.5-D and 3-D models, based on quasi-dynamic simulations. (a) Shallow fault. (b) Slightly buried fault. (c) Deeply buried fault. Parameters used are $V_{PL} = 10^{-9}$ m s$^{-1}$, $\sigma = 27.5$ MPa, $L = 0.06$ m, $a = 0.01$, $b_{VW} = 0.015$, $b_{VS} = 0.005$, $V_0 = 10^{-6}$ m s$^{-1}$, $\mu = 32$ GPa, $c_s = 3000$ m s$^{-1}$ and $W_{3D} = 100$ km.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>**Figure A1.** Calibration of *W* for 2.5-D and 3-D models, based on quasi-dynamic simulations. (a) Shallow fault. (b) Slightly buried fault. (c) Deeply buried fault. Parameters used are $V_{PL} = 10^{-9}$ m s$^{-1}$, $\sigma = 27.5$ MPa, $L = 0.06$ m, $a = 0.01$, $b_{VW} = 0.015$, $b_{VS} = 0.005$, $V_0 = 10^{-6}$ m s$^{-1}$, $\mu = 32$ GPa, $c_s = 3000$ m s$^{-1}$ and $W_{3D} = 100$ km.</figcaption></figure>