of 55
manuscript submitted to
JGR: Solid Earth
Scale dependence of earthquake rupture prestress in
1
models with enhanced weakening: implications for
2
event statistics and inferences of fault stress
3
Val`ere Lambert
1
Nadia Lapusta
1
,
2
Daniel Faulkner
3
4
1
Seismological Laboratory, California Institute of Technology, Pasadena, California, USA
5
2
Department of Mechanical and Civil Engineering, California Institute of Technology, Pasadena,
6
California, USA
7
3
School of Environmental Sciences, University of Liverpool, Liverpool, UK
8
Key Points:
9
Local shear prestress varies significantly within and among ruptures, being close
10
to the quasi-static fault strength in nucleation regions.
11
Efficient weakening allows rupture propagation over areas of lower prestress, lead-
12
ing to lower average prestress over larger rupture areas.
13
Fault models with more efficient dynamic weakening produce fewer smaller events
14
and result in systematically lower b-values.
15
Corresponding author: Val`ere Lambert,
vlambert@caltech.edu
–1–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
Abstract
16
Determining conditions for earthquake slip on faults is a key goal of fault mechanics highly
17
relevant to seismic hazard. Previous studies have demonstrated that enhanced dynamic
18
weakening (EDW) can lead to dynamic rupture of faults with much lower shear stress
19
than required for rupture nucleation. We study the stress conditions before earthquake
20
ruptures of different sizes that spontaneously evolve in numerical simulations of earth-
21
quake sequences on rate-and-state faults with EDW due to thermal pressurization of pore
22
fluids. We find that average shear stress right before dynamic rupture (aka shear pre-
23
stress) systematically varies with the rupture size. The smallest ruptures have prestress
24
comparable to the local shear stress required for nucleation. Larger ruptures weaken the
25
fault more, propagate over increasingly under-stressed areas due to dynamic stress con-
26
centration, and result in progressively lower average prestress over the entire rupture.
27
The effect is more significant in fault models with more efficient EDW. We find that, as
28
a result, fault models with more efficient weakening produce fewer small events and re-
29
sult in systematically lower b-values of the frequency-magnitude event distributions. The
30
findings 1) illustrate that large earthquakes can occur on faults that appear not to be
31
critically stressed compared to stresses required for slip nucleation; 2) highlight the im-
32
portance of finite-fault modeling in relating the local friction behavior determined in the
33
lab to the field scale; and 3) suggest that paucity of small events or seismic quiescence
34
may be the observational indication of mature faults that operate under low shear stress
35
due to EDW.
36
1 Introduction
37
Determining the absolute level and controlling factors of the stress state on faults
38
has profound implications for earthquake physics, seismic hazard assessment, and the
39
role of faulting in plate tectonics and geodynamics. Numerous lines of field evidence sug-
40
gest that the average shear stress acting on mature faults must be low, 20 MPa or less,
41
in comparison to the expected shear resistance of 100 - 200 MPa averaged over the seis-
42
–2–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
mogenic depth, given rock overburden and hydrostatic pore fluid pressure, along with
43
typical quasi-static friction coefficients of 0.6 - 0.85 (aka ”Byerlee friction”) measured
44
in laboratory experiments (Brune et al., 1969; Henyey & Wasserburg, 1971; Sibson, 1975;
45
Byerlee, 1978; Lachenbruch & Sass, 1980; Townend & Zoback, 2004; Rice, 2006; Suppe,
46
2007; Tanikawa & Shimamoto, 2009; Nankali, 2011; Fulton et al., 2013; Gao & Wang,
47
2014). Such evidence includes the lack of a substantial heat flow anomaly around ma-
48
ture faults that would be expected for fault slip at 100 MPa or more (Brune et al., 1969;
49
Henyey & Wasserburg, 1971; Lachenbruch & Sass, 1980; Nankali, 2011; Gao & Wang,
50
2014), inferences of steep angles between the principal stress direction and fault plane
51
(Townend & Zoback, 2004), analyses of the fault core obtained by drilling through shal-
52
low parts of faults that have experienced major recent events, including the 2011
M
w
53
9.0 Tohoku-oki event (Tanikawa & Shimamoto, 2009; Fulton et al., 2013), the geome-
54
try of thrust-belt wedges (Suppe, 2007), and the existence of long-lived narrow shear zones
55
that do not exhibit any evidence of melting (Sibson, 1975; Rice, 2006). Note that such
56
evidence for apparent fault weakness pertains predominantly to mature faults, whereas
57
some studies suggest that smaller, less mature faults may sustain the expected high shear
58
stresses given Byerlee friction values and overburden minus hydrostatic pore pressure (e.g.
59
Townend & Zoback, 2000).
60
A relatively straightforward explanation for the low-stress operation of mature faults
61
is that they may be persistently weak (Figure 1), due to the presence of anomalously low
62
quasi-static friction coefficients and/or low effective normal stress from pervasive fluid
63
overpressure (Brown et al., 2003; Faulkner et al., 2006; Bangs et al., 2009; Collettini et
64
al., 2009; Lockner et al., 2011). However, most materials with low quasi-static friction
65
coefficients (less than 0.5) under laboratory conditions tend to exhibit velocity-strengthening
66
behavior (Ikari et al., 2011), which would preclude spontaneous nucleation of dynamic
67
ruptures. Moreover, while evidence of substantial fluid overpressure has been documented
68
for many subduction zones (Brown et al., 2003; Bangs et al., 2009), there remains much
69
debate over the ubiquity of chronic near-lithostatic fluid overpressurization along faults
70
–3–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
in other tectonic settings, such as continental faults, with some borehole measurements
71
suggesting fluid pressure levels more consistent with hydrostatic conditions (Townend
72
& Zoback, 2000; Zoback et al., 2010).
73
An alternative hypothesis for explaining such low-stress, low-heat operation is that
74
mature faults are indeed strong at slow, quasi-static sliding rates but undergo consid-
75
erable enhanced dynamic weakening at seismic slip rates, which has been widely hypoth-
76
esized in theoretical studies and documented in laboratory experiments (Figure 1, dashed
77
black line; Sibson, 1973; Tsutsumi & Shimamoto, 1997; Rice, 2006; Wibberley et al., 2008;
78
Di Toro et al., 2011; Noda et al., 2009; Acosta et al., 2018). The presence of enhanced
79
dynamic weakening on natural faults has been questioned by the expectation that en-
80
hanced dynamic weakening would produce much larger static stress drops than typical
81
values of 1 to 10 MPa inferred from earthquakes on natural faults (Allmann & Shearer,
82
2009; Ye et al., 2016b). The expectation is based on a common assumption that the shear
83
prestress over the entire rupture area should be near the static strength of the fault while
84
the final shear stress should be near the dynamic strength of the fault, resulting in a large
85
static stress change. However, a number of numerical and laboratory studies have demon-
86
strated that, once nucleated, dynamic ruptures can propagate under regions with pre-
87
stress conditions that are well below the expected static strength, based on prescribed
88
or measured quasi-static friction coefficients and confining conditions (Zheng & Rice, 1998;
89
Noda et al., 2009; Lu et al., 2010; Dunham et al., 2011; Gabriel et al., 2012; Fineberg
90
& Bouchbinder, 2015) while the final shear stress could be higher then dynamic shear
91
stress for pulse-like ruptures, with both inferences promoting reasonable stress drops. Such
92
studies have often considered a single dynamic rupture nucleated artificially and prop-
93
agating over uniform prestress conditions.
94
Recent numerical studies of earthquake sequences have shown that fault models
95
with a combination of both hypotheses for low-stress operation, including some chronic
96
fluid overpressure as well as mild-to-moderate enhanced dynamic weakening due to the
97
–4–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
thermal pressurization of pore fluids, work well for reproducing a range of observations
98
(Perry et al., 2020; Lambert et al., in press). These include reasonable static stress drops
99
between 1 - 10 MPa nearly independent of earthquake magnitude, the seismologically
100
inferred increase in average breakdown energy with rupture size, the radiation ratios be-
101
tween 0.1 and 1 inferred for natural events, and the heat flow constraints. The simula-
102
tions produce mainly crack-like or mild pulse-like ruptures, with no significant under-
103
shoot. The near magnitude-invariance of average static stress drop arises in these fault
104
models because enhanced dynamic weakening results in both lower average prestress and
105
lower average final shear stress for larger ruptures with larger slip, with the average static
106
stress drops being nearly magnitude-independent. These studies suggest that distinguish-
107
ing between the conditions required for rupture nucleation and propagation is important
108
for assessing the relationship between laboratory friction measurements, seismological
109
observations and the absolute stress conditions on faults.
110
Here, we use and expand upon the set of numerical models from Perry et al. (2020)
111
and Lambert et al. (in press) to document the variability of prestress on a fault that arises
112
from the history of previous ruptures, and to study the relation between the size of dy-
113
namic rupture events and the average shear prestress over the rupture area. We also ex-
114
amine how the complexity of earthquake sequences, in terms of the variability of rup-
115
ture size, differs with the efficiency of dynamic weakening. We study these behaviors in
116
the context of simulations of sequences of earthquakes and slow slip, which allow the pre-
117
stress conditions before earthquakes to be set by the loading conditions, evolving fault
118
shear resistance (including weakening and healing), and stress redistribution by prior slip,
119
as would occur on natural faults. Moreover, our simulations resolve the spontaneous nu-
120
cleation process with the natural acceleration of slow unsteady slip prior to dynamic rup-
121
ture. The constitutive relations for the evolving fault resistance and healing adopted in
122
our models have been formulated as a result of a large body of laboratory, field and the-
123
oretical work (e.g Sibson, 1973; Dieterich, 1979; Ruina, 1983; Rice, 2006; Wibberley et
124
al., 2008; Di Toro et al., 2011). Indeed, laboratory experiments of fault shear resistance
125
–5–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
at both slow and fast slip rates have been indispensible for our understanding of fault
126
behavior and for formulating fault models such as the ones used in this study. The mod-
127
elling allows us to examine the implications of the laboratory-derived constitutive be-
128
haviors for the larger-scale behavior of faults, and we compare our inferences of average
129
shear prestress from relatively large-scale finite-fault modeling to field measurements of
130
crustal stresses acting on mature faults and small-scale laboratory measurements of the
131
shear resistance of typical fault materials.
132
2 Building on laboratory constraints to model larger-scale fault be-
133
havior
134
Laboratory experiments have been instrumental for exploring aspects of fault re-
135
sistance during both slow and fast sliding (10
9
m/s - 1 m/s, Figure 1). Experiments
136
with slow sliding velocities (
<
10
3
m/s) are critical for formulating fault constitutive
137
laws that form the basis for understanding the nucleation of earthquake ruptures. High-
138
velocity laboratory friction experiments have demonstrated enhanced dynamic weaken-
139
ing of faults and elucidated a range of mechanisms by which this dynamic weakening can
140
occur (e.g. Han et al., 2007; Wibberley et al., 2008; Goldsby & Tullis, 2011; Di Toro et
141
al., 2011; Faulkner et al., 2011; De Paola et al., 2015; Acosta et al., 2018). Most slow-
142
and high-velocity experiments measure or infer the relevant quantities - slip, slip rate,
143
shear stress etc - averaged over the sample and examine the evolution of shear resistance
144
corresponding to a particular history of loading, such as imposed variations in the dis-
145
placement rate of the loading piston, and the particular fault conditions (normal stress,
146
temperature, pore fluid pressure, etc.). Some experimental studies imposed the expected
147
sliding motion during earthquakes in order to directly relate laboratory stress measure-
148
ments to seismological quantities, such as static stress drop and breakdown energy (e.g.
149
Sone & Shimamoto, 2009; Fukuyama & Mizoguchi, 2010; Nielsen et al., 2016).
150
–6–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
To understand the full implications of the evolution of shear resistance measured
151
in small-scale experiments for slip at larger scales along natural faults, they are synthe-
152
sized into mathematical formulations and used in numerical modeling, for the following
153
reasons. During slipping events on a finite fault over scales of tens of meters to kilome-
154
tres - much larger than the experimental scale - the fault does not slip uniformly with
155
a predetermined slip-rate history. Rather, the slip event initiates on a portion of the fault
156
and then spreads along the fault, with varying slip-rate histories and final slips at dif-
157
ferent points along the fault. This is captured in inversions of large earthquakes (e.g. Heaton,
158
1990; Simons et al., 2011; Ye et al., 2016a; Tinti et al., 2016) and, to a degree, in larger-
159
scale experiments, sometimes involving analog materials (Lu et al., 2010; McLaskey et
160
al., 2014; Svetlizky & Fineberg, 2014; Yamashita et al., 2015; Rubino et al., 2017). In
161
the process, the slip (1) transfers stress to the more locked portions of the fault and (2)
162
enters portions of the fault with different conditions - such as levels of shear pre-stress,
163
pore fluid pressure, etc - and potentially different friction and hydraulic properties. Hence
164
the resulting coupled evolution of shear resistance and slip rate at different locations on
165
the fault is often quite different and, through stress transfer, strongly dependent on the
166
entire slip process at all locations throughout the rupture. These nonlinear and often dy-
167
namic feedback processes on the scales of tens of meters to kilometers can currently be
168
only captured through numerical modeling.
169
Many numerical models of earthquake source processes utilize insight from labo-
170
ratory experiments that indicate that the resistance to shear
τ
along a fault depends on
171
the sliding rate
V
and the quality and/or lifetime of the local contacts, typically param-
172
eterized by a state variable
θ
with units of time, as well as on the effective normal stress
173
σ
=
σ
p
acting on the fault, with
σ
being the normal stress and
p
being the pore fluid
174
pressure localized within the shearing layer (e.g. Dieterich, 1979; Marone, 1998). For con-
175
tinuum problems involving frictional sliding, the motion within the continuum is gov-
176
erned by the balance of linear momentum, subject to the boundary condition that trac-
177
tions are given by the constitutive law of the interface. For frictional sliding without changes
178
–7–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
in the elastodynamic normal stress, which is the case considered in this work, the bound-
179
ary condition reduces to the shear stress being equal to the shear resistance on the in-
180
terface (
y
= 0):
181
τ
stress
(
x,y
= 0
,z
;
t
) =
τ
resistance
(
x,y
= 0
,z
;
t
)
=
f
(
V,θ
)(
σ
p
)
.
(1)
An important concept in the rate-and-state formulation of the friction coefficient
f
(
V,θ
)
182
is that the friction coefficient is not a fixed property of the interface but evolves over time,
183
facilitating the time-dependent changes of shear resistance and hence shear stress along
184
the fault during shear.
185
The most commonly used formulation of rate-and-state laws is the Dieterich-Ruina
186
formulation (Dieterich, 1979; Ruina, 1983):
187
f
(
V,θ
) =
[
f
+
a
ln
V
V
+
b
ln
θV
L
]
,
(2)
where
f
is a reference steady-state friction coefficient at the reference sliding rate
V
,
188
L
is the characteristic slip distance, and
a
and
b
are the direct effect and evolution ef-
189
fect parameters, respectively. Our fault models are governed by a form of the laboratory-
190
derived Dieterich-Ruina rate-and-state friction law regularized for zero and negative slip
191
rates (Lapusta et al., 2000; Noda & Lapusta, 2010). The evolution of the state variable
192
can be described by various evolution laws; we employ the aging law (Ruina, 1983):
193
̇
θ
= 1
V θ
L
,
(3)
–8–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
which describes evolution during sliding as well as time-dependent healing in near-stationary
194
contact. In our models, the shear resistance and shear stress also change due to the evo-
195
lution of pore fluid pressure
p
.
196
We conduct numerical simulations following the methodological developments of
197
Lapusta et al. (2000), Noda and Lapusta (2010) and Lambert et al. (in press) in order
198
to solve the elastodynamic equations of motion with the fault boundary conditions, in-
199
cluding the evolution of pore fluid pressure and temperature on the fault coupled with
200
off-fault diffusion. The simulations solve for mode III slip on a 1-D fault embedded into
201
a 2-D uniform, isotropic, elastic medium (Figure 2). The potential types of slip on the
202
fault include sequences of earthquakes and aseismic slip (SEAS) and they are simulated
203
in their entirety, including the nucleation process, dynamic rupture propagation, post-
204
seismic slip that follows the event, and interseismic period between seismic events that
205
can last up to tens or hundreds of years and host steady and transient slow slip (Fig-
206
ure 2).
207
The simulated fault in our models contains a 24-km-long segment with velocity-
208
weakening (VW) frictional properties where earthquake ruptures may nucleate and prop-
209
agate, surrounded by velocity-strengthening (VS) segments that inhibit rupture nucle-
210
ation and propagation. Our simulations include enhanced dynamic weakening due to the
211
thermal pressurization of pore fluids, which occurs when pore fluids within the fault shear-
212
ing layer heat up and pressurize during dynamic rupture, reducing the effective normal
213
stress and shear resistance (Sibson, 1973; Rice, 2006; Noda & Lapusta, 2010). Thermal
214
pressurization is one potential mechanism for enhanced weakening; qualitatively simi-
215
lar results should hold for models with other types of enhanced dynamic weakening. We
216
follow the thermal pressurization formulation of Noda and Lapusta (2010) (Supplemen-
217
tary Materials).
218
For the purpose of comparing local frictional behavior with the average prestress
219
for dynamic ruptures of varying sizes, we focus this study on simulated ruptures that ar-
220
–9–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
rest within the VW region, where the friction properties are uniform with a quasi-static
221
reference friction of 0.6, consistent with many materials exhibiting VW behavior in lab-
222
oratory experiments (Ikari et al., 2011). We examine the evolution of the apparent fric-
223
tion coefficient, or the ratio of the current shear stress
τ
to the interseismic drained ef-
224
fective normal stress (
σ
p
int
), where
p
int
is the interseismic drained value of the pore
225
pressure. The ”drained” refers to the effective stress with ambient pore pressure unaf-
226
fected by slip processes such as dilatancy, compaction, or thermal pressurization.
227
We examine fault models with varying levels of ambient fluid overpressure in terms
228
of the effective normal stress, as well as varying degrees of efficiency in enhanced weak-
229
ening due to thermal pressurization. The parameter values we have chosen (Tables 1-
230
3) are motivated by prior studies that have reproduced a range of seismological obser-
231
vations as well as low-stress, low-heat operation of mature faults (Perry et al., 2020; Lam-
232
bert et al., in press). The parameter values also facilitate our goal of examining ruptures
233
in fault models with a range of efficiency in enhanced dynamic weakening. We define the
234
beginning and end of dynamic rupture,
t
ini
and
t
fin
respectively, as well as the ruptured
235
area Ω, using a slip velocity threshold (
V
thresh
= 0
.
01 m/s) for seismic slip, based on
236
previous studies (Perry et al., 2020; Lambert et al., in press). Note that
t
ini
and
t
fin
re-
237
fer to the beginning and end of the entire rupture event, which starts when one location
238
on the fault reaches the threshold velocity and ends when all points on the fault drop
239
below the threshold velocity. In the following, we use ”rupture” to refer to such dynamic
240
slip events, unless noted otherwise. Further description of the numerical methodology
241
can be found in the Supplementary Materials.
242
3 Evolution of local slip and shear resistance and notions of failure
243
Our simulations capture the evolution of motion and shear stress across the fault
244
over sequences of earthquakes spanning several thousands of years (Figure 2C). The ini-
245
tial distributions of shear stress and other quantities such as the slip rate are assumed
246
–10–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
to be uniform along most of the VW region of the fault at the start of our simulations,
247
other than a small region of initially high prestress near the VW-VS boundary to nu-
248
cleate the first rupture in the earthquake sequence. The distributions of shear stress and
249
slip along the fault evolve to become highly variable throughout periods of fast earthquake-
250
producing slip as well as slow aseismic slip and fault locking. Below we review how the
251
rate-and-state friction framework allows the model to represent both creeping, locked,
252
and seismically slipping fault areas as well as transitions between these different styles
253
of slip.
254
During dynamic rupture, the evolution of slip rate and shear stress can be partic-
255
ularly complex and variable along the fault. At points where individual ruptures nucle-
256
ate, the slip rate gradually accelerates towards seismic slip rates and shear stress at the
257
beginning of rupture,
t
ini
, is relatively high, with the apparent friction coefficient
τ/
(
σ
258
p
int
) close to the quasi-static reference friction of 0.6. As seismic slip rates are reached,
259
τ/
(
σ
p
int
) drops substantially due to thermal pressurization of pore fluids in a man-
260
ner qualitatively consistent with the enhanced dynamic weakening observed in high-velocity
261
laboratory friction experiments (Figure 2H). The evolution of slip rate and shear stress
262
outside of the nucleation region is even more complicated: The shear stress at
t
ini
, prior
263
to the arrival of the rupture front, can be much lower than the shear stress levels where
264
the rupture nucleates, then increases to a higher peak shear stress that reflects the in-
265
terseismic fault healing and rate-and-state direct effect and is achieved due to the dy-
266
namic stress concentration at the rupture front, and then decreases due to weakening
267
with seismic slip (Figure 2H vs. I). Consistently, the slip rate rapidly increases to seis-
268
mic values at the beginning of slip and then decreases, as in a typical Yoffe-like behav-
269
ior for dynamic ruptures (Figure 2G; e.g Tinti et al., 2005). Thus, even with the uni-
270
form normal stress and uniform parameters of the assumed friction and pore pressure
271
equations within the seismogenic VW region, the prestress conditions throughout the
272
rupture area can be highly variable and, in part, substantially different between regions
273
of rupture nucleation and rupture propagation.
274
–11–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
Note that the peak shear stress during dynamic rupture of fault locations outside
275
the nucleation zone can correspond to much higher apparent friction coefficient (e.g., 0.95
276
in Figure 2I) than the reference friction coefficient (
f
= 0
.
6 in this study). This is due
277
to both the direct effect at the rupture tip and the high, interseismically ”healed” value
278
of the state variable
θ
, as discussed in Lambert and Lapusta (2020) and the Supplemen-
279
tary Materials (equation S3). As follows from the first line of equation S3, the difference
280
between the peak friction coefficient and
f
due to the direct effect of
a
ln(
V
peak
/V
) would
281
be 0.14 to 0.16 for
V
peak
= 1 to 10 m/s and other parameters of our model, with the
282
rest due to the much larger value of the ”healed” state variable than that for sliding at
283
the reference sliding rate.
284
The local evolution of shear stress throughout the VW seismogenic zone differs among
285
points based on the long-term history of motion, including both local slip as well as slip
286
across the entire fault. For example, a point at the center of the VW region (z = 0 km)
287
of one of our simulations (fault model TP 3 in Table 2, as shown in Figure 2C) experi-
288
ences substantial slip only during the largest earthquake ruptures that span the entire
289
VW domain, resulting in a relatively simple and quasi-repetitive pattern of stress accu-
290
mulation and weakening over sequences of earthquakes (Figure 3A & C). In contrast, an-
291
other point in the VW region closer to the VS boundary (z = -9.6 km) experiences dif-
292
ferent amounts of slip during dynamic ruptures of varying size, resulting in a more com-
293
plicated evolution of shear stress with accumulating slip (Figure 3B & D).
294
In between individual earthquakes, the VS regions of the fault creep (i.e., slowly
295
slip) with the slip rate close to the prescribed tectonic plate rate, due to that rate be-
296
ing imposed on the fault areas nearby, with occasional quasi-static accelerations due to
297
post-seismic slip (Figure 4, left column). The creep penetrates into the VW regions nearby,
298
creating fault areas prone to earthquake nucleation (Jiang & Lapusta, 2016; Michel et
299
al., 2017) (Figure 4, right column). These points of the VW region close to the VS re-
300
gion (within one or so nucleation length) are reloaded due to creep and post-seismic slip
301
–12–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
from previous rupture within the VS regions. The loading rate at these points near the
302
VS-VW boundary varies over time depending on the rate of motion in the VS region,
303
which in turn depends on the previous history of co-seismic slip during dynamic ruptures
304
in the VW region.
305
The slip rate and apparent friction at points close to the VW-VS boundary are typ-
306
ically brought to near steady conditions around the loading plate rate, however both ex-
307
hibit small oscillations as these points continue to be loaded by creep in the VS region,
308
resulting in further acceleration, slip and weakening, and thus the transmission of stress
309
further into the VW region until a sufficiently large area is loaded to sustain rupture nu-
310
cleation and acceleration to seismic slip rates (Figure 4E-G). This oscillatory behavior
311
is consistent with predictions from the stability analysis of a single degree-of-freedom spring-
312
slider undergoing frictional slip, where the amplitude of the oscillations is expected to
313
grow as the spring stiffness decreases below a critical stiffness value until (Gu et al., 1984).
314
The effective stiffness of the slipping fault zone in a continuum model is inversely pro-
315
portional to the slipping zone side (Rice & Ruina, 1983), decreasing with the increas-
316
ing slipping region. Note that this rate-and-state nucleation process has been used to
317
explain the period-dependent response of microseismicity to periodic stress perturbations
318
in Nepal, where seismicity shows significant variations in response to annual monsoon-
319
induced stress variations but not to semidiurnal tidal stresses of the same magnitude (Ader
320
et al., 2014).
321
In contrast, much of the VW region further away from the VS regions is essentially
322
locked, which is expressed in the rate-and-state formulation as sliding at very low, but
323
still non-zero, slip rates that are many orders of magnitude smaller than the loading rate
324
(Figure 5A-B). This differential motion between the VS and VW regions loads points
325
in the VW region (Figure 5C-D), gradually increasing shear stress there (e.g., between
326
700 and 800 years in Figure 5C). Note that the interseismic stressing rate is higher at
327
locations closer to the creeping regions than further away from it (Figures 5C vs. 5D vs.
328
–13–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
4F), as one would expect. At the same time, the essentially locked points within the VW
329
region experience time-dependent healing of the local shear resistance encapsulated in
330
the increase of the state variable
θ
(Figure 5E-F). One of the manifestations of this heal-
331
ing is that larger interseismic increases in the state variable generally lead to higher peak
332
shear stress during dynamic rupture propagation (Lambert & Lapusta, 2020). Despite
333
the increase in the state variable, its value is far below the steady-state one for the very
334
low interseismic slip rates, consistent with continuing healing prior to dynamic rupture
335
(Figure 5G-H). Depending on whether the local shear stressing rate (which increases the
336
shear stress
τ
on the left of equation 1) is larger or smaller than the rate of healing (ex-
337
pressed by the last,
θ
term on the right hand side of equation 2), the local slip rate (that
338
enters the second term of equation 2) increases (as between 700 and 800 years in Fig-
339
ure 5A) or decreases, i.e., the fault is accelerating towards failure or becomes even more
340
locked. However, most of the locked points of the fault never accelerate close to failure
341
interseismically; rather, they fail due to stress concentrations from dynamic events, seen
342
as vertical lines in Figure 5C-D.
343
We note that healing on natural faults, in the presence of fluids and depth-dependent
344
elevated temperatures, can be affected by a number of mechanisms that are not captured
345
by the basic state evolution equation (Yasuhara et al., 2005; Tenthorey & Cox, 2006; Chen
346
et al., 2015a, 2015b). Incorporating more realistic healing into shear resistance formu-
347
lations and numerical modelling is an important goal for future work. This can be done
348
by modifying the evolution of the state variable
θ
or adding other state variables that
349
would encode healing. Yet, qualitatively, additional healing mechanisms would have sim-
350
ilar effects on the simulations as the current rate-and-state healing, in that the healing
351
would modify the peak shear resistance and the subsequent evolution of the resistance
352
based on the interseismic fault state, potentially further amplifying differences in shear
353
resistance evolution for different points along the fault (e.g., nucleation points vs. locked
354
points) that our simulations already highlight.
355
–14–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
The presence of time-dependent healing as well as persistent, potentially unper-
356
ceivable, slow (quasi-static) motion and its acceleration under variable levels of shear stress
357
illustrate how the concepts of failure, and hence strength, are not easily defined for fric-
358
tional sliding. For realistic frictional interfaces, the precise value of a static friction co-
359
efficient is ill-defined, since no interface loaded in shear is perfectly static; rather creep
360
processes occur at slow, unperceivable slip rates at any level of shear loading (Dieterich
361
& Kilgore, 1994; Bhattacharya et al., 2017) and/or over parts of the contacting inter-
362
faces (Rubinstein et al., 2004, 2006; Ben-David et al., 2010). Hence the transition from
363
locked interfaces to detectable slip is always a gradual process (although it may be oc-
364
curring faster than the time scales of interest/observation in many applications). This
365
reality is reflected in lab-derived fault constitutive relations such as rate-and-state fric-
366
tion. Since failure typically refers to the presence of irreversible or inelastic deformation,
367
frictional interfaces may be considered failing under any style or rate of motion, be it dur-
368
ing slow steady sliding, transient slow slip, or dynamic rupture. Therefore, any mean-
369
ingful notion of strength first requires definition of the failure of interest, e.g., reaching
370
seismic slip rates of the order of 1 m/s. Without such explicit definition, failure is then
371
implicitly defined as transition from locked to slipping and corresponds to sliding with
372
a detectable velocity; for laboratory experiments or observational studies, this would im-
373
ply that whether the interface is locked or slipping depends on the instrumental preci-
374
sion for detectable motion.
375
In this study, we would like to compare the shear stress values required for aseis-
376
mic slip nucleation and for dynamic rupture propagation. During spontaneous aseismic
377
slip nucleation, the slip rates evolve from very low to seismic, passing in the process through
378
the slip rate equal to the tectonic loading rate
V
pl
. In the standard rate-and-state fric-
379
tion, at each fixed sliding rate
V
, the friction coefficient eventually evolves to a steady-
380
state value
f
ss
(
V
) (equation S2; for very small slip rates, the regularized formulation of
381
equation S5 needs to be considered). Under slow loading, aseismic earthquake nucleation
382
on a finite fault is typically a gradual process, with many points within the nucleation
383
–15–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
zone being close to the steady state (Figure 4; Rubin & Ampuero, 2005; Kaneko & La-
384
pusta, 2008). While the steady-state values of friction depend on the sliding rate, the
385
dependence is relatively minor at the low, quasi-static slip rates between the plate rate
386
of approximately 10
9
m/s and sub-seismic slip rates of
<
10
3
m/s (Figure 1) which
387
are relevant for fault creep and earthquake nucleation, and for which the standard rate-
388
and-state formulation is (approximately) valid. The product of this collection of steady-
389
state quasi-static friction coefficients and the interseismic drained effective stress gives
390
the shear resistance of faults at sustained slow sliding rates, which we call the
steady-
391
state quasi-static fault shear resistance
(referred to in short as local SSQS shear resis-
392
tance). As the representative value of such local SSQS shear resistance, we choose the
393
shear resistance of the fault steadily creeping at the prescribed long-term tectonic plate
394
rate
V
pl
(which the fault would have long-term if it were slipping stably), with the in-
395
terseismic drained value of the pore pressure
p
int
:
396
τ
V
pl
ss
(
z,t
) = (
σ
p
int
)
f
ss
(
V
pl
)
(4)
In our models,
τ
V
pl
ss
/
(
σ
p
int
) = 0
.
63 within the VW region. Note that choosing
V
in-
397
stead of
V
pl
would result in a similar value of
τ
V
pl
ss
/
(
σ
p
int
) =
f
= 0
.
6.
398
In the following section, we compare this representative value of local SSQS shear
399
resistance to the spatial distribution of shear stress prior to dynamic ruptures in our sim-
400
ulations. Note that the local SSQS shear resistance is similar to what is typically viewed
401
as ”frictional fault strength” in the sense of Byerlee (1978), i.e., this is the resistance that
402
needs to be met for noticeable quasi-static slip with the loading rate or another refer-
403
ence rate.
404
–16–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
4 Larger ruptures associated with lower shear prestress over the rup-
405
ture scale but higher prestress over smaller scales near nucleation
406
The interseismic periods in between individual earthquake ruptures in our simu-
407
lations vary from months to decades, depending on the size of the rupture and the stress
408
state resulting from the history of prior slip along the fault. Our earthquake sequence
409
simulations produce a wide variety of rupture sizes due to heterogeneous prestress con-
410
ditions along the fault that spontaneously arise in our models.
411
Let us consider the evolution of slip and shear stress in representative simulated
412
spontaneous ruptures of increasing sizes within the same simulation (Figure 6). Over se-
413
quences of rupture events, the shear stress conditions prior to and after individual dy-
414
namic ruptures become spatially heterogeneous. This stress heterogeneity is due in part
415
to the history of spatially variable slip and local static stress drop produced in previous
416
ruptures, as well as stress relaxation and redistribution due to aseismic slip. In addition,
417
while our simulated fault models are loaded by a constant long-term loading rate of
V
pl
,
418
the effective loading conditions along the fault interface vary in space and time due to
419
differences in slip rate along the fault. Ruptures nucleate preferentially in regions with
420
the highest shear prestress, which in our models occur near the creeping regions as dis-
421
cussed in section 3 (Figure 6). The ruptures then propagate into the less stressed areas
422
of the fault. Put another way, the average prestress over the nucleation region is higher
423
than the average prestress over the entire ruptured region (Figure 7A vs. B), as we quan-
424
tify in the following.
425
We compute the average shear prestress right before a dynamic rupture event over
426
the entire future rupture area (which we do as post-processing of data in our simulation).
427
We also compute the average shear prestress over the slow-slip nucleation zone, which
428
we call the
nucleation stress
. We compare these average shear stress measures with the
429
local steady-state quasi-static (SSQS) fault shear resistance
τ
V
pl
ss
, which is related to the
430
local fault constitutive properties during slow slip and given by equation 4.
431
–17–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
Averaging of spatially variable stress fields can be done in several different ways
432
(Noda & Lapusta, 2012; Noda et al., 2013). The simplest definition of the average shear
433
prestress over the rupture region Ω is the spatially averaged prestress
τ
A
ini
acting in the
434
overall slip direction at the beginning of the rupture
t
ini
, given by:
435
τ
A
ini
=
τ
(
z,t
ini
)
dz
dz
.
(5)
We can similarly define the spatially averaged nucleation stress
τ
A
nucl
within the nu-
436
cleation region. We define the nucleation region to be the fault segment between the ex-
437
panding stress fronts at the initiation of dynamic rupture; the size of the nucleation re-
438
gions in our simulations is comparable to the theoretical nucleation size estimate
h
RA
439
of Rubin and Ampuero (2005) (equation S6, Figure S1).
440
Not surprisingly and consistent with prior studies, we find that the spatially av-
441
eraged nucleation stress
τ
A
nucl
for our simulated ruptures is comparable to the local SSQS
442
shear resistance
τ
V
pl
ss
(Figure 7A). As a consequence, it does not significantly depend on
443
the ultimate rupture size or slip. Since the nucleation stress here is computed at the be-
444
ginning of dynamic rupture, it is then the shear stress within the nucleation zone at the
445
end of the nucleation, when parts of the zone slip with near-dynamic slip rates approach-
446
ing 10
2
m/s. That is why the nucleation stress is systematically slightly lower than the
447
local SSQS shear resistance defined as the steady-state shear resistance to slip with the
448
(lower) plate rate. The difference between the nucleation stress and local SSQS shear
449
resistance could be more substantial if dynamic weakening were efficient enough to af-
450
fect some portion of the earthquake nucleation region (Segall & Rice, 2006).
451
In contrast, the spatially averaged prestress over the entire ruptured area
τ
A
ini
tends
452
to decrease with the rupture size and increasingly deviate from the local SSQS shear re-
453
sistance and nucleation stress for increasingly efficient dynamic weakening (Figures 6 &
454
7B). Such behavior is also true for another average prestress measure, the energy-based
455
–18–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
average prestress
τ
E
ini
(Noda & Lapusta, 2012), which is the average shear prestress weighted
456
by the final slip of the rupture, and hence represents the average prestress associated with
457
the potency of the impending rupture:
458
τ
E
ini
=
τ
(
z,t
ini
)
δ
fin
(
z
)
dz
δ
fin
(
z
)
dz
(6)
where
δ
fin
(
z
) =
δ
(
z,t
fin
)
δ
(
z,t
ini
) is the final local slip accrued in the rupture. We
459
denote
τ
E
with a bar as it not only represents an average over space but also requires
460
knowledge of the final slip of the rupture.
τ
E
ini
differs from the spatially-averaged pre-
461
stress
τ
A
ini
over the rupture area when the resulting slip distribution is not uniform. We
462
find that
τ
E
ini
and
τ
A
ini
for our simulated ruptures are comparable and vary similarly with
463
the rupture size and efficiency of dynamic weakening, with the values of
τ
E
ini
being slightly
464
larger (Figure S2).
465
The finding that larger ruptures are associated with smaller average shear prestress
466
over the ruptured area may appear counterintuitive. Why do smaller ruptures not be-
467
come larger if they are more favorably prestressed? To understand this behavior, let us
468
consider the prestress averaged over several fixed scales around the nucleation region for
469
ruptures of different sizes. We locate the VW-VS boundary next to which each of our
470
simulated ruptures nucleate and average the prestress along the VW region over fixed
471
distances (1, 2, 4, 8, 12 and 16 km) from the corresponding VW-VS boundary (Figure
472
8; shown for fault model TP4 from Table 2). While the spatially-averaged prestress over
473
the entire rupture length decreases with increasing rupture size, we see that the prestress
474
spatially-averaged over smaller fixed scales is generally higher for larger ruptures than
475
for smaller ruptures (Figure 8 warmer vs cooler colored triangles). For smaller ruptures,
476
the average shear stress over scales just larger than their total rupture length is lower
477
than the average prestress of larger ruptures with comparable length to the fixed aver-
478
aging scales (Figure 8, triangles below the circles). This confirms that the smaller rup-
479
tures arrest because the prestress conditions ahead of the rupture are too low to sustain
480
–19–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
further rupture propagation. For larger ruptures, the average prestress levels at scales
481
smaller than their total rupture length are generally higher or comparable to the aver-
482
age prestress over smaller ruptures with the length comparable to the fixed averaging
483
scales (Figure 8, triangles above the circles). This finding suggests that larger ruptures
484
have higher, more favorable average prestress conditions at smaller scales compared to
485
smaller ruptures, which facilitates continued rupture propagation. Hence we find that
486
the shear prestress prior to our simulated ruptures of varying sizes self-organizes into a
487
spatial distribution of scale-dependent average shear stress that governs the rupture oc-
488
currence.
489
5 Role of dynamic stress transfers and motion-dependent local shear
490
resistance
491
Such scale- and motion-dependent average fault shear prestress before ruptures re-
492
sults from two related and interacting factors. First, as dynamic rupture propagates, some
493
of the released energy is carried by waves along the fault, creating a substantial stress
494
concentration near the rupture tip that is a well-known feature of dynamic rupture (e.g.,
495
Freund, 1990). The stress concentration enables rupture propagation over regions where
496
the prestress is lower than the local SSQS shear resistance, drawing the local shear stress
497
up to the peak stress before the subsequent stress drop due to local weakening (black
498
lines in Figure 6). The dynamic stress concentration increases with the rupture dimen-
499
sion and/or slip and thus allows larger ruptures to continue propagating over regions with
500
lower, and hence less favorable, prestress conditions (Figure 6). This is illustrated in this
501
work for largely crack-like ruptures that occur in the presented models with mild to mod-
502
erate enhanced dynamic weakening (Lambert et al., in press), but similar conclusions
503
would be reached for pulse-like ruptures provided that they satisfy the observational con-
504
straint of magnitude-independent stress drops, which implies that ruptures with larger
505
magnitudes would have larger average slip and hence larger stress concentrations. Note
506
that a pulse-like rupture with the same or similar spatial distribution of the slip rate (and
507
–20–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
hence the same local slip) propagating along the fault would result in a similar stress con-
508
centration at the rupture tip regardless of the rupture length; however, in that scenario,
509
pulses with larger rupture propagation lengths would have systematically lower static
510
stress drops, as the stress drops would be proportional to the (uniform) pulse slip divided
511
by ever increasing propagation lengths.
512
Second, the evolving local shear resistance substantially depends on both the prior
513
history of slip events on the fault through fault prestress and on the motion during the
514
current rupture event through dynamic stress transfers that add substantial time-dependent
515
loading. This pronounced dependence is due to strong coupling between the evolving mo-
516
tion, the resulting shear heating, and the evolving shear resistance. As a result, the evo-
517
lution of local slip rate and local shear resistance (1) significantly differs at different fault
518
locations of each rupture (despite uniform constitutive properties) and (2) significantly
519
differs at the same fault location for different ruptures (Figures 2D-I and 6D-E).
520
These two factors create a substantial positive feedback, in which larger ruptures
521
with more slip generate larger stress concentrations, leading to faster and larger slip, which
522
dynamically causes more fault weakening, which in turn promotes more/faster slip, more
523
energy release, larger stress concentrations, and increasing rupture sizes.
524
The result that larger ruptures are associated with lower average prestress indicates
525
the need for increasingly less favorable stress conditions to arrest growing ruptures. For
526
a given rupture size, if the prestress ahead of the rupture is favorable, then the rupture
527
would continue to grow until it experiences sufficiently unfavorable prestress conditions,
528
thus lowering the overall average prestress. Alternatively, the rupture may be forcibly
529
arrested by other means such as strong geometric or rheological barriers. For example,
530
ruptures propagating over higher prestress conditions within the VW region can be ar-
531
rested by fault regions with VS properties; in those cases, the overall average prestress
532
conditions would depend on the properties of the VS regions (Perry et al., 2020). De-
533
–21–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
tailed study of the implications of fault geometry and heterogeneity for rupture arrest
534
and the average stress conditions prior to rupture is an important topic for future work.
535
6 Comparison of finite-fault modeling to single-degree-of-freedom rep-
536
resentations
537
As captured in field observations of natural earthquakes and reflected in our sim-
538
ulations, sufficiently large earthquake ruptures nucleate on a subsection of the fault and
539
then propagate through other sections of the fault. Capturing such space-dependent be-
540
havior is typically called ”finite-fault” modeling, in contrast to the point source that con-
541
siders a spatially average representation of an event, as if it occurs at one ”point”. A typ-
542
ical numerical model of a point source is the single-degree-of-freedom system (SDOF)
543
of a slider with friction pulled by a spring (e.g. Dieterich, 1979; Ruina, 1983; Rice & Ru-
544
ina, 1983). Small-scale laboratory experiments often measure properties averaged over
545
a sample and are typically modeled as a SDOF spring-slider systems.
546
The significant role of spatially varying prestress conditions and dynamic stress trans-
547
fers during rupture propagation in determining the rupture behavior implies that cap-
548
turing the finite-fault nature of the process is essential for determining the stress evo-
549
lution characteristic of dynamic rupture. For example, several laboratory studies applied
550
variable slip rates histories inferred from natural earthquakes to rock samples, measured
551
the resulting shear resistance, and then related laboratory stress measurements to seis-
552
mological source properties such as breakdown energy and stress drops (e.g. Sone & Shi-
553
mamoto, 2009; Fukuyama & Mizoguchi, 2010; Nielsen et al., 2016). Such experiments
554
have provided invaluable data about the local shear resistance of faults, specifically en-
555
hanced dynamic weakening, that have informed theoretical and numerical modeling of
556
finite faults (e.g. Zheng & Rice, 1998; Rice, 2006; Noda et al., 2009; Noda & Lapusta,
557
2010; Dunham et al., 2011; Gabriel et al., 2012; Perry et al., 2020; Lambert et al., in press),
558
including the current study. However, the interpretation of such experiments needs to
559
–22–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
take into account their SDOF nature. For example, to improve alignment etc, the ex-
560
periments often impose pre-sliding at slow slip rates (of the order of micron/s) prior to
561
imitating seismic motion. That procedure results in the shear prestress before seismic
562
slip comparable to the local SSQS shear resistance (equation 4) and near steady-state
563
values of the state variable, as appropriate for a location within a nucleation zone. In
564
contrast, our simulations show that most points on a fault through which the rupture
565
propagates have much lower shear prestress and much larger values of the state variable
566
corresponding to well-healed fault (Figures 6 and 9B). Furthermore, the experiments of-
567
ten apply smoothened slip-rate histories obtained from finite-fault inversions, while the
568
stress concentration at the tip of dynamic rupture makes the slip rate variation much
569
more dramatic.
570
To illustrate the differences for the shear resistance evolution obtained with such
571
experimental procedures versus the one from our simulated finite-fault models, let us com-
572
pare the local fault behavior during one of our dynamic ruptures with a SDOF calcu-
573
lation. In the SDOF calculation, we use the same fault properties (equations 3, S4 and
574
S7-8) and same parameter values as in the finite-fault VW regions but apply quasi-static
575
presliding and modified, smoothened slip rates motivated by the laboratory procedures
576
of Fukuyama and Mizoguchi (2010) (further details in Supplementary Materials). We
577
conduct the comparison for two fault locations, one in the nucleation region and one within
578
dynamic rupture propagation region (Figure 9). These SDOF calculations are success-
579
ful at reproducing the presence of the enhanced dynamic weakening with slip as occurs
580
during dynamic ruptures and generally capture the more moderate slip evolution and
581
behavior of points within the nucleation region of our simulated dynamic ruptures. At
582
the same time, the overall shear stress evolution during typical propagation of the dy-
583
namic rupture substantially differs from that of the SDOF calculation, with notable fea-
584
tures including the low initial stress (which depends on prior slip history) relative to the
585
SSQS shear resistance, the much more dramatic increase in shear stress associated with
586
the dynamic rupture front (which arises due to the more healed fault coupled with the
587
–23–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
dynamic stress concentration), and the shear stress evolution at the end of slip (which
588
depends on the final slip distribution over the entire finite fault) (Figure 9).
589
7 Implications for earthquake statistics
590
A notable feature of the scale dependence of average prestress before dynamic rup-
591
ture is that, as an earthquake grows larger, the prestress needed for further propagation
592
decreases (Figure 7B). In addition, the higher the weakening rate, the easier it should
593
be for a rupture to have favorable prestress conditions to continue growing, rather than
594
arresting as a smaller earthquake. Hence one could hypothesize that the more efficient
595
the enhanced dynamic weakening, the smaller the complexity of the resulting earthquake
596
sequences, with increasing representation of larger events at the expense of smaller events.
597
This is exactly what our modeling shows (Figure 10). The fault models with in-
598
creasingly more efficient weakening produce earthquake sequences with increasingly fewer
599
small events and decreasing b-values of the cumulative size distribution (Figure 10). Fault
600
models with even more efficient dynamic weakening than considered in this study, such
601
as those that produce sharp self-healing pulses, result in relatively simple earthquake se-
602
quences consisting of only large events (Lambert et al., in press). The fault models gov-
603
erned by relatively mild to more moderate weakening as considered in this work develop
604
a wider range of earthquake sizes, due to a feedback loop of more likely rupture arrest
605
due to milder weakening creating stress heterogeneity that in turn makes rupture arrest
606
more likely. This result is consistent with those of previous quasi-dynamic earthquake
607
sequence simulations demonstrating complex earthquake sequences with b-values around
608
0.75 on faults with standard rate-and-state friction only and milder quasi-dynamic stress
609
transfer (Cattania, 2019). Our study shows that the b-values decrease to 0.5 for fully
610
dynamic simulations without enhanced dynamic weakening, and further decrease to 0.25
611
or so for the most efficient weakening considered in this study.
612
–24–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
While the frequency-magnitude distribution of seismicity over relatively large re-
613
gions, such as Northern or Southern California, is generally well-described by Gutenberg-
614
Richter scaling with typical b-values near unity (E. Field et al., 2013), whether such scal-
615
ing applies to individual fault segments and/or their immediate surroundings is a topic
616
of active research (Wesnousky, 1994; Ishibe & Shimazaki, 2012; Kagan et al., 2012; Page
617
& Felzer, 2015; Page & van der Elst, 2018; E. H. Field et al., 2017). Estimates of b-values
618
associated with individual fault segments can exhibit considerable variability (e.g. be-
619
tween 0.5 and 1.5 along faults in California; Tormann et al., 2014), and are sensitive to
620
a number of factors, including the magnitude of completeness of the relevant earthquake
621
catalog and the choice of observation region and time window (Tormann et al., 2014; Page
622
& Felzer, 2015; Ishibe & Shimazaki, 2012; Page & van der Elst, 2018). A number of stud-
623
ies suggest that the rate of large earthquakes on major faults, such as the San Andreas
624
Fault, is elevated above what would be expected given typical Gutenberg-Richter scal-
625
ing from smaller magnitude events (Schwartz & Coppersmith, 1984; E. H. Field et al.,
626
2017). In particular, some mature fault segments that have historically hosted large earth-
627
quakes, such as the Cholame and Carrizo segments of the San Andreas Fault, exhibit sub-
628
stantial deviations from typical Gutenberg-Richter scaling, being nearly absent of small
629
earthquakes (Sieh, 1978; Wesnousky, 1994; Bouchon & Karabulut, 2008; Hauksson et al.,
630
2012; Jiang & Lapusta, 2016; Michailos et al., 2019). Our findings suggest that the paucity
631
of microseismicity on such mature fault segments may indicate that they undergo sub-
632
stantial dynamic weakening during earthquakes ruptures.
633
8 Discussion
634
Our simulations demonstrate that the average shear prestress required for rupture
635
propagation can be considerably lower than the average shear stress required for the rup-
636
ture nucleation. This is because the quasi-static nucleation process is governed by rel-
637
atively small stress changes and hence requires favorable prestress conditions - close to
638
the local steady-state quasi-static shear resistance - to proceed. In contrast, during dy-
639
–25–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
namic rupture, the rupture front is driven by larger wave-mediated dynamic stress con-
640
centrations, which are more substantial for larger ruptures and facilitate rupture prop-
641
agation over less favorably stressed regions, resulting in the spatially-averaged prestress
642
over the ruptured area being much lower than the average local SSQS shear resistance.
643
More efficient weakening facilitates larger dynamic stress changes at the rupture front,
644
allowing propagation over even less favorable prestress conditions. Our results highlight
645
the significance of heterogeneity in prestress, or shear resistance, for the nucleation and
646
ultimate arrest of finite ruptures, even in fault models that have otherwise uniform ma-
647
terial and confining properties.
648
The decrease in averaged prestress with rupture length can be interpreted as a de-
649
crease in the average quasi-static friction coefficient
τ
A
ini
/
(
σ
p
int
) with rupture size (Fig-
650
ure 7). The average quasi-static friction coefficients for ruptures on the scale of the nu-
651
cleation size are consistent with the prescribed quasi-static reference friction coefficient
652
near typical Byerlee values. However, as we average the prestress over larger rupture lengths,
653
the average quasi-static friction coefficient can considerably decrease depending on the
654
efficiency in weakening.
655
The presence of enhanced dynamic weakening draws the average shear stress along
656
larger regions of the fault below the local SSQS consistent with earthquake nucleation,
657
resulting in lower average shear stress conditions in terms of both the average prestress
658
for larger ruptures and the average dynamic resistance associated with shear heating dur-
659
ing ruptures (Figure 11). The models presented in this study with mild-to-moderate en-
660
hanced weakening include considerable persistent fluid overpressurization to maintain
661
low-heat, low-stress conditions with average dynamic shear resistance during seismic slip
662
rates below 10 MPa; however the degree of fluid overpressure required to maintain low-
663
heat conditions is less than that with comparable rate-and-state properties but no en-
664
hanced weakening. The presence of some enhanced dynamic weakening is also needed
665
for persistently weak fault models due to chronic fluid overpressure in order to ensure
666
–26–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
that static stress drops are not too small, as they would otherwise be with low effective
667
stress and small changes in the friction coefficient due to standard rate-and-state laws
668
(Figures 11 and S3; Lambert et al., in press). Fault models with more efficient dynamic
669
weakening have been shown to be able to reproduce low-stress operation and reasonable
670
static stress drops with quasi-static friction coefficients around Byerlee values and higher
671
effective normal stress (e.g.
100 MPa; Noda et al., 2009; Dunham et al., 2011; Lam-
672
bert et al., in press). Earthquake sequence simulations of such fault models typically con-
673
sist of only large ruptures (Lambert et al., in press), consistent with the notion that large
674
fault areas governed by efficient weakening maintain substantially lower average shear
675
stresses than that required for nucleation. These findings further strengthen the conclu-
676
sion of prior studies that enhanced dynamic weakening can help explain the discrepancy
677
between laboratory values of (quasi-static) friction coefficients around 0.6 and geophys-
678
ical inferences of low effective coefficients of friction (
<
0.2), along with mild average static
679
stress drops of 1 to 10 MPa, over fault areas that host large earthquakes (e.g Marone,
680
1998; Suppe, 2007; Allmann & Shearer, 2009; Noda et al., 2009; Dunham et al., 2011;
681
Ikari et al., 2011; Gao & Wang, 2014; Ye et al., 2016b; Perry et al., 2020; Lambert et al.,
682
in press).
683
The scale dependence of average prestress before ruptures can also be interpreted
684
as a scale dependence of
average fault strength
, since the average prestress represents a
685
measure of how much shear stress that fault region can hold before failing in a rupture.
686
Given this interpretation, our simulations suggest that faults maintain lower average shear
687
stresses, and hence appear weaker, at larger scales than at smaller scales. This interpre-
688
tation is conceptually consistent with laboratory measurements of scale-dependent yield
689
stress for rocks and a number of engineering materials, which demonstrate decreasing
690
material strength with increasing scale (Jaeger & Cook, 1976; Bandis et al., 1981; Greer
691
et al., 2005; Pharr et al., 2010; Uchic et al., 2004; Yamashita et al., 2015; Thom et al.,
692
2017). Note that our larger simulated ruptures, even with more efficient weakening, still
693
require higher average shear stresses over smaller scales in order to nucleate and grow.
694
–27–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
Thus the lower average prestress levels that allow continued failure in dynamic ruptures
695
at larger scales only become relevant once the rupture event has already nucleated and
696
sufficiently grown over smaller scales. This consideration suggests that the critical stress
697
conditions for rupture occurrence are governed not by a single stress quantity but by a
698
distribution of scale-dependent stress criteria for rupture nucleation and continued prop-
699
agation. An important implication of our findings is that the critical stress for earthquake
700
occurrence may not be governed by a simple condition such as a certain level of Coulomb
701
stress. Given our findings, in order to reason about the stress conditions critical for a
702
rupture to occur, it is important to consider both the size of the rupture and the weak-
703
ening behavior, and hence the style of motion, that may occur throughout rupture prop-
704
agation.
705
The scale dependence of fault material strength has also been hypothesized to ex-
706
plain the measured scaling of roughness on natural fault surfaces (Brodsky et al., 2016).
707
Dynamic rupture simulations on geometrically irregular faults motivated by such rough-
708
ness measurements have indicated an additional contribution to fault shear resistance
709
arising from roughness drag during rupture propagation (Fang & Dunham, 2013). Fur-
710
ther examination of the scale dependence of average shear resistance across faults includ-
711
ing realistic fault geometry is an important topic for future work.
712
A common assumption is that the shear prestress over the entire ruptured area must
713
be near the local static (or quasi-static) strength, comparable to the SSQS shear resis-
714
tance discussed in this study. We demonstrate that the assumption is not necessarily valid
715
and that faults with enhanced dynamic weakening and history of large earthquake rup-
716
tures would, in fact, be expected to have low average shear stress over large enough scales.
717
At the same time, the state of stress needs to be heterogeneous, with the average stresses
718
over small scales (comparable to earthquake nucleation) being close to the (much higher)
719
local SSQS shear resistance in some places. Thus, while individual measurements of low
720
resolved shear stress onto a fault may suggest that those locations appear to not be crit-
721
–28–
ESSOAr | https://doi.org/10.1002/essoar.10506240.1 | Non-exclusive | First posted online: Tue, 16 Feb 2021 10:01:59 | This content has not been peer reviewed.