JOURNAL OF GEOPHYSICAL RESEARCH
Supporting Information for ”Scale dependence of
1
earthquake rupture prestress in models with
2
enhanced weakening: Implications for event statistics
3
and inferences of fault stress”
4
Val`ere Lambert
1
Nadia Lapusta
1
,
2
Daniel Faulkner
3
1
Seismological Laboratory, California Institute of Technology, Pasadena, California, USA
5
1
,
2
Department of Mechanical and Civil Engineering, California Institute of Technology, Pasadena, California, USA
6
3
School of Environmental Sciences, University of Liverpool, Liverpool, UK
7
Contents of this file
8
1. Text S1 to S2
9
2. Figures S1 to S6
10
3. Table S1
11
S1. Methodology for simulations of sequences of earthquakes and aseismic slip
12
with and without the thermal pressurization of pore fluids
13
In order to conduct numerical simulations of sequences of spontaneous earthquakes and
14
aseismic slip, we utilize the spectral boundary integral method to solve the elastodynamic
15
equations of motion with the friction boundary conditions, including the evolution of pore
16
fluid pressure and temperature on the fault coupled with off-fault diffusion (Lapusta et al.,
17
2000; Noda & Lapusta, 2010). Our fault models are governed by a form of the laboratory-
18
July 21, 2021, 10:21am
X - 2
:
derived Dieterich-Ruina rate-and-state friction law regularized for zero and negative slip
19
rates, with the state evolution governed by the aging law (Rice & Ben-Zion, 1996; Noda
20
& Lapusta, 2010). The most commonly used formulation of rate-and-state laws is the
21
Dieterich-Ruina formulation (Dieterich, 1979; Ruina, 1983):
22
τ
=
σf
(
V,θ
) = (
σ
−
p
)
[
f
∗
+
a
ln
V
V
∗
+
b
ln
θV
∗
D
RS
]
,
(S1)
where
f
∗
is a reference steady-state friction coefficient at reference sliding rate
V
∗
,
D
RS
23
is the characteristic slip distance, and
a
and
b
are the direct effect and evolution effect
24
parameters, respectively. During steady-state sliding (
̇
θ
= 0), the friction coefficient is
25
expressed as:
26
f
ss
(
V
) =
f
∗
+ (
a
−
b
) ln
V
V
∗
,
(S2)
where the combination of frictional properties (
a
−
b
)
>
0 results in steady-state velocity-
27
strengthening (VS) behavior, where stable slip is expected, and properties resulting in
28
(
a
−
b
)
<
0 lead to steady-state velocity-weakening (VW) behavior, where accelerating
29
slip and hence stick-slip occur for sufficiently large regions.
30
31
The peak shear stress during dynamic rupture propagation can correspond to a much
32
higher apparent friction coefficient than the reference friction coefficient
f
∗
or the similar
33
steady-state friction coefficient at seismic slip rates of the order of 1 m/s. Assuming that
34
the fault has been locked interseismically with the state variable healing to a value
θ
int
35
and the slip rate rapidly accelerates to the peak slip rate
V
peak
upon arrival of the rupture
36
July 21, 2021, 10:21am
:
X - 3
front with negligible evolution of the state variable
θ
≈
θ
int
, the peak friction can be
37
approximately given as:
38
τ
peak
/
(
σ
−
p
int
) =
f
∗
+
a
ln
V
peak
V
∗
+
b
ln
θ
int
θ
ss
(
V
∗
)
(S3)
=
τ
ss
(
V
peak
)
(
σ
−
p
int
)
+
b
ln
θ
int
θ
ss
(
V
peak
)
=
τ
ss
(
V
pl
)
(
σ
−
p
int
)
+ (
a
−
b
) ln
V
peak
V
pl
+
b
ln
θ
int
θ
ss
(
V
peak
)
Note that
V
peak
V
∗
V
pl
and
θ
int
θ
ss
(
V
∗
)
θ
ss
(
V
peak
) for typical seismic slip rates
39
and interseismic durations of healing. The last two terms on the third line gives the dif-
40
ference between the local SSQS shear resistance described in the main text and the peak
41
shear resistance, where the last term typically dominates for periods of extending healing
42
and higher values of
θ
int
. Consequently, for a given dynamic slip rate
V
peak
, the better
43
healed the interface with higher
θ
ini
, the higher the peak friction during dynamic rupture
44
(Lambert & Lapusta, 2020).
45
46
The standard Dieterich-Ruina formulation (equation S1) has been empirically-
47
determined from laboratory experiments at sliding rates between 10
−
9
m/s to around
48
10
−
3
m/s. Under the standard logarithmic formulation, friction becomes negative as the
49
slip rate
V
approaches zero and is undefined for zero or negative slip rates (Figure S5).
50
The standard formulation may be regularized near
V
= 0 such that the shear resistance
51
remains positive for all positive values of
V
(Rice & Ben-Zion, 1996):
52
τ
(
V,θ
) =
a
σ
sinh
−
1
[
V
2
V
∗
exp
(
f
∗
+
b
log(
θV
∗
/D
RS
)
a
)]
,
(S4)
July 21, 2021, 10:21am
X - 4
:
with the steady-state shear resistance given by:
53
τ
ss
(
V
) =
a
σ
sinh
−
1
[
V
2
V
∗
exp
(
f
∗
+
b
log(
V
∗
/V
)
a
)]
.
(S5)
Theoretical justification for such regularization has been provided by drawing analogy
54
between the direct velocity effect and the exponential formulation of thermally-activated
55
creep at contact junctions, where the contact shear stress acts as a biasing factor (Rice
56
et al., 2001). The standard logarithmic rate-dependent formulation is derived when only
57
considering forward activated jumps, which may be dominant under significant shear
58
stress and conditions relevant to most laboratory experiments. The regularized formu-
59
lation (equation S4) arises when including the presence of backward jumps, which are
60
equally probable as forward jumps for
τ
= 0, as in the full thermally-activated creep the-
61
ory. The logarithmic and regularized formulations are equivalent for conditions consistent
62
with laboratory experiments, and differ only for very low slip rates (Figure S5).
63
64
Earthquakes may nucleate only if the VW region is larger than the nucleation size
h
∗
.
65
For 2D problems, two theoretical estimates of the nucleation size in mode III are (Rice &
66
Ruina, 1983; Rubin & Ampuero, 2005):
67
h
∗
RR
=
π
4
μD
RS
(
b
−
a
)(
σ
−
p
)
;
h
∗
RA
=
2
π
μD
RS
b
(
b
−
a
)
2
(
σ
−
p
)
,
(S6)
where
μ
is the shear modulus. The simulated fault in our models contains a 24-km region
68
with VW frictional properties surrounded by VS regions to create a 72-km frictional re-
69
gion. Outside of this frictional regions, the fault moves with a prescribed plate rate
V
pl
70
July 21, 2021, 10:21am
:
X - 5
to provide tectonic-like loading (Figure 2A of main text).
71
72
The thermal pressurization of pore fluids is governed in our simulations by the follow-
73
ing coupled differential equations for temperature and pore pressure evolution (Noda &
74
Lapusta, 2010):
75
∂T
(
y,z
;
t
)
∂t
=
α
th
∂
2
T
(
y,z
;
t
)
∂y
2
+
τ
(
z
;
t
)
V
(
z
;
t
)
ρc
exp(
−
y
2
/
2
w
2
)
√
2
πw
,
(S7)
∂p
(
y,z
;
t
)
∂t
=
α
hy
∂
2
p
(
y,z
;
t
)
∂y
2
+ Λ
∂T
(
y,z
;
t
)
∂t
,
(S8)
where
T
is the temperature of the pore fluid,
α
th
is the thermal diffusivity,
τV
is the shear
76
heating source distributed over a Gaussian shear layer of half-width
w
,
ρc
is the specific
77
heat,
y
is the distance normal to the fault plane,
α
hy
is the hydraulic diffusivity, and Λ
78
is the coupling coefficient that gives pore pressure change per unit temperature change
79
under undrained conditions. To approximate the effects of off-fault yielding we employ
80
a velocity limit of
V
max
= 15 m/s, as discussed in detail in Lambert et al. (2021). This
81
approximation is motivated by detailed dynamic rupture simulations with off-fault yield-
82
ing (Andrews, 2004), with the value of velocity limited corresponding to a representative
83
seismogenic depth of 10 km.
84
85
Our simulations include fault models with varying levels of ambient fluid overpressure
86
in terms of effective normal stress and as well as degrees of efficiency due to enhanced
87
weakening due to thermal pressurization. Parameters for the simulations are given in
88
Tables 1-3. Note that the stress changes associated with standard rate-and-state friction
89
have a relatively mild logarithmic dependence on slip rate and are directly proportional
90
July 21, 2021, 10:21am
X - 6
:
to the effective confining stress. As such, persistently weak rate-and-state fault models
91
with low effective normal stress and no enhanced weakening result in generally mild static
92
stress drops (
≤
2 MPa) for typical frictional parameters measured in the laboratory
93
(Figure 2 of main text). Thus, the inclusion of at least mild enhanced dynamic weakening
94
is required for fault models with low effective normal stress, such as due to substantial
95
fluid overpressurization, to produce average static stress drops between 1 - 10 MPa, as
96
typically inferred for natural earthquakes (Figures 11 of main text and S3; Lambert et
97
al., 2021).
98
99
In order to examine the prestress at the beginning of dynamic ruptures, we define the
100
beginning and end of dynamic rupture, as well as the ruptured area, based on a slip veloc-
101
ity threshold (
V
thresh
= 1 cm/s) for seismic slip. We have found in previous studies that
102
varying
V
thresh
between by 10
−
3
to 10
−
1
m/s results in minor variations of the determined
103
rupture timing and area, within 1% (Perry et al., 2020; Lambert et al., 2021).
104
105
We find that both the spatially averaged and energy-based prestress over the entire
106
length of our simulated ruptures decreases with increasing rupture size in a manner con-
107
sistent with a power-law relationship between the average prestress and the ratio of the
108
rupture length to nucleation size in our 1-D fault models (Figure S3):
109
τ
A
ini
/
(
σ
−
p
int
) =
d
2
−
d
1
log
10
[
λ
rupt
h
∗
RA
]
,
(S9)
τ
E
ini
/
(
σ
−
p
int
) =
e
2
−
e
1
log
10
[
λ
rupt
h
∗
RA
]
,
(S10)
July 21, 2021, 10:21am
:
X - 7
where
d
1
,
d
2
,
e
1
and
e
2
are fit parameters (Table S6). The intercept parameters
d
2
and
110
e
2
are comparable to the prescribed average local SSQS shear resistance divided by the
111
interseismic effective normal stress
τ
V
pl
ss
/
(
σ
−
p
int
)
≈
0
.
63. (Table S6), and
d
1
and
e
1
in-
112
crease by about a factor of 5 between our models with no enhanced weakening (RS1 and
113
2) and moderately efficient enhanced dynamic weakening (TP3 and 4).
114
115
Our fault models with more efficient enhanced dynamic weakening produce fewer smaller
116
events than those with mild to moderate enhanced weakening, as can be observed in
117
the frequency-magnitude statistics (Figure 10 of the main text). To create frequency-
118
magnitude histograms we compute the seismic moment
M
0
=
μA
δ
for ruptures, where
119
μ
is the shear modulus,
A
is the rupture area and
δ
is the average slip in the rupture.
120
As our simulations are 2-D, we compute the moment by assuming a circular rupture area
121
A
=
π
(
λ
rupt
/
2)
2
, where
λ
rupt
is the rupture length. We calculate b-values by fitting a lin-
122
ear trend to the cumulative frequency-magnitude distributions of each simulation (Figure
123
10E-H of the main text) between moment magnitude bins of 4.5 and 7.5.
124
125
S2. Single-degree-of-freedom representation of laboratory experiments
126
We compare the evolution of local slip rate and shear stress in our simulated dynamic
127
ruptures with single-degree-of-freedom (SDOF) calculations motivated by high-velocity
128
laboratory experiments that impose variable seismic slip rates to infer shear resistance
129
evolution and often compare their findings with seismological observations (Sone & Shi-
130
mamoto, 2009; Fukuyama & Mizoguchi, 2010). The SDOF calculations are governed by
131
the same rate-and-state friction with enhanced dynamic weakening due to thermal pres-
132
July 21, 2021, 10:21am
X - 8
:
surization as in our fault model TP4. Our SDOF calculations impose a slip-rate history, as
133
typically done in laboratory experiments, and solve for the evolution of shear stress, state
134
variable, temperature and pore pressure using equation 3 of the main text and equations
135
S4 and S7-8 given the initial state. We assume initial conditions where sliding has been
136
maintained until steady-state conditions at the slip rate of
V
= 0
.
1 mm/s, comparable to
137
the initial conditions of Fukuyama and Mizoguchi (2010). We then impose two different
138
slip rate functions characterized by regularized Yoffe functions (Tinti et al., 2005), with
139
total slip of 1.95 m (comparable to our simulated slip) and maximum slip rate of 2 m/s.
140
Tinti et al. (2005) regularized the stress singularity in the analytical Yoffe function by
141
convolving it with a triangular function of half-width
t
s
. The regularized Yoffe functions
142
are characterized by two time-scales, the half-width
t
s
and the rise time
t
r
. For the two
143
examples shown in Figure 9 of the main text, we choose values of
t
r
= 3s with
t
s
= 0
.
1
t
r
144
for RYF1 and
t
r
= 1
.
4s with
t
s
= 0
.
4
t
r
for RYF2, in order to compare pulses with more
145
pronounced and gradual accelerations that produce the same slip and peak slip rate.
146
July 21, 2021, 10:21am
:
X - 9
References
Andrews, D. J. (2004, 06). Rupture Models with Dynamically Determined Breakdown
147
Displacement.
Bulletin of the Seismological Society of America
,
94
(3), 769-775. doi:
148
10.1785/0120030142
149
Dieterich, J. H. (1979). Modeling of rock friction 1. Experimental results and constitutive
150
equations.
Journal of Geophysical Research
,
84
(B5), 2161-2168.
151
Fukuyama, E., & Mizoguchi, K. (2010). Constitutive parameters for earthquake rupture
152
dynamics based on high-velocity friction tests with variable sliprate.
International
153
Journal of Fracture
,
163
(1), 15–26. doi: 10.1007/s10704-009-9417-5
154
Lambert, V., & Lapusta, N. (2020). Rupture-dependent breakdown energy in fault
155
models with thermo-hydro-mechanical processes.
Solid Earth
,
11
(6), 2283-2302. doi:
156
10.5194/se-11-2283-2020
157
Lambert, V., Lapusta, N., & Perry, S. (2021). Propagation of large earthquakes as self-
158
healing pulses or mild cracks.
Nature
,
591
(7849), 252–258. doi: 10.1038/s41586-021
159
-03248-1
160
Lapusta, N., Rice, J. R., Ben-Zion, Y., & Zheng, G. (2000). Elastodynamic analysis
161
for slow tectonic loading with spontaneous rupture episodes on faults with rate-
162
and state-dependent friction.
Journal of Geophysical Research
,
105
, 765-789. doi:
163
10.1029/2000JB900250
164
Noda, H., & Lapusta, N. (2010). Three-dimensional earthquake sequence simulations
165
with evolving temperature and pore pressure due to shear heating: Effect of hetero-
166
geneous hydraulic diffusivity.
Journal of Geophysical Research
,
115
, B123414. doi:
167
10.1029/2010JB007780
168
July 21, 2021, 10:21am
X - 10
:
Perry, S. M., Lambert, V., & Lapusta, N. (2020). Nearly magnitude-invariant stress drops
169
in simulated crack-like earthquake sequences on rate-and-state faults with thermal
170
pressurization of pore fluids.
Journal of Geophysical Research: Solid Earth
,
125
(3),
171
e2019JB018597. doi: 10.1029/2019JB018597
172
Rice, J. R., & Ben-Zion, Y. (1996). Slip complexity in earthquake fault models.
Proceed-
173
ings of the National Academy of Sciences of the United States of America
,
93
(9),
174
3811–3818. doi: 10.1073/pnas.93.9.3811
175
Rice, J. R., Lapusta, N., & Ranjith, K. (2001). Rate and state dependent friction and the
176
stability of sliding between elastically deformable solids.
Journal of the Mechanics
177
and Physics of Solids
,
49
(9), 1865-1898.
178
Rice, J. R., & Ruina, A. L. (1983). Stability of steady frictional slipping.
Journal of
179
Applied Mechanics
,
50
(2), 343-349.
180
Rubin, A., & Ampuero, J.-P. (2005). Earthquake nucleation on (aging) rate and state
181
faults.
Journal of Geophysical Research: Solid Earth
,
110
(B11).
182
Ruina, A. (1983). Slip instability and state variable friction laws.
Journal of Geophysical
183
Research
,
88
(B12), 10359-10370.
184
Sone, H., & Shimamoto, T. (2009). Frictional resistance of faults during accelerating
185
and decelerating earthquake slip.
Nature Geoscience
,
2
(10), 705–708. doi: 10.1038/
186
ngeo637
187
Tinti, E., Spudich, P., & Cocco, M. (2005). Earthquake fracture energy inferred from
188
kinematic rupture models on extended faults.
Journal of Geophysical Research: Solid
189
Earth
,
110
(B12). doi: 10.1029/2005JB003644
190
July 21, 2021, 10:21am
:
X - 11
Rupture length / Theoretical nucleation size
0
10
20
30
40
Measured nucleation size
theoretical nucleation size
0.0
1.2
2.0
0.4
0.8
1.6
TP 3
TP 2
TP 1
RS 2
RS 1
TP 4
Figure S1.
The measured nucleation sizes of the simulated ruptures are comparable to
the theoretical estimate
h
∗
RA
, within a factor of 2.
July 21, 2021, 10:21am
X - 12
:
0.0
0.6
0.4
0.2
0
10
20
30
40
A)
Average local steady-state quasi-static friction in VW region
More efficient
enhanced weakening
TP 3
TP 2
TP 1
RS 2
RS 1
B)
Spa.-Averaged prestress
effective normal stress
Rupture length / Theoretical nucleation size
TP 4
Energy-based prestress
effective normal stress
0.0
0.6
0.4
0.2
0
10
20
30
40
Figure S2.
The spatially-averaged prestress
τ
A
ini
and energy-averaged prestress
τ
E
ini
are
generally comparable and decrease with increasing rupture size and efficiency of weaken-
ing.
July 21, 2021, 10:21am
:
X - 13
0.0
0.6
0.4
0.2
1
10
A)
Average local steady-state quasi-static friction in VW region
More efficient
enhanced weakening
TP 3
TP 2
TP 1
RS 2
RS 1
B)
Spa.-Averaged prestress
effective normal stress
Rupture length / Theoretical nucleation size
TP 4
Energy-based prestress
effective normal stress
0.0
0.6
0.4
0.2
20
30
40
1
10
20
30
40
Figure S3.
The spatially-averaged prestress
τ
A
ini
and energy-averaged prestress
τ
E
ini
are
generally comparable and decrease with increasing rupture size and efficiency of weaken-
ing.
July 21, 2021, 10:21am
X - 14
:
Rupture length / Theoretical nucleation size
Energy-based static stress drop
effective normal stress
0
10
20
30
40
0.0
0.6
0.4
0.2
Spa.- avg. static stress drop
effective normal stress
0
10
20
30
40
0.0
0.6
0.4
0.2
A)
B)
TP 3
TP 2
TP 1
RS 2
RS 1
TP 4
Figure S4.
The (A) spatially-averaged and (B) energy-based average static stress drops
for ruptures represent relatively mild decreases in average shear stress with respect to the
effective normal stress. Persistently weak fault models with low effective normal stress
≤
20 MPa and relatively mild weakening, such as from standard rate-and-state friction
(RS1 and RS2) produce potentially too small average static stress drops
≤
2 MPa, whereas
models with mild to moderate enhanced weakening (TP1-4) produce realistic average
static stress drops of 1 - 10 MPa.
July 21, 2021, 10:21am