# The role of gouge production in the seismic behavior of rough faults: A
numerical study

###### Abstract

Fault zones mature through the accumulation of earthquakes and the wearing of
contact asperities at multiple scales. This study examines how wear-induced
gouge production affects the evolution of fault seismicity, focusing on
earthquake nucleation, recurrence, and moment partitioning. Using 2D quasi-
dynamic simulations integrating rate-and-state friction with Archard’s law of
wear, we model the space-time distribution of gouge and its effect on the
critical slip distance. The study reveals a shift from single to multi-rupture
nucleation, marked by increased foreshock activity. The recurrence interval
undergoes two separate phases: an initial phase of steady increase followed by
a secondary phase of unpredictable behavior. Finally, we observe a transition
in the moment partitioning from faster to slower slip rates and a decrease in
the moment released per cycle relative to the case where no gouge formation is
simulated. This research sheds light on wear-driven mechanisms affecting fault
slip behavior, offering valuable insights into how the evolution of gouge
along a fault affects its seismic potential.

Institute for Building Materials, ETH Zürich, Switzerland Department of Earth,
Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology,
USA

=name==email address=

We consider the effect of gouge production in earthquake cycle simulations on
rough faults.

Gouge production causes heterogeneity in fracture energy, favoring partial
ruptures and more variable recurrence intervals.

The fault releases less seismic energy overall, favoring aseismic slip and
creep-dominated nucleation.

## Plain Language Summary

Our research explores the evolution of faults – cracks in the Earth’s crust
where earthquakes originate – as they undergo wear from friction. Through
simplified computer models, we examine how gouge (the small rock particles
created by wear) influences the initiation, frequency, and magnitude of
earthquakes. Our findings indicate that as more gouge accumulates on the fault
over time, the initiation process of earthquakes shifts, leading to an
increase in precursory activity before a major earthquake. Additionally, the
time between earthquakes changes, initially becoming longer before turning
erratic. We also observed a gradual decrease in the total energy released by
faults as they wear down, that is, faults with more gouge host earthquakes of
relatively smaller magnitudes.

##  1 Introduction

Rough faults are known to produce large volumes of gouge during frictional
slip [Power . (1988), Wang Scholz (1994), McBeck . (2021), Mollon . (2023)].
This is due to their steep topography with asperities that interlock and
produce gouge through wear [Power . (1988), Aghababaei . (2016), Frérot .
(2018), Milanese . (2020)]. Sheared by tectonic stresses in the Earth’s crust,
natural faults often exhibit far-from-equilibrium behavior, withstanding large
elastic energy gradients that they dissipate through the nucleation of
frictional instabilities (i.e. earthquakes). As these instabilities unfold,
some of the stored elastic energy is emitted in the form of seismic waves
while most of it is dissipated through a set of different phenomena, such as
frictional heating, crack propagation, inelastic deformation, off-fault damage
[Nielsen . (2016), Reches . (2019), Cocco . (2023)] and other non-linear
dissipation processes. These unstable fronts originate when the interlocking
contact asperities are loaded beyond their stability threshold, triggering a
positive feedback mechanism known as dynamic weakening. In this process, the
more they slip, the weaker they become, leading to even larger slip in a self-
reinforcing cycle [Sammis Sornette (2002)] that relaxes the fault by re-
distributing the frictional stresses all along it.

The re-distribution of stress on the fault following instability is a well-
known mechanism that has been modeled extensively using numerical methods
[Lapusta Rice (2003), Aagaard Heaton (2008), Radiguet . (2013), Bruhat .
(2020), Cattania Segall (2021), Albertini . (2021), Schär . (2021)]. However,
gouge production by wear is also known to play an important role in the self-
organization and maturation of natural faults [Duarte . (2009), Milanese .
(2019), Ostapchuk . (2021), Mollon . (2023)], yet most models of fault
dynamics ignore its contribution. Both observational and experimental studies
have pointed out the stabilizing effect of gouge on fault slip [Marone .
(1990), Marone (1998), B. Carpenter . (2012), BM. Carpenter . (2014)],
enabling velocity strengthening mechanisms that favor aseismic creep over
dynamic weakening. Still, given the scale and complexity of the problem, most
studies cannot afford to capture the long-term evolution of gouge accumulation
and its effect on the overall seismic cycle. Although several numerical works
have studied the effect of gouge on the frictional behavior of simulated
faults [Mora Place (1998), Mair Abe (2008), Abe Mair (2009), Ferdowsi Rubin
(2021), Mollon . (2023)], few of them incorporate a mechanism for dynamic
production of gouge coupled to frictional slip. Gouge is likely to play a key
role in rupture dynamics by modulating the slip weakening distance
$D_{\mathrm{c}}$ [Marone Kilgore (1993)], which in turn affects the fracture
energy and can determine rupture propagation and arrest.

Bizzarri (2010) simulated dynamic gouge production in a spring-slider system,
modelling its effect on $D_{\mathrm{c}}$, and found the recurrence interval to
decrease with the amount of gouge, but they neglected the effect of spatial
heterogeneity. Moreover, Tal (2023) used geometric considerations to account
for gouge creation through the removal of material at the interface and found
that the magnitude of events would remain higher with wear than without it,
but they did not model its effect on the critical slip distance. The effect of
gouge on the critical slip distance is important because it shapes the
fracture energy profile along the fault, with notorious consequences on fault-
slip dynamics. Yet, a thorough analysis of the complex feedback mechanics that
underlie such coupling between friction (through $D_{\mathrm{c}}$) and wear-
induced heterogeneous gouge production is missing.

Here, we use a well-established framework for simulating cycles of earthquakes
Segall Bradley (2012); Cattania Segall (2021) to explore fault dynamics
induced by frictional wear on rough faults and shed light on the role that
frictional wear plays in the seismicity of a fault over cycles of earthquakes.
In particular, we analyze the evolution of seismicity patterns on rough faults
as a result of wear processes that account for the generation and accumulation
of gouge material. With this approach, we investigate the effect of
progressive gouge accumulation on the recurrence interval, the nucleation
dynamics, and partition of the moment between seismic and aseismic components.

##  2 Numerical model

###  2.1 Governing and constitutive equations

We employ a quasi-dynamic boundary element code Segall Bradley (2012) to
conduct 2-D plane-strain simulations of a rough fault (see Fig.1a). Fault slip
is governed by the following equation:

$\tau(x)-\tau_{\mathrm{f}}(x)=\frac{\mu}{2c_{\mathrm{s}}}V(x)$ (1)

where $x$ is the distance along the fault, $\tau_{\mathrm{f}}$ is the
frictional resistance, and $\tau$ is the shear stress due to remote loading
and elastic transfer of stress between fault elements. The right-hand side is
the radiation damping term, which represents the stress change due to
radiation of plane S-waves, with $c_{\mathrm{s}}=3700~{}\mathrm{m/s}$ the
shear wave speed, $\mu=30~{}\mathrm{GPa}$ the shear modulus, and $V(x)$ the
slip rate.

Frictional resistance evolves according to rate-and-state friction Dieterich
(1979); Ruina (1983) in its regularized form:

$\tau_{\mathrm{f}}(V,\theta)=a~{}\sigma_{\mathrm{n}}~{}\text{asinh}\Bigl{[}\frac{V}{2V_{\mathrm{0}}}\text{exp}\Bigl{(}\frac{\Psi}{a}\Bigr{)}\Bigr{]}$
(2)

where $\Psi=f_{0}+b\ln(V_{\mathrm{0}}\theta/D_{\mathrm{c}})$,
$\sigma_{\mathrm{n}}$ is the normal contact pressure, $f_{0}=0.6$ and
$\mathrm{V}_{0}=10^{-6}~{}\mathrm{m/s}$ are reference values for the friction
coefficient and the slip rate, $D_{\mathrm{c}}$ is the critical slip distance,
$a=0.015$ and $b=0.020$ two constitutive parameters and $\theta$ is the state
variable, which represents the local contact area of asperities, expressed in
time-units as a proxy for the asperity ‘life-span’. For the state evolution,
we employ the aging law Dieterich (1979); Ruina (1983):

$\dot{\theta}=1-\frac{V\theta}{D_{\mathrm{c}}}.$ (3)

This constitutive behavior causes stress to increase from the background
stress $\tau_{0}$ with fault load reaching its peak strength at
$\Delta\tau_{\mathrm{0p}}$ and then unload down to
$\tau_{0}-\Delta\tau_{\mathrm{0r}}$ over a slip distance proportional to
$D_{\mathrm{c}}$ (see constitutive behavior of a fault-point in Fig. 1b).

Figure 1: Simulation setup and sample data. a) Simulation setup. We impose a
pure-shear fixed loading rate $\dot{\tau}$ on a rough fault of length $L$ with
a rate ($V$) and state ($\theta$) constitutive law. b) Frictional stress
$\tau-\tau_{0}$ at a given point on the fault (marked in c with a white dashed
line) as a function of the local slip $\delta-\delta_{0}$ for a given cycle,
where $\tau_{0}$ and $\delta_{0}$ are the shear and slip at the beginning of
the cycle ($t=t_{0}$). $\Delta\tau_{\mathrm{0p}}$ stands for initial-to-peak
stress, $\Delta\tau_{\mathrm{pr}}$ peak-to-residual, and
$\Delta\tau_{\mathrm{0r}}$ initial-to-residual stress changes. c) Space-time
map of the slip rate during a given cycle. d) Example of one of the roughness
profiles along with the associated normal stress distribution
$\sigma_{\mathrm{n}}(x)$ and its evolution throughout the cycle (red shows
increase, blue shows decrease). The horizontal dashed line marks the average
background normal stress on the fault $\sigma_{\mathrm{n}}^{0}=1\mathrm{MPa}$.
e) Slip profile at the end of the cycle (blue) and evolution of
$D_{\mathrm{c}}(x)$ during the cycle as a result of the wear coupling (red).
The shaded area corresponds to the difference between the start (brighter red)
and the end (darker red) of the cycle. All panels refer to the same cycle.
Vertical dotted lines in panels c-e indicate stress peaks.

###  2.2 Roughness

To model fault roughness, we use the method from Cattania Segall (2021) by
adopting a fractal surface geometry (see Fig. 1d), characterized by power
spectral density

$P_{\mathrm{h}}=C_{\mathrm{h}}|k|^{-\beta}$ (4)

where we assume the shaping of the fault’s geometry to be well described by a
non-stationary process with stationary increments Milanese . (2019). In such
case, $\beta=2H+1$, where $H$ is the Hurst exponent. For natural faults, this
is typically between 0.6 and 0.8 Renard . (2013); Beeler (2023); here, we set
$H=0.7$. To keep the simulations feasible in terms of computation time, we
only include wave numbers up to $k_{\mathrm{max}}=2\pi(N_{\mathrm{max}}-1)/L$,
where $L=15.26$ km is the size of the fault (see Fig. 1a) and
$N_{\mathrm{max}}$ is set to $L/L_{\mathrm{n}}=15$, $L_{\mathrm{n}}=1.02$ km
being the nucleation length, calculated from Rubin Ampuero (2005) as

$L_{\mathrm{n}}=\frac{1}{\pi}\frac{\mu^{\prime}D_{\mathrm{c}}}{b\sigma_{\mathrm{n}}}\Bigl{(}\frac{b}{b-a}\Bigr{)}^{2}~{},$
(5)

where $\mu^{\prime}=\mu/(1-\nu)$ and $\mathrm{\nu}=0.25$ is the Poisson’s
ratio. Moreover, we define the minimum wavelength of the heterogeneous stress
profile on the fault by
$\lambda_{\mathrm{min}}=2\pi/k_{\mathrm{max}}=1/(N_{\mathrm{max}}-1)L$, where
$N_{\mathrm{max}}$ is always fixed at 15. Therefore,
$\lambda_{\mathrm{min}}=1.09\;\mathrm{km}$, which is very close to the
nucleation length $L_{\mathrm{n}}=1.02\;\mathrm{km}$ (Eq. 5). This allows us
to capture more explicitly the interaction between nucleating fronts and
strength peaks (asperities) along the fault.

###  2.3 Loading

In order to simplify the loading process, we restrict the stress-rate tensor
to a pure shear component, which is the one responsible for slip initiation.
The resulting shear loading rate is given by

$\dot{\tau}=\frac{\dot{\sigma}_{\mathrm{D}}}{2}\sin(2\alpha_{1}+2\alpha_{2}),$
(6)

where $\mathrm{\dot{\sigma}_{D}}=0.25~{}\mathrm{MPa/year}$ is the differential
stress rate, $\alpha_{1}$ is the loading angle, and $\alpha_{2}$ the slope of
roughness. By setting the loading angle to $45^{\circ}$, we can fix the
spatially averaged effective normal stress $\sigma_{\mathrm{n}}^{0}$, which we
set to a uniform value of $1~{}\mathrm{MPa}$. The choice of a low value of
$\sigma_{\mathrm{n}}^{0}$ is dictated by computational costs. Slip on a rough
fault causes normal stress perturbations, and the time step is controlled by
the smallest normal stress in the domain. Since negative stress changes could
potentially induce fault opening if we assume a purely elastic response, we
introduce a minimum value of $\sigma_{\mathrm{min}}$ below which the normal
stress cannot fall, where
$\sigma_{\mathrm{min}}=1~{}\mathrm{kPa}\ll\sigma_{\mathrm{n}}^{0}$. This
threshold ultimately determines the simulation time step, introducing a
limitation to the maximum amount of time that we can simulate in a reasonable
run-time. In our case, the simulations model a period of around 30 years. By
lowering the value of $\sigma_{\mathrm{n}}^{0}$ to $1~{}\mathrm{MPa}$ instead
of $10~{}\mathrm{MPa}$, which would be more realistic, we are able to simulate
around 35 to 40 earthquake cycles, while for
$\sigma_{\mathrm{n}}^{0}=10~{}\mathrm{MPa}$, we only get around 4 in those
same 30 years (see Supplementary Material, case HNS).

###  2.4 Gouge production

To model gouge production on the fault, we employ Archard’s wear law Archard
(1953). First developed in the context of adhesive wear for metal-on-metal
contact, this linear model also captures well the formation of gouge in rock
interfaces Wang Scholz (1994), either by adhesive, abrasive or comminution
mechanisms Hughes . (2020). Archard’s model describes the total volume of wear
produced as a linear function of the normal stress and the total amount of
slip. To account for spatial heterogeneity, we apply it to every point in
space ($x$) of our system so that the local wear rate in units of gouge
thickness per unit time ($t$) is given by

$S_{g}(x,t)=K\frac{\sigma_{\mathrm{n}}(x)V(x,t)}{\eta},$ (7)

where $K=10$ is a dimensionless wear coefficient of the same order of
magnitude as the one in Boneh Reches (2018), $V(x,t)$ is the slip rate, and
$\eta=1.49~{}\mathrm{GPa}$ is the hardness, which is representative of soft
materials like gypsum or calcite Broz . (2006) that produce larger volumes of
gouge as compared to harder materials like quartz. This allows us to shorten
the time of the simulations, producing more gouge per sliding distance. In
particular, we obtain an average gouge production rate of $6.7~{}\mathrm{mm}$
of gouge per meter of accumulated frictional slip, which can double for
regions under higher normal stress, where $\sigma_{\mathrm{n}}>1\mathrm{MPa}$
(see Fig. 1d).

To quantify the effect of wear on the frictional response of the fault, we
coupled it linearly to our rate-and-state friction model through the
$D_{\mathrm{c}}$ parameter, in agreement with experimental observations Marone
Kilgore (1993). As a result, $D_{\mathrm{c}}$ evolves in space and time
according to

$D_{\mathrm{c}}(x,t)=D_{\mathrm{c}}^{0}+\mathrm{\gamma}g(x,t),$ (8)

where $D_{\mathrm{c}}^{0}=10^{-4}~{}\mathrm{m}$ is the initial uniform value
for $D_{\mathrm{c}}$, obtained from Scholz (1988) for deep faults,
$\gamma=0.01$ is a dimensionless coupling constant whose value is derived from
the experimental work by Marone Kilgore (1993), and $g(x,t)$ is the gouge
field representing the thickness of the gouge layer. The latter evolves
according to the classical transport or continuity equation

$\dot{g}(x,t)=-\nabla(V(x,t)g(x,t))+S_{g}(x,t)$ (9)

where the gradient term represents the outflow of gouge as a result of
advection by the embedding velocity field, and $S_{g}$ is the source term,
given in Eq. 7. A representative evolution of $D_{\mathrm{c}}(x,t)$ during a
cycle is characterized by a pronounced increase over regions with high normal
pressure (see Eq. 7) and accumulated slip (see elliptical slip profile in
Fig.1e). The advection of gouge is only responsible for a small percentage of
the layer’s thickness, but it does help to smooth out the distribution of
$D_{\mathrm{c}}$, especially for higher normal pressures (See Fig. S8 in the
Supplementary Material).

##  3 Results and discussion

To isolate the effect of frictional wear and separate it from that of
roughness, we ran two different sets of simulations, the reference simulations
without wear contributions to the fault friction (i.e. $\gamma=0$) and
simulations implementing our proposed model (i.e. $\gamma\neq 0$). Each set
consists of 10 simulations, that differ only in the randomness of the fault’s
fractal geometry, introduced through random phase shifts in the frequency
components of the power spectrum. That is, simulations within a given set all
have different roughness profiles, but these profiles are the same across
different sets to be able to compare them coherently and isolate the effect of
wear.

Two extra sets of 10 simulations were run, where $\gamma=0$ and
$D_{\mathrm{c}}=3\cdot\mathrm{D_{c}^{0}}$, for
$\lambda_{\mathrm{min}}/L_{\mathrm{n}}=1$ and
$\lambda_{\mathrm{min}}/L_{\mathrm{n}}=1/3$. This is done to study the
differences between a situation where $D_{\mathrm{c}}$ is higher but uniform
and another situation where $D_{\mathrm{c}}$ is increasingly heterogeneous due
to wear. Furthermore, another 2 additional sets were run, with
$\sigma_{\mathrm{n}}^{0}=10\mathrm{MPa}$, to study the effect of normal stress
on gouge formation (see Supplementary Material). In total, we ran 60
simulations, corresponding to six different cases (see Table 1).

The resolution is maintained at $\Lambda_{0}/\Delta x>2$ over the entire
domain at all times, in agreement with Cattania Segall (2021), where
$\Lambda_{0}\sim\mu D_{\mathrm{c}}/b\sigma_{\mathrm{n}}$ is the size of the
quasi-static process zone at zero rupture speed Andrews (1976); Rubin Ampuero
(2005); Perfettini Ampuero (2008). Finally, we define seismic events as slip
events where the slip velocity anywhere on the fault exceeds the threshold
velocity $V_{\mathrm{dyn}}$, here $\sim 0.4~{}\mathrm{cm/s}$, which depends on
the normal stress on the fault as
$V_{\mathrm{dyn}}=2a\sigma_{\mathrm{n}}^{0}c_{\mathrm{s}}/\mu$ Rubin Ampuero
(2005).

Model |  $\lambda_{\mathrm{min}}$ |  $\eta$ |  $\sigma_{\mathrm{n}}^{0}$ |  $L$ |  $L_{\mathrm{n}}$ |  $D_{\mathrm{c}}^{0}$
---|---|---|---|---|---|---
NG |  1.09 km |  0 |  1 MPa |  15.26 km |  1.02 km |  $10^{-4}$ m/s
G |  1.09 km |  0.01 |  1 MPa |  15.26 km |  1.02 km |  $10^{-4}$ m/s
NG-X1 |  3.27 km |  0 |  1 MPa |  45.79 km |  3.06 km |  $3\cdot 10^{-4}$ m/s
NG-X2 |  1.09 km |  0 |  1 MPa |  15.26 km |  3.06 km |  $3\cdot 10^{-4}$ m/s
NG-HNS |  109 m |  0 |  10 MPa |  1.52 km |  102 m |  $10^{-4}$ m/s
G-HNS |  109 m |  0.01 |  10 MPa |  1.52 km |  102 m |  $10^{-4}$ m/s
Table 1: Different sets of simulations with corresponding parameters.

###  3.1 Nucleation dynamics: Single-rupture vs. Multiple-rupture

We analyze the nucleation dynamics by considering the space-time distribution
of the slip rate ($V(x,t)$) throughout different cycles, as illustrated in
Figs. 2a-d for an intermediate (Cycle 20) and a late cycle (Cycle 30) in both
NG and G cases.

One of the main results of this analysis is that the presence of gouge (case
G) causes a progressive increase of slow slip events preceding the mainshock,
together with a notable intensification of foreshocks; in contrast to the case
without gouge (NG), where the earthquake initiation process features a system-
spanning rupture. This does not imply that the NG case cannot exhibit a
similar nucleation behavior to G. Indeed, Cattania Segall (2021) shows how,
for sufficiently high normal stress perturbations, a fault can produce slow
slip and foreshocks. However, through the positive feedback mechanism that
links the stress and accumulated slip to the critical slip distance $D_{c}(x)$
by the formation of gouge, this process is enhanced. The situation is reversed
in Tal (2023), where they observe a reduction of aseismic slip when wear is
accounted for. We believe this is due to the fact that they do not model the
link between gouge formation and frictional traction at the interface,
effectively ignoring the effect of wear on the fracture energy.

As the amplitude of both $D_{\mathrm{c}}(x)$ and $\sigma_{\mathrm{n}}(x)$
grows from cycle to cycle, enhanced by the accumulation of slip (see slip
profile in Fig. 1e), G goes from displaying only one small partial rupture
followed by a full rupture in cycle 20 (Fig. 2b) to a more complex process in
cycle 30 (Fig. 2d), involving two partial ruptures that break 60% of the fault
before triggering the full rupture. The magnitude of these partial ruptures
can be deduced from the slip rate (red line in Figs. 2i-iv) and the moment
rate per unit length on the fault (blue line in Figs. 2i-iv), which is
computed as

$\dot{M}_{\mathrm{l}}=\mu\int V(x,t)dx\quad.$ (10)

Prior to the mainshock, all simulations show an aseismic peak in the slip rate
(Figs. 2i-iii), while G-Cycle 30 shows two seismic peaks (Fig. 2iv). As for
the full rupture, the moment rate is much higher in cycle 20 than it is in
cycle 30, indicating a decrease in magnitude (see blue lines in Figs. 2ii and
iv). Both of these ruptures, partial and full, are prompted by the loading
effect of expanding aseismic fronts that precede them. These fronts (indicated
by white arrows in Figs. 2b and d) seem to originate from small local events,
which nucleate at the edge of weakening patches growing in between asperities.
Very short-lived, these small events are self-arrested events, unable to
penetrate the strength barriers ahead of them. Instead, they evolve into the
aforementioned, stable and slowly growing aseismic fronts, before becoming
unstable again as soon as they overcome such constraining barriers.

Figure 2: Earthquake nucleation process for representative cycles in
simulations without (NG) and with (G) gouge production. Data for both cases is
shown for an intermediate cycle (cycle 20) and a late cycle (cycle 30) of a
given set of simulations. i-iv) Evolution of the maximum fault slip rate (red)
and the moment rate per unit length (dark blue), where earthquakes are
signaled with dotted lines as the peaks in the maximum slip rate going above
the dynamic threshold and separated at least by one minute. a-d) Color maps of
the slip rate distribution over the entire fault from the beginning to the end
of the cycle. White ’x’ signs mark the location of every asperity failure on
the fault, which can trigger either a partial or a full rupture. e,f)
Approximations of $G_{\mathrm{s}}(x)$ and $G_{\mathrm{c}}(x)$ (Eqs. 11 and 12)
for both cases, along with the slip distribution $\delta(x)$ (dotted green
lines). $G_{\mathrm{c}}(x)$ is plotted for a range of observed co-seismic slip
rates going from $10^{-3}$ m/s to $10^{-4}$ m/s. Data from cycle 20 is shown
in dark color, while that from cycle 30 is displayed in lighter color.

The overall process of slip nucleation and propagation is described here in
detail. As the fault is uniformly loaded, slip starts nucleating at the
weakest parts, where the peak frictional strength ($\Delta\tau_{\mathrm{0p}}$)
is lower. When they reach their peak strength, these weaker spots abruptly
change their behavior, from strengthening to weakening as they start to slip,
transferring all the excess stress outwards to their surroundings before
reaching a plateau of strength (see Fig. 1b) as they approach the residual
level. The stress accumulating on the borders of the slipping area as a result
of the redistribution of the external loading pushes the neighboring parts to
start weakening too, creating a slip-weakening front around the patch that
drives its own growth, while the central parts of the patch stay close to the
residual strength level. After they reach a certain critical length (analogous
to Eq. 5), which depends on the local heterogeneities, these patches become
unstable and start growing spontaneously, potentially propagating across the
whole interface. This is what happens in NG (see Fig. 2a, c), where the
instability generated is strong enough to break through all asperities,
triggering a propagating front that ruptures the entire fault. It is worth
noting here that, as the instability unfolds, the parts that first started
weakening immediately reach their weakening limit and go into approximately
constant frictional strength. This becomes particularly important for the
study of partial ruptures in the fault exhibited by case G (right column of
Fig. 2), which shows different nucleation dynamics that favor aseismic slip. A
recent study Castellano . (2023) reported simulations in which frictional
instabilities that arrested after hitting a barrier, led to a continuous
stable expansion of the containing slipping patch with slow slip inside that
loads the neighboring strong asperities until they break. This behavior has
been reported in other numerical studies Rubin Ampuero (2005); Garagash
Germanovich (2012); Dublanchet (2018); Cattania Segall (2021); Sáez Lecampion
(2024), and it is what we observe too (see aseismic fronts indicated by white
arrows on Figs. 2b,d), suggesting it could be a leading mechanism for slow
slip events and accompanying tremor (resulting from the self-arrested
ruptures).

To understand the nucleation mechanics at play, we can refer to the energy
balance of nucleating patches, illustrated in Figs. 2e-f. According to
classical nucleation theory, a phase grows spontaneously whenever the change
in free energy ($\Delta\mathcal{G}$) relative to its change in size becomes
negative. Since $\Delta\mathcal{G}=G_{\mathrm{c}}-G_{\mathrm{s}}$, where
$G_{\mathrm{c}}$ is the fracture energy and $G_{\mathrm{s}}$ is the static
energy release rate, the patch will spontaneously grow whenever
$G_{\mathrm{s}}>G_{\mathrm{c}}$. These two terms are defined Lawn (1993);
Rubin Ampuero (2005) as

$G_{\mathrm{c}}(x)=\frac{\Delta\tau_{\mathrm{pr}}^{2}}{2W^{*}}=\frac{1}{2}b\sigma_{\mathrm{n}}(x)D_{c}(x)\Big{[}\ln\Big{(}\frac{V_{\mathrm{co}}\theta_{\mathrm{i}}}{D_{\mathrm{c}}\Omega}\Big{)}\Big{]}^{2}$
(11)
$G_{\mathrm{s}}(x)=\frac{K^{2}}{\mu^{\prime}}=\frac{\pi\ell}{2\mu^{\prime}}\Delta\tau_{\mathrm{0r}}^{2}=\frac{\pi\ell}{2\mu^{\prime}}\Big{[}\dot{\tau}T\Big{]}^{2},$
(12)

where $W^{*}$ is the slip weakening rate, $K$ is the stress intensity factor,
$\ell$ is the crack length, $T$ is the recurrence interval, $V_{\mathrm{co}}$
is the co-seismic slip rate, $\Omega=V\theta/D_{\mathrm{c}}\approx 1$ for
steady-state and $\theta_{\mathrm{i}}$ is the state variable prior to the
arrival of the crack-tip, which can be approximated by $T$ assuming
$\dot{\theta}\approx 1$ interseismically. Note that the stress drop in eq. 12
is not the coseismic stress drop during a single rupture, but the total shear
stress imposed by loading since the previous full rupture. We assume that this
is entirely released by the combination of the ongoing rupture and all
previous earthquakes and creep within it, all of which contribute to the
stress concentration outside the crack tip quantified by $K$.

As gouge builds up around the strongest points (high $\sigma_{n}$) surrounding
nucleating patches, modulated by the amount of accumulated slip, and makes
$D_{\mathrm{c}}$ larger (see Eq.8), the fracture energy $G_{\mathrm{c}}$
increases accordingly (see Fig. 2f), making it harder for these fronts to
propagate further. This increase in fracture energy is proportional to the
total slip accumulated on the fault (dotted lines in Figs. 2e-f), in agreement
with Archard’s law (Eq. 7). As a result, the energy release rate of nucleating
patches is not enough to overcome these barriers, as observed in Figs. 2e-f,
where arrest events tend to coincides with points on the fault where
$G_{\mathrm{c}}>G_{\mathrm{s}}$, within uncertainties derived from the choice
of co-seismic slip rate. For example, the rupture in the 20${}^{t}h$ cycle
(NG) does not arrest, as predicted by the fracture mechanics criterion (Fig.
2i,a,e); during the $30th$ cycle, a higher energy barrier causes a brief
deceleration (Fig. 2iii,e), which may then allow the rupture to overcome the
barrier by lowering the fracture energy. In contrast, simulations with gouge
are characterized by higher energy barriers which coincide with rupture arrest
in the simulations (Fig. 2b,d,f), with some minor discrepancies.

These arrested fronts then continue to slip aseismically ($\sim
10^{-6}\mathrm{m/s}$) within the arrest boundaries, as previously mentioned,
expanding quasi-statically as the fault is loaded until they finally overcome
the energy barriers, become unstable, and propagate seismically. This process
would explain the transition over time towards more stable mechanisms for
stress release, such as slow slip and micro-seismicity, akin to those observed
in nature Scognamiglio . (2018); Tape . (2018); Mouchon . (2023) and in the
lab Cebry . (2022); Selvadurai . (2023). In contrast, the extra case NG-X1
features system-spanning events (as in NG) due to the spatial stretching of
the fracture energy profile, while NG-X2 also exhibits a complex nucleation
pattern with numerous partial ruptures like G (see Supplementary Material).
This means that it is the shape of the fracture energy profile relative to the
energy release rate that dictates nucleation dynamics on the fault.

To summarize, the introduced coupling between friction and wear opens up a new
pathway for the fault’s build-up of heterogeneity through
$D_{\mathrm{c}}~{}(x)$, on top of the re-distribution of stress. This adds
further complexity to the nucleation of events, which transition from a single
full rupture to a complex mix of slow slip fronts and foreshocks, resulting in
richer fault dynamics.

###  3.2 Recurrence interval

We examine the evolution of the recurrence interval ($T_{i}$) where $i$ is the
cycle number. To measure it systematically, we identified the timestamps when
the maximum velocity along the fault drops below the creeping velocity
threshold $V_{\mathrm{cr}}$. This approach captures the entire cycle by
detecting fault locking following a single rupture, or a sequence of closely
spaced partial ruptures that collectively rupture the entire fault. Therefore,
it spans the time from one post-seismic phase to the next, encompassing all
intermediate phases of the seismic cycle, and is illustrated in Fig. 3a for a
given cycle in a simulation series.

The time evolution of $T_{i}$ is seen to be affected considerably by the
formation of gouge. In NG, it starts at 290 days, then drops and stabilizes at
about 270 days after 10 cycles. Its behavior is quite predictable, with a
notable increase in the standard deviation with time (see Fig. 3d). On the
other hand, simulations with wear-friction coupling ($\gamma=0.01$) show two
distinct phases. Initially, $T$ diverges from the non-wear case (NG),
stabilizing at a 10-day longer average. Subsequently (after around 25 cycles),
it becomes shorter and more chaotic, with large variances observed in specific
simulation sets. By the $30^{\mathrm{th}}$ cycle, it becomes quite erratic and
unpredictable, showing pronounced fluctuations. This is in agreement with our
previous reasoning about partial ruptures, which become more numerous with the
number of cycles (see Fig. 3b) and can alter the recurrence interval by
changing the stress distribution of the fault mid-cycle.

Figure 3: Evolution of the recurrence interval over multiple cycles. a)
Maximum slip rate over time of intermediate cycle (20) for both cases.
Vertical dashed lines show the end of the cycle, while horizontal lines
indicate the creeping (bottom) and dynamic (top) threshold velocities. b)
Evolution of rupture size distribution averaged over simulations in 4-sized
groups, together with the $75^{\mathrm{th}}$ percentile. Rupture lengths are
estimated as contiguous parts of the fault slipping above $V_{\mathrm{dyn}}$.
c) Average number of ruptures per cycle. Shaded area shows the standard
deviation. d) Evolution of the average recurrence interval across cycles for
both NG and G. Dashed lines show specific cases corresponding to the set of
simulations shown in a.

The main explanation for these results comes from the growing heterogeneity in
the fault as a result of the wear coupling. In spring-slider models that
ignore the effect of spatial variability Bizzarri (2010, 2012), the recurrence
interval is seen to decrease with gouge due to the effect of $D_{\mathrm{c}}$
on the critical stiffness of the spring ($k_{\mathrm{cr}}$). We observe a
similar reduction of the recurrence interval in NG-X1 (see Supplementary
Material). However, our results for G reveal the importance of heterogeneity.
With the increasing amplitude of $D_{\mathrm{c}}~{}(x)$ and
$\sigma_{\mathrm{n}}~{}(x)$ following the distribution of accumulated slip
(Fig. 2f) generated by the self-reinforced wear and stress transfer mechanisms
fed by frictional instabilities, it becomes more difficult for slip fronts to
propagate across the fault. This is because they have to overcome larger
energy barriers. As a consequence, instead of a full rupture breaking the
whole interface, partial ruptures become more frequent capturing part of the
energy budget and saving it for later ruptures. We can clearly see this in
Figs. 3b,c), where we count ruptures by taking the peaks of the maximum
rupture length that coincide with peaks in the maximum slip rate going over
the dynamic threshold. Hence, in the beginning, $T^{\mathrm{G}}$ becomes
larger than $T^{\mathrm{NG}}$, since asperities are tougher, due to the
accumulation of gouge and its effect on the fracture energy through
$D_{\mathrm{c}}~{}(x)$ (see Eq.8), but the differences in fracture energy are
not too large, so failure can propagate with ease. Basically, the failure
mechanism of G is initially the same as that of NG but slightly delayed by the
toughening of individual asperities. As fracture energy differences become
larger, partial ruptures start to proliferate (Fig. 3b,c), delaying the main
shock of the current cycle Turner . (2024) but also creating stress
concentrations where the main shock will eventually nucleate.

When the energy released by the failure of the weakest asperities in the fault
becomes too small compared to the fracture energy of the strongest and
toughest ones, these ruptures fail to propagate throughout the entire
interface (see Fig. 2f). Instead, they are arrested by the tougher asperities,
loading them for the next rupture to nucleate. In principle, one would expect
that having multiple ruptures instead of a single full rupture would delay the
end of the cycle, but in reality, it could also anticipate it, depending on
the stress state of the fault. This would explain why the recurrence interval
becomes so unpredictable. Initially, when the differences in fracture energy
are not too large, $T$ depends exclusively on the fracture energy of the first
asperity to fail, since this one unchains the instability of the whole
interface (single full rupture). However, as the spread in fracture energy
between asperities grows (see Fig. 2f), the first failing asperity is not
enough to break the entire interface anymore. This adds complexity to the
recurrence interval, which starts to depend on the completion of multiple
ruptures (see Fig. 3c), further contributing to the unpredictability. This
behavior can also be observed in NG-X2 (see Table 1), where the recurrence
interval is generally longer from the beginning, with a downward tendency, and
then turns erratic after the $20^{th}$ cycle (see Fig.S5c in the Supplementary
Material).

In conclusion, the evolution of the recurrence interval can be split into two
phases: A growing phase (with respect to NG) and a chaotic phase. Initially,
$T$ becomes longer than in the NG case as a result of individual asperities
becoming tougher by increasing their critical slip distance, and full ruptures
dominate. When the spread in fracture energy goes over a certain threshold,
not all instabilities are able to break through the entire interface, making
the recurrence interval dependent on multiple ruptures instead of one and,
therefore, more unpredictable.

###  3.3 Partitioning of the moment budget

We also analyze the distribution of the moment budget between seismic and
aseismic components and how it evolves with time. We compute the moment per
unit length ($M_{\mathrm{l}}$) from the slip rate using

$M_{\mathrm{l}}=\mu\iint V(x,t)dxdt~{},$ (13)

and compute its distribution over the sliding velocity for both cases at
different cycles averaged over the set of simulations (see Figs. 4a,b).

Both NG and G show two peaks (or modes) that coincide with the creeping
velocity $V_{\mathrm{cr}}$ and the dynamic velocity $V_{\mathrm{dyn}}$. The
latter is used to split the moment in two ranges: aseismic (ASM) for
$V<V_{\mathrm{dyn}}$, and seismic (SM) for $V\geq V_{\mathrm{dyn}}$. G shows,
as expected, a clear shift in the slip-rate distribution of released moment
from seismic to aseismic components between earlier (Fig. 4a) and later cycles
(Fig. 4b). This is again explained by the rise of heterogeneity and complexity
in the system, favoring more stable mechanisms for stress release, such as
slow slip and micro-seismicity.

Figure 4: Comparative partitioning of the moment budget per unit length for
both cases. a, b) Average distribution of the moment over a whole range of
slip rates: aseismic (ASM) and seismic (SM), for cycles 12 to 34 and cases NG
(a) and G (b). c) Evolution of both slip rate components (Aseismic [ASM] and
Seismic [SM]) across cycles (left axis), as well as the moment per cycle per
unit length (right axis), for both NG and G cases.

Fig. 4c shows, on its left axis, the partitioning of the fault’s moment per
unit length among aseismic (ASM) and seismic (SM) velocity ranges for 12
different cycles along the sequence. As expected, G experiences an increase of
its aseismic component at the expense of the seismic one, which decreases by
10% of the total moment. Conversely, the percentage of each component for NG
remains rather constant, with the aseismic oscillating around 32 %, as in
NG-X1 (See Table 1), which has the same relative shape of $G_{\mathrm{c}}$ and
$G_{\mathrm{s}}$. In NG-X2, however, the partitioning also remains constant
but the aseismic component takes over with over 60% of the total released
moment. This is due to its squeezed fracture energy profile, which generates
much more aseismic activity resulting from multiple rupture arrests.

Finally, the right axis of Fig. 4c shows the seismic moment per cycle and unit
length ($M_{l}^{i}$), which is halved in G with respect to NG after the
$24^{\mathrm{th}}$ cycle, corresponding to a 0.2 reduction in magnitude
($\Delta M_{\mathrm{w}}=-0.2$) over 12 cycles, according to Kanamori’s formula
for the moment magnitude scale Kanamori (1977). This is in contrast to what is
reported by Tal (2023), where they see that wear prevents the magnitude of
events from going down. A possible explanation for this is that, contrary to
how we model wear, where heterogeneity in the fracture energy always increases
over time, their implementation has a homogenizing effect on the fracture
energy, through the polishing of roughness and the redistribution of normal
stresses at the interface.

In summary, the enhancement of spatial variations in fracture energy as a
result of gouge accumulation promotes rupture arrest on the fault, which in
turn favors aseismic slip over seismic, reducing the amount of energy
dissipated during the seismic cycle.

##  4 Conclusion

In this study, we explored how gouge production influences the mechanical
behavior and seismogenic potential of rough faults. By integrating rate-and-
state friction with Archard’s model of wear, we examined the long-term impact
on several key aspects of fault behavior. Our findings indicate that the
recurrence interval of seismic events undergoes two distinct phases: an
initial phase of steady increase followed by a secondary phase of
unpredictability and chaos. We also observed a significant transition in the
nucleation mechanics, shifting from a single full rupture mechanism to a more
intricate, creep-dominated process involving increased microseismic activity.
Additionally, this shift towards slower slip processes results in a notable
reduction of the seismic moment per cycle released by the fault. This
highlights the crucial role of frictional wear in altering fault dynamics over
time.

##  5 Open Research

The simulation data generated in this study has been deposited in the ETH
Research Collection database under accession code ethz-b-XXXXXX

###### Acknowledgements.

This work was supported by an ETH Zurich Doc.Mobility Fellowship. E.M.
acknowledges funding from the Swiss National Science Foundation under grant
P500PN 202863.

## References

  * Aagaard Heaton (2008) aagaardConstrainingFaultConstitutive2008Aagaard, BT. Heaton, TH.  2008\.  Constraining Fault Constitutive Behavior with Slip and Stress Heterogeneity Constraining fault constitutive behavior with slip and stress heterogeneity.  Journal of Geophysical Research: Solid Earth113B4.  10.1029/2006JB004793
  * Abe Mair (2009) abeEffectsGougeFragment2009Abe, S. Mair, K.  2009\.  Effects of Gouge Fragment Shape on Fault Friction: New 3D Modelling Results Effects of gouge fragment shape on fault friction: New 3D modelling results.  Geophysical Research Letters3623.  10.1029/2009GL040684
  * Aghababaei . (2016) aghababaeiCriticalLengthScale2016Aghababaei, R., Warner, DH. Molinari, JF.  201606\.  Critical Length Scale Controls Adhesive Wear Mechanisms Critical length scale controls adhesive wear mechanisms.  Nature Communications7111816.  10.1038/ncomms11816
  * Albertini . (2021) albertiniStochasticPropertiesStatic2021Albertini, G., Karrer, S., Grigoriu, MD. Kammer, DS.  202102\.  Stochastic Properties of Static Friction Stochastic properties of static friction.  Journal of the Mechanics and Physics of Solids147104242.  10.1016/j.jmps.2020.104242
  * Andrews (1976) andrewsRuptureVelocityPlane1976Andrews, DJ.  1976\.  Rupture Velocity of Plane Strain Shear Cracks Rupture velocity of plane strain shear cracks.  Journal of Geophysical Research (1896-1977)81325679–5687.  10.1029/JB081i032p05679
  * Archard (1953) archardContactRubbingFlat1953Archard, JF.  195308\.  Contact and Rubbing of Flat Surfaces Contact and Rubbing of Flat Surfaces.  Journal of Applied Physics248981–988.  10.1063/1.1721448
  * Beeler (2023) beelerScaleDependenceFaultSurface2023Beeler, NM.  2023\.  On the Scale-Dependence of Fault Surface Roughness On the Scale-Dependence of Fault Surface Roughness.  Journal of Geophysical Research: Solid Earth1282e2022JB024856.  10.1029/2022JB024856
  * Bizzarri (2010) bizzarriRecurrenceEarthquakesRole2010Bizzarri, A.  2010\.  On the Recurrence of Earthquakes: Role of Wear in Brittle Faulting On the recurrence of earthquakes: Role of wear in brittle faulting.  Geophysical Research Letters3720.  10.1029/2010GL045480
  * Bizzarri (2012) bizzarriWhatCanPhysical2012Bizzarri, A.  201212\.  What Can Physical Source Models Tell Us about the Recurrence Time of Earthquakes? What can physical source models tell us about the recurrence time of earthquakes?  Earth-Science Reviews1154304–318.  10.1016/j.earscirev.2012.10.004
  * Boneh Reches (2018) bonehGeotribologyFrictionWear2018Boneh, Y. Reches, Z.  201805\.  Geotribology - Friction, Wear, and Lubrication of Faults Geotribology - Friction, wear, and lubrication of faults.  Tectonophysics733171–181.  10.1016/j.tecto.2017.11.022
  * Broz . (2006) brozMicrohardnessToughnessModulus2006Broz, ME., Cook, RF. Whitney, DL.  200601\.  Microhardness, Toughness, and Modulus of Mohs Scale Minerals Microhardness, toughness, and modulus of Mohs scale minerals.  American Mineralogist911135–142.  10.2138/am.2006.1844
  * Bruhat . (2020) bruhatInfluenceFaultRoughness2020Bruhat, L., Klinger, Y., Vallage, A. Dunham, EM.  202003\.  Influence of Fault Roughness on Surface Displacement: From Numerical Simulations to Coseismic Slip Distributions Influence of fault roughness on surface displacement: From numerical simulations to coseismic slip distributions.  Geophysical Journal International22031857–1877.  10.1093/gji/ggz545
  * B. Carpenter . (2012) carpenterFrictionalPropertiesSliding2012Carpenter, B., Saffer, D. Marone, C.  201208\.  Frictional Properties and Sliding Stability of the San Andreas Fault from Deep Drill Core Frictional properties and sliding stability of the San Andreas fault from deep drill core.  Geology408759–762.  10.1130/G33007.1
  * BM. Carpenter . (2014) carpenterFrictionalHeterogeneitiesCarbonatebearing2014Carpenter, BM., Scuderi, MM., Collettini, C. Marone, C.  2014\.  Frictional Heterogeneities on Carbonate-Bearing Normal Faults: Insights from the Monte Maggio Fault, Italy Frictional heterogeneities on carbonate-bearing normal faults: Insights from the Monte Maggio Fault, Italy.  Journal of Geophysical Research: Solid Earth119129062–9076.  10.1002/2014JB011337
  * Castellano . (2023) castellanoNucleationFrictionalSlip2023Castellano, M., Lorez, F. Kammer, DS.  202304\.  Nucleation of Frictional Slip: A Yielding or a Fracture Process? Nucleation of frictional slip: A yielding or a fracture process?  Journal of the Mechanics and Physics of Solids173105193.  10.1016/j.jmps.2022.105193
  * Cattania Segall (2021) cattaniaPrecursorySlowSlip2021Cattania, C. Segall, P.  2021\.  Precursory Slow Slip and Foreshocks on Rough Faults Precursory Slow Slip and Foreshocks on Rough Faults.  Journal of Geophysical Research: Solid Earth1264e2020JB020430.  10.1029/2020JB020430
  * Cebry . (2022) cebryCreepFrontsComplexity2022Cebry, SBL., Ke, CY., Shreedharan, S., Marone, C., Kammer, DS. McLaskey, GC.  202211\.  Creep Fronts and Complexity in Laboratory Earthquake Sequences Illuminate Delayed Earthquake Triggering Creep fronts and complexity in laboratory earthquake sequences illuminate delayed earthquake triggering.  Nature Communications1316839.  10.1038/s41467-022-34397-0
  * Cocco . (2023) coccoFractureEnergyBreakdown2023Cocco, M., Aretusini, S., Cornelio, C., Nielsen, SB., Spagnuolo, E., Tinti, E. Di Toro, G.  2023\.  Fracture Energy and Breakdown Work During Earthquakes Fracture Energy and Breakdown Work During Earthquakes.  Annual Review of Earth and Planetary Sciences511217–252.  10.1146/annurev-earth-071822-100304
  * Dieterich (1979) dieterichModelingRockFriction1979Dieterich, JH.  1979\.  Modeling of Rock Friction: 1. Experimental Results and Constitutive Equations Modeling of rock friction: 1. Experimental results and constitutive equations.  Journal of Geophysical Research: Solid Earth84B52161–2168.  10.1029/JB084iB05p02161
  * Duarte . (2009) duarteNoiseSlidingFriction2009Duarte, M., Vragovic, I., Molina, JM., Prieto, R., Narciso, J. Louis, E.  200901\.  $1/F$ Noise in Sliding Friction under Wear Conditions: The Role of Debris $1/f$ Noise in Sliding Friction under Wear Conditions: The Role of Debris.  Physical Review Letters1024045501.  10.1103/PhysRevLett.102.045501
  * Dublanchet (2018) dublanchetDynamicsEarthquakePrecursors2018Dublanchet, P.  201802\.  The Dynamics of Earthquake Precursors Controlled by Effective Friction The dynamics of earthquake precursors controlled by effective friction.  Geophysical Journal International2122853–871.  10.1093/gji/ggx438
  * Ferdowsi Rubin (2021) ferdowsiSlideHoldSlideProtocolsFrictional2021Ferdowsi, B. Rubin, AM.  2021\.  Slide-Hold-Slide Protocols and Frictional Healing in Discrete Element Method (DEM) Simulations of Granular Fault Gouge Slide-Hold-Slide Protocols and Frictional Healing in Discrete Element Method (DEM) Simulations of Granular Fault Gouge.  Journal of Geophysical Research: Solid Earth12612e2021JB022125.  10.1029/2021JB022125
  * Frérot . (2018) frerotMechanisticUnderstandingWear2018Frérot, L., Aghababaei, R. Molinari, JF.  201805\.  A Mechanistic Understanding of the Wear Coefficient: From Single to Multiple Asperities Contact A mechanistic understanding of the wear coefficient: From single to multiple asperities contact.  Journal of the Mechanics and Physics of Solids114172–184.  10.1016/j.jmps.2018.02.015
  * Garagash Germanovich (2012) garagashNucleationArrestDynamic2012Garagash, DI. Germanovich, LN.  2012\.  Nucleation and Arrest of Dynamic Slip on a Pressurized Fault Nucleation and arrest of dynamic slip on a pressurized fault.  Journal of Geophysical Research: Solid Earth117B10.  10.1029/2012JB009209
  * Hughes . (2020) hughesFrictionalBehaviourWear2020Hughes, A., Kendrick, JE., Lamur, A., Wadsworth, FB., Wallace, PA., Di Toro, G. Lavallée, Y.  202011\.  Frictional Behaviour, Wear and Comminution of Synthetic Porous Geomaterials Frictional Behaviour, Wear and Comminution of Synthetic Porous Geomaterials.  Frontiers in Earth Science8.  10.3389/feart.2020.562548
  * Kanamori (1977) kanamoriEnergyReleaseGreat1977Kanamori, H.  1977\.  The Energy Release in Great Earthquakes The energy release in great earthquakes.  Journal of Geophysical Research (1896-1977)82202981–2987.  10.1029/JB082i020p02981
  * Lapusta Rice (2003) lapustaNucleationEarlySeismic2003Lapusta, N. Rice, JR.  2003\.  Nucleation and Early Seismic Propagation of Small and Large Events in a Crustal Earthquake Model Nucleation and early seismic propagation of small and large events in a crustal earthquake model.  Journal of Geophysical Research: Solid Earth108B4.  10.1029/2001JB000793
  * Lawn (1993) lawnFractureBrittleSolids1993Lawn, B.  1993\.  Fracture of Brittle Solids Fracture of Brittle Solids (2 ).  CambridgeCambridge University Press.  10.1017/CBO9780511623127
  * Mair Abe (2008) mair3DNumericalSimulations2008Mair, K. Abe, S.  200809\.  3D Numerical Simulations of Fault Gouge Evolution during Shear: Grain Size Reduction and Strain Localization 3D numerical simulations of fault gouge evolution during shear: Grain size reduction and strain localization.  Earth and Planetary Science Letters274172–81.  10.1016/j.epsl.2008.07.010
  * Marone (1998) maroneLaboratoryDerivedFrictionLaws1998Marone, C.  1998\.  Laboratory-Derived Friction Laws and Their Application to Seismic Faulting Laboratory-Derived Friction Laws and Their Application to Seismic Faulting.  Annual Review of Earth and Planetary Sciences261643–696.  10.1146/annurev.earth.26.1.643
  * Marone Kilgore (1993) maroneScalingCriticalSlip1993Marone, C. Kilgore, B.  199304\.  Scaling of the Critical Slip Distance for Seismic Faulting with Shear Strain in Fault Zones Scaling of the critical slip distance for seismic faulting with shear strain in fault zones.  Nature3626421618–621.  10.1038/362618a0
  * Marone . (1990) maroneFrictionalBehaviorConstitutive1990Marone, C., Raleigh, CB. Scholz, CH.  1990\.  Frictional Behavior and Constitutive Modeling of Simulated Fault Gouge Frictional behavior and constitutive modeling of simulated fault gouge.  Journal of Geophysical Research: Solid Earth95B57007–7025.  10.1029/JB095iB05p07007
  * McBeck . (2021) mcbeckDecryptingHealedFault2021McBeck, J., Mair, K. Renard, F.  202104\.  Decrypting Healed Fault Zones: How Gouge Production Reduces the Influence of Fault Roughness Decrypting healed fault zones: How gouge production reduces the influence of fault roughness.  Geophysical Journal International2252759–774.  10.1093/gji/ggab003
  * Milanese . (2019) milaneseEmergenceSelfaffineSurfaces2019Milanese, E., Brink, T., Aghababaei, R. Molinari, JF.  201903\.  Emergence of Self-Affine Surfaces during Adhesive Wear Emergence of self-affine surfaces during adhesive wear.  Nature Communications1011116.  10.1038/s41467-019-09127-8
  * Milanese . (2020) milaneseRoleInterfacialAdhesion2020Milanese, E., Brink, T., Aghababaei, R. Molinari, JF.  202010\.  Role of Interfacial Adhesion on Minimum Wear Particle Size and Roughness Evolution Role of interfacial adhesion on minimum wear particle size and roughness evolution.  Physical Review E1024043001.  10.1103/PhysRevE.102.043001
  * Mollon . (2023) mollonLaboratoryEarthquakesSimulations2023Mollon, G., Aubry, J. Schubnel, A.  2023\.  Laboratory Earthquakes Simulations—Typical Events, Fault Damage, and Gouge Production Laboratory Earthquakes Simulations—Typical Events, Fault Damage, and Gouge Production.  Journal of Geophysical Research: Solid Earth1282e2022JB025429.  10.1029/2022JB025429
  * Mora Place (1998) moraNumericalSimulationEarthquake1998Mora, P. Place, D.  1998\.  Numerical Simulation of Earthquake Faults with Gouge: Toward a Comprehensive Explanation for the Heat Flow Paradox Numerical simulation of earthquake faults with gouge: Toward a comprehensive explanation for the heat flow paradox.  Journal of Geophysical Research: Solid Earth103B921067–21089.  10.1029/98JB01490
  * Mouchon . (2023) mouchonSubdailySlowFault2023Mouchon, C., Frank, WB., Radiguet, M., Poli, P. Cotte, N.  2023\.  Subdaily Slow Fault Slip Dynamics Captured by Low-Frequency Earthquakes Subdaily Slow Fault Slip Dynamics Captured by Low-Frequency Earthquakes.  AGU Advances44e2022AV000848.  10.1029/2022AV000848
  * Nielsen . (2016) nielsenFractureEnergyFriction2016Nielsen, S., Spagnuolo, E., Violay, M., Smith, S., Di Toro, G. Bistacchi, A.  201610\.  G: Fracture Energy, Friction and Dissipation in Earthquakes G: Fracture energy, friction and dissipation in earthquakes.  Journal of Seismology2041187–1205.  10.1007/s10950-016-9560-1
  * Ostapchuk . (2021) ostapchukAcousticEmissionReveals2021Ostapchuk, A., Morozova, K., Markov, V., Pavlov, D. Popov, M.  2021\.  Acoustic Emission Reveals Multiple Slip Modes on a Frictional Fault Acoustic Emission Reveals Multiple Slip Modes on a Frictional Fault.  Frontiers in Earth Science9.
  * Perfettini Ampuero (2008) perfettiniDynamicsVelocityStrengthening2008Perfettini, H. Ampuero, JP.  2008\.  Dynamics of a Velocity Strengthening Fault Region: Implications for Slow Earthquakes and Postseismic Slip Dynamics of a velocity strengthening fault region: Implications for slow earthquakes and postseismic slip.  Journal of Geophysical Research: Solid Earth113B9.  10.1029/2007JB005398
  * Power . (1988) powerRoughnessWearBrittle1988Power, WL., Tullis, TE. Weeks, JD.  1988\.  Roughness and Wear during Brittle Faulting Roughness and wear during brittle faulting.  Journal of Geophysical Research: Solid Earth93B1215268–15278.  10.1029/JB093iB12p15268
  * Radiguet . (2013) radiguetSurvivalHeterogeneousStress2013Radiguet, M., Kammer, DS., Gillet, P. Molinari, JF.  201310\.  Survival of Heterogeneous Stress Distributions Created by Precursory Slip at Frictional Interfaces Survival of Heterogeneous Stress Distributions Created by Precursory Slip at Frictional Interfaces.  Physical Review Letters11116164302.  10.1103/PhysRevLett.111.164302
  * Reches . (2019) rechesEnergyfluxControlSteadystate2019Reches, Z., Zu, X. Carpenter, BM.  201907\.  Energy-Flux Control of the Steady-State, Creep, and Dynamic Slip Modes of Faults Energy-flux control of the steady-state, creep, and dynamic slip modes of faults.  Scientific Reports9110627.  10.1038/s41598-019-46922-1
  * Renard . (2013) renardConstantDimensionalityFault2013Renard, F., Candela, T. Bouchaud, E.  2013\.  Constant Dimensionality of Fault Roughness from the Scale of Micro-Fractures to the Scale of Continents Constant dimensionality of fault roughness from the scale of micro-fractures to the scale of continents.  Geophysical Research Letters40183–87.  10.1029/2012GL054143
  * Rubin Ampuero (2005) rubinEarthquakeNucleationAging2005Rubin, AM. Ampuero, JP.  2005\.  Earthquake Nucleation on (Aging) Rate and State Faults Earthquake nucleation on (aging) rate and state faults.  Journal of Geophysical Research: Solid Earth110B11.  10.1029/2005JB003686
  * Ruina (1983) ruinaSlipInstabilityState1983Ruina, A.  1983\.  Slip Instability and State Variable Friction Laws Slip instability and state variable friction laws.  Journal of Geophysical Research: Solid Earth88B1210359–10370.  10.1029/JB088iB12p10359
  * Sáez Lecampion (2024) saezFluiddrivenSlowSlip2024Sáez, A. Lecampion, B.  202402\.  Fluid-Driven Slow Slip and Earthquake Nucleation on a Slip-Weakening Circular Fault Fluid-driven slow slip and earthquake nucleation on a slip-weakening circular fault.  Journal of the Mechanics and Physics of Solids183105506.  10.1016/j.jmps.2023.105506
  * Sammis Sornette (2002) sammisPositiveFeedbackMemory2002Sammis, CG. Sornette, D.  200202\.  Positive Feedback, Memory, and the Predictability of Earthquakes Positive feedback, memory, and the predictability of earthquakes.  Proceedings of the National Academy of Sciences99suppl_12501–2508.  10.1073/pnas.012580999
  * Schär . (2021) scharNucleationFrictionalSliding2021Schär, S., Albertini, G. Kammer, DS.  202108\.  Nucleation of Frictional Sliding by Coalescence of Microslip Nucleation of frictional sliding by coalescence of microslip.  International Journal of Solids and Structures225111059.  10.1016/j.ijsolstr.2021.111059
  * Scholz (1988) scholzCriticalSlipDistance1988Scholz, CH.  198812\.  The Critical Slip Distance for Seismic Faulting The critical slip distance for seismic faulting.  Nature3366201761–763.  10.1038/336761a0
  * Scognamiglio . (2018) scognamiglioComplexFaultGeometry2018Scognamiglio, L., Tinti, E., Casarotti, E., Pucci, S., Villani, F., Cocco, M.Dreger, D.  2018\.  Complex Fault Geometry and Rupture Dynamics of the MW 6.5, 30 October 2016, Central Italy Earthquake Complex Fault Geometry and Rupture Dynamics of the MW 6.5, 30 October 2016, Central Italy Earthquake.  Journal of Geophysical Research: Solid Earth12342943–2964.  10.1002/2018JB015603
  * Segall Bradley (2012) segallRoleThermalPressurization2012Segall, P. Bradley, AM.  201205\.  The Role of Thermal Pressurization and Dilatancy in Controlling the Rate of Fault Slip The Role of Thermal Pressurization and Dilatancy in Controlling the Rate of Fault Slip.  Journal of Applied Mechanics793031013.  10.1115/1.4005896
  * Selvadurai . (2023) selvaduraiModelingFrictionalPrecursory2023Selvadurai, PA., Galvez, P., Mai, PM. Glaser, SD.  202301\.  Modeling Frictional Precursory Phenomena Using a Wear-Based Rate- and State-Dependent Friction Model in the Laboratory Modeling frictional precursory phenomena using a wear-based rate- and state-dependent friction model in the laboratory.  Tectonophysics847229689.  10.1016/j.tecto.2022.229689
  * Tal (2023) talModelingEarthquakeCycles2023Tal, Y.  202307\.  Modeling Earthquake Cycles and Wear on Rough Faults with the Mortar Finite Element Method Modeling earthquake cycles and wear on rough faults with the mortar finite element method.  Geophysical Journal International2341190–209.  10.1093/gji/ggad063
  * Tape . (2018) tapeEarthquakeNucleationFault2018Tape, C., Holtkamp, S., Silwal, V., Hawthorne, J., Kaneko, Y., Ampuero, JP.West, ME.  201807\.  Earthquake Nucleation and Fault Slip Complexity in the Lower Crust of Central Alaska Earthquake nucleation and fault slip complexity in the lower crust of central Alaska.  Nature Geoscience117536–541.  10.1038/s41561-018-0144-2
  * Turner . (2024) turnerPartialRupturesCannot2024Turner, AR., Hawthorne, JC. Cattania, C.  2024\.  Partial Ruptures Cannot Explain the Long Recurrence Intervals of Repeating Earthquakes Partial Ruptures Cannot Explain the Long Recurrence Intervals of Repeating Earthquakes.  Journal of Geophysical Research: Solid Earth1291e2023JB027870.  10.1029/2023JB027870
  * Wang Scholz (1994) wangWearProcessesFrictional1994Wang, W. Scholz, CH.  1994\.  Wear Processes during Frictional Sliding of Rock: A Theoretical and Experimental Study Wear processes during frictional sliding of rock: A theoretical and experimental study.  Journal of Geophysical Research: Solid Earth99B46789–6799.  10.1029/93JB02875