manuscript submitted to
JGR-Solid Earth
Resolving simulated sequences of earthquakes and fault
1
interactions: implications for physics-based seismic
2
hazard assessment
3
Val`ere Lambert
1
and Nadia Lapusta
1
,
2
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
Key Points:
8
•
Long-term interactions of fault segments qualitatively differ among fully versus
9
quasi-dynamic simulations, and when using oversized cells.
10
•
Reproducing frequency-magnitude distributions and static stress drops is not suf-
11
ficient to constrain the rate of multi-segment ruptures.
12
•
Simulated earthquake sequences can differ due to compounded effects of numer-
13
ical errors, even when individual ruptures are well-resolved.
14
Corresponding author: Val`ere Lambert,
vlambert@caltech.edu
–1–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
Abstract
15
Physics-based numerical modeling of earthquake source processes strives to predict quan-
16
tities of interest for seismic hazard, such as the probability of an earthquake rupture jump-
17
ing between fault segments. How to assess the predictive power of numerical models re-
18
mains a topic of ongoing debate. Here, we investigate how sensitive are the outcomes of
19
numerical simulations of sequences of earthquakes and aseismic slip to choices in numer-
20
ical discretization and treatment of inertial effects, using a simplified 2-D crustal fault
21
model with two co-planar segments separated by a creeping barrier. Our simulations demon-
22
strate that simplifying inertial effects and using oversized cells significantly affects the
23
resulting earthquake sequences, including the rate of two-segment ruptures. We find that
24
a number of fault models with different properties and modeling assumptions can pro-
25
duce comparable frequency-magnitude statistics and static stress drops but have rates
26
of two-segment ruptures ranging from 0 (single-segment ruptures only) to 1 (two-segment
27
ruptures only). For sufficiently long faults, we find that long-term sequences of events
28
can substantially differ even among simulations that are well-resolved by standard con-
29
siderations. In such simulations, some outcomes, such as static stress drops, are stable
30
among adequately-resolved simulations, whereas others, such as the rate of two-segment
31
ruptures, can be highly sensitive to numerical procedures and physical assumptions, and
32
hence cannot be reliably inferred. Our results emphasize the need to examine the po-
33
tential dependence of simulation outcomes on the modeling procedures and resolution,
34
particularly when assessing their predictive value for seismic hazard assessment.
35
1 Introduction
36
Earthquakes occur in the context of fault networks and many large earthquakes span
37
several fault segments. This reality brings about the issue of fault interaction and high-
38
lights the need for simulating earthquake source processes over several fault segments
39
and regional-scale fault networks. How dynamic ruptures navigate fault segmentation
40
has strong implications for seismic hazard analysis (Field, 2019). Earthquakes are ca-
41
–2–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
pable of jumping fault segments. For example, the 1992 Landers earthquake succeeded
42
in rupturing across at least 4 fault segments, amounting to a Mw 7.3 event (Sieh et al.,
43
1993). The 2016 Mw 7.8 Kaikoura earthquake ruptured at least 21 segments of the Marl-
44
borough fault system (Ulrich et al., 2019). Increasingly, seismological observations show
45
that it is not uncommon to see ruptures navigating and triggering subsequent ruptures
46
within fault networks, including the recent 2019 Mw 6.4 and 7.1 Ridgecrest earthquakes
47
(Ross et al., 2019), and the 2012 Mw 8.6 and 8.2 Indian Ocean earthquakes (Wei et al.,
48
2013), the largest and second-largest recorded strike-slip earthquakes to date. Yet, in any
49
given seismogenic region, the record of past large events is not long enough to forecast
50
the behavior of ruptures with respect to the existing fault segments, specifically how likely
51
would the rupture be to jump between nearby segments, prompting the discussion on
52
whether and how physics-based models may inform this and other questions important
53
for seismic hazard assessment (Field, 2019).
54
Determining what conditions allow a dynamic rupture to propagate or arrest are
55
key to understanding the maximum potential magnitude of an earthquake. Previous mod-
56
eling of single fully dynamic ruptures have shown great success in investigating earth-
57
quake propagation in nonplanar and multi-segment fault models, including step-overs
58
and branched geometries (Harris et al., 1991; Harris & Day, 1993, 1999; Kame et al., 2003;
59
Duan & Oglesby, 2006; Dunham et al., 2011a; Galvez et al., 2014; Douilly et al., 2015;
60
Lozos et al., 2015; Hu et al., 2016; Withers et al., 2018; Ando & Kaneko, 2018; Wollherr
61
et al., 2019; Ulrich et al., 2019). In particular, such modeling has shown that the abil-
62
ity of a rupture to propagate across segments depends on the stresses before the rupture
63
and shear resistance assumptions, as well as the geometry of the fault system. However,
64
single-rupture simulations need to select initial conditions and need additional assump-
65
tions to incorporate the effect of previous seismic and aseismic slip.
66
Fault processes involve both sequences of dynamic events and complex patterns of
67
quasi-static slip. Simulating this behavior in its entirety is a fascinating scientific prob-
68
–3–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
lem. However, even for the more pragmatic goal of physics-based predictive modeling
69
of destructive large dynamic events, it is still important to consider sequences of earth-
70
quakes and aseismic slip (SEAS), since prior slip events, including aseismic slip, may de-
71
termine where earthquakes would nucleate as well as modify stress and other initial con-
72
ditions before dynamic rupture. Furthermore, such simulations provide a framework for
73
determining physical properties consistent with a range of observations including geode-
74
tically recorded surface motions, microseismicity, past (including paleoseismic) events,
75
and thermal constraints, and hence may inform us about the current state of a fault seg-
76
ment or system and potential future rupture scenarios (e.g. Lapusta et al., 2000; Lapusta
77
& Rice, 2003; Liu & Rice, 2005; Ben-Zion & Rice, 1997; Chen & Lapusta, 2009; Kaneko
78
et al., 2010; Segall et al., 2010; Barbot et al., 2012; Noda & Lapusta, 2013; Erickson &
79
Dunham, 2014; Erickson & Day, 2016; Jiang & Lapusta, 2016; Lambert & Barbot, 2016;
80
Allison & Dunham, 2018; Lin & Lapusta, 2018; Cattania, 2019; Perry et al., 2020; Lam-
81
bert et al., 2021). However, simulating long-term slip histories is quite challenging be-
82
cause of the variety of temporal and spatial scales involved.
83
Recently, several earthquake simulators have been developed with the goal of sim-
84
ulating millions of earthquake ruptures over regional fault networks for tens of thousands
85
of years (Tullis et al., 2012; Richards-Dinger & Dieterich, 2012; Shaw et al., 2018). The
86
term ”simulators” typically refers to approaches that employ significant simplifications,
87
compared to most SEAS simulations, in solution procedures and physical processes, in
88
order to simulate earthquake sequences on complex, regional scale 3-D fault networks
89
for long periods of time. For example, earthquake simulators typically account only for
90
the quasi-static stress transfer due to earthquake events, ignoring wave-mediated stress
91
changes, aseismic slip/deformation, and fluid effects; employ approximate rule-based up-
92
date schemes for earthquake progression instead of solutions of the governing continuum
93
mechanics equations; and use oversized numerical cells. Such simplifications are currently
94
necessary to permit simulations of hundreds of thousands of events over hundreds of fault
95
segments that comprise the regional networks (Shaw et al., 2018). Earthquake simula-
96
–4–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
tors have matched a number of regional-scale statistical relations, including the Gutenberg-
97
Richter frequency-magnitude scaling (Shaw et al., 2018), and highlighted the importance
98
of large-scale fault and rupture interactions.
99
Here, we examine the sensitivity of the long-term interaction of fault segments to
100
choices in numerical discretization and representations of inertial effects in simulated se-
101
quences of earthquakes and aseismic slip, using a relatively simple 2-D model of two co-
102
planar strike-slip fault segments separated by a velocity-strengthening (VS) barrier. We
103
explore how considerations for adequate numerical resolution and convergence depend
104
on the physical assumptions and complexity of earthquake sequences as well as on the
105
modeling outcome of interest. We especially focus on the rate of earthquake ruptures jump-
106
ing across the VS barrier and examine whether reproducing comparable earthquake frequency-
107
magnitude statistics and static stress drops provides sufficient predictive power for the
108
jump rate, a quantity of interest to seismic hazard studies (Field, 2019).
109
2 Model setup and numerical resolution
110
Our simulations are conducted following the methodological developments of Lapusta
et al. (2000), Noda and Lapusta (2010) and Lambert et al. (2021). We consider a one-
dimensional (1-D) fault embedded into a 2-D uniform, isotropic, elastic medium (Fig-
ure 1). The 2-D model approximates a faulted crustal plate coupled to a moving sub-
strate using the idea of a constrained continuum (Lehner et al., 1981; Johnson, 1992).
Fault slip may vary spatially along-strike but it is depth-averaged through a prescribed
seismogenic thickness
λ
S
= 15 km, beneath which the elastic domain is coupled to a
substrate moving at the prescribed loading rate (
V
pl
= 10
−
9
m/s). The elastodynamic
equation for the depth-averaged displacement along-strike
u
(
x,y,t
) is given by (Lehner
et al., 1981; Kaneko & Lapusta, 2008):
Z
2
∂
2
u
∂x
2
+
∂
2
u
∂y
+
1
λ
2
eff
(
1
2
sign(
y
)
V
pl
t
−
u
)
=
1
c
s
∂
2
u
∂t
2
,
(1)
–5–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
where
λ
eff
= (
π/
4)
λ
S
and
Z
= 1
/
(1
−
ν
), with
ν
being the Poisson’s ratio. The effec-
111
tive wave speed along-strike for the crustal plane model is
c
L
=
Zc
s
, where
c
s
is the
112
shear wave speed. The along-strike slip is then given by
δ
(
x,t
) =
u
(
x,y
= 0
+
,t
)
−
113
u
(
x,y
−
,t
).
114
Our simulations resolve sequences of earthquakes and aseismic slip (SEAS) in their
115
entirety, including the gradual development of frictional instability and spontaneous nu-
116
cleation, dynamic rupture propagation, post-seismic slip that follows the event, and the
117
interseismic period between events (Figure 2). In all models, frictional resistance along
118
the fault interface is governed by the standard laboratory-derived rate-and-state friction
119
law with the state evolution described by the aging law (Dieterich, 1979; Ruina, 1983):
120
τ
=
σf
= (
σ
−
p
)
[
f
∗
+
a
ln
V
V
∗
+
b
ln
V
∗
θ
D
RS
]
,
(2)
dθ
dt
= 1
−
V θ
D
RS
,
(3)
where
σ
= (
σ
−
p
) is the effective normal stress,
σ
is the normal stress,
p
is the pore
121
pressure,
τ
is the shear stress,
f
is the friction coefficient,
V
is the slip velocity,
θ
is the
122
state variable,
D
RS
is the characteristic slip for the evolution of the state variable,
f
∗
is
123
the reference steady-state friction coefficient corresponding to a reference slip rate
V
∗
,
124
and
a
and
b
are the direct and evolution effect constitutive parameters, respectively.
125
At steady-state (constant slip velocity), the shear stress and state variable evolve
126
to their steady-state values
τ
ss
and
θ
ss
given by:
127
τ
ss
(
V
) = (
σ
−
p
)
[
f
∗
+ (
a
−
b
) ln
V
V
∗
]
,
(4)
θ
ss
(
V
) =
D
RS
V
.
(5)
–6–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
The combination of frictional properties such that (
a
−
b
)
>
0 results in steady-state
128
velocity-strengthening (VS) behavior, where the shear resistance increases with an in-
129
crease in slip velocity and where stable slip is expected. If (
a
−
b
)
<
0 then the fault
130
exhibits velocity-weakening (VW) behavior, in which case an increase in slip velocity leads
131
to a decrease in shear resistance, making these regions of the fault potentially seismo-
132
genic if their size exceeds a critical nucleation size.
133
Two theoretical estimates of the nucleation size in mode II are (Rice & Ruina, 1983;
134
Rubin & Ampuero, 2005):
135
h
∗
RR
=
π
4
μL
(1
−
ν
)(
b
−
a
)(
σ
−
p
)
;
h
∗
RA
=
2
π
μLb
(1
−
ν
)(
b
−
a
)
2
(
σ
−
p
)
,
(6)
where
μ
is the shear modulus. The estimate
h
∗
RR
was derived from the linear stability
136
analysis of steady frictional sliding by Rice and Ruina (1983). It also represents the crit-
137
ical cell size for steady-state quasi-static sliding such that larger cells can become un-
138
stable on their own. Thus
h
∗
RR
represents a key length scale to resolve for slow interseis-
139
mic processes and earthquake nucleation (Rice & Ruina, 1983; Lapusta et al., 2000). The
140
estimate
h
∗
RA
was determined in the parameter regime
a/b >
0
.
5 using the energy bal-
141
ance of a quasi-statically expanding crack (Rubin & Ampuero, 2005), and provides an
142
estimate of the minimum size for a slipping region that releases enough stored energy
143
to result in the radiation of waves.
144
We aim to explore the impact of numerical resolution on the long-term simulated
145
slip behavior of sequences of earthquakes and aseismic slip. The nucleation size,
h
∗
, es-
146
timated by either
h
∗
RR
or
h
∗
RA
from equation (6), is one length-scale that clearly needs
147
to be well resolved. Early resolution studies for sequences of events showed that reso-
148
lution of the nucleation scale
h
∗
RR
by 20 to 40 cells is required for stable numerical re-
149
sults (Lapusta et al., 2000). Later, the need to resolve the nucleation size by at least 20
150
cells was shown to be due to the more stringent criterion of resolving the region where
151
shear resistance breaks down at the rupture front, often referred to as the cohesive zone.
152
–7–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
The cohesive zone can be an order of magnitude smaller than the nucleation size, depend-
153
ing on the constitutive description (Day et al., 2005; Lapusta & Liu, 2009). The size of
154
the cohesive zone depends on the weakening rate
W
of shear stress with slip associated
155
with the constitutive law. The quasi-static estimate Λ
0
of the cohesive zone size at near-
156
zero rupture speed and constant
W
is given by:
157
Λ
0
=
C
1
μ
′
W
,
(7)
where
C
1
is a constant,
μ
′
=
μ
for mode III, and
μ
′
=
μ/
(1
−
ν
) for mode II (Rice,
158
1980). For standard rate-and-state friction with the aging form of the state variable evo-
159
lution, the weakening rate is given by
W
=
D
RS
/
(
b
σ
) (Lapusta & Liu, 2009) and:
160
Λ
0
=
C
1
μ
′
D
RS
b
σ
.
(8)
If one assumes that the traction distribution within the cohesive zone is linear, then the
161
constant
C
1
can be approximated as
C
1
= 9
π/
32 (Rice, 1980).
162
For fully dynamic rupture simulations, continuously resolving the breakdown pro-
163
cess at the rupture front becomes even more challenging as the cohesive zone size Λ ex-
164
hibits a contraction with increasing rupture speed
v
R
(e.g Rice, 1980):
165
Λ = Λ
0
A
−
1
(
v
R
);
A
−
1
II
=
(1
−
ν
)
c
2
s
D
v
2
R
(1
−
v
2
R
/c
2
s
)
1
/
2
;
A
−
1
III
= (1
−
v
2
R
/c
2
s
)
1
/
2
,
(9)
where
D
= 4(1
−
v
2
R
/c
2
s
)
1
/
2
(1
−
v
2
R
/c
2
p
)
1
/
2
−
(2
−
v
2
R
/c
2
s
)
1
/
2
with
c
p
=
√
2(1
−
ν
)
/
(1
−
2
ν
)
c
s
.
166
Note that
A
−
1
(0
+
) = 1, giving the quasi-static cohesive zone estimate Λ
0
when
v
R
=
167
0
+
. As the rupture speed approaches the limiting wave speed,
v
R
→
c
R
(Rayleigh wave
168
speed) for mode II and
v
R
→
c
s
(shear wave speed) for mode III, one has
A
−
1
(
v
R
)
→
169
0 and the width of the breakdown region approaches zero. Hence it becomes increasingly
170
more challenging to resolve the rupture front during fully dynamic simulations if the rup-
171
–8–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
ture accelerates towards the limiting speeds. Such acceleration typically occurs during
172
long enough propagation of dynamic rupture over favorable prestress, unless impeded
173
by additional factors such as unfavorable prestress or situations with increasing effec-
174
tive breakdown energy, e.g., due to off-fault inelasticity, thermal pressurization of pore
175
fluids, or navigating fault roughness (Poliakov et al., 2002; Andrews, 2005; Rice, 2006;
176
Okubo et al., 2019; Dunham et al., 2011a; Perry et al., 2020; Lambert & Lapusta, 2020).
177
Simulations of faults with rate-and-state friction and the aging form of the state vari-
178
able evolution embedded in elastic bulk result in ruptures with near-constant breakdown
179
energy (Perry et al., 2020) and this holds for most cases considered in this study. In sec-
180
tion 7, we show that adding an approximation of off-fault inelasticity to our simulations
181
that reduces the rupture speeds does not alter our conclusions.
182
In our model, the fault contains a frictional domain consisting of two VW regions
183
of length
λ
V W
= 32 km that are separated by a 2-km-long VS region that impedes rup-
184
ture propagation. We select large enough values of the velocity strengthening in the cen-
185
tral VS region so that the region acts like a barrier, requiring ruptures to jump/renucleate
186
on the other side of the barrier to propagate over the second segment. This region is a
187
proxy for what would be a gap in the fault connectivity, at least at the surface, requir-
188
ing the ruptures to jump across. The remainder of the frictional region surrounding the
189
VW segments has more mild VS properties (Figure 1). At the edges of the model, out-
190
side of the frictional domain, fault slip is prescribed at the loading plate rate. Values for
191
the model parameters used in our simulations are provided in Tables 1 and 2. We first
192
examine models with lower instability ratio
λ
V W
/h
∗
RR
that result in quasi-periodic se-
193
quences of events, and then consider models with higher instability ratios that result in
194
more complex earthquake sequences and qualitatively different convergence behavior.
195
–9–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
3 Resolving quasi-periodic fully dynamic sequences of earthquakes and
196
aseismic slip (SEAS)
197
Let us consider the simulated slip behavior of fault model M1 with instability ra-
198
tio
λ
V W
/h
∗
RR
= 21 (Table 2). Its quasi-static cohesive zone (Λ
0
= 1
.
1 km) should be
199
well-resolved by cell sizes of 12.5 and 25 m, with 88 and 44 cells over Λ
0
, respectively;
200
the nucleation size is even larger and hence also well-resolved. Consistently with these
201
considerations, these two well-discretized simulations produce the same relatively sim-
202
ple quasi-periodic sequences of earthquake events that periodically jump across the VS
203
barrier (Figure 2A & B). We clearly see that the results are the same for the two sim-
204
ulations with different resolutions, including the local evolution of slip rate and shear stress
205
during ruptures late in the earthquake sequence (Figure 2D-E). Note that the cohesive
206
zone evolves throughout the rupture process, shrinking with the increasing rupture speed
207
by 3-4 times in these simulations(Figure 2F-H) and the spatial discretization is fine enough
208
to adequately characterize the rupture front throughout the entire dynamic process. Both
209
simulations have the jump rate is 0.54; we define this rate of ruptures jumping across
210
the VS barrier within a given time period as the total number of ruptures that propa-
211
gate towards the barrier and result in seismic slip on both fault segments divided by the
212
total number of ruptures that propagate towards the barrier.
213
The variability between different ruptures in fault model M1 is generally mild, as
214
shown by their frequency-magnitude histograms (Figure 2C). To create the frequency-
215
magnitude histograms, we compute the moment for each simulated event in our 2-D mod-
216
els as
M
=
μ
δA
where the rupture area is defined with respect to the rupture length
217
L
R
and seismogenic depth
λ
S
, as
A
= (
π/
4)
L
2
R
when
L
R
≤
λ
S
and
A
=
L
R
λ
S
when
218
L
R
> λ
S
.
219
The quasi-periodic nature of events observed over the first 4000 years in well-resolved
220
simulations of fault model M1 persists in longer-duration simulations over 20,000 years,
221
resulting in similar long-term jump rates of 0.48 to 0.54 depending on the time interval
222
–10–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
considered (Figure 3). We also examine simulations of fault model M1 with different ini-
223
tial shear stress conditions and find that the long-term sequences of events converge to
224
the same quasi-periodic behavior upon adequate discretization, despite the initial few
225
events being different (Figure 3A vs B; details of initial shear stress distributions S1 and
226
S2 are provided in the Supplementary Materials). Simulations of fault model M1 thus
227
exhibit
long-term numerical convergence
upon adequate discretization, producing vir-
228
tually indistinguishable long-term slip behavior and a consistent rate of two-segment rup-
229
tures among simulations with differing initial conditions, after a sufficiently large initial
230
sequence of events.
231
Let us now consider simulations that use larger computational cells. The cell sizes
232
of 250 m and 125 m resolve the quasi-static cohesive zone Λ
0
with 4.5-9 cells (Figure 4).
233
While this resolution seems adequate (Day et al., 2005), one can anticipate that the dy-
234
namic shrinking of the cohesive zone size by 3-4 times would result in a more marginal
235
resolution of 1-3 cells. Indeed, we see that the simulated long-term sequences of events
236
and jump rates differ substantially from those of the well-resolved simulations (Figures
237
2A & B vs. 4A & B). Considering even larger cell sizes of 500 m and 1000 m brings fur-
238
ther differences in the event sequences and jump rates (Figure 5), with the earthquake
239
sequences that look plausible and not obviously numerically compromised even for the
240
largest cell sizes (Supplementary Figure S1). Note that the jump rate in simulations with
241
marginal and oversized cells is neither systematically larger nor smaller than the range
242
0.48-0.54 from the well-resolved cases, but varies from 0.25 to 0.95 depending on the choice
243
of numerical discretization.
244
Increasingly poor resolution of the dynamic cohesive zone at the rupture front and,
245
for the largest cell sizes, of the nucleation zone results in an increasing abundance of small
246
events (Figure 5), as had been shown in previous studies (Rice, 1993; Rice & Ben-Zion,
247
1996; Lapusta & Liu, 2009). Inadequate resolution of the dynamic rupture front prevents
248
simulating the actual stress concentration and promotes event arrest. Inadequate res-
249
–11–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
olution of the nucleation size enables individual cells or small number of cells to fail in-
250
dependently due to the inadequate resolution of the stress interactions (Rice, 1993; Rice
251
& Ben-Zion, 1996; Lapusta & Liu, 2009). Using sufficiently oversized cells can result in
252
power-law statistics in terms of the frequency-magnitude distribution of simulated earth-
253
quake ruptures (Figure 5E-J; Rice, 1993; Rice & Ben-Zion, 1996).
254
Note that the suggested minimum average resolution of 3 cells of the (variable) co-
255
hesive zone from the dynamic rupture study by Day et al. (2005) is not adequate for con-
256
vergent results in these earthquake sequence simulations. That criterion would be achieved
257
in this model for a cell size between the 250 m and 125 m. Yet the simulated long-term
258
behavior for those cell sizes is clearly different from the better-resolved and convergent
259
results with the cell sizes of 25 m and 12.5 m. At the same time, the criterion by Day
260
et al. (2005) works well for a single dynamic rupture as intended, since the first dynamic
261
events in simulations with cell sizes 12.5 m, 25 m, 125 m, and 250 m are quite similar
262
to each other (Supplementary Figure S2). The events are not identical, however; for ex-
263
ample, the average slip with the resolution of 12.5 m and 125 m differs by 0.7%. Clearly,
264
these differences - acceptable for a single event - accumulate in these highly nonlinear
265
solutions, resulting in different event statistics and jump rate (Figure 5).
266
We find that our fully dynamic 2-D simulations of fault model M1, which include
267
uniform VW properties with relatively mild weakening due to standard rate-and-state
268
friction, converge when the quasi-static cohesive zone estimate Λ
0
is discretized by at
269
least 22 cells, which translates to the average resolution of the dynamically variable co-
270
hesive zone size of 10-15 cells. Fault models with additional or different ingredients, such
271
as fault heterogeneity/roughness, more efficient weakening, 3D elastodynamics with 3D
272
faults, or different instability ratio, would require further considerations for resolution
273
requirements that result in convergent simulations. For example, as we discuss in sec-
274
tion 6, the convergence and resolution properties of models with higher instability ra-
275
tios, which result in more complex earthquake sequences, are qualitatively different.
276
–12–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
In the more complicated earthquake sequences observed in under-resolved simu-
277
lations of fault model M1, some statistics, such as the rate of two-segment ruptures, de-
278
pends on the specific period that one considers throughout the simulation. To explore
279
the variability in the event statistics and jump rate across the VS barrier in models with
280
different numerical resolution, we examine the jump rate over different 2000-year peri-
281
ods throughout longer term simulations of 20,000 years, using a sliding window of 1000
282
years starting at the beginning of the simulation (19 periods total; Figure 5). The choice
283
of a 2000-year period allows us to have a sufficient number (
∼
20) of large earthquakes
284
within a period to estimate jump rates. We also consider the outcomes for two differ-
285
ent initial conditions S1 and S2, as before. For the well-resolved simulations exhibiting
286
long-term convergence, the frequency-magnitude and 2000-year jump rate statistics for
287
simulations with different initiation conditions are comparable, with the jump rate for
288
all 2000-year periods being consistent with the overall 20,000 year jump rate (Figure 5A-
289
B). As the numerical resolution worsens, the sequences of events become more complex
290
with greater variability in rupture sizes and increased production of smaller events (Fig-
291
ure 5C-J). The jump rate during any 2000-year period also becomes more variable in marginally-
292
resolved simulations and can considerably differ from both the 20,000-year jump rate of
293
the same simulation as well as from the true jump rate in the well-resolved simulations.
294
Note that, despite being clearly affected by numerical resolution, the frequency-magnitude
295
and jump-rate distributions of inadequately resolved simulations can appear generally
296
consistent among simulations with similar cell sizes and different initial conditions (Fig-
297
ure 5 left vs. right columns). In other words, even if simulations using marginal or over-
298
sized cells produce comparable statistical properties for different initial conditions, these
299
characteristics do not necessarily represent robust features of the physical system but
300
rather may still be numerical artifacts.
301
–13–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
4 Interaction of fault segments in simulations with quasi-dynamic ap-
302
proximation for inertial effects
303
Many numerical studies of long-term fault behavior utilize quasi-dynamic solutions
304
to the equations of motion, in which the wave-mediated stress transfers during the co-
305
seismic phase are replaced with a radiation damping approximation (Rice, 1993). The
306
quasi-dynamic approximation substantially reduces the computational expense of the sim-
307
ulation, as the consideration of stress redistribution by waves requires substantial ad-
308
ditional storage and computational expense. Considerable insight into fault mechanics
309
has been derived from studies using quasi-dynamic formulations, particularly when such
310
approximations are used to incorporate new physical effects that may otherwise result
311
in prohibitive computational expense, as well as in scenarios where it may be argued that
312
inertial effects are relatively mild, such as during earthquake nucleation or during aseis-
313
mic slip transients (Rice, 1993; Segall & Rice, 1995; Liu & Rice, 2005, 2007; Rubin &
314
Ampuero, 2005; Segall et al., 2010; Liu, 2014; Lambert & Barbot, 2016; Erickson et al.,
315
2017; Allison & Dunham, 2018). However, as with all approximations, it is important
316
to be aware of how such simplifications modify the outcome of study (Thomas et al., 2014).
317
Let us review the quasi-dynamic approximation for inertial effects during sliding
318
and study their implications for the long-term interaction of two fault segments. In the
319
2D boundary integral formulation, the elastodynamic shear stress along a 1D fault plane,
320
can be expressed as (Cochard & Madariaga, 1994; Perrin et al., 1995):
321
τ
(
x,t
) =
τ
0
(
x,t
) +
φ
static
(
x,t
) +
φ
dynamic
(
x,t
)
−
ηV
(
x,t
)
,
(10)
where
τ
0
(
x,t
) are the ”loading” tractions (i.e. the stress induced on the fault plane if
322
it were constrained against any slip),
φ
static
(
x,t
) and
φ
dynamic
(
x,t
) represent the static
323
and dynamic contributions to the stress transfer along the fault, respectively, and the
324
last term represents radiation damping (
η
=
μ/
(2
c
s
) for mode III).
325
–14–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
The static solution for the equations of motion would only contain
φ
static
, which
326
depends only on the current values of slip along the fault. However, the static solution
327
does not exist during dynamic rupture when inertial effects becomes important.
φ
dynamic
328
and
ηV
both arise due to the inertial effects.
φ
dynamic
represents the wave-mediated stress
329
interactions along the interface and this term is challenging to compute as it requires cal-
330
culating convolutions on time and storing the history of deformation. Radiation damp-
331
ing
ηV
is much easier to incorporate as it depends on the current slip rate, and repre-
332
sents part of the radiated energy (Rice, 1993). The quasi-dynamic approximation, in which
333
φ
dynamic
is ignored and only
ηV
is included, allows the solution to exist during inertially-
334
controlled dynamic rupture. However, the solution is altered from the true elastodynamic
335
representation.
336
Let us consider the long-term behavior of fault model M1, as examined in section
337
3, but now using the quasi-dynamic approximation. For well-resolved quasi-dynamic sim-
338
ulations of fault model M1, we find that the long-term slip behavior of the two fault seg-
339
ment system is even simpler than for the fully dynamic case, with ruptures being exclu-
340
sively isolated to individual segments and the jump rate being zero (Figure 6A). For sim-
341
ulations with the increasing cell size, and thus decreasing spatial resolution, we see in-
342
creased variability in the size of the individual ruptures, to the point where some marginally-
343
resolved simulations produce ruptures that jump across the VS barrier, whereas well-
344
resolved simulations of the same fault model never do (Figure 6B-C). The increasing cell
345
size also leads to increased production of smaller events and more complicated fault be-
346
havior, similarly to the fully dynamic simulations (Figure 6D-F).
347
In addition to substantially reducing the computational expense associated with
348
calculating the wave-mediated stress transfers, quasi-dynamic simulations place milder
349
constraints on the spatial resolution since the cohesive zone always remains near the quasi-
350
static estimate, Λ
≈
Λ
0
= Λ(
v
R
= 0
+
) (Figure 6G-H). This is because the stress trans-
351
fer calculated for the ruptures is always quasi-static, and the much stronger stress trans-
352
–15–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
fer due to waves is ignored (Figure 7E). As a result, the quasi-dynamic simulations pro-
353
duce significantly smaller slip velocities and rupture speeds than the fully dynamic ones
354
7A-C, consistent with previous studies (Lapusta & Liu, 2009; Thomas et al., 2014).
355
One can attempt to enhance the slip rates and rupture speeds in the quasi-dynamic
356
simulations by reducing the radiation damping term
η
; this can be interpreted as increas-
357
ing the effective shear wave speed in the radiation damping term
c
enh.
s
=
βc
s
, thus al-
358
lowing for higher slip rates (Lapusta & Liu, 2009). We compare the enhanced quasi-dynamic
359
simulations (
β
= 3) with the standard quasi-dynamic (
β
= 1) and fully dynamic sim-
360
ulations of fault model M1 (Figure 8). Decreasing the radiation damping increases the
361
effective rupture speed and slip rate (Figure 7A -C) in comparison to the standard quasi-
362
dynamic simulation, however, for the parameters considered, it does not substantially
363
alter the long-term interactions of the two fault segments, nor match the rate of ruptures
364
jumping across the VS barrier in the fully dynamic case (Figure 8).
365
In comparing the three simulations with different treatment of the inertial effects,
366
it is clear that the fully dynamic ruptures result in higher slip rates and narrowing of
367
the cohesive zone (Figure 7). For simulations with standard rate-and-state friction, the
368
peak shear stresses vary mildly from fully dynamic versus quasi-dynamic representations,
369
as they are limited by the shear resistance of the fault, which has a relatively mild log-
370
arithmic dependence on slip rate. However, the stress transfer along the fault substan-
371
tially differs for fully dynamic versus quasi-dynamic representations (Figure 7E). The
372
difference between the stress transfer term and the shear stress is accommodated by the
373
radiation damping
ηV
, which results in higher slip rates
V
to balance the larger dynamic
374
stress stranfers (Figure 7C - E). Hence while the resolved peak shear stresses along the
375
fault may be comparable due to the specific choice of the constitutive relationship, the
376
rupture dynamics and kinematics, as seen through the stress transfer, slip rate, and rup-
377
ture speed along the fault, differ considerably with and without the inclusion of full in-
378
ertial effects.
379
–16–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
These larger dynamic stress transfers facilitate the triggering and continued prop-
380
agation of slip on the neighboring fault segment, rather than leaving the rupture to al-
381
ways be arrested by the creeping barrier, as in the well-resolved quasi-dynamic simula-
382
tions (Figure 8). Decreasing the radiation damping term allows for somewhat higher slip
383
rates and arbitrarily higher rupture speeds, but it does not mimic the full effects in the
384
dynamic stress transfer, particularly at the rupture front. As the result, the fully dynamic
385
simulations have higher jump rates. The differences between fully dynamic and quasi-
386
dynamic approximations can be even more substantial for models with enhanced weak-
387
ening at seismic slip rates from the flash heating of contact asperities or the thermal pres-
388
surization of pore fluids (Thomas et al., 2014).
389
5 Constraining rupture jump rates using earthquake frequency-magnitude
390
statistics
391
Two common observations about natural earthquakes and regional seismicity are
392
the average static stress drops between 1 to 10 MPa independently of the event magni-
393
tude (e.g Allmann & Shearer, 2009; Ye et al., 2016) as well as the frequency-magnitude
394
statistics of earthquakes within a region, which commonly follow the Gutenberg-Richter
395
power law relation (Field et al., 2013). Earthquake simulators are capable of matching
396
these observations (Shaw et al., 2018). An important question is whether matching these
397
constraints endows simulators with predictive power for other quantities of interest to
398
seismic hazard assessment, such as the probability of multiple fault-segment ruptures,
399
despite using approximations for inertial effects and oversized computational cells.
400
Let us consider this question using simulations of earthquake sequences in five fault
401
models with the same fault geometry but different friction properties and different as-
402
sumptions about inertial effects, and one additional model in which the effective seismo-
403
genic depth
λ
S
is slightly reduced from 15 to 14 km (Figure 9, Table 2). All six mod-
404
els have comparable nucleation and quasi-static cohesive zone sizes (Table 2) and use over-
405
–17–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
sized cells of ∆
x
= 1000 m (An example of well-resolved simulations with similar con-
406
clusions is given in section 7). The six simulations produce comparable frequency-magnitude
407
distributions, characterized by a b-value of 0.3-0.4 for 4000 years of the simulated time.
408
All six simulations also produce ruptures with comparable average static stress drops (Sup-
409
plementary Figure S3), with values typically between 1 and 10 MPa, as commonly in-
410
ferred for natural earthquakes (Allmann & Shearer, 2009; Ye et al., 2016).
411
However, the probability of a rupture jumping across the VS barrier varies dramat-
412
ically among the six simulations, ranging from 0 to near 100%. This substantial variabil-
413
ity in jump rate for simulations with comparable frequency-magnitude statistics persists
414
in longer-duration simulations over 20,000 years, where both the 20,000-year jump rate
415
and distributions of jump rates within individual 2000-year periods can substantially dif-
416
fer (Figure 10). In particular, fault model M1 results in a jump rate of 0 for the quasi-
417
dynamic simulation and near 1 for the fully dynamic simulation (Figures 9 and 10A vs.
418
D), despite having similar frequency-magnitude statistics. This case illustrates how us-
419
ing approximations for inertial effects may considerably bias estimates of the actual rate
420
of multi-segment ruptures, even if the frequency-magnitude statistics and static stress
421
drops are comparable. In addition, the suite of simulations suggest that the probabil-
422
ity of ruptures jumping across the VS barrier is sensitive to variations in the frictional
423
parameters, effective normal stress, as well as minor changes in the seismogenic depth.
424
The results from our simple 2-D modeling suggest that reproducing static stress
425
drops and frequency-magnitude statistics does not provide substantial predictive power
426
for the long-term interaction of fault segments. These results are perhaps not surpris-
427
ing given that many combinations of rate-and-state properties and effective normal stress
428
may produce ruptures with comparable static stress changes (Supplementary Figure S3),
429
but different overall levels of shear resistance. Moreover, numerical studies have shown
430
that fault models including enhanced dynamic weakening may also produce nearly magnitude-
431
invariant static stress drops with reasonable values between 1 - 10 MPa (Perry et al., 2020).
432
–18–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
Such models with enhanced weakening result in larger dynamic stress variations which
433
may mediate longer-range interactions among faults. However, enhanced weakening can
434
also draw the average stress along the fault further away from nucleation conditions (Jiang
435
& Lapusta, 2016; Lambert et al., 2021), which may produce less favorable conditions for
436
dynamic triggering (Ulrich et al., 2019).
437
Similarly, a number of studies have demonstrated that power-law frequency-magnitude
438
statistics can be reproduced in many models, including discrete fault models (Burridge
439
& Knopoff, 1967; Bak & Tang, 1989; Olami et al., 1992), continuum fault models that
440
are inadequately resolved and therefore numerically discrete (Ben-Zion & Rice, 1995),
441
and continuum models with larger instability ratio (Wu & Chen, 2014; Cattania, 2019).
442
In other words, Gutenberg-Richter statistics is consistent with a model having many po-
443
tential rupture sizes, such as many individual faults of varying size or even a single fault
444
that can host earthquakes of many sizes, between the nucleation size and fault dimen-
445
sions. Therefore Gutenberg-Richter statistics may be compatible with a range of fault
446
and/or fault network properties, and may not pose a considerable physical constraint
447
on its own.
448
6 Resolution and convergence of SEAS simulations of faults with higher
449
instability ratios
450
As discussed in section 3, we find that the discretization required to achieve long-
451
term numerical convergence in simulations of fault model M1, with instability ratio of
452
λ
V W
/h
∗
RR
≈
21, is more stringent than the current standards based on simulations of
453
single dynamic ruptures and shorter SEAS simulations with lower instability ratios (Day
454
et al., 2005; Lapusta & Liu, 2009). It has been demonstrated that fault models with rel-
455
atively low instability ratios can result in quasi-periodic behavior, as seen in fault model
456
M1 (Figure 2), whereas increasing the instability ratio can lead to more variable sequences
457
of events with partial-segment ruptures of different rupture size, potentially consistent
458
–19–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
with Gutenberg-Richter scaling (e.g Lapusta et al., 2000; Lapusta & Rice, 2003; Wu &
459
Chen, 2014; Michel et al., 2017; Cattania, 2019). As simulations with higher instabil-
460
ity ratios can produce ruptures with a wider variety of rupture sizes, with the rupture
461
size depending on the prestress conditions before rupture nucleation, one could hypoth-
462
esize that simulations of fault models with higher instability ratios may be more sensi-
463
tive to how the evolution of shear stress is resolved over long-term fault behavior.
464
To test that, let us consider sequences of events in fault model M5 (Table 2), which
465
has smaller characteristic slip distance, hence smaller nucleation size (
h
∗
RR
≈
603 m),
466
and larger instability ratio (
λ
V W
/h
∗
RR
= 53 vs. 21 in M1). Interestingly, we find that
467
the long-term sequence of simulated events in this model is not the same for finely-discretized
468
simulations with cell sizes of 25, 12.5 and 6.25 m (Figure 11), in which the quasi-static
469
cohesive zone Λ
0
is resolved by 18, 36 and 72 cells, respectively. The simulations pro-
470
duce nearly identical fault behavior for the first several hundred years of simulated time,
471
but then eventually begin to differ (Figure 11A-C).
472
Let us consider the first event in the three simulations of model M5 with fine dis-
473
cretization (Figure 11A-C), which all have the same initial conditions. If we examine the
474
local evolution of shear stress vs. slip at two spatial points in the simulations, the results
475
are virtually identical (Figure 12A-B), suggesting that a single dynamic rupture in these
476
finely-discretized simulations is adequately resolved. The evolution of shear stress and
477
slip rate at the rupture front with time is also well-resolved for each individual simula-
478
tion. While the different spatial resolutions result in small variations in the timing and
479
magnitude of the resolved properties at specified locations (Figure 12C-F), these differ-
480
ences are well within of what is considered well-resolved and convergent in prior stud-
481
ies (e.g Day et al., 2005). Early in the rupture, shortly after nucleation (near
x
= 30
482
km), the rupture front is almost identical in the three simulations. (Figure 12C & E).
483
As the rupture continues, small numerical differences for different resolutions result in
484
minor differences in the rupture, such as less than 0.08% difference in the rupture ar-
485
–20–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
rival time and 2% difference in the peak slip rate between the two best-resolved simu-
486
lations at the location close to the end of the rupture (Figure 12D & F). Such minor dif-
487
ferences arise even for fine resolutions due to cumulative effects of slightly different rep-
488
resentations of the solution by the discrete cells; for example, the fixed computational
489
cells sample slightly different portions of the passing rupture front, leading to small ac-
490
cumulating differences in the magnitude of the shear stress and slip rate.
491
These small differences - that do not substantially alter the resulting rupture char-
492
acteristics of individual events - do eventually alter the resulting earthquake sequences.
493
For several ruptures early on in finely-discretized simulations, the slip and shear stress
494
distributions before and after individual events are virtually indistinguishable (Figure
495
13A-B). However, eventually the small variations accumulate, resulting in enough dif-
496
ferences in prestress conditions to cause more substantial differences in rupture lengths
497
and amounts of slip within individual events, as well as changes in timing and location
498
of earthquake nucleations (Figure 13C-E). As a result, the long-term history of sequences
499
of slip events is altered (Figure 13F), including the rate of ruptures that jump across the
500
VS barrier. We hypothesize that this alteration occurs for higher but not lower insta-
501
bility ratios due to more complex earthquake sequences in the latter case, although this
502
issue requires further study.
503
Despite the specific sequences of events being different in the finely-discretized sim-
504
ulations shown in Figure 11A-C, we do find that certain outcomes are quite similar be-
505
tween these simulations, such as relationships between average static stress drop and seis-
506
mic moment, average slip and rupture length, and breakdown energy and average slip,
507
as well as general characteristics of the evolution of average shear stress and shear heat-
508
ing with time (Figure 14). Other parameters, such as the rate of ruptures jumping from
509
one fault segment to another, are sensitive to numerical resolution even in these finely-
510
resolved simulations, although they have relatively similar values (from 0.64 to 0.78). This
511
highlights how the criteria for adequate discretization in numerical simulations can de-
512
–21–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
pend on both the physical problem being considered and the outcome of interest. Note
513
that while it is plausible that further discretization of fault model M5 would result in
514
eventual convergence, and thus potentially a true rate of two-segment ruptures, the spa-
515
tial discretization considered in this study is already much finer than those considered
516
in most numerical SEAS studies, especially in more realistic models of 2D faults in 3D
517
media which are often challenged to resolve Λ
0
by even 3 cells.
518
While the specific rate of ruptures jumping across the VS barrier varies among these
519
finely-discretized simulations of fault model M5, it is possible that some broader statis-
520
tical features of the jump rate are more robust. We examine the frequency-magnitude
521
and 2000-year jump rate statistics for the long-term sequences of events in simulations
522
of model M5 with different discretization. While the distributions mildly vary among finely-
523
discretized simulations with differing cell sizes (12.5 m and 25 m), they are comparable
524
(Figure 15 and Supplementary Figure S4). Thus, one can ascertain information about
525
the probability distribution for the rate of multi-segment ruptures, even if specific re-
526
sults vary due to numerical discretization. Such small numerical perturbations could po-
527
tentially be considered representative of various sources of physical perturbations on nat-
528
ural faults, and the statistical consistency of the distributions could be explored by pro-
529
ducing ensembles of simulations with varying initial conditions. However, our results sug-
530
gest that it is still important to sufficiently resolve the rupture process as the statisti-
531
cal distributions for rupture properties in simulations using oversized cells can be more
532
substantially impacted by numerical artifacts and considerably vary from simulations
533
with finer discretization (Figures 15 and 16).
534
7 Resolution and convergence in SEAS simulations with moderate rup-
535
ture speeds due to an approximation for off-fault plasticity
536
While the 2-D fault models discussed in this study can be considered relatively sim-
537
ple, in some ways they can be particularly challenging to resolve. In fault models with
538
–22–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.
manuscript submitted to
JGR-Solid Earth
purely elastic bulk, dynamic ruptures are able to accelerate to rupture speeds close to
539
the limiting values
c
L
(e.g. Figure 7 for fault model M1), making it difficult to resolve
540
the significantly shrinking cohesive zone Λ. For example, during fully dynamic ruptures
541
in simulations of fault model M5, the rupture speed approaches 0
.
99
c
L
and the cohesive
542
zone shrinks more than 7 times to about 63.5 m. In real rocks, high slip rates and hence
543
high strain rates associated with dynamic rupture would be mitigated by off-fault inelas-
544
tic behavior around the rupture front, which would contribute to limiting the rupture
545
speed (Andrews, 2004; Dunham et al., 2011b).
546
In order to examine how conditions for resolution and convergence may differ in
547
long-term SEAS simulations with more moderate rupture speeds, we approximate the
548
effects of off-fault yielding by employing a limit on the slip velocity, as suggested by Andrews
549
(2004) and discussed in detail in Lambert et al. (2021). We consider long-term fully dy-
550
namic simulations of fault model M5 with the slip velocity limited to 2 m/s in order to
551
maintain rupture speeds around 0.8
c
L
, consistent with the cohesive zone shrinking by
552
about a factor of 2 from the quasi-static estimate.
553
Surprisingly, the finely-discretized simulations of fault model M5 with limited rup-
554
ture speed still produce differing sequences of events, despite the rupture front and lo-
555
cal behavior being well-resolved and nearly identical for cell sizes of 6.25 to 25 m (Fig-
556
ure 17 and Supplementary Figure S5). As with the standard fully dynamic simulations
557
without the plasticity approximation, well-resolved simulations of fault model M5 with
558
the velocity limit are nearly identical for the initial few sequences (Figure 18A-B). How-
559
ever, the sequences of events begin to differ due to slight differences in how the evolu-
560
tion of shear stress is resolved during a slow-slip transient within the nucleation region
561
of an impending rupture, resulting in a 3-year delay between the nucleation of the sub-
562
sequent rupture in each simulation (Figure 18C-D). As discussed earlier for the standard
563
fully dynamic simulations, the small differences in prestress lead to mild differences in
564
slip and rupture size in subsequent events, which eventually compound to produce more
565
–23–
ESSOAr | https://doi.org/10.1002/essoar.10506727.1 | Non-exclusive | First posted online: Fri, 9 Apr 2021 06:49:59 | This content has not been peer reviewed.