# Shear localisation controls the dynamics of earthquakes

Fabian Barras1 and Nicolas Brantut2
1 The Njord Centre, University of Oslo, Oslo, Norway
2 Department of Earth Sciences, University College London, London, UK

( )

###### Abstract

Earthquakes are produced by the propagation of rapid slip along tectonic
faults. The propagation dynamics is governed by a balance between elastic
stored energy in the surrounding rock, and dissipated energy at the
propagating tip of the slipping patch. Energy dissipation is dictated by the
mechanical behaviour of the fault, which is itself the result of feedbacks
between thermo-hydro-mechanical processes acting at the mm to sub-mm scale.
Here, we numerically simulate shear ruptures using a dual scale approach,
allowing us to couple a sub-mm description of inner fault processes and km-
scale elastodynamics, and show that the sudden localisation of shear strain
within a shear zone leads to the emergence of classical cracks driven by a
constant fracture energy. The fracture energy associated to strain
localisation is substantially smaller than that predicted assuming uniform
shearing. We show the existence of a unique scaling law between the localised
shearing width and the rupture speed. Our results indicate that earthquakes
are likely to be systematically associated to extreme strain localisation.

##  1 Introduction

Earthquake sources correspond to slip events dynamically propagating along
faults. At crustal scale, faults can be viewed as two-dimensional surfaces,
across which the displacement field is discontinuous. However, geological and
experimental observations show that “slip” across faults is the result of
shear deformation across narrow layers of highly comminuted, transformed or
partially melted rocks. In the shallow continental crust, fault core materials
are often made of fine-grained siliclastic and clay gouges, with a porosity
filled with pressurised water (e.g. Scholz, 1988; Rice, 2006). The dynamics of
ruptures in crustal faults is controlled by the rheology of these water-
saturated fault gouges.

During earthquakes, faults slide at elevated slip rates of the order of metres
per second, which leads to dramatic weakening of fault gouge materials
(Scholz, 2019, chap. 2). In dry materials, weakening is most likely controlled
by the local temperature rise arising from dissipation of frictional work,
combined with thermally activated rheology of the rock-forming minerals (e.g
Rice, 2006; Beeler et al., 2008; Proctor et al., 2014; De Paola et al., 2015;
Yao et al., 2015; Pozzi et al., 2021; Harbord et al., 2021). In the presence
of fluids, an additional weakening mechanism is expected, due to the
differential thermal expansion of the pore fluids and the solid pore space:
upon heating, the fluid pressure rises, effective normal stress decreases and
the frictional strength drops. This so-called “thermal pressurisation”
mechanism, initially proposed by Sibson (1975) as a temperature-limiting
process in wet rocks, has been shown to produce realistic predictions for
thermal evolution and energy dissipation during earthquakes (e.g. Rice, 2006;
Viesca and Garagash, 2015), and is a potential candidate to explain some of
the complexity observed in natural earthquakes (e.g. Noda and Lapusta, 2013)
and the operation of plate boundary faults at low ambient stress (e.g. Noda et
al., 2009; Lambert et al., 2021).

The thickness of the actively deforming zone determines the shear heating rate
and how easily fluids and heat diffuse away from the fault plane, and thus has
a tremendous influence on the resulting rupture dynamics (Andrews, 2002; Noda
et al., 2009; Noda and Lapusta, 2010; Viesca and Garagash, 2015; Lambert et
al., 2021). While geological and experimental observations can be used to
constrain the thickness of actively deforming fault gouge material, the range
of acceptable values spans more than 3 orders of magnitude, from fractions of
millimetres to centimetres (Rice, 2006), and it is one of the key unknown that
limits our ability to determine the efficiency of thermal weakening mechanisms
in nature.

The influence of shear zone width on earthquake propagation is further
complicated by the fact that this parameter is likely evolving during seismic
slip: strain localisation is expected to be part of the fault weakening
process. Several mechanisms might be responsible for strain localisation
during earthquake slip, including granular rearrangements and grain size
reduction (e.g., Mair and Abe, 2008; Hermundstad et al., 2010), shear heating
coupled to thermal weakening (e.g., Braeck and Podladchikov, 2007), thermal
pressurisation (e.g., Sulem et al., 2011; Rice et al., 2014; Platt et al.,
2014) and thermal decomposition (e.g., Veveakis et al., 2012; Platt et al.,
2015). In all cases, the strain localisation process is associated to a rapid
reduction in shear strength, and we therefore expect strain localisation to
exert a strong control on the overall dynamics of rupture.

Here, we demonstrate and quantify how strain localisation impacts rupture
dynamics: we run dynamic rupture simulations and compute the fracture energy
associated with the localisation process, and find a relationship between the
rupture speed and the degree of strain localisation within the fault gouge. We
use the case of thermal pressurisation as a representative thermal weakening
process that is compatible with seismological observations (Rice, 2006; Viesca
and Garagash, 2015; Lambert et al., 2021), and is known to spontaneously lead
to strain localisation (Rice et al., 2014; Platt et al., 2014). We argue that
the interplay between rupture dynamics and strain localisation analysed here
applies to most thermal weakening processes in rocks and engineering
materials.

Figure 1: Schematic of the dual scale setup governing the propagation of
localised shear bands. The dynamic rupture extends over large (kilometric)
scales along the fault ($x-z$ plane), whereas the frictional strength is
determined by solving a coupled diffusion problem, such as thermal
pressurisation in this paper, where strain rate spontaneously evolves over
submillimetre scales across the fault gouge (in the $y$ direction).

##  2 Shear localisation and faulting

As a general starting point for our analysis, let us consider slip on
geological fault as the deformation distributed across a narrow pre-existing
shear zone (e.g., a fault gouge made of unconsolidated grains), and examine
the conditions under which shear strain can be further localised within this
pre-existing zone. Let us assume that the shear stress $\tau$ across the shear
zone is function of a set of variables that includes the shear strain rate
$\dot{\gamma}$ and another diffusive quantity $\vartheta$:

$\tau\sim f(\dot{\gamma},\vartheta),$ (1)

and that the rate of work produced by shearing acts as a source term in the
diffusion of $\vartheta$:

$\dot{\vartheta}=\beta\tau\dot{\gamma}+\alpha\nabla^{2}\vartheta,$ (2)

where $\nabla^{2}$ denotes the Laplace operator, $\alpha$ is a diffusivity and
$\beta$ is analogous to the Taylor-Quiney coefficient (Taylor and Quinney,
1934) for the cases when $\vartheta$ corresponds to temperature. In relation
to the first condition, one can define

$g^{\prime}(\dot{\gamma},\vartheta)=\frac{\partial f}{\partial\dot{\gamma}},$
(3)

which describes the rate-dependent rheology of the material. Natural examples
include viscous creep of rocks at elevated temperature, granular material
rheology (e.g. Jop et al., 2006) and rate-and-state friction. Similarly, one
can define

$h^{\prime}(\dot{\gamma},\vartheta)=\frac{\partial f}{\partial\vartheta},$ (4)

to describe the effect of $\vartheta$ on the material rheology. In practice,
this diffusive quantity $\vartheta$ will often correspond to temperature and
$h^{\prime}$ describe thermal weakening effects. $\vartheta$ could also
correspond to fluid pressure in a porous material whose strength is reduced by
an increase in pore fluid pressure (following the concept of effective stress
discussed later in Equation 7). It can also account for the combined effect of
pressure and temperature as in the case of thermal pressurisation that will be
discussed later in this manuscript. If conditions (1) and (2) are met, a
linear stability analysis (detailed in Appendix A) demonstrates that uniform
shearing at a given time $t=t_{0}$ becomes unstable if:

$\frac{h^{\prime}_{0}}{g^{\prime}_{0}}<0,$ (5)

with
$\\{f_{0},g^{\prime}_{0},h^{\prime}_{0}\\}=\\{f,g^{\prime},h^{\prime}\\}|_{t=t_{0}}$.
Moreover, the analysis also shows that only perturbation wavelengths $\lambda$
greater than a critical wavelength are unstable:

$\lambda>2\pi\sqrt{-\frac{\alpha}{\beta
f_{0}}\frac{g^{\prime}_{0}}{h^{\prime}_{0}}}\equiv\lambda_{\rm c},$ (6)

which indicates that such instability leads to the localisation of shear
strain down to some thickness $W_{\rm loc}\sim\lambda_{\rm c}/2$. Remarkably,
this type of localisation instability can also arise within rate-strengthening
materials ($g^{\prime}_{0}>0$) providing that $h^{\prime}_{0}<0$, as it is
often the case with thermal weakening mechanisms. As a result, shear flow
concentrates over a thickness much smaller than the initial width of the shear
zone and leads to a substantive drop of the associated shear stress. Examples
of this strain localisation instability have been described by Rice et al.
(2014); Platt et al. (2014) in the context of thermal pressurisation in
crustal faults, as well as Bai (1982) in the context of adiabatic shear
banding in metals. In this work, we quantitatively investigate how strain
localisation across the shear zone drives the rapid acceleration of slip and
the propagation of rupture front along the fault plane during earthquakes.

##  3 Model

Here, we analyse the process of thermal pressurisation, which has been shown
to be a realistic dynamic weakening mechanism (e.g. Rice, 2006; Viesca and
Garagash, 2015; Lambert et al., 2021) and that undergoes the localisation
instability outlined above (Rice et al., 2014; Platt et al., 2014).

In this case, the diffusive variable $\vartheta$ corresponds to pore fluid
pressure $p$ that affects the effective normal stress $(\sigma_{\rm n}-p)$ in
the shear zone and, thereby, its shear strength together with a rate-dependent
friction coefficient:

$\tau(\dot{\gamma},p)=f_{\rm rsf}(\dot{\gamma})(\sigma_{\rm n}-p).$ (7)

In Equation (7), we adopt the rate-strengthening rheology $f_{\rm
rsf}(\dot{\gamma})$ of Platt et al. (2014) detailed in Appendix B.2. Moreover,
fluid diffusion across the shear zone is governed by a coupled system of
thermal and hydraulic diffusion equations (see Equations 31 in Appendix). This
thermo-hydraulic coupling is caused by the different compressibilities and
thermal expansivities of the solid matrix and the pore fluid, and describes
the change in pore pressure produced by local temperature variations in the
gouge.

Exploring the interplay between strain localisation and rupture dynamics is a
challenging dual-scale problem: it requires solving for heat and fluid
diffusion at the scale of the fault core (from millimeters to centimeters in
natural fault zones) together with the elastodynamics governing the
propagation of the earthquake rupture along the fault (elastic waves moving at
kilometers per second in crustal rocks). We follow Noda et al. (2009) and take
advantage of the separation of scale to solve thermal pressurisation only
across the fault (along the $y$ axis in Figure 1). In this one-dimensional
configuration, Platt et al. (2014) found numerically that shear localisation,
under imposed constant slip rate $V$ across the gouge, stabilises to a finite
width that is well approximated by

$W_{\rm rsf}(V)\simeq 6.9\frac{A\rho c}{\Lambda f_{\rm
c}}\frac{(\sqrt{\alpha_{\rm hy}}+\sqrt{\alpha_{\rm th}})^{2}}{V(f_{\rm
c}+2A)},$ (8)

where $\alpha_{\rm hy,th}$ correspond respectively to the hydraulic and
thermal diffusivities, $\rho c$ is the heat capacity of the gouge material,
$\Lambda$ is the thermal pressurisation parameter that describes the change of
pore fluid pressure caused by an increase in temperature in the gouge. The
characteristic shear strength

$\tau_{\rm c}=f_{\rm c}(\sigma_{\rm n}-p_{0})=f_{\rm rsf}(V/h)(\sigma_{\rm
n}-p_{0})$ (9)

is a function of the initial uniform strain rate $\dot{\gamma}=V/h$ and
background pore pressure $p_{0}$. In Equation (8), the constant $A$
corresponds to the rate-strengthening coefficient that describes the ”direct”
response of the shear zone to a change in strain rate similar to standard
rate-and-state models (see Equation (32) in Appendix). Once the localisation
width $W_{\rm rsf}$ is much smaller that the diffusion length scale, thermal
pressurisation can be well approximated by a slip-on-a-plane solution assuming
$\dot{\gamma}(y)=V\bar{\delta}(y)$, with $\bar{\delta}$ being the Dirac delta
distribution. In this situation, the shear stress across the shear zone is
only function of the imposed slip rate $V$ and the accumulated slip
$\delta=Vt$ (Rice, 2006; Mase and Smith, 1987):

$\tau_{\rm sp}(\delta;V)=\tau_{\rm
c}\exp\Big{(}\frac{\delta}{L^{*}}\Big{)}\mathrm{erfc}\Big{(}\sqrt{\frac{\delta}{L^{*}}}\Big{)},$
(10)

where

$L^{*}(V)=\frac{4}{f^{2}_{\rm c}}\Big{(}\frac{\rho
c}{\Lambda}\Big{)}^{2}\frac{(\sqrt{\alpha_{\rm hy}}+\sqrt{\alpha_{\rm
th}})^{2}}{V}.$ (11)

During earthquakes, slip rate across the shear zone is far from being constant
and evolves rapidly near the tip of the propagating rupture. Next, we aim to
analyse the coupling between strain localisation, slip acceleration and
rupture dynamics in a simple faulting geometry that is sufficient to capture
its key physical aspects. We consider a planar fault within an infinite linear
elastic medium sliding in anti-plane shear (mode-III). In this configuration
shown in Figure 1, the shear traction at the fault $\tau(x,t)$ can be related
to the slip rate across the fault $V(x,t)$ and, thereby, the strain rate in
the shear zone, following

$\tau(x,t)=\tau_{\rm b}-\frac{\mu}{2c_{\rm s}}V(x,t)+\phi(x,t)=\tau_{\rm
b}-\frac{\mu}{2c_{\rm
s}}\int_{-h/2}^{h/2}\dot{\gamma}(x,y,t)\>\mathrm{d}y+\phi(x,t).$ (12)

In the equation above, $\mu$ is the shear modulus, $c_{\rm s}$ the shear wave
speed of the linear elastic medium surrounding the shear zone, $\phi$ is the
non-local dynamic contribution accounting for the history of slip along the
interface, and $\tau_{\rm b}$ represents the far-field background stress.
Equation (12) allows us to define the characteristic seismic slip rate $V_{\rm
c}$ and associated uniform strain rate $\dot{\gamma}_{\rm c}$ as

$\dot{\gamma}_{\rm c}=\frac{V_{\rm c}}{h}=\frac{2c_{\rm s}\tau_{\rm
b}}{h\mu},$ (13)

which are used in the remainder of the paper together with the related
characteristic shear strength $\tau_{\rm c}$ (9). The elastodynamic equation
(12) couples the strain rate in the shear zone $\dot{\gamma}$ to the shear
stress $\tau$ and allows us to implement a dual-scale coupled numerical scheme
that solves the rupture elastodynamics along the shear zone together with
pressure and temperature diffusion across the shear zone. The details of our
coupled numerical scheme are given in Appendix D.

##  4 Results

Figure 2: Dynamic rupture driven by shear localisation simulated with the
coupled model. The top panels a) and b) respectively present snapshots at
different times of the longitudinal profile of slip rate and shear stress
during which the rupture accelerates from sixty to about ninety percent of the
shear wave velocity. Note that the simulated domain is symmetric with respect
to the nucleation position $x=0$ such that another rupture tip moves toward
the negative positions. The bottom panels c) and d) present the profile of
strain rate $\dot{\gamma}$, pressure $p$ and temperature $T$ at the successive
positions of the rupture tip highlighted by black dots in panel a) and b). See
Appendix C and Table 1 for further details on the dimensional analysis behind
this coupled problem and the dimensionless scaling used to plot the data in
the different panels.

In our simulations, the shear zone is initially creeping at aseismic slip
velocity and, at time $t=0$, failure is nucleated by rising the pore pressure
near the center of the fault $x=0$ (further details of nucleation procedure
and parameter values are given in Appendix D). Initially, acceleration of slip
is mostly concentrated in the nucleation region, followed by a rapid lateral
rupture propagation whereby the slip rate increases in an expanding region
beyond the initial nucleation patch, concomitantly with a shear stress drop
linked with thermal pressurisation of pore fluids inside the gouge and intense
strain localisation (Figure 2). Rupture acceleration coincides with larger
slip velocities and stress drop at the tip (Figure 3a-c) and more intense
localisation of shear deformation across the gouge where up to four orders of
magnitude larger strain rate concentrates on less than five percent of the
thickness of the shear zone (Figure 3b). Interestingly, the peak slip rate and
drop of shear stress measured at different positions along the fault arise for
the same characteristic slip value $\delta_{\rm loc}$ and coincides with
intense strain localisation. The observed value of $\delta_{\rm loc}$ is
identical to the one reported from one-dimensional simulation under imposed
velocity and is in the order of magnitude of $\delta_{\rm c}$ (see Figures 9
and 10 of Platt et al. (2014)). Remarkably, this observation enables us to
apply the one-dimensional theory discussed in the previous section to derive
predictions of the shear zone dynamics after strain localisation. For
instance, the slip-on-a-plane solution described in Equation (10) can be used
to capture the magnitude of the residual shear stress reached immediately
after strain localisation $\tau_{\rm res}\approx\tau_{\rm
sp}(\delta=\delta_{\rm c};V=V_{\rm tip})$, with $V_{\rm tip}$ being the slip
rate observed at the rupture tip (see Figure 3c and related caption).
Moreover, once the localisation instability arises, the thickness of actively
strained material at various positions along the interface collapses on a
single $W_{\rm loc}(V)$ curve, which follows the prediction given in Equation
(8).

Figure 3: Time evolution of the elastic variables at different locations along
the interface during the dynamic rupture shown in Figure 2. Line colors relate
to the positions along the interface $x/L_{\rm c}$ and the associated
propagation speeds of the rupture $v_{\rm r}/c_{\rm s}$, whereas the arrows
point to the direction of forward time evolution. Slip rate (a) and shear
stress (c) versus slip revealing how the peak slip rate is associated to
abrupt stress drop and arise at the same amount of cumulated slip $\delta_{\rm
loc}\approx 0.3\delta_{\rm c}$. (b) Slip rate versus width of strain rate
localisation $W_{\rm loc}$ measured from the $\dot{\gamma}(y)$ profiles
following the procedure shown in Figure 3 of Platt et al. (2014). The
different post-peak delocalisation trajectories collapse along a single
prediction given in Equation (8). The dashed lines in panel (c) correspond to
the prediction $\tau_{\rm sp}(\delta;V_{\rm tip})$ and gives a good prediction
of the residual shear stress reached after strain localisation. The slip rate
at the rupture tip $V_{\rm tip}$ is approximated by $V$ at the mid-time
between the peaks in shear stress and in slip rate. (A more precise definition
of the tip position is discussed and computed later in the context of Figure
4.)

###  4.1 Rupture dynamics driven by shear localisation

Next, we quantitatively demonstrate how strain localisation is the driving
mechanism of the propagating rupture. To do so, we analyze snapshots of the
propagating rupture and the near-tip evolution of the macroscopic and
microscopic mechanical variables (Figure 4). Ahead of the propagating tip
(point A), the shear zone is creeping with uniform shear strain rate. As the
rupture approaches, the strain rate builds up uniformly across the gouge
(point B) until the localisation instability arises (point C) together with a
rapid increase in macroscopic slip rate $V$ and abrupt drop of shear stress
$\tau$. In the wake of the rupture (point D), the profile of strain rate
across the gouge progressively delocalises, following the decay of the
macroscopic slip rate given by the prediction $W_{\rm rsf}(V)$ shown in Figure
3b. The near-tip evolution of $V$ and $\tau$ is reminiscent to the singular
solutions at the tip of a dynamic fracture (Freund, 1990). Defining
$\alpha^{2}_{\rm s}=1-v^{2}_{\rm r}/c^{2}_{\rm s}$, the analogy to linear
elastic fracture mechanics (LEFM) can be quantitatively tested by rescaling
the slip rate and stress according to

$V\frac{\mu\alpha_{\rm s}}{2v_{\rm r}}=\tau-\tau_{\rm
res}=\Delta\tau=\frac{K}{\sqrt{2\pi(x-x_{\rm tip})}}$ (14)

and fitting the dynamic fracture solution following the procedure of Barras et
al. (2020). The stress intensity factor $K$, residual stress $\tau_{\rm res}$
and position of the tip $x_{\rm tip}$ are the free parameters that are fitted
simultaneously to match the near-tip decrease of $V$ behind the rupture tip
and increase of $\tau$ ahead of the rupture tip.

Figure 4: Snapshot near the tip of the propagating rupture shown in Figure 2.
Bottom panel presents the spacial evolution of the shear stress and slip rate,
which are simultaneously fitted by the fracture mechanics prediction shown by
the dashed blue curve. (See the main text for details on the fitting
procedure). Top panels show the strain rate profile across the shear zone
observed at the instants A,B,C and D corresponding to the black dots in the
bottom panel.

The good agreement with dynamic fracture solution (dashed blue curves in
Figure 4) confirms the crack-like nature of the simulated rupture process near
the tip of the slipping patch. Such agreement allows us to use the inverted
value of $K$ and invoke the crack-tip energy balance to compute the rupture
energy

$G_{\rm c}=\frac{K^{2}}{2\mu\alpha_{\rm s}},$ (15)

which corresponds to the part of dissipated energy that governs the
propagation of the rupture. In seismology, extracting the fracture energy of
natural earthquakes still eludes the resolution of seismic inversions, such
that the breakdown work is often used as a proxy for $G_{\rm c}$ and
integrates the excess of work on top of residual friction (Tinti et al., 2005;
Cocco et al., 2023). For systems where frictional weakening does not
necessarily reach a well-defined residual value, the breakdown work is defined
as (Abercrombie and Rice, 2005):

$E_{\rm
BD}(\delta)=\int_{0}^{\delta}\Big{(}\tau(\delta^{\prime})-\tau(\delta)\Big{)}\>\mathrm{d}\delta^{\prime}.$
(16)

In our numerical simulations, the integration of $E_{\rm BD}$ at different
locations along the interface reveals a clear plateau over an order of
magnitude in slip (Figure 5), which indicates the portion of $E_{\rm BD}$ that
corresponds to $G_{\rm c}$ following Brener and Bouchbinder (2021).
Remarkably, we can then quantitatively verify that the two independent
estimates of the rupture energy (from the near-tip singularity and the
integration of the breakdown work) are in excellent agreement (gray horizontal
line in Figure 5) as another proof of the crack-like nature of the rupture
dynamics. Furthermore, the observed plateau in $E_{\rm BD}$ is clearly
associated to the rapid stress drop caused by localisation instability (see
$\tau(\delta)$ profile in Figure 3c) and confirms that rapid strain
localisation is the driving mechanism of the propagating rupture. In addition,
the magnitude of $G_{\rm c}$ associated to strain localisation is more than
five times smaller than that expected from uniform shearing under adiabatic
undrained conditions ($\sim\tau_{\rm c}\delta_{\rm c}$).

The interplay between strain localisation and rupture dynamics can be further
established by relating the thinnest localisation width observed at a given
location along the interface to the local speed of the rupture (Figure 6).
Following the behavior reported in Figure 3a, the dynamic stress drop caused
by the localisation instability can be well estimated by the slip-on-a-plane
solution (10). Together with the elastodynamic relation (14), it relates the
slip rate near the rupture tip $V_{\rm tip}$ to the rupture speed $v_{\rm r}$,
which can further be combined to the solution $W_{\rm rsf}(V_{\rm tip})$ of
Equation (8):

$\begin{cases}V_{\rm tip}=\displaystyle\frac{2v_{\rm r}}{\mu\alpha_{\rm
s}}\Big{(}\tau_{\rm c}-\tau_{\rm sp}(\delta_{\rm c};V_{\rm tip})\Big{)},\\\
W_{\rm loc}=W_{\rm rsf}(V_{\rm tip}).\end{cases}$ (17)

The implicit relation above provides a universal relationship between $v_{\rm
r}$ and the degree of localisation $W_{\rm loc}$ observed at the scale of the
gouge. Its good agreement with different simulations data (Figure 6) unveil
the dynamical feedback loop at play between strain localisation and rupture
dynamics: the dynamic stress drop caused by strain localisation drives the
acceleration of slip, which further amplifies the dynamic stress drop and
promotes rupture acceleration.

Figure 5: Breakdown energy integrated from the $\tau$ versus $\delta$ profiles
at different positions along the fault, for the simulation shown in Figure 3
(labelled “strong strain localisation”, blue dots) and for another simulation
with the same parameters but larger hydraulic diffusivity (labelled “weak
strain localisation”, brown dots). The rupture simulated with larger hydraulic
diffusivity shows very weak strain localisation. (Additional plots of the
dynamics observed during that simulation are given in Figures 7 and 8 of the
Appendix). The rapid loss of stress caused by strain localisation creates an
horizontal plateau whose associated magnitude is well predicted by the rupture
energy inverted from the dynamic fracture fit shown in Figure 4 and
highlighted here by the horizontal gray line.

##  5 Discussion

Figure 6: Minimum strain rate localisation width $W_{\rm loc}$ versus
instantaneous rupture speed $v_{\rm r}$ computed during rupture propagation
for several simulations using the same parameters but different background
stresses. Heterogeneous simulations with steps in background stress were
conducted with an initial plateau of
$\tau_{\mathrm{b}}/\tau_{\mathrm{c}}=0.59$ around the nucleation zone, and a
sharp drop down to a smaller value at position $x/L_{\mathrm{c}}=\pm 11.52$
away from the center of the nucleation region. In heterogeneous simulations,
rupture speeds may vary nonmonotonically during propagation, initially
increasing and subsequently decreasing when encountering a large downstep in
stress. Regardless of the details of the dynamics, the relationship between
peak localised width and rupture speed is well approximated by the theoretical
prediction proposed in Equation (17).

Our simulations demonstrate that strain localisation produces a loss of
stress-bearing capacity of shear zones that can create and sustain earthquake
rupture. The abrupt drop of shear stress produces an accelerating crack-like
rupture in agreement with the predictions of LEFM. Notably, the rupture is
driven by a well-defined fracture energy that corresponds to the edge-
localised dissipation during the localisation process.

Such a behaviour is in contrast with that of ruptures driven by thermal
pressurisation without internal strain localisation (Viesca and Garagash,
2015), for which breakdown work uniformly increases with increasing slip,
i.e., without a well-defined residual strength at any point within the
propagating rupture. Similarly, the integrated breakdown work for simulated
ruptures that feature weak strain localisation lack a well-defined edge-
localised rupture energy (Figure 5). In this case, the rupture tip singularity
significantly deviates from fracture mechanics predictions (14) and the
rupture dynamics is no longer governed by local near-tip energy balance
(Brener and Bouchbinder, 2021), as further discussed in Appendix D.6. Far from
the rupture tip, our simulations show further shear weakening driven by the
diffusion process (2) that continues at a slower pace as strain delocalises,
and we approach again the slip-on-a-plane asymptotics described by Viesca and
Garagash (2015).

Strain localisation within a preexisting gouge material is strongly correlated
to the dynamics of fault slip, and specifically to the rupture speed (Figure
6). The degree of strain localisation increases with increasing rupture speed,
with a narrowing of the deformed, heated and pressurised region, approaching
1/1000th of the initial shear zone width. Despite the complexity of the
problem, quantitative estimates can be obtained by a simple analytical
approximation (Equation 8) adapted from Platt et al. (2014), so that the
original predictions for peak localised width listed in Rice et al. (2014);
Platt et al. (2014) still apply. Ideally, we could use the relationship
between width and rupture speed depicted in Figure 6 to interpret the
localisation features observed in the geological record in terms of rupture
dynamics. However, strain localisation in rocks is not exclusively associated
with dynamic ruptures and fast slip rates (e.g. Evans and Kohlstedt, 1995),
and only careful micro- and nano-structural studies can be relied upon to
determine the seismic nature of geological structures, notably via detection
of features characteristic of frictional heating (Rowe and Griffith, 2015).
Keeping this caveat in mind, our results highlight that the degree of strain
localisation may be used as a complementary indicator of seismic slip: indeed,
simulations leading to dynamic ruptures are always associated with strong
localisation, with typical width in the sub-millimeter range (see also Daub et
al., 2008).

In this paper, we chose thermal pressurisation as the driving mechanism for
localisation and implemented a numerical scheme coupling small-scale diffusive
processes across the shear zone to long-range elastodynamic coupling along the
shear zone. Our results can however be generalized to other type of
localisation instability arising within shear zones where (1) the shear stress
in the shear zone is function of a set of variables that includes the shear
strain rate $\dot{\gamma}$ and another diffusive quantity $\vartheta$
(Equation 1), and (2) the rate of work produced by shearing acts as a source
term in the diffusion of $\vartheta$ (Equation 2). Importantly, shear
localisation can produce and sustain rupture in shear zones having a rate-
strengthening rheology ($g^{\prime}_{0}>0$) often interpreted as a token of
stability and aseismic slip.

If the conditions (5) and (6) are fulfilled, a localisation instability can
develop and lead to an abrupt drop of shear stress, which leads to the
emergence of a well-defined edge-localised fracture energy and LEFM-like
rupture. Far from the tip, any diffusion-driven weakening leads to $E_{\rm
BD}\sim\delta^{2/3}$ at large slip (Brantut and Viesca, 2017). Therefore, the
behavior summarized in Figure 5 is expected to arise for any type of
localisation-driven rupture, including those where the rheology is controlled
by temperature, such as superplasticity (e.g. Green et al., 2015; De Paola et
al., 2015; Pozzi et al., 2021; Harbord et al., 2021). Indeed, simulations of
high speed deformation in metals, which are also rate-hardening and
temperature-sensitive, tend to exhibit similar characteristics, with the
emergence of a localisation-driven dissipation at the edge of propagating
shear bands (Bonnet-Lebouvier et al., 2002).

Our work demonstrates how localisation instabilities arising across a creeping
shear zone create an abrupt drop of shear stress that promotes the propagation
of classical dynamic ruptures over large distances along the shear zone.
Whether frictional systems are governed by classical fracture mechanics or by
nonlinear friction is an important and debated question in geophysics (e.g.
Svetlizky et al., 2019; Lambert et al., 2021; Paglialunga et al., 2022; Cocco
et al., 2023). Strain localisation is an abrupt structural weakening mechanism
that provides a clear separation between the cohesive zone and the interior of
the slipping patch, hence justifying the small-scale yielding hypothesis.
However, the relative simplicity of the rupture tip behaviour does not
preclude any complexity of the overall rupture style. Away from the rupture
tip, thermal and hydraulic diffusion and strain delocalisation maintain a slow
decay of the shear stress, which is prone to impact how earthquake ruptures
stop (Paglialunga et al., 2022).

## References

  * Abercrombie and Rice (2005) Abercrombie, R. E., and J. R. Rice, Can observations of earthquake scaling constrain slip weakening ?, Geophys. J. Int., 162, 406–424, 2005.
  * Abramowitz and Stegun (1972) Abramowitz, M., and I. A. Stegun, Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables, Dover, New York, 1972.
  * Andrews (2002) Andrews, D. J., A fault constitutive relation accounting for thermal pressurization of pore fluid, J. Geophys. Res., 107(B12), 2363, doi:10.1029/2002JB001942, 2002.
  * Bai (1982) Bai, Y. L., Thermo-plastic instability in simple shear, J. Mech. Phys. Solids, 30(4), 195–207, 1982.
  * Barras et al. (2020) Barras, F., M. Aldam, T. Roch, E. A. Brener, E. Bouchbinder, and J.-F. Molinari, The emergence of crack-like behavior of frictional rupture: Edge singularity and energy balance, Earth Planet. Sci. Lett., 531, 115,978, 2020.
  * Beeler et al. (2008) Beeler, N. M., T. E. Tullis, and D. L. Goldsby, Constitutive relationships and physical basis of fault strength due to flash heating, J. Geophys. Res., 113, B01401, doi:10.1029/2007JB004988, 2008.
  * Bonnet-Lebouvier et al. (2002) Bonnet-Lebouvier, A.-S., A. Molinari, and P. Lipinski, Analysis of the dynamic propagation of adiabatic shear bands, Int. J. Solids Struct., 39, 4249–4269, 2002.
  * Braeck and Podladchikov (2007) Braeck, S., and Y. Y. Podladchikov, Spontaneous thermal runaway as an ultimate failure mechanism of materials, Phys. Res. Lett., 98(9), 095,504, doi:10.1103/PhysRevLett.98.095504, 2007.
  * Brantut and Viesca (2017) Brantut, N., and R. C. Viesca, The fracture energy of ruptures driven by flash heating, Geophys. Res. Lett., 44, 6718–6725, doi:10.1002/2017GL074110, 2017.
  * Brener and Bouchbinder (2021) Brener, E. A., and E. Bouchbinder, Theory of unconventional singularities of frictional shear cracks, J. Mech. Phys. Solids, 153, 104466, 2021.
  * Cocco et al. (2023) Cocco, M., S. Aretusini, C. Cornelio, S. B. Nielsen, E. Spagnuolo, E. Tinti, and G. Di Toro, Fracture energy and breakdown work during earthquakes, Annu. Rev. Earth Planet. Sci., 51, 217–252, 2023.
  * Daub et al. (2008) Daub, E. G., M. L. Manning, and J. M. Carlson, Shear strain localization in elastodynamic rupture simulations, Geophys. Res. Lett., 35, L12310, doi:10.1029/2008GL033835, 2008.
  * De Paola et al. (2015) De Paola, N., R. E. Holdsworth, C. Viti, C. Collettini, and R. Bullock, Can grain size sensitive flow lubricate faults during the initial stages of earthquake propagation?, Earth Planet. Sci. Lett., 431, 48–58, 2015.
  * Evans and Kohlstedt (1995) Evans, B., and D. L. Kohlstedt, Rheology of rocks, in Rock Physics and Phase Relations: A Handbook of Physical Constants, edited by T. J. Ahrens, AGU Reference Shelf, pp. 148–165, American Geophysical Union, Washington, DC, 1995.
  * Freund (1990) Freund, L. B., Dynamic Fracture Mechanics, Cambridge University Press, Cambridge, UK, 1990.
  * Green et al. (2015) Green, H. W., F. Shi, K. Bozhilov, G. Xia, and Z. Reches, Phase transformation and nanometric flow cause extreme weakening during fault slip, Nat. Geosci., 8, 484–489, 2015.
  * Harbord et al. (2021) Harbord, C., N. Brantut, E. Spagnuolo, and G. Di Toro, Fault friction during simulated seismic slip pulses, J. Geophys. Res., 126, e2021JB022149, doi:10.1029/2021JB022149, 2021.
  * Hermundstad et al. (2010) Hermundstad, A. M., E. G. Daub, and J. M. Carlson, Energetics of strain localization in a model of seismic slip, J. Geophys. Res., 115, B06320, doi:10.1029/2009JB006960, 2010.
  * Jop et al. (2006) Jop, P., Y. Forterre, and O. Pouliquen, A constitutive law for dense granular flows, Nature, 441(7094), 727–730, 2006.
  * Lambert et al. (2021) Lambert, V., N. Lapusta, and S. Perry, Propagation of large earthquakes as self-healing pulses or mild cracks, Nature, 591, 252–258, 2021\.
  * Mair and Abe (2008) Mair, K., and S. Abe, 3D numerical simulations of fault gouge evolution during shear: Grain size reduction and strain localization, Earth Planet. Sci. Lett., 274, 72–81, 2008.
  * Mase and Smith (1987) Mase, C. W., and L. Smith, Effects of frictional heating on the thermal, hydrologic, and mechanical response of a fault, J. Geophys. Res., 92(B7), 6249–6272, 1987.
  * Morrissey and Geubelle (1997) Morrissey, J. W., and P. H. Geubelle, A numerical scheme for mode III dynamic fracture problems, Int. J. Num. Meth. Eng., 40, 1181–1196, 1997\.
  * Noda and Lapusta (2010) Noda, H., and N. Lapusta, Three-dimensional earthquake sequence simulations with evolving temperature and pore pressure due to shear heating: Effect of heterogeneous hydraulic diffusivity, J. Geophys. Res., 115, B12314, doi:10.1029/2010JB007780, 2010.
  * Noda and Lapusta (2013) Noda, H., and N. Lapusta, Stable creeping fault segments can become destructive as a result of dynamic weakening, Nature, 493, 518–523, 2013\.
  * Noda et al. (2009) Noda, H., E. M. Dunham, and J. R. Rice, Earthquake ruptures with thermal weakening and the operation of major faults at low overall stress levels, J. Geophys. Res., 114, B07302, doi:10.1029/2008JB006143, 2009\.
  * Paglialunga et al. (2022) Paglialunga, F., F. X. Passelègue, N. Brantut, F. BArras, M. Lebihain, and M. Violay, On the scale dependence in the dynamics of frictional ruptre: Constant fracture energy versus size-dependent breakdown work, Earth Planet. Sci. Lett., 584, 117,442, 2022.
  * Platt et al. (2014) Platt, J. D., J. W. Rudnicki, and J. R. Rice, Stability and localization of rapid shear in fluid-saturated fault gouge, 2. Localized zone width and strength evolution, J. Geophys. Res., 119, 4334–4359, doi:10.1002/2013JB010711, 2014.
  * Platt et al. (2015) Platt, J. D., N. Brantut, and J. R. Rice, Strain localization driven by thermal decomposition during seismic shear, J. Geophys. Res., 120, doi:10.1002/2014JB011493, 2015.
  * Pozzi et al. (2021) Pozzi, G., N. De Paola, S. Nielsen, R. E. Holdsworth, T. Tesei, M. Thieme, and S. Demouchy, Coseismic fault lubrication by viscous deformation, Nat. Geosci., 14, 437–442, 2021.
  * Proctor et al. (2014) Proctor, B. P., T. M. Mitchell, G. Hirth, D. Goldsby, F. Zorzi, J. D. Platt, and G. Di Toro, Dynamic weakening of serpentinite gouge and bare-surfaces at seismic slip rates, J. Geophys. Res., 119, doi:10.1002/2014JB011057, 2014.
  * Rice (2006) Rice, J. R., Heating and weakening of faults during earthquake slip, J. Geophys. Res., 111, B05311, doi:10.1029/2005JB004006, 2006.
  * Rice et al. (2014) Rice, J. R., J. W. Rudnicki, and J. D. Platt, Stability and localization of rapid shear in fluid-saturated fault gouge, 1. Linearized stability analysis, J. Geophys. Res., 119, 4311–4333, doi:10.1002/2013JB01071, 2014.
  * Rowe and Griffith (2015) Rowe, C. D., and W. A. Griffith, Do faults preserve a record of seismic slip: A second opinion, J. Struct. Geol., 78, 1–26, 2015.
  * Scholz (1988) Scholz, C. H., The critical slip distance for seismic faulting, Nature, 336, 761–763, 1988.
  * Scholz (2019) Scholz, C. H., The mechanics of earthquakes and faulting, Cambridge University Press, New York, NY, 2019.
  * Sibson (1975) Sibson, R. H., Generation of pseudotachylyte by ancient seismic faulting, Geophys. J. Int., 43(3), 775–794, doi:10.1111/j.1365-246X.1975.tb06195.x, 1975.
  * Sulem et al. (2011) Sulem, J., I. Stefanou, and E. Veveakis, Stability analysis of undrained adiabatic shearing of a rock layer with Cosserat microstructure, Granular Matter, 13, 261–268, doi:10.1007/s10035-010-0244-1, 2011.
  * Svetlizky et al. (2019) Svetlizky, I., E. Bayart, and J. Fineberg, Brittle fracture theory describes the onset of frictional motion, Annu. Rev. Cond. Mat. Phys., 10, 253–273, 2019.
  * Taylor and Quinney (1934) Taylor, G. I., and H. Quinney, The latent energy remaining in a metal after cold working, Proc. Roy. Soc. Lond. A, 143, 307–326, 1934.
  * Tinti et al. (2005) Tinti, E., P. Spudich, and M. Cocco, Earthquake fracture energy inferred from kinematic rupture models on extended faults, J. Geophys. Res., 110, B12303, doi:10.1029/2005JB003644, 2005.
  * Veveakis et al. (2012) Veveakis, E., J. Sulem, and I. Stefanou, Modeling of fault gouges with Cosserat Continuum Mechanics: Influence of thermal pressurization and chemical decomposition as coseismic weakening mechanisms, J. Struct. Geol., 38, 254–264, 2012.
  * Viesca and Garagash (2015) Viesca, R. C., and D. I. Garagash, Ubiquitous weakening of faults due to thermal pressurization, Nat. Geosci., doi:10.1038/ngeo2554, 2015.
  * Yao et al. (2015) Yao, L., S. Ma, J. D. Platt, A. Niemeijer, and T. Shimamoto, The crucial role of temperature in high-velocity weakening of faults: Experiments on gouge using host blocks of different thermal conductivities, Geology, doi:10.1130/G37310.1, 2015.

##  Appendix A Linear perturbation analysis

Let us consider a shear zone as in Figure 1 initially creeping under imposed
shear stress $\tau_{0}=f_{0}$ and uniform strain rate $\dot{\gamma}_{0}(t)$
and field $\vartheta_{0}(t)$ conditions across the shear zone following the
constitutive Equations (1) and (2):

$\partial_{t}\vartheta_{0}=\beta f_{0}\dot{\gamma}_{0},$ (18)

with $\partial_{t}$ denoting a partial time derivative. At $t=t_{0}$, small
perturbations to the uniform configuration are introduced such that the
evolution of the three variables of interest can be written as

$\left\\{\begin{array}[]{c}\tau(y,t)=f_{0}+\tau_{1}(y,t)\\\
\dot{\gamma}(y,t)=\dot{\gamma}_{0}(t)+\dot{\gamma}_{1}(y,t)\\\
\vartheta(y,t)=\vartheta_{0}(t)+\vartheta_{1}(y,t)\end{array}\right.,$ (19)

with
$\\{\tau_{1}/f_{0},\dot{\gamma}_{1}/\dot{\gamma}_{0},\vartheta_{1}/\vartheta_{0}\\}\ll
1$. Using the definitions of $g^{\prime}_{0}$ and $h^{\prime}_{0}$ given in
the main text (3),(4), the constitutive law governing the shear zone rheology
can be expanded as

$f_{0}+\tau_{1}=f_{0}+g^{\prime}_{0}\dot{\gamma}_{1}+h^{\prime}_{0}\vartheta_{1}+\mathcal{O}(\dot{\gamma}^{2}_{1},\vartheta^{2}_{1}),$
(20)

and further simplified as

$\tau_{1}=g^{\prime}_{0}\dot{\gamma}_{1}+h^{\prime}_{0}\vartheta_{1},$ (21)

neglecting the higher-order terms. The diffusion equation can be similarly
linearised

$\partial_{t}\vartheta_{0}+\partial_{t}\vartheta_{1}=\beta(f_{0}\dot{\gamma}_{0}+f_{0}\dot{\gamma}_{1}+\tau_{1}\dot{\gamma}_{0})+\alpha\partial^{2}_{y}\vartheta_{1}$
(22)

and simplified further using the uniform shear solution (18). The remaining
terms write

$\partial_{t}\hat{\vartheta}_{1}=\beta(f_{0}\dot{\hat{\gamma}}_{1}+\hat{\tau}_{1}\dot{\gamma}_{0})-4\pi^{2}k^{2}\alpha\hat{\vartheta}_{1},$
(23)

using the following spectral decomposition

$\\{\hat{\tau}_{1},\dot{\hat{\gamma}}_{1},\hat{\vartheta}_{1}\\}(k,t)=\int^{\infty}_{-\infty}\\{\tau_{1},\dot{\gamma}_{1},\vartheta_{1}\\}(y,t)e^{-2\pi
iy}\>\mathrm{d}y.$ (24)

From the conservation of momentum and neglecting inertia through the thickness
of the shear zone ($\partial_{y}\tau=0$), one has $\hat{\tau}_{1}=0$ such that
Equation (21) writes

$\dot{\hat{\gamma}}_{1}=-\frac{h^{\prime}_{0}}{g_{0}^{\prime}}\hat{\vartheta}_{1}.$
(25)

Equation (23) becomes

$\partial_{t}\hat{\vartheta}_{1}=-\hat{\vartheta}_{1}\Big{(}\beta
f_{0}\frac{h_{0}^{\prime}}{g_{0}^{\prime}}+4\pi^{2}k^{2}\alpha\Big{)}$ (26)

and allows the following solution

$\hat{\vartheta}_{1}=\hat{\vartheta}_{0}\exp{\Big{(}(t-t_{0})(-4\pi^{2}k^{2}\alpha-\beta
f_{0}\frac{h_{0}^{\prime}}{g_{0}^{\prime}})\Big{)}},$ (27)

with $\hat{\vartheta}_{0}$ being the initial value of $\hat{\vartheta}_{1}$ at
$t=t_{0}$ Consequently, uniform creeping across the shear zone is stable if

$k^{2}>-\frac{\beta
f_{0}}{4\pi^{2}\alpha}\frac{h^{\prime}_{0}}{g^{\prime}_{0}},$ (28)

which leads to the conditions (5) and (6) in the main text.

##  Appendix B Model

###  B.1 Elastodynamics (longitudinal problem)

The longitudinal problem is solved assuming that the shear zone lays between
two semi-infinite half-spaces. The elastodynamics is then solved using a
boundary integral formulation that relates the traction $\sigma_{yi}(x,y=\pm
h/2,t)$ acting within the shear zone to the respective displacements of each
surfaces of the surrounding elastic wall rock $u_{i}(x,y=\pm h/2,t)$. We
consider homogeneous elastic properties and anti-plane shear conditions, such
that the shear traction in the shear zone $\tau(x,t)=\sigma_{yz}$ can be
related to the differential slip
$\delta(x,t)=u_{z}(x,y=h/2,t)-u_{z}(x,y=-h/2,t)$ and slip rate
$V(x,t)=\dot{u}_{z}(x,y=h/2,t)-\dot{u}_{z}(x,y=-h/2,t)$ following Equation
(12), in which $\phi$ accounts for the non-local elastodynamic interactions
and is evaluated in the Fourier domain as:

$\Phi(k,t)=-\frac{1}{2}\mu|k|D(k,t)+\frac{1}{2}\mu|k|\int_{-\infty}^{t}W_{\rm
III}\Big{(}|k|c_{\rm
s}(t-t^{\prime})\Big{)}\dot{D}(k,t^{\prime})\>\mathrm{d}t^{\prime},$ (29)

with $\Phi(k,t)$, $D(k,t)$ and $\dot{D}(k,t)$ being respectively the Fourier
transform pairs of $\phi(x,t)$, $\delta(x,t)$ and $V(x,t)$. The mode III
elastodynamic kernel $W_{\rm III}$ is defined in Morrissey and Geubelle (1997)
as a function of the Bessel function of the first kind $J_{1}$:

$W_{\rm III}(T)=\int_{T}^{\infty}\frac{J_{1}(\zeta)}{\zeta}d\zeta.$ (30)

###  B.2 Thermal pressurisation (transverse problem)

The transverse problem accounts for the thermo-hydro-mechanical shear response
of a fluid-saturated fault gouge. Taking advantage of the time and length
scale separation between the $x$ and $y$ directions existing in the problem of
interest, heat and fluid flows through the granular material are only solved
along the $y$ direction. Next, the assumption of fluid flow dominated by
viscosity (Darcy’s law with hydraulic diffusivity $\alpha_{\rm hy}$) and heat
transfer by conduction (Fourier’s law with thermal diffusivity $\alpha_{\rm
th}$) are invoked to obtained a coupled set of diffusion equations that
describes change of temperature $T$ caused by frictional shear and the
associated change in fluid pressure $p$ due to thermal expansion of the porous
medium. The thermal pressurisation equations read (e.g. Rice, 2006):

$\begin{aligned} \frac{\partial T}{\partial t}&=\frac{\tau\dot{\gamma}}{\rho
c}+\alpha_{\rm th}\frac{\partial^{2}T}{\partial y^{2}},\\\ \frac{\partial
p}{\partial t}&=\Lambda\frac{\partial T}{\partial t}+\alpha_{\rm
hy}\frac{\partial^{2}p}{\partial y^{2}}.\end{aligned}.$ (31)

In the equations above, $\rho$ and $c$ are respectively the density and
specific heat capacity of the gouge and $\Lambda$ is a mean-field parameter
that describes the change of pore fluid pressure caused by an increase in
temperature in the gouge. $\tau$ and $\dot{\gamma}$ are respectively the shear
stress and strain rate in the gouge. Following Rice et al. (2014); Platt et
al. (2014), we assume that the gouge follows a weak rate-hardening behavior:

$\tau=f_{\rm rsf}(\dot{\gamma})(\sigma_{\rm
n}-p)=\Big{(}f_{0}+A\ln(\dot{\gamma}/\dot{\gamma}_{0})\Big{)}(\sigma_{\rm
n}-p),$ (32)

defined by the friction parameters $f_{0}$, $A$ and $\dot{\gamma}_{0}$ and the
normal stress $\sigma_{\rm n}$. By analogy to rate-and-state friction, $A$
controls the ”direct” strengthening response of the gouge whereas the ”long
term” evolution of the shear stress is given by the thermo-hydro-mechanical
response of the gouge through the evolution of $p(y,t)$.

###  B.3 Coupling conditions

Throughout the rupture, the separation of scales between the longitudinal and
transverse direction allows us to neglect 1) fluid and heat flows in the $x$
direction and 2) inertial effects through the thickness of the gouge.
Consequently, the shear stress is invariant across the gouge and corresponds
to

$\tau\Big{(}x,-h/2\leq y\leq h/2,t\Big{)}\equiv\tau(x,t).$ (33)

Conversely, the strain rate varies across the gouge and is related to the
macroscopic slip velocity via the following integration:

$\int_{-h/2}^{h/2}\dot{\gamma}\>\mathrm{d}y=\int_{-h/2}^{h/2}\frac{\partial\dot{u}}{\partial
y}\>\mathrm{d}y=\dot{u}(x,y=h/2,t)-\dot{u}(x,y=-h/2,t)\equiv V(x,t).$ (34)

The stress (33) and kinematic (34) conditions above allow for coupling the
longitudinal and transverse problems.

##  Appendix C Dimensional analysis

Defining $V_{\rm c}=2c_{\rm s}\tau_{\rm b}/\mu$, $t_{\rm c}=\delta_{\rm
c}/V_{\rm c}$ and $L_{\rm c}=c_{\rm s}t_{\rm c}$, the elastodynamic equation
(12) can be rewritten in dimensionless form as

$\tilde{\tau}(\tilde{x},\tilde{t})=1+\tilde{V}(\tilde{x},\tilde{t})+\tilde{\phi}(\tilde{x},\tilde{t}),$
(35)

with $\tilde{\tau}=\tau/\tau_{\rm b}$, $\tilde{x}=x/L_{\rm c}$,
$\tilde{t}=t/t_{\rm c}$, $\tilde{V}=V/V_{\rm c}$ and

$\Phi(\tilde{k},\tilde{t})=-\tilde{D}(\tilde{k},\tilde{t})+|\tilde{k}|\int_{-\infty}^{\tilde{t}}W_{\rm
III}\Big{(}|\tilde{k}|(\tilde{t}-\tilde{t}^{\prime})\Big{)}\dot{\tilde{D}}(\tilde{k},\tilde{t}^{\prime})\>\mathrm{d}\tilde{t}^{\prime}.$
(36)

The characteristic slip distance $\delta_{\rm c}$ will be set later by the
rheology of the interface. Defining $\dot{\gamma}_{\rm c}=V_{\rm c}/h$ and
$f_{\rm c}=f_{\rm rsf}(\dot{\gamma}_{\rm c})$, the effective friction
relationship can be rewritten as

$\displaystyle\tilde{\tau}$ $\displaystyle=\Big{(}1+\frac{A}{f_{\rm
c}}\ln(\dot{\gamma}/\dot{\gamma}_{\rm
c})\Big{)}\Big{(}1-\frac{p-p_{0}}{\sigma_{\rm n}-p_{0}}\Big{)}\frac{f_{\rm
c}(\sigma_{\rm n}-p_{0})}{\tau_{\rm b}}$ (37)
$\displaystyle\equiv\Big{(}1+\tilde{z}^{-1}\ln(\dot{\tilde{\gamma}})\Big{)}\Big{(}1-\tilde{p}\Big{)}\tilde{\eta}^{-1}$

In the equation above, $\tilde{z}+\ln(\dot{\tilde{\gamma}})$ is further
approximated as $\operatorname{arcsinh}(e^{\tilde{z}}\dot{\tilde{\gamma}}/2)$
in order to rewrite the constitutive relationship as

$\tilde{\tau}=\frac{1-\tilde{p}}{\tilde{z}\tilde{\eta}}\operatorname{arcsinh}\Big{(}\frac{\dot{\tilde{\gamma}}}{2}e^{\tilde{z}}\Big{)},$
(38)

which has the advantage of regularizing the stick-to-slip transition as
$\dot{\tilde{\gamma}}\to 0$. In dimensionless form, the kinematic condition
(34) writes then

$\tilde{V}(\tilde{x},\tilde{t})=\int_{-1/2}^{1/2}\dot{\tilde{\gamma}}(\tilde{x},\tilde{y},\tilde{t})\;\>\mathrm{d}\tilde{y}=\int_{-1/2}^{1/2}2e^{-\tilde{z}}\sinh\Big{(}\frac{\tilde{z}\tilde{\eta}\tilde{\tau}(\tilde{x},\tilde{t})}{1-\tilde{p}(\tilde{x},\tilde{y},\tilde{t})}\Big{)}\;\>\mathrm{d}\tilde{y}.$
(39)

The set of diffusion equations governing thermal pressurisation (31) can then
be written in the following dimensionless form

$\displaystyle\frac{\partial\tilde{T}}{\partial\tilde{t}}$
$\displaystyle=\tilde{\eta}\tilde{\tau}\dot{\tilde{\gamma}}+\tilde{\alpha}_{\rm
th}\frac{\partial^{2}\tilde{T}}{\partial\tilde{y}^{2}}$ (40)
$\displaystyle\frac{\partial\tilde{p}}{\partial\tilde{t}}$
$\displaystyle=\frac{\partial\tilde{T}}{\partial\tilde{t}}+\tilde{\alpha}_{\rm
hy}\frac{\partial^{2}\tilde{p}}{\partial\tilde{y}^{2}}$

and provides the definition of the characteristic slip

$\delta_{\rm c}=\frac{\rho ch}{\Lambda f_{\rm c}}.$ (41)

Remarkably, in this dimensionless framework summarized in Table 1, the
behavior of the coupled system is controlled only by the following four
dimensionless quantities:

$\left\\{\begin{aligned} \tilde{\alpha}_{\rm th}&=\alpha_{\rm th}\frac{t_{\rm
c}}{h^{2}}=\alpha_{\rm th}\frac{\mu\rho c}{2c_{\rm s}\tau_{\rm b}f_{\rm
c}\Lambda h},\\\ \tilde{\alpha}_{\rm hy}&=\alpha_{\rm hy}\frac{t_{\rm
c}}{h^{2}}=\alpha_{\rm hy}\frac{\mu\rho c}{2c_{\rm s}\tau_{\rm b}f_{\rm
c}\Lambda h},\\\ \tilde{z}&=\frac{f_{\rm c}}{A},\\\
\tilde{\eta}&=\frac{\tau_{\rm b}}{f_{\rm c}(\sigma_{\rm n}-p_{0})},\\\
\end{aligned}\right.$ (42)

recalling that $f_{\rm c}=f_{\rm rsf}\Big{(}\dot{\gamma}=\frac{2c_{\rm
s}\tau_{\rm b}}{h\mu}\Big{)}$. The reference set of parameters used for the
simulations are representative to a fault zone at 7 km depth (Platt et al.,
2014) and corresponds to $\tilde{\alpha}_{\rm th}=0.0061$,
$\tilde{\alpha}_{\rm hy}=0.075$, $\tilde{z}=15$ and $\tilde{\eta}=0.64$.

Quantities | Variables | Definition
---|---|---
Pressure | $\displaystyle\tilde{p}$ | $\displaystyle(p-p_{0})/(\sigma-p_{0})$
Temperature | $\displaystyle\tilde{T}$ | $\displaystyle\Lambda(T-T_{0})/(\sigma-p_{0})$
Distance across gouge | $\displaystyle\tilde{y}$ | $\displaystyle y/h$
Shear stress | $\displaystyle\tilde{\tau}$ | $\displaystyle\tau/\tau_{\rm b}$
Slip velocity | $\displaystyle\tilde{v}=V/V_{\rm c}$ | $\displaystyle v\times\mu/(2c_{\rm s}\tau_{\rm b})$
Slip | $\displaystyle\tilde{\delta}=\delta/\delta_{\rm c}$ | $\displaystyle\delta\times\Lambda f_{\rm c}/(\rho ch)$
Distance along fault | $\displaystyle\tilde{x}=x/L_{\rm c}$ | $\displaystyle x\times(2\tau_{\rm b})/(\mu\delta_{\rm c})$
Time | $\displaystyle\tilde{t}=t/t_{\rm c}$ | $\displaystyle t\times(2c_{\rm s}\tau_{\rm b})/(\mu\delta_{\rm c})$
Strain rate | $\displaystyle\dot{\tilde{\gamma}}=\dot{\gamma}/\dot{\gamma}_{\rm c}$ | $\displaystyle\dot{\gamma}\times(\mu h)/(2c_{\rm s}\tau_{\rm b})$
Table 1: Table of dimensionless variables

##  Appendix D Numerical methods

###  D.1 System of diffusion equations

The two diffusion Equations (40) solved in the transverse direction can be
generically written as

$\frac{\partial\tilde{\Psi}}{\partial\tilde{t}}=\tilde{S}+\tilde{D}\frac{\partial^{2}\tilde{\Psi}}{\partial\tilde{y}^{2}},$
(43)

and are numerically integrated using Crank-Nicholson method with an explicit
source term:

$\displaystyle\frac{\tilde{\Psi}_{j}^{m+1}-\tilde{\Psi}_{j}^{m}}{\Delta\tilde{t}}=\tilde{S}^{m}_{j}$
$\displaystyle+\frac{\tilde{D}}{\Delta\tilde{y}_{j+1}+\Delta\tilde{y}_{j}}\Big{(}\frac{\tilde{\Psi}_{j+1}^{m+1}-\tilde{\Psi}_{j}^{m+1}}{\Delta\tilde{y}_{j+1}}+\frac{\tilde{\Psi}_{j+1}^{m}-\tilde{\Psi}_{j}^{m}}{\Delta\tilde{y}_{j+1}}-\frac{\tilde{\Psi}_{j}^{m+1}-\tilde{\Psi}_{j-1}^{m+1}}{\Delta\tilde{y}_{j}}-\frac{\tilde{\Psi}_{j}^{m}-\tilde{\Psi}_{j-1}^{m}}{\Delta\tilde{y}_{j}}\Big{)},$

with $\tilde{t}=m\Delta\tilde{t}$ and an irregular grid with
$\Delta\tilde{y}_{j}$ the spacing between the node $j$ and $j-1$. Regrouping
the unknowns at the next time step on the left hand side, the equation above
can be expressed as the following linear system:

$\tilde{\Psi}^{m+1}_{j}=(\mathbf{P}_{lj})^{-1}(\tilde{S}^{m}_{l}+\mathbf{Q}_{ln}\>\tilde{\Psi}^{m}_{n}).$
(44)

The matrix $\mathbf{P}_{lj}$ is tridiagonal and is efficiently inverted using
Thomas algorithm.

###  D.2 Elastodynamics

The elastodynamics along the fault is solved through Equations (12) and (29).
Along a longitudinal grid of equally spaced sampling points $\tilde{x}_{i}$,
the non-local elastodynamic contribution is then integrated as

$\tilde{\Phi}^{m}_{k}=-k\tilde{q}_{0}\tilde{\Omega}^{m}_{k}+k\tilde{q}_{0}\sum_{n=0}^{m}W_{{\rm
III},k}^{m-n}\dot{\tilde{\Omega}}_{k}^{n}\Delta\tilde{t},$ (45)

with $W_{{\rm III},k}^{m}=W_{\rm III}(k\tilde{q}_{0}m\Delta\tilde{t})$ being
the convolution kernel whose integration according to Equation (30) is
computed through a polynomials approximation (Morrissey and Geubelle, 1997;
Abramowitz and Stegun, 1972). $\tilde{\Omega}_{k}$, $\dot{\tilde{\Omega}}_{k}$
and $\tilde{\Phi}_{k}$ are the discrete Fourier mode of respectively
$\tilde{\delta}_{i}$, $\tilde{V}_{i}$ and $\tilde{\phi}_{i}$. The discrete
Fourier transform introduces a longitudinal periodic boundary condition with
period $\tilde{X}$ such that the fundamental wave number
$\tilde{q}_{0}=2\pi/\tilde{X}$. Combining Equations (35) and (39), the
continuity of shear stress through the gouge at time
$\tilde{t}=m\Delta\tilde{t}$ requires that:

$\zeta(\tilde{\tau}^{m}_{i})\equiv\tilde{\tau}^{m}_{i}-1+\sum_{j|\;y_{j}\in[-\frac{1}{2},\frac{1}{2}]}2e^{-\tilde{z}}\sinh\Big{(}\frac{\tilde{\eta}\tilde{z}\tilde{\tau}^{m}_{i}}{1-\tilde{p}^{m}_{i,j}}\Big{)}\;\Delta\tilde{y}-\tilde{\phi}^{m}_{i}=0.$
(46)

$\zeta$ is a monotonic and differentiable function of $\tilde{\tau}^{m}_{i}$
whose root is found iteratively using a Newton-Raphson scheme. Finally, the
slip velocity is integrated from the solution $\tilde{\tau}^{m}_{i}$ as

$\tilde{V}^{m}_{i}=\sum_{j|\;y_{j}\in[\frac{1}{2},\frac{1}{2}]}2e^{-\tilde{z}}\sinh\Big{(}\frac{\tilde{\eta}\tilde{z}\tilde{\tau}^{m}_{i}}{1-\tilde{p}^{m}_{i,j}}\Big{)}\;\Delta\tilde{y}.$
(47)

###  D.3 Coupled scheme

At time $\tilde{t}=m\Delta\tilde{t}$, the coupled numerical scheme is
integrated in time by the following predictor-corrector scheme:

  1. 1.

Compute $\tilde{p}^{m+1}_{j}$ and $\tilde{T}^{m+1}_{j}$ following Crank-
Nicholson integration scheme (Equation 44)

  2. 2.

Integrate predictor slip:
$\tilde{\delta}^{*}_{i}=\tilde{\delta}^{m}_{i}+\Delta\tilde{t}\;\tilde{V}^{m}_{i}$

  3. 3.

Using the predictor slip, compute the predictor velocity $\tilde{V}^{*}_{i}$
by solving Equations (45-46-47)

  4. 4.

Correct slip integration:
$\tilde{\delta}^{m+1}_{i}=\tilde{\delta}^{m}_{i}+\Delta\tilde{t}(\tilde{V}^{m}_{i}+\tilde{V}^{*}_{i})/2$

  5. 5.

Integrate the non-local dynamic contribution $\tilde{\phi}^{m+1}_{i}$ using
Equation (45)

  6. 6.

Compute frictional stress $\tilde{\tau}^{m+1}_{i}$ using Equation (46)

  7. 7.

Compute slip velocity $\tilde{V}^{m+1}_{i}$ using Equation (47)

###  D.4 Initial and boundary conditions

The system of interest consists of a shear zone lying between two semi-
infinite linear elastic wall rocks as shown in Figure 1. The system is
initially at rest under homogeneous pressure $p_{0}$, temperature $T_{0}$,
shear $\tau_{b}$ and normal $\sigma_{\rm n}$ stress. At time $t=0$, the
rupture event is nucleated by introducing a Gaussian perturbation of the pore
pressure at the center of the fault such that

$\tilde{p}(\tilde{x},\tilde{y})=\tilde{\Pi}\exp\Big{(}\frac{\tilde{x}^{2}}{2\tilde{\xi}^{2}}\Big{)},\;\;\;\mathrm{for}\;\tilde{y}\in[-1/2,1/2].$
(48)

The amplitude and standard deviation are typically set respectively to
$\tilde{\Pi}=0.6$ and $\tilde{\xi}=5\%\tilde{X}$. The system is assumed to be
infinite in the transverse directions, whereas periodic boundary conditions in
the longitudinal direction (with period $\tilde{X}$) are assumed to compute
the discrete Fourier modes $\tilde{\Omega}_{k}$, $\dot{\tilde{\Omega}}_{k}$
and $\tilde{\Phi}_{k}$ used in Equation (45).

###  D.5 Stability and convergence

Crank-Nicholson time integration defined in Equation (44) is unconditionally
stable and the definition of the stable time step is set by time integration
of the elastodynamic model following the Courant-Friedrichs-Lewy condition:

$\Delta\tilde{t}=\beta_{\rm CFL}\Delta\tilde{x}.$ (49)

A value of $\beta_{\rm CFL}=0.1$ has been used in the simulations reported in
this paper. The longitudinal direction is regularly discretized with $N_{x}$
sampling points, such that $\Delta\tilde{x}=\tilde{X}/N_{x}$. In the
transverse direction, the gouge thickness is regularly sampled by $N_{y}+1$
points, such that $\Delta\tilde{y}=1/N_{y}$. To capture the exponential decay
of pressure and temperature fields outside the gouge layer, the numerical
domain symmetrically stretches into the surrounding wall rock, which is
sampled with $N_{y}$ additional points that are logarithmically-spaced such
that $\Delta\tilde{y}_{j=1}=\Delta\tilde{y}_{j=2N_{y}}=1$. Background pore
pressure and temperature $\tilde{p}=\tilde{T}=0$ are imposed at the two edges
($y_{j=0}$ and $y_{j=2N_{y}}$) of the numerical grid, whose distance from the
center of the gouge should exceed the diffusion length scale:

$\tilde{y}_{j=2N_{y}}>2\sqrt{\tilde{\alpha}N_{t}\Delta\tilde{t}},$ (50)

with $\tilde{\alpha}=\sqrt{\tilde{\alpha}^{2}_{\rm hy}+\tilde{\alpha}^{2}_{\rm
th}}$ and $N_{t}$ being the total number of time steps.
Typical simulations use $N_{x}=2048$ and $N_{y}=200$, which are sufficient to
reach numerical convergence.

###  D.6 Rupture dynamics with weak strain localisation

The brown data set shown in Figures 5, 7 and 8 correspond to a shear zone with
(4$\times$) larger diffusivities, which lead to weaker strain localisation and
larger localisation width following the predictions of Equations (6) and (8).
As a result, the shear stress does not feature a sharp drop but rather a
progressive weakening driven by the increase in pressure and temperature with
no well-defined residual shear stress. Dynamic rupture can still be produce by
using (4$\times$) larger domain and nucleation size but significantly differ
from the dynamics of shear crack. In absence of a well-defined residual shear
stress, the near-tip singular solution becomes (e.g. Brener and Bouchbinder,
2021):

$V\sim(x-x_{\rm{tip}})^{\xi-1}\sin(\xi\pi),$ (51)

with the singularity $\xi$ being one half at the vicinity of a standard LEFM
shear crack (as in Equation 14). Strong shear localisation produces a well-
defined residual stress and a near tip singularity similar to standard shear
crack, as demonstrated in the main text (Figure 4). This is no longer true for
failure with weak strain localisation that features significantly different
near-tip singularity and rupture dynamics (Figure 8).

Figure 7: Same as Figure 3 but for a simulation with four-times larger thermal
and hydraulic diffusivities. Consequently, strain localisation remains very
weak during the initial stage of the rupture when $v_{\rm r}<0.9c_{\rm s}$.
Figure 8: Snapshot of the slip velocity concentration near the rupture tip
moving at $v_{\rm r}=0.85c_{\rm s}$ for the simulation with large thermal and
hydraulic diffusivities and weak strain localisation shown in Figure 7. The
near-tip singularity significantly deviates from LEFM prediction of Equation
(14) but is well described by the asymptotic solution (51) with an exponent
$\xi=0.63$.