of 29
1.
Introduction
Determining the absolute level and controlling factors of the stress state on faults has profound implica
-
tions for earthquake physics, seismic hazard assessment, and the role of faulting in plate tectonics and
geodynamics. Numerous lines of field evidence suggest that the average shear stress acting on mature faults
must be low, 20 MPa or less, in comparison to the expected shear resistance of 100–200 MPa averaged over
the seismogenic depth, given rock overburden and hydrostatic pore fluid pressure, along with typical qua
-
si-static friction coefficients of 0.6–0.85 (aka “Byerlee friction”) measured in laboratory experiments (Brune
Abstract
Determining conditions for earthquake slip on faults is a key goal of fault mechanics highly
relevant to seismic hazard. Previous studies have demonstrated that enhanced dynamic weakening (EDW)
can lead to dynamic rupture of faults with much lower shear stress than required for rupture nucleation.
We study the stress conditions before earthquake ruptures of different sizes that spontaneously evolve
in numerical simulations of earthquake sequences on rate-and-state faults with EDW due to thermal
pressurization of pore fluids. We find that average shear stress right before dynamic rupture (aka shear
prestress) systematically varies with the rupture size. The smallest ruptures have prestress comparable
to the local shear stress required for nucleation. Larger ruptures weaken the fault more, propagate over
increasingly under-stressed areas due to dynamic stress concentration, and result in progressively lower
average prestress over the entire rupture. The effect is more significant in fault models with more efficient
EDW. We find that, as a result, fault models with more efficient weakening produce fewer small events
and result in systematically lower b-values of the frequency-magnitude event distributions. The findings
(a) illustrate that large earthquakes can occur on faults that appear not to be critically stressed compared
to stresses required for slip nucleation; (b) highlight the importance of finite-fault modeling in relating the
local friction behavior determined in the lab to the field scale; and (c) suggest that paucity of small events
or seismic quiescence may be the observational indication of mature faults that operate under low shear
stress due to EDW.
Plain Language Summary
Understanding the evolution of stress on faults over periods
of slow and fast motion is crucial for assessing how earthquakes start, grow, and ultimately stop. Here
we use computer models to explore the stress conditions required for simulated earthquake ruptures to
occur. We find that the critical stress conditions for rupture propagation depend on the size of the rupture
and how efficiently the fault shear resistance weakens during fast slip. In particular, larger earthquakes
on faults that experience increasingly more efficient weakening during rupture can propagate under
systematically lower stress conditions than those required for rupture nucleation. As a result, we find
that faults that exhibit efficient weakening can host predominantly large earthquakes at the expense of
smaller earthquakes. Our findings illustrate how large earthquakes can occur on faults that appear to be
understressed compared to expected conditions for rupture nucleation. Moreover, our findings support
a body of work suggesting that the scarcity of small earthquakes on some major mature fault segments,
like the central section of the San Andreas Fault, may indicate that they experience substantial weakening
during large earthquakes, a consideration that may be particularly useful for earthquake early warning
systems.
LAMBERT ET AL.
© 2021. The Authors.
This is an open access article under
the terms of the
Creative Commons
Attribution
License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.
Scale Dependence of Earthquake Rupture Prestress in
Models With Enhanced Weakening: Implications for
Event Statistics and Inferences of Fault Stress
Valère Lambert
1
, Nadia Lapusta
1,2
, and Daniel Faulkner
3
1
Seismological Laboratory, California Institute of Technology, Pasadena, CA, USA,
2
Department of Mechanical
and Civil Engineering, California Institute of Technology, Pasadena, CA, USA,
3
School of Environmental Sciences,
University of Liverpool, Liverpool, UK
Key Points:
Local shear prestress varies
significantly within and among
ruptures, being close to the quasi-
static fault strength in nucleation
regions
Efficient weakening allows rupture
propagation over areas of lower
prestress, leading to lower average
prestress over larger rupture areas
Fault models with more efficient
dynamic weakening produce
fewer smaller events and result in
systematically lower
b
-values
Supporting Information:
Supporting Information may be found
in the online version of this article.
Correspondence to:
V. Lambert,
vlambert@caltech.edu
Citation:
Lambert, V., Lapusta, N., & Faulkner, D.
(2021). Scale dependence of earthquake
rupture prestress in models with
enhanced weakening: Implications for
event statistics and inferences of fault
stress.
Journal of Geophysical Research:
Solid Earth
,
126
, e2021JB021886.
https://doi.org/10.1029/2021JB021886
Received 11 FEB 2021
Accepted 12 SEP 2021
10.1029/2021JB021886
RESEARCH ARTICLE
1 of 29
Journal of Geophysical Research: Solid Earth
LAMBERT ET AL.
10.1029/2021JB021886
2 of 29
et al.,
1969
; Byerlee,
1978
; Fulton et al.,
2013
; Gao & Wang,
2014
; Henyey & Wasserburg,
1971
; Lachenbruch
& Sass,
1980
; Nankali,
2011
; Rice,
2006
; Sibson,
1975
; Suppe,
2007
; Tanikawa & Shimamoto,
2009
; Townend
& Zoback,
2004
). Such evidence includes the lack of a substantial heat flow anomaly around mature faults
that would be expected for fault slip at 100 MPa or more (Brune et al.,
1969
; Gao & Wang,
2014
; Henyey
& Wasserburg,
1971
; Lachenbruch & Sass,
1980
; Nankali,
2011
), inferences of steep angles between the
principal stress direction and fault plane (Townend & Zoback,
2004
), analyses of the fault core obtained
by drilling through shallow parts of faults that have experienced major recent events, including the 2011
w
EM
9.0 Tohoku-oki event (Fulton et al.,
2013
; Tanikawa & Shimamoto,
2009
), the geometry of thrust-belt
wedges (Suppe,
2007
), and the existence of long-lived narrow shear zones that do not exhibit any evidence
of melting (Rice,
2006
; Sibson,
1975
). Note that such evidence for apparent fault weakness pertains pre
-
dominantly to mature faults, whereas some studies suggest that smaller, less mature faults may sustain the
expected high shear stresses given Byerlee friction values and overburden minus hydrostatic pore pressure
(e.g., Townend & Zoback,
2000
).
A relatively straightforward explanation for the low-stress operation of mature faults is that they may be
persistently weak (Figure
1
), due to the presence of anomalously low quasi-static friction coefficients and/
or low effective normal stress from pervasive fluid overpressure (Bangs et al.,
2009
; Brown et al.,
2003
;
Collettini et al.,
2009
; Faulkner et al.,
2006
; Lockner et al.,
2011
). Mature fault zones typically have a layer of
fine gouge material at their core (e.g., Faulkner et al.,
2006
). However, most materials with low quasi-static
friction coefficients (less than 0.5) under laboratory conditions tend to exhibit velocity-strengthening be
-
havior (Ikari et al.,
2011
), which would preclude spontaneous nucleation of dynamic ruptures. Moreover,
while evidence of substantial fluid overpressure has been documented for many subduction zones (Bangs
et al.,
2009
; Brown et al.,
2003
), there remains much debate over the ubiquity of chronic near-lithostatic flu
-
id overpressurization along faults in other tectonic settings, such as continental faults, with some borehole
Figure 1.
Field observations suggest that the average effective friction on mature faults must be low (
E
0.1). One explanation for this inferred low effective
friction would be that mature faults are persistently weak, such as from the presence of fault materials with persistently low friction coefficients
/()
Ep

(red). Faults may also be persistently weak while having actual friction coefficients that are persistently high (
E
0.2, blue), but require substantial chronic fluid
overpressure in order to maintain low effective fault friction. A number of laboratory experiments indicate that the coefficient of friction for many materials
relevant to seismogenic faults is around 0.6–0.8 at low sliding rates, but drops dramatically to lower values (
E
0.2) at higher slip rates relevant to seismic slip,
consistent with the notion of quasi-statically strong, but dynamically weak behavior (dashed black line).
Journal of Geophysical Research: Solid Earth
LAMBERT ET AL.
10.1029/2021JB021886
3 of 29
measurements suggesting fluid pressure levels more consistent with hydrostatic conditions (Townend &
Zoback,
2000
; Zoback et al.,
2010
).
An alternative hypothesis for explaining such low-stress, low-heat operation is that mature faults are indeed
strong at slow, quasi-static sliding rates but undergo considerable enhanced dynamic weakening at seismic
slip rates, which has been widely hypothesized in theoretical studies and documented in laboratory exper
-
iments (Figure
1
, dashed black line; Acosta et al.,
2018
; Di Toro et al.,
2011
; Noda et al.,
2009
; Rice,
2006
;
Sibson,
1973
; Tsutsumi & Shimamoto,
1997
; Wibberley et al.,
2008
). The presence of enhanced dynamic
weakening on natural faults can be questioned by the expectation that enhanced dynamic weakening would
result in magnitude-dependent static stress drops (Beeler et al.,
2012
; Gao & Wang,
2014
), with larger rup
-
tures resulting in higher stress drops than typical values of 1–10 MPa inferred from earthquakes on natural
faults (Allmann & Shearer,
2009
; Ye et al.,
2016b
). The expectation is based on a common assumption that
the shear prestress over the entire rupture area should be near the static strength of the fault while the final
shear stress should be near the dynamic resistance of the fault, resulting in a large static stress change for
more efficient weakening. However, a number of numerical and laboratory studies have demonstrated that,
once nucleated, dynamic ruptures can propagate under regions with prestress conditions that are well be
-
low the expected static strength, based on prescribed or measured quasi-static friction coefficients and con
-
fining conditions (Dunham et al.,
2011
; Fineberg & Bouchbinder,
2015
; Gabriel et al.,
2012
; Lu et al.,
2010
;
Noda et al.,
2009
; Schmitt et al.,
2015
; Zheng & Rice,
1998
) while the final shear stress could be higher than
dynamic shear stress for pulse-like ruptures, with both inferences promoting reasonable stress drops. Such
studies have often considered a single dynamic rupture nucleated artificially and propagating over uniform
prestress conditions.
Recent numerical studies of earthquake sequences have shown that fault models with a combination of both
hypotheses for low-stress operation, including some chronic fluid overpressure as well as mild-to-moderate
enhanced dynamic weakening due to the thermal pressurization of pore fluids, work well for reproducing a
range of observations (Lambert et al.,
2021
; Perry et al.,
2020
). These include reasonable static stress drops
between 1 and 10 MPa nearly independent of earthquake magnitude, the seismologically inferred increase
in average breakdown energy with rupture size, the radiation ratios between 0.1 and 1 inferred for natural
events, and the heat flow constraints. The simulations produce mainly crack-like or mild pulse-like rup
-
tures, with no significant undershoot. The near magnitude-invariance of average static stress drop arises in
these fault models because enhanced dynamic weakening results in both lower average prestress and lower
average final shear stress for larger ruptures with larger slip, with the average static stress drops being near
-
ly magnitude-independent. These studies suggest that distinguishing between the conditions required for
rupture nucleation and propagation is important for assessing the relationship between laboratory friction
measurements, seismological observations and the absolute stress conditions on faults.
Here, we use and expand upon the set of numerical models from Perry et al. (
2020
) and Lambert et al. (
2021
)
to document the variability of prestress on a fault that arises from the history of previous ruptures, and to
study the relation between the size of dynamic rupture events and the average shear prestress over the rup
-
ture area. We also examine how the complexity of earthquake sequences, in terms of the variability of rup
-
ture size, differs with the efficiency of dynamic weakening. We study these behaviors in the context of sim
-
ulations of sequences of earthquakes and slow slip, which allow the prestress conditions before earthquakes
to be set by the loading conditions, evolving fault shear resistance (including weakening and healing), and
stress redistribution by prior slip, as would occur on natural faults. Moreover, our simulations resolve the
spontaneous nucleation process with the natural acceleration of slow unsteady slip prior to dynamic rup
-
ture. The constitutive relations for the evolving fault resistance and healing adopted in our models have
been formulated as a result of a large body of laboratory, field and theoretical work (e.g., Dieterich,
1979
; Di
Toro et al.,
2011
; Rice,
2006
; Ruina,
1983
; Sibson,
1973
; Wibberley et al.,
2008
). Indeed, laboratory experi
-
ments of fault shear resistance at both slow and fast slip rates have been indispensible for our understanding
of fault behavior and for formulating fault models such as the ones used in this study. The modeling allows
us to examine the implications of the laboratory-derived constitutive behaviors for the larger-scale behavior
of faults, and we compare our inferences of average shear prestress from relatively large-scale finite-fault
modeling to field measurements of crustal stresses acting on mature faults and small-scale laboratory meas
-
urements of the shear resistance of typical fault materials.
Journal of Geophysical Research: Solid Earth
LAMBERT ET AL.
10.1029/2021JB021886
4 of 29
2.
Building on Laboratory Constraints to Model Larger-Scale Fault Behavior
Laboratory experiments have been instrumental for exploring aspects of fault resistance during both slow
and fast sliding (
9
10
E
–1 m/s, Figure
1
). Experiments with slow sliding velocities (
3
10
E
m/s) are critical
for formulating fault constitutive laws that form the basis for understanding the nucleation of earthquake
ruptures. High-velocity laboratory friction experiments have demonstrated enhanced dynamic weakening
of faults and elucidated a range of mechanisms by which this dynamic weakening can occur (e.g., Acosta
et al.,
2018
; De Paola et al.,
2015
; Di Toro et al.,
2011
; Faulkner et al.,
2011
; Goldsby & Tullis,
2011
; Han
et al.,
2007
; Wibberley et al.,
2008
). Most slow- and high-velocity experiments measure or infer the relevant
quantities—slip, slip rate, shear stress etc.—averaged over the sample and examine the evolution of shear
resistance corresponding to a particular history of loading, such as imposed variations in the displacement
rate of the loading piston, and the particular fault conditions (normal stress, temperature, pore fluid pres
-
sure, etc.). Some experimental studies imposed the expected sliding motion during earthquakes in order
to directly relate laboratory stress measurements to seismological quantities, such as static stress drop and
breakdown energy (e.g., Fukuyama & Mizoguchi,
2010
; Nielsen et al.,
2016
; Sone & Shimamoto,
2009
).
To understand the full implications of the evolution of shear resistance measured in small-scale experi
-
ments for slip at larger scales along natural faults, they are synthesized into mathematical formulations and
used in numerical modeling, for the following reasons. During slipping events on a finite fault over scales of
tens of meters to kilometres—much larger than the experimental scale—the fault does not slip uniformly
with a predetermined slip-rate history. Rather, the slip event initiates on a portion of the fault and then
spreads along the fault, with varying slip-rate histories and final slips at different points along the fault.
This is captured in inversions of large earthquakes (e.g., Heaton,
1990
; Simons et al.,
2011
; Tinti et al.,
2016
;
Ye et al.,
2016a
) and, to a degree, in larger-scale experiments, sometimes involving analog materials (Lu
et al.,
2010
; McLaskey et al.,
2014
; Rubino et al.,
2017
; Svetlizky & Fineberg,
2014
; Yamashita et al.,
2015
).
In the process, the slip (a) transfers stress to the more locked portions of the fault and (b) enters portions of
the fault with different conditions—such as levels of shear pre-stress, pore fluid pressure, etc.—and poten
-
tially different friction and hydraulic properties. Hence the resulting coupled evolution of shear resistance
and slip rate at different locations on the fault is often quite different and, through stress transfer, strongly
dependent on the entire slip process at all locations throughout the rupture. These nonlinear and often
dynamic feedback processes on the scales of tens of meters to kilometers can currently be only captured
through numerical modeling.
Many numerical models of earthquake source processes utilize insight from laboratory experiments that
indicate that the resistance to shear
E
along a fault depends on the sliding rate
EV
and the quality and/or
lifetime of the local contacts, typically parameterized by a state variable
E
with units of time, as well as on
the effective normal stress
Ep

 
acting on the fault, with
E
being the normal stress and
Ep
being the
pore fluid pressure localized within the shearing layer (e.g., Dieterich,
1979
; Marone,
1998
). For continuum
problems involving frictional sliding, the motion within the continuum is governed by the balance of linear
momentum, subject to the boundary condition that tractions are given by the constitutive law of the inter
-
face. For frictional sliding without changes in the elastodynamic normal stress
E
, which is the case consid
-
ered in this work, the boundary condition reduces to the shear stress being equal to the shear resistance on
the interface (
0
Ey
):
stress
resistance
(, 0,;)
(, 0,;)
( , )(
).
xyzt
xyzt
fVp


(1)
An important concept in the rate-and-state formulation of the friction coefficient
(,)
EfV
is that the friction
coefficient is not a fixed property of the interface but evolves over time, facilitating the time-dependent
changes of shear resistance and hence shear stress along the fault during shear.
The most commonly used formulation of rate-and-state laws is the Dieterich-Ruina formulation (Dieter
-
ich,
1979
; Ruina,
1983
):
*
*
*
RS
( , )
ln
ln
,
VV
fVfab
VD




(2)
where
*
Ef
is a reference steady-state friction coefficient at the reference sliding rate
*
EV
,
RS
ED
is the character
-
istic slip distance, and
Ea
and
Eb
are the direct effect and evolution effect parameters, respectively. Our fault
Journal of Geophysical Research: Solid Earth
LAMBERT ET AL.
10.1029/2021JB021886
5 of 29
models are governed by a form of the laboratory-derived Dieterich-Ruina rate-and-state friction law regu
-
larized for zero and negative slip rates (Lapusta et al.,
2000
; Noda & Lapusta,
2010
). The evolution of the
state variable can be described by various evolution laws; we employ the aging law (Ruina,
1983
):
RS
1,
V
D
 
(3)
which describes evolution during sliding as well as time-dependent healing in near-stationary contact. In
our models, the shear resistance and shear stress also change due to the evolution of pore fluid pressure
Ep
.
We conduct numerical simulations following the methodological developments of Lapusta et al. (
2000
),
Noda and Lapusta (
2010
), and Lambert et al. (
2021
) in order to solve the elastodynamic equations of mo
-
tion with the fault boundary conditions, including the evolution of pore fluid pressure and temperature on
the fault coupled with off-fault diffusion. The simulations solve for mode III slip on a 1-D fault embedded
into a 2-D uniform, isotropic, elastic medium (Figure
2
). The potential types of slip on the fault include
sequences of earthquakes and aseismic slip (SEAS) and they are simulated in their entirety, including the
nucleation process, dynamic rupture propagation, postseismic slip that follows the event, and interseismic
period between seismic events that can last up to tens or hundreds of years and host steady and transient
slow slip (Figure
2
).
The simulated fault in our models contains a 24-km-long segment with velocity-weakening (VW) frictional
properties where earthquake ruptures may nucleate and propagate, surrounded by velocity-strengthening
(VS) segments that inhibit rupture nucleation and propagation. Our simulations include enhanced dynamic
weakening due to the thermal pressurization of pore fluids, which occurs when pore fluids within the fault
shearing layer heat up and pressurize during dynamic rupture, reducing the effective normal stress and
shear resistance (Noda & Lapusta,
2010
; Rice,
2006
; Sibson,
1973
). Thermal pressurization is one potential
mechanism for enhanced weakening; qualitatively similar results should hold for models with other types
of enhanced dynamic weakening. We follow the thermal pressurization formulation of Noda and Lapus
-
ta (
2010
) (Supporting Information
S1
). We approximate the effects of off-fault yielding by employing a limit
on the slip velocity
max
15
EV
m/s, which is motivated by detailed dynamic rupture simulations with off-
fault yielding (Andrews,
2004
) and discussed in detail in Lambert et al. (
2021
).
For the purpose of comparing local frictional behavior with the average prestress for dynamic ruptures
of varying sizes, we focus this study on simulated ruptures that arrest within the VW region, where the
friction properties are uniform with a quasi-static reference friction of 0.6, consistent with many materials
exhibiting VW behavior in laboratory experiments (Ikari et al.,
2011
). We examine the evolution of the
apparent friction coefficient, or the ratio of the current shear stress
E
to the interseismic drained effective
normal stress
int
()
Ep
, where
int
Ep
is the interseismic drained value of the pore pressure. The term “drained”
refers to the effective stress with ambient pore pressure unaffected by slip processes such as dilatancy, com
-
paction, or thermal pressurization. The values used for the direct and evolution effect parameters,
Ea
and
Eb
respectively, within the VW region are consistent with typical laboratory values (Tables
1
3
; Blanpied
et al.,
1991
,
1995
; Marone,
1998
). Within the VS regions. We use higher values of the direct effect
Ea
in order
to more efficiently stop VW-region-spanning ruptures, which helps to maintain a reasonable domain size
for computational expense. The properties of the VS regions should not appreciably alter the conclusions of
this work, as we focus on the average stress measures for “partial” ruptures that are arrested predominantly
in the VW region. The properties of larger ruptures that span the entire VW region and that arrest in the
fault regions with VS properties are more sensitive to the VS properties, as discussed in Perry et al. (
2020
).
We examine fault models with varying levels of ambient fluid overpressure in terms of the effective nor
-
mal stress, as well as varying degrees of efficiency in enhanced weakening due to thermal pressurization.
The parameter values we have chosen (Tables
1
3
) are motivated by prior studies that have reproduced a
range of seismological observations as well as low-stress, low-heat operation of mature faults (Lambert
et al.,
2021
; Perry et al.,
2020
). The parameter values also facilitate our goal of examining ruptures in fault
models with a range of efficiency in enhanced dynamic weakening. We refer to the weakening behavior
of fault models as being mild, moderate or efficient in comparison to the underlying weakening of the
standard rate-and-state friction. This can be approximately quantified as the difference between the steady-
state frictional resistance at the plate rate and seismic slip rates, equal to
()()ln()

pbaVV
int
seispl
/
, which
corresponds to about
10%
E
times the interseismic effective normal stress for
seis
1
EV
m/s. Mild weakening