of 24
Supplementary Material: Absolute stress levels in
models of low-heat faults: Links to geophysical
observables and differences for crack-like ruptures and
self-healing pulses
Val`ere Lambert
a,
, Nadia Lapusta
b,c
a
Department of Earth and Marine Sciences, University of California, Santa Cruz,
California, USA
b
Seismological Laboratory, California Institute of Technology, Pasadena, California, USA
c
Department of Mechanical and Civil Engineering, California Institute of Technology,
Pasadena, California, USA
1. Fault models with rate-and-state friction laws and thermal pres-
1
surization of pore fluids
2
Our fault models are governed by a form of the laboratory-derived Dieterich-
3
Ruina rate-and-state friction law (Dieterich, 1979; Ruina, 1983):
4
τ
=
σf
(
V,θ
) = (
σ
p
)

f
+
a
log
V
V
+
b
log
θV
D
RS

,
(1)
where
σ
is the effective normal stress given the normal stress
σ
and pore
5
fluid pressure
p
,
f
is a reference steady-state friction coefficient at reference
6
Corresponding author
Email address:
valerelambert@gmail.com
(Val`ere Lambert)
Preprint submitted to Earth and Planetary Science Letters
June 22, 2023
sliding rate
V
,
D
RS
is the characteristic slip distance, and
a
and
b
are the
7
direct effect and evolution effect parameters, respectively. In our simulations,
8
we use a form of the Dieterich-Ruina law regularized for zero and negative
9
slip rates (Rice and Ben-Zion, 1996; Lapusta et al., 2000; Noda and Lapusta,
10
2010).
11
We use the aging form of the state variable evolution (Ruina, 1983):
12
13
̇
θ
= 1
V θ
D
RS
,
(2)
which aims to incorporate the evolution of contacts during sliding as well as
14
time-dependent healing during sufficiently slow motion. During steady-state
15
sliding (
̇
θ
= 0), the friction coefficient is expressed as:
16
17
f
ss
(
V
) =
f
+ (
a
b
) log
V
V
,
(3)
where the combination of frictional properties (
a
b
)
>
0 results in steady-
18
state velocity-strengthening (VS) behavior that leads to stable slip and (
a
19
b
)
<
0 results in steady-state velocity-weakening (VW) behavior that can
20
leads to spontaneously accelerating slip and hence seismic events for suffi-
21
ciently large regions.
22
2
The local steady-state quasi-static (SSQS) shear resistance (equation 1 of
23
main text), which we find is consistent with the average stress levels within
24
the nucleation region of simulated ruptures (Supplementary Figure 2), is
25
defined based on the steady-state friction coefficient (equation 3) determined
26
at the prescribed tectonic plate loading rate
V
pl
:
27
f
ss
(
V
pl
) =
f
+ (
a
b
) log
V
pl
V
.
(4)
For the VW properties used in our models (Supplementary Table 1) and plate
28
rate of 10
9
m/s, the steady-state friction coefficient for our fault models is
29
0.63 which is slightly higher than the reference friction coefficient
f
of 0.6.
30
Earthquake ruptures can nucleate only if the VW region is larger than the
31
nucleation size
h
. For 2-D problems, two theoretical estimates of the nucle-
32
ation size in mode III are (Rice and Ruina, 1983; Rubin and Ampuero, 2005):
33
34
h
RR
=
π
4
μD
RS
(
b
a
)(
σ
p
)
;
h
RA
=
2
π
μD
RS
b
(
b
a
)
2
(
σ
p
)
,
(5)
where
μ
is the shear modulus.
h
RR
is the estimate of spatial scale that can
35
host accelerating slip (Rice and Ruina, 1983).
h
RA
is the upper bound, for
36
3
the parameter regime 1
/
2
< b/a <
1, on the size of the quasi-statically
37
expanding nucleation zone before it transitions to dynamic rupture (Rubin
38
and Ampuero, 2005).
39
In our models, the VW region is roughly 44 times the estimated nucle-
40
ation size
h
RA
, ensuring that the regions would be seismogenic. We define
41
the nucleation region for our simulated ruptures as the region between ex-
42
panding stress fronts at the initiation of dynamic rupture, which is generally
43
comparable to the
h
RA
for our models. Simulated earthquakes in our models
44
nucleate preferably in regions of higher prestress at the creeping-to-locked
45
boundary between the VW and VS regions where the stressing rate for the
46
locked VW region is highest, as observed in previous modeling studies (e.g.
47
Lambert et al., 2021a, and references therein).
48
Such rate-and-state laws are the most established representation of fault
49
shear resistance for relatively low slip rates from tectonic plate rates of 10
9
50
m/s to around 10
3
m/s. At higher slip rates, further enhanced dynamic
51
weakening of fault frictional resistance can occur, for example, due to rapid
52
shear heating, as widely documented in high-velocity laboratory experiments
53
(Tullis, 2007; Wibberley et al., 2008; Di Toro et al., 2011; Acosta et al., 2018)
54
and supported by theoretical studies (Rice, 2006). One particular mechanism
55
4
for such enhanced weakening is the thermal pressurization of pore fluids,
56
which occurs when pore fluids within the fault shearing layer heat up, expand,
57
and pressurize during dynamic rupture, reducing the effective normal stress,
58
and therefore shear resistance (Sibson, 1973; Rice, 2006; Noda and Lapusta,
59
2010).
60
Our simulations also incorporate enhanced dynamic weakening due to
61
the thermal pressurization of pore fluids, which is governed by the following
62
coupled differential equations for pore pressure and temperature evolution
63
(Figure S1; Noda and Lapusta, 2010):
64
65
∂p
(
y,z
;
t
)
∂t
=
α
hy
2
p
(
y,z
;
t
)
∂y
2
+ Λ
∂T
(
y,z
;
t
)
∂t
,
(6)
66
∂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
,
(7)
where
T
is the temperature of the pore fluid,
α
hy
is the hydraulic diffusivity,
67
α
th
is the thermal diffusivity,
τV
is the source of shear heating distributed
68
over a Gaussian shear layer of half-width
w
,
ρc
is the specific heat,
y
is the
69
distance normal to the fault plane, and Λ is the coupling coefficient that
70
gives pore pressure change per unit temperature change under undrained
71
conditions.
72
5
We also approximate the effects of off-fault yielding during dynamic rup-
73
ture in our models by employing a velocity limit of
V
max
= 15 m/s, as dis-
74
cussed in detail in Lambert et al. (2021b). This approximation is motivated
75
by detailed dynamic rupture simulations with off-fault yielding (Andrews,
76
2004), with the value of the velocity limit corresponding to a representative
77
seismogenic depth of 10 km.
78
2. Measures of average fault shear stress
79
Here, we review different methods for averaging the shear stress along a
80
fault over space and time. We denote methods for averaging shear stress over
81
space with supercripts (e.g.
τ
A
), whereas quantities that represent temporal
82
averages or require information about the full time-dependent slip process
83
are denoted by a bar (e.g.
τ
). The formulations for average stress measures
84
are presented in this work for a 1-D fault; the generalization of the average
85
stress measures to 2-D interfaces in 3-D media are presented in Noda and
86
Lapusta (2012).
87
88
In calculating average stress measures for individual dynamic ruptures,
89
we define the rupture domain as the collection of points over the entire fault
90
6
domain Ω where the local slip velocity exceeds the seismic velocity threshold:
91
92
rupt
=
{
z
|
V
(
z
)
V
thresh
}
.
(8)
One may also compute the averages over the entire VW domain Ω
vw
with
93
velocity-weakening frictional properties:
94
95
vw
=
{
z
|
(
a
(
z
)
b
(
z
)
<
0
}
=
{|
z
|≤
12km
}
.
(9)
For our 2-D models with 1-D faults, the spatial averaging is performed over
96
the fault length, such as the rupture length
λ
rupt
and length of the velocity-
97
weakening region
λ
vw
:
98
99
λ
rupt
=
Z
rupt
dz
;
λ
vw
=
Z
vw
dz
= 24 km
(10)
For model-spanning ruptures in our earthquake sequence simulations, the
100
rupture length
λ
rupt
is comparable to the length of the VW region
λ
vw
. Thus
101
average shear stress measures taken over the scale of the model-spanning
102
earthquakes are expected to be comparable to averages over the total VW
103
domain.
104
7
2.1. Spatially-averaged shear stress
105
The simplest definition of average shear stress along an interface is the
106
spatially-averaged shear stress
τ
A
(
t
), which is just the spatial average of shear
107
stress
τ
(
z,t
) acting in the overall slip direction. For example, the spatially-
108
averaged shear stress over the entire VW region in our models at time
t
is
109
given by:
110
111
τ
A
vw
(
t
) =
R
vw
τ
(
z,t
)
dz
λ
vw
.
(11)
For a given rupture, the spatially-averaged prestress
τ
A
ini
and final stress
τ
A
fin
112
are given by:
113
114
τ
A
ini
τ
A
rupt
(
t
ini
) =
R
rupt
τ
(
z,t
ini
)
dz
λ
rupt
and
τ
A
fin
τ
A
rupt
(
t
fin
) =
R
rupt
τ
(
z,t
fin
)
dz
λ
rupt
,
(12)
The spatially-averaged rupture stress
τ
A
rupt
, averaged over both the rupture
115
area and duration, is given by:
116
117
τ
A
rupt
=
R
t
fin
t
ini
R
rupt
τ
(
z,t
)
dz
λ
rupt
T
rupt
,
(13)
8
where
T
rupt
=
t
fin
t
ini
is the rupture duration.
118
119
While simple to define, the spatially-averaged shear stress is not generally
120
the most useful definition. It is impossible to directly compute on natural
121
faults, as generally we cannot infer the spatial distribution of absolute stress
122
from remote observations. But this measure is also difficult to put constraints
123
on or interpret, as it does not directly characterize the dissipation or describe
124
the local evolution of shear resistance (or frictional behavior) during rupture
125
(Noda and Lapusta, 2012).
126
2.2. Dissipation-based average shear stress
127
A more physically meaningful definition of average fault shear stress is
128
the dissipation-based average shear stress
τ
D
(
t
), which is the spatial average
129
of shear stress
τ
(
z,t
) acting in the overall slip direction weighted by the slip
130
rate distribution. The dissipation-based average shear stress at time
t
over
131
the entire VW region and a given rupture area Ω
rupt
is given by:
132
133
τ
D
vw
(
t
) =
R
vw
τ
(
z,t
)
V
(
z,t
)
dz
R
vw
V
(
z,t
)
dz
and
τ
D
rupt
(
t
) =
R
rupt
τ
(
z,t
)
V
(
z,t
)
dz
R
rupt
V
(
z,t
)
dz
,
(14)
respectively.
134
9
τ
D
is equal to the increment in dissipation per increment in potency. The
135
potency per unit length
P
over a given domain Ω in our 2-D simulations with
136
1-D faults is related to the average slip
δ
in the principal slip direction and
137
fault length
λ
. The potency over the VW region Ω
vw
and a given ruptured
138
region Ω
rupt
are given by:
139
140
P
vw
=
Z
vw
Z
t
0
V
(
z,t
)
dt
dz
=
δ
vw
λ
vw
and
P
rupt
=
Z
rupt
Z
t
fin
t
ini
V
(
z,t
)
dt
dz
=
δ
rupt
λ
rupt
,
(15)
respectively. Here,
δ
vw
is the cumulative average slip across the VW region
141
up to time
t
and
δ
rupt
is the average final slip accrued in the rupture between
142
the initiation and termination of the rupture,
t
ini
and
t
fin
, respectively.
143
One can define the dissipation-based average rupture stress
τ
D
rupt
as the
144
average shear stress associated with dissipation throughout a rupture.
τ
D
rupt
145
is then related to the total dissipated energy per unit potency within the
146
rupture,
E
Diss
/P
rupt
, by integrating
τ
D
over the average slip of the rupture
147
process:
148
τ
D
rupt
=
E
Diss
P
rupt
=
R
t
fin
t
ini
R
rupt
τ
(
z,t
)
V
(
z,t
)
dzdt
R
t
fin
t
ini
R
rupt
V
(
z,t
)
dzdt
=
R
δ
rupt
0
τ
D
rupt
(
t
(
δ
))
δ
rupt
.
(16)
10
One can see that
τ
D
and
τ
A
are equivalent when the slip rate
V
is uniform
149
over the interface (
dV/dz
= 0), as may be the case in well-controlled labora-
150
tory experiments where loading displacements may be practically uniformly
151
applied over the sample. In contrast, for interfaces experiencing spatially-
152
varying slip rates, such as during ruptures or larger meter-scale laboratory
153
experiments which experience spatially variable loading and slip (McLaskey
154
and Lockner, 2014; McLaskey et al., 2015; Yamashita et al., 2015; Mclaskey
155
and Yamashita, 2017), the value of
τ
D
(
t
) is controlled by the local shear resis-
156
tance opposing regions of fastest motion, hence the regions dissipating most
157
of the work imposed on the system. Thus,
τ
D
is a more useful definition for
158
the
average shear resistance on the fault
, which can potentially be inferred
159
from heat measurements. However, similarly to
τ
A
, the dissipation-based
160
average shear stress
τ
D
is not directly related to the average local frictional
161
behavior on the interface (Noda and Lapusta, 2012).
162
2.3. Energy-based average stress
163
The energy-based average rupture stress
τ
E
rupt
is related to the average
164
shear stress associated with the work, or strain energy change ∆
W
, given
165
the potency of the rupture:
166
11
167
τ
E
rupt
=
W
P
rupt
=
1
2

τ
E
ini
+
τ
E
fin

.
(17)
Here,
τ
E
ini
and
τ
E
fin
refer to average shear stresses at the beginning and end of
168
rupture, weighted by the final slip:
169
170
τ
E
ini
=
R
rupt
τ
(
z,t
ini
)
δ
rupt
(
z
)
dz
R
rupt
δ
rupt
(
z
)
dz
;
τ
E
fin
=
R
rupt
τ
(
z,t
fin
)
δ
rupt
(
z
)
dz
R
rupt
δ
rupt
(
z
)
dz
.
(18)
The energy-based average shear stress can be related to the shear tractions
171
and slip along the interface by integrating the work rate along a virtual quasi-
172
static process connecting the beginning and termination of a rupture (Noda
173
and Lapusta, 2012). One can construct a version of this process that results
174
in an average curve with slip that aims to preserve the local shear resistance
175
behavior as well as serve as an energy-partitioning diagram. By considering
176
the average slip
δ
rupt
in the overall slip direction after a rupture, one can scale
177
the final slip throughout the rupture area
̃
δ
(
z
) =
δ
rupt
(
z
)
/
δ
rupt
to provide a
178
weighting based on the relative contribution of the local final slip at each
179
point to the average slip or potency. One can then introduce the integration
180
variable
δ
which represents the increment in average slip from 0 to
δ
rupt
.
181
If the slip and slip rate on the fault are always parallel to the overall slip
182
12
direction, as they are in our 2-D models with 1-D faults, then the energy-
183
based average stress can be written as a function of
δ
:
184
τ
E
(
δ
) =
R
rupt
τ
(
z,δ
̃
δ
(
z
))
δ
rupt
(
z
)
dz
R
rupt
δ
rupt
(
z
)
dz
.
(19)
τ
E
has been shown to capture the average local frictional behavior on the
185
interface (Noda and Lapusta, 2012). Its plot versus the accumulating slip
186
can represent the energy partitioning during ruptures, with the area under
187
the curve equal to the dissipated energy per unit fault area and the initial
188
and final points contributing to the trapezoid illustrating the strain energy
189
released (Figure 4 of main text).
190
Note that the energy-based average stress cannot be determined from the
191
distribution of shear stress and slip along the interface at individual instants
192
in time during dynamic rupture, but requires knowledge of the entire rupture
193
event after the dynamic process is over. Therefore all energy-based average
194
stress quantities require knowledge of the final slip, and hence are denoted
195
by a bar, as in
τ
E
(
δ
). Note also that the average of
τ
E
(
δ
) as
δ
varies from 0
196
to the average rupture slip
δ
rupt
is equal to
τ
D
rupt
as the construction of
τ
E
(
δ
)
197
preserves the dissipated energy.
198
.
199
13
Parameter
Symbol
Value
Loading slip rate
V
pl
10
9
m/s
Shear wave speed
c
s
3299 m/s
Shear modulus
μ
36 GPa
Thermal diffusivity
α
th
10
6
m
2
/s
Specific heat
ρc
2.7 MPa/K
Shear zone half-width
w
10 mm
Rate-and-state parameters
Reference slip velocity
V
10
6
m/s
Reference friction coefficient
f
0.6
Rate-and-state direct effect (VW)
a
0.010
Rate-and-state evolution effect (VW)
b
0.015
Rate-and-state direct effect (VS)
a
0.050
Rate-and-state evolution effect (VS)
b
0.003
Length scales for simulations
Fault length
λ
96 km
Frictional domain
λ
fr
72 km
Velocity-weakening region
λ
vw
24 km
Cell size
z
3.3 m
Models with thermal pressurization of pore fluids, TP 1-6 in Table 2
Quasi-static cohesive zone
Λ
0
84 m
2D Nucleation size (Rice & Ruina, 1983)
h
RR
226 m
2D Nucleation size (Rubin & Ampuero, 2005)
h
RA
550 m
Models with only rate-and-state friction, RS 1-2 in Table 2
Quasi-static cohesive zone
Λ
0
106 m
2D Nucleation size (Rice & Ruina, 1983)
h
RR
282 m
2D Nucleation size (Rubin & Ampuero, 2005)
h
RA
688 m
Table 1: Model parameters used in simulations unless otherwise specified.
14
Parameter
Symbol
TP 1
TP 2
Interseismic drained effective stress ̄
σ
= (
σ
p
int
)
25 MPa
25 MPa
Shear zone half-width
w
10 mm
10 mm
Characteristic slip
D
RS
1 mm
1 mm
Rate-and-state direct effect (VS)
a
0.05
0.05
TP coupling coefficient
Λ
0.1 MPa/K
0.34 MPa/K
Hydraulic diffusivity
α
hy
10
3
m
2
/s
10
3
m
2
/s
Parameter
Symbol
TP 3
TP 4
Interseismic drained effective stress ̄
σ
= (
σ
p
int
)
25 MPa
50 MPa
Shear zone half-width
w
10 mm
10 mm
Characteristic slip
D
RS
1 mm
2 mm
Rate-and-state direct effect (VS)
a
0.025
0.05
TP coupling coefficient
Λ
0.34 MPa/K 0.34 MPa/K
Hydraulic diffusivity
α
hy
10
4
m
2
/s
10
3
m
2
/s
Parameter
Symbol
TP 5
TP 6
Interseismic drained effective stress ̄
σ
= (
σ
p
int
)
25 MPa
100 MPa
Shear zone half-width
w
1 mm
10 mm
Characteristic slip
D
RS
1 mm
4 mm
Rate-and-state direct effect (VS)
a
0.025
0.05
TP coupling coefficient
Λ
0.34 MPa/K 0.34 MPa/K
Hydraulic diffusivity
α
hy
10
4
m
2
/s
10
5
m
2
/s
Parameter
Symbol
RS 1
RS 2
Interseismic drained effective stress ̄
σ
= (
σ
p
int
)
10 MPa
20 MPa
Characteristic slip
D
RS
0.5 mm
1 mm
Table 2: Parameters for fault models including the thermal pressurization of pore fluids
(TP 1-6) and models including only standard rate-and-state friction (RS 1-2).
15
Seismogenic
VW
a-b < 0
VS
a-b > 0
VS
a-b > 0
z
y
Slip rate, V
pl
= 10
-9
m/s is prescribed
x
Linear elastic infinite space treated
by spectral BIE method
A)
B)
z
y
x
Frictional heating
Fluid pressure,
p
Temperature,
T
Figure 1: Modeling of sequences of earthquakes and aseismic slip on a 1-D rate-and-state
fault in a 2-D elastic medium with (a) a velocity-weakening (VW) seismogenic region
surrounded by two velocity-strengthening (VS) sections and (b) enhanced dynamic weak-
ening due to the thermal pressurization of pore fluids. The evolution of temperature and
pore fluid pressure due to shear heating and off-fault diffusion is computed throughout
our simulations.
16
0.0
0.6
0.4
0.2
0
10
20
30
40
A)
50
Average local SSQS friction in VW region
Ruptures forcibly
arrested by
VS regions
More efficient
enhanced weakening
TP 3
TP 2
TP 1
RS 2
RS 1
Spa.-averaged prestress
effective normal stress
B)
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
50
Figure 2: Scaling of measures of average rupture prestress with rupture length. The
(A) spatially-averaged prestress
τ
A
ini
and (B) energy-averaged prestress
τ
E
ini
are generally
comparable and decrease with increasing rupture size for ruptures arresting within the
VW region. For ruptures that are forcibly arrested in the VS region, the average prestress
depends on the VS properties and can increase as the average effective weakening decreases
(Perry et al., 2020)
17