of 59
manuscript submitted to
JGR: Solid Earth
Community-Driven Code Comparisons for Three-Dimensional
1
Dynamic Modeling of Sequences of Earthquakes and Aseismic Slip
2
(SEAS)
3
Junle Jiang
1
, Brittany A. Erickson
2
, Valère R. Lambert
3
,
4
Jean-Paul Ampuero
4
, Ryosuke Ando
5
, Sylvain D. Barbot
6
, Camilla Cattania
7
,
5
Luca Dal Zilio
8
–
9
, Benchun Duan
10
, Eric M. Dunham
11
, Alice-Agnes Gabriel
12
–
13
,
6
Nadia Lapusta
8
, Duo Li
12
, Meng Li
14
, Dunyu Liu
15
, Yajing Liu
16
,
7
So Ozawa
5
, Casper Pranger
9
–
12
, Ylona van Dinther
9
–
14
8
1
School of Geosciences, University of Oklahoma, Norman, OK, USA
9
2
Department of Computer and Information Science, University of Oregon, Eugene, OR, USA
10
3
Department of Earth and Planetary Sciences, University of California, Santa Cruz, CA, USA
11
4
Université Côte d’Azur, IRD, CNRS, Observatoire de la Côte d’Azur, Géoazur, France
12
5
Department of Earth and Planetary Science, University of Tokyo, Japan
13
6
Department of Earth Sciences, University of Southern California, Los Angeles, CA, USA
14
7
Department of Earth, Atmospheric, and Planetary Sciences, Massachusetts Institute of Technology, Cambridge, MA, USA
15
8
Seismological Laboratory and Department of Mechanical and Civil Engineering,
16
California Institute of Technology, Pasadena, CA, USA
17
9
Institute of Geophysics, Department of Earth Sciences, ETH Zurich, Zurich, Switzerland
18
10
Department of Geology and Geophysics, Texas A&M University, College Station, TX, USA
19
11
Department of Geophysics, Stanford University, Stanford, CA, USA
20
12
Department of Earth and Environmental Sciences, Ludwig-Maximilians-Universität München, Munich, Germany
21
13
Institute of Geophysics and Planetary Physics, Scripps Institution of Oceanography,
22
University of California, San Diego, CA, USA
23
14
Department of Earth Sciences, Utrecht University, Utrecht, Netherlands
24
15
Institute for Geophysics, University of Texas at Austin, TX, USA
25
16
Department of Earth and Planetary Sciences, McGill University, Montréal, QC, Canada
26
Key Points:
27

We pursue community efforts to develop code verification benchmarks for three-dimensional
28
earthquake rupture and crustal faulting problems
29

We assess the agreement and discrepancies of seismic and aseismic fault behavior among
30
simulations based on different numerical methods
31

Our comparisons lend confidence to numerical codes and reveal sensitivities of model observables
32
to major computational and physical factors
33
Corresponding author: Junle Jiang,
jiang@ou.edu
–1–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
Abstract
34
Dynamic modeling of sequences of earthquakes and aseismic slip (SEAS) provides a self-consistent,
35
physics-based framework to connect, interpret, and predict diverse geophysical observations across
36
spatial and temporal scales. Amid growing applications of SEAS models, numerical code verification
37
is essential to ensure reliable simulation results but is often infeasible due to the lack of analytical
38
solutions. Here, we develop two benchmarks for three-dimensional (3D) SEAS problems to compare
39
and verify numerical codes based on boundary-element, finite-element, and finite-difference methods,
40
in a community initiative. Our benchmarks consider a planar vertical strike-slip fault obeying a
41
rate- and state-dependent friction law, in a 3D homogeneous, linear elastic whole-space or half-space,
42
where spontaneous earthquakes and slow slip arise due to tectonic-like loading. We use a suite of
43
quasi-dynamic simulations from 10 modeling groups to assess the agreement during all phases
44
of multiple seismic cycles. We find excellent quantitative agreement among simulated outputs
45
for sufficiently large model domains and small grid spacings. However, discrepancies in rupture
46
fronts of the initial event are influenced by the free surface and various computational factors.
47
The recurrence intervals and nucleation phase of later earthquakes are particularly sensitive to
48
numerical resolution and domain-size-dependent loading. Despite such variability, key properties
49
of individual earthquakes, including rupture style, duration, total slip, peak slip rate, and stress
50
drop, are comparable among even marginally resolved simulations. Our benchmark efforts offer
51
a community-based example to improve numerical simulations and reveal sensitivities of model
52
observables, which are important for advancing SEAS models to better understand earthquake
53
system dynamics.
54
Plain Language Summary
55
Earthquakes and fault zone processes occur over time scales ranging from milliseconds to millennia
56
and longer. Computational models are increasingly used to simulate sequences of earthquakes and
57
aseismic slip (SEAS). These simulations can be connected to diverse geophysical observations,
58
offering insights into earthquake system dynamics. To improve these simulations, we pursue community
59
efforts to design benchmarks for 3D SEAS problems. We involve earthquake researchers around
60
the globe to compare simulation results using different numerical codes. We identify major factors
61
that contribute to the discrepancies among simulations. For example, the spatial dimension and
62
resolution of the computational model can affect how earthquakes start and grow, as well as how
63
frequently they recur. Code comparisons are more challenging when we consider the Earth’s surface
64
in the simulations. Fortunately, we find that several key characteristics of earthquakes are accurately
65
reproduced in simulations, such as the duration, total movement, maximum speed, and stress change
66
on the fault, even when model resolutions are not ideal. These exercises are important for promoting
67
a new generation of advanced models for earthquakes. Understanding the sensitivity of simulation
68
–2–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
outputs will help test models against real-world observations. Our community efforts can serve as a
69
useful example to other geoscience communities.
70
1 Introduction
71
Physics-based computational models of dynamic processes in the Earth are increasingly used
72
to understand and predict observations from the lab and field across spatial and temporal scales,
73
addressing fundamental questions in various branches of solid Earth research. In earthquake science,
74
models of earthquake source processes are aimed at capturing dynamic earthquake ruptures from
75
seconds to minutes and slow slip processes subject to short-term anthropogenic or environmental
76
forcing, or tectonic loading over timescales of years and longer. For individual earthquakes,
dynamic
77
rupture simulations
have emerged as powerful tools to reveal the influence of fault structure, geometry,
78
constitutive laws, and prestress on earthquake rupture propagation and associated ground motion
79
(e.g.,
Andrews
, 1976a,b;
Ben-Zion
, 2001;
Bhat et al.
, 2007;
Bizzarri and Cocco
, 2003, 2006;
Day
,
80
1982;
Das and Aki
, 1977;
Duan and Day
, 2008;
Dunham et al.
, 2011a,b;
Gabriel et al.
, 2012;
81
Harris et al.
, 1991, 2021;
Kozdon and Dunham
, 2013;
Lozos et al.
, 2011;
Ma and Beroza
, 2008;
82
Madariaga et al.
, 1998;
Mikumo and Miyatake
, 1978, 1993;
Nielsen et al.
, 2000;
Olsen et al.
, 1997;
83
Ripperger et al.
, 2007;
Shi and Day
, 2013;
Tinti et al.
, 2021;
Wollherr et al.
, 2019;
Xu et al.
, 2015).
84
These simulations are limited to single-event scenarios and subject to imposed artificial prestress
85
conditions and ad hoc nucleation procedures. For larger-scale fault network systems,
earthquake
86
simulators
aim to produce complex spatiotemporal characteristics of seismicity over millennial
87
time scales (
Richards-Dinger and Dieterich
, 2012;
Robinson and Benites
, 1995, 1996, 2001;
Shaw
88
et al.
, 2018;
Tullis et al.
, 2012). The formidable computational demand inevitably requires simplification
89
and approximation of some key physical features that could influence or dominate earthquake and
90
fault interactions, such as seismic waves, slow slip, tectonic loading, and inelastic response.
91
To understand earthquake system dynamics, it has been widely recognized that we need models
92
that simulate fault behavior over multiple seismic events and the intervening periods of aseismic
93
deformation. To address this need, numerical simulations of Sequences of Earthquakes and Aseismic
94
Slip (SEAS) are developed to consider all phases of earthquake faulting, from slow loading to
95
earthquake nucleation, propagation and termination over time scales of milliseconds to millennia
96
in a unified, self-consistent framework (Figure 1;
Ben-Zion and Rice
, 1995;
Lapusta et al.
, 2000;
97
Rice
, 1993). While retaining computational rigor, SEAS models incorporate the structure, rock
98
properties, friction, and rheology of a fault zone, and produce the pre-, inter-, and post-seismic slip
99
and the resulting stress redistribution that ultimately lead to spontaneous earthquake nucleation
100
and dynamic ruptures. SEAS models can include many physical processes relevant to long-term
101
slip, such as evolving shear resistance of the fault zone affected by shear heating, fluid effects, and
102
interseismic healing, wave-mediated inertial effects during dynamic rupture, folding, viscoelasticity,
103
–3–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
and fluid flow (e.g.,
Allison and Dunham
, 2018;
Barbot
, 2018;
Lambert and Barbot
, 2016;
Noda
104
and Lapusta
, 2010;
Sathiakumar et al.
, 2020;
Thomas et al.
, 2014;
Zhu et al.
, 2020). This modeling
105
framework can help determine and quantify which physical factors control diverse observables such
106
as ground deformation and shaking, and the frequency, size, and rupture style of microseismicity
107
and large earthquakes. SEAS models also bridge the domains of dynamic rupture simulations and
108
earthquake simulators, providing physically justified approximations and self-consistent choices for
109
initial conditions and earthquake nucleation procedures.
110
Developments in SEAS models over the past two decades have led to increased diversity and
111
complexity of models and closer connections between simulations and observations from the lab
112
and field. For example, numerical models have been combined with seismic and geodetic observations
113
to study fault frictional properties (e.g.,
Barbot et al.
, 2009;
Dublanchet et al.
, 2013;
Floyd et al.
,
114
2016;
Hori et al.
, 2004;
Jiang and Fialko
, 2016;
Johnson et al.
, 2006;
Mitsui and Iio
, 2011;
Tymofyeyeva
115
et al.
, 2019), tremor and slow slip (e.g.,
Dal Zilio et al.
, 2020;
Dublanchet
, 2018;
Hawthorne and
116
Rubin
, 2013;
Luo and Ampuero
, 2018;
Mele Veedu and Barbot
, 2016;
Shibazaki and Iio
, 2003;
117
Wang and Barbot
, 2020), foreshock and aftershock sequences (e.g.,
Cattania and Segall
, 2021;
118
Kaneko and Lapusta
, 2008;
Perfettini and Avouac
, 2007;
Noda et al.
, 2013), and characteristics
119
of small and large earthquake ruptures (e.g.,
Barbot et al.
, 2012;
Cattania and Segall
, 2019;
Chen
120
and Lapusta
, 2009;
Jiang and Lapusta
, 2016, 2017;
Lambert and Lapusta
, 2021). The framework
121
of earthquake sequence modeling is also adopted in diverse settings, which include subduction
122
zones (e.g.,
Hori et al.
, 2004;
Liu and Rice
, 2005, 2007;
Li and Liu
, 2016, 2017;
Shi et al.
, 2020;
123
Van Dinther et al.
, 2013), collision zones (e.g.,
Dal Zilio et al.
, 2018;
Michel et al.
, 2017;
Qiu et al.
,
124
2016), and induced seismicity phenomena (e.g.,
Dieterich et al.
, 2015;
Kroll and Cochran
, 2021;
125
McClure and Horne
, 2011), among many applications.
126
While researchers continue to build more advanced and detailed SEAS models, verification
127
of different numerical codes is essential to ensure credible and reproducible results, and sustain
128
scientific progress. In practice, analytical solutions are generally not available, even for simple
129
SEAS problems, and convergence of simulations to a high-resolution reference case may not always
130
detect systematic issues in complex numerical codes. An alternative means for verifying model
131
results are comparisons of independent numerical codes from different research groups. As an
132
example, the SCEC/USGS Spontaneous Rupture Code Verification Project pioneered the code
133
comparison exercise and improved confidence in the outcomes of dynamic rupture simulations
134
(
Barall and Harris
, 2015;
Day et al.
, 2005;
Harris et al.
, 2009, 2018).
135
Verification of SEAS models is confronted with distinct challenges, due to the wide range
136
of spatial and temporal scales that characterize the earthquake source behavior and the diversity
137
of numerical algorithms and codes. For example, codes based on the spectral boundary element
138
method (SBEM) (
Barbot
, 2021;
Lapusta and Rice
, 2003;
Lapusta and Liu
, 2009) are highly efficient
139
–4–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
in solving for fully dynamic earthquake ruptures, albeit with relatively simple fault geometry and
140
bulk. Codes based on the boundary element method (BEM) (e.g.,
Barbot
, 2019;
Kato
, 2016;
Liu
,
141
2013;
Luo et al.
, 2017;
Nakata et al.
, 2012;
Rice and Tse
, 1986;
Segall and Bradley
, 2012;
Tse
142
and Rice
, 1986) can efficiently simulate earthquake ruptures in problems with more complex fault
143
geometry, often with the approximation of inertia (i.e.,
quasi-dynamic
earthquakes). Codes based
144
on the finite difference method (FDM) (e.g.,
Allison and Dunham
, 2018;
Erickson and Dunham
,
145
2014;
Erickson et al.
, 2017;
Herrendörfer et al.
, 2018;
Mckay et al.
, 2019;
Pranger
, 2020), finite
146
element method (FEM) (e.g.,
Liu et al.
, 2020;
Luo et al.
, 2020;
Tal and Hager
, 2018), and spectral
147
element method (SEM) (e.g.,
Kaneko et al.
, 2011;
Thakur et al.
, 2020) can flexibly incorporate
148
geometrical and structural complexity in earthquake simulations, usually at a greater computational
149
cost than BEM. For all these codes, common challenges lie in the interaction between the highly
150
nonlinear nature of the SEAS problems and numerical round-off errors, which can lead to the divergence
151
of model behaviors with increasing simulated time (
Lambert and Lapusta
, 2021). Simulation techniques
152
are further complicated when additional physical factors, e.g., fault roughness, material heterogeneities,
153
and bulk inelastic responses, are incorporated or approximated (e.g.,
Abdelmeguid et al.
, 2019;
154
Dal Zilio et al.
, 2022;
Romanet and Ozawa
, 2021). However, considering such complexity may
155
be crucial in our efforts to understand earthquakes and predict seismic hazards.
156
This study represents ongoing community efforts in the SEAS working group, supported by
157
the Southern California Earthquake Center (SCEC) to perform code verification exercises for SEAS
158
models. We reported the community initiative and results from our first two benchmarks, BP1-QD
159
and BP2-QD, for two-dimensional (2D) SEAS problems in
Erickson et al.
(2020). We gather 11
160
independent modeling groups using different numerical codes to participate and compare 2D SEAS
161
simulations. Through code comparisons, we identify how various computational factors, such
162
as the numerical resolution, domain size, and boundary conditions, influence simulation results
163
in 2D antiplane problems. Our exercises demonstrated excellent agreement in simulations with
164
a sufficiently small grid spacing and large domain size, lending confidence to the participating
165
numerical codes. We also found that artificial complexity in earthquake patterns can arise due
166
to insufficient numerical resolution for key physical length scales, although ensemble-averaged
167
measures, such as earthquake recurrence times, are more robust than observables from individual
168
simulations, even at poor numerical resolutions.
169
As our community and code capabilities grow, we have made substantial progress in benchmark
170
efforts for three-dimensional (3D) SEAS problems. Here, we present our recent development of
171
two new 3D benchmarks, BP4 and BP5. The dramatically increased computational demand for 3D
172
problems requires us to balance the simplicity and realism of the benchmark problems (Section 2).
173
Although we present the complete benchmark descriptions that include both fully dynamic (FD;
174
including inertia) and quasi-dynamic (QD; approximating inertia) formulations of earthquake ruptures,
175
–5–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
our code comparison results are limited to the quasi-dynamic problems. We examine choices of
176
numerical implementations among the modeling groups to ensure consistent comparisons of a large
177
set of 3D simulations (Section 3). We also design new strategies and metrics for code verification
178
for complex 3D simulations that are often done at the upper limit of numerical resolutions (Section 4).
179
In particular, we explore the sensitivity of diverse model outputs and observables to major computational
180
and physical factors. Through these efforts, we aim to improve and promote a new generation of
181
rigorous, robust numerical codes for SEAS problems, and to inform and interact with other communities
182
that are tackling similar computational challenges in nonlinear, multiscale, multi-physics problems
183
(e.g.
Buiter et al.
, 2016;
Matsui et al.
, 2016;
Maxwell et al.
, 2014;
Nearing et al.
, 2018).
184
2 Community Benchmark Development
185
2.1 Strategy for Benchmark Design
186
We follow the principle of starting simple and incrementally adding complexity in the design
187
process of SEAS benchmarks. For 2D benchmark problems (BP1-QD and BP2-QD), a 1D fault in
188
a 2D antiplane setting was considered to explore how the computational domain size and boundary
189
conditions affect simulation results and how numerical resolution (grid spacing or cell size) influences
190
earthquake patterns and statistics (
Erickson et al.
, 2020). Overall, we aim to verify different numerical
191
codes through a detailed comparison of simulated fault behavior over multiple time scales. These
192
efforts require a better understanding of the dependence of fault slip history on fault properties,
193
friction laws, initial conditions, model spin-up, and other factors.
194
Our findings and experience from 2D benchmark exercises prepare us for more complicated
195
3D benchmark problems. We need to design 3D benchmarks that are tractable for the widest suite
196
of numerical codes and thereby maximize participation of modelers, especially considering the
197
higher computational cost of 3D simulations and distinct capabilities of different codes in the community.
198
For example, codes based on the spectral boundary element method, e.g., BICyclE (
Lapusta and
199
Liu
, 2009), are efficient in solving for quasi-dynamic or fully dynamic earthquake ruptures, but
200
rely on periodic boundary conditions and free surface approximations. Methods based on the finite
201
element method, e.g., EQsimu (
Liu et al.
, 2020), can incorporate more complicated fault geometries
202
and bulk, including a rigorous treatment of the free surface, but need to balance the domain size
203
with a reasonable computational cost.
204
While we can in principle compare the full spectrum of fault behavior in SEAS models, the
205
focus of our exercise here is on reproducing earthquake nucleation, rupture, and recurrence. With
206
the computational cost in mind, we design benchmark problems where a direct comparison of
207
individual earthquakes is feasible (hence a consistent nucleation location is desirable). We then
208
assess the agreement of important model observables and their sensitivity to computational and
209
–6–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
physical factors. A better understanding of the roles of various inputs and outputs in SEAS models
210
will guide us in developing more complicated benchmarks and validating SEAS models in future.
211
Since the participation of many modelers is essential to the success of the code verification
212
exercise, we seek to build a consensus in the community at the outset of our benchmark design
213
process. We conducted surveys among the interested modelers to decide on the most preferred
214
benchmark problems. For instance, we have chosen to focus on quasi-dynamic problems for our
215
initial 3D benchmarks, BP4 and BP5, given that many numerical codes cannot yet incorporate
216
full inertial effects but adopt the radiation damping approximation (
Rice
, 1993). While we assess a
217
myriad of simulation outputs and develop metrics for model comparisons, we are flexible about the
218
submitted simulation data, given that sometimes substantial code development is needed. During
219
the subsequent development following initial comparisons of benchmark BP4, we learned lessons
220
about the computational cost and have accordingly revised the model parameters and output types
221
for benchmark BP5, hence some minor differences exist between the two benchmarks.
222
2.2 Benchmark Problem Setup
223
We have developed two benchmarks, BP4 and BP5, for 3D SEAS simulations (Figure 2). Our
224
first 3D benchmark problem, BP4, considers a 3D homogeneous, isotropic, linear elastic whole
225
space in
R
3
, defined by
=
¹
1
–푥
2
–푥
3
º 2 ¹1
–
3
, where
1
,
2
, and
3
refer to the coordinates in
226
the fault-normal, along-strike, and along-dip directions, respectively. A vertical strike-slip fault
227
is embedded at
1
=
0
. We use the notation “
̧
” and “
” to refer to the side of the fault with
1
228
positive and negative, respectively. We assume 3D motion, denoting components of the displacement
229
vector
as
=
¹
–푡
º
–푖
=
1
–
2
–
3
, in the
-direction. The second 3D benchmark problem, BP5,
230
involves a fault with half the vertical dimension in a 3D half-space, defined by
=
¹
1
–푥
2
–푥
3
º 2 ¹1
–
1º¹1
–
1º¹
0
–
,
231
with a free surface at
3
=
0
and
3
as positive downward. Several model parameters in BP5 are
232
adjusted to allow for reduced computational demand compared with BP4.
233
Each benchmark problem branches into two versions, depending on the treatment of the inertial
234
effect, i.e., quasi-dynamic (QD) or fully dynamic (FD) earthquake ruptures, which are assigned
235
with different suffixes in benchmark names (e.g., BP4-QD or BP4-FD). Full descriptions of these
236
benchmarks are available online on the SEAS code comparison platform (
https://strike.
237
scec.org/cvws/seas/
) and also included as supplementary materials. We summarize below
238
the governing equations, constitutive laws, and initial and interface conditions that are important for
239
understanding SEAS simulations for both QD and FD problems, and related numerical resolution
240
issues. For consistency and clarity, we have changed a few notations from the original benchmark
241
descriptions.
242
–7–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
The 3D fault zone motion is governed by the momentum balance equation, or the equilibrium
243
equation if inertia is neglected:
244
2
휕푡
2
=
r
for FD problems;
(1a)
245
0
=
r
for QD problems,
(1b)
246
247
where
is the displacement vector,
is the stress tensor, and
is the material density. Hooke’s
248
law relates the stress tensor
to strain tensor
by
249
푖푗
=
퐾휖
푘푘
푖푗
̧
2

푖푗
1
3
푘푘
푖푗

– 푖– 푗
=
1
–
2
–
3
–
(2)
250
where
and
are the bulk and shear moduli, respectively, and the use of subscript
follows the
251
Einstein summation convention. The strain-displacement relations are given by
252
푖푗
=
1
2

휕푢
휕푥
̧
휕푢
휕푥

– 푖– 푗
=
1
–
2
–
3
•
(3)
253
2.2.1 Boundary and Interface Conditions
254
We have a boundary condition at the surface (
3
=
0
) (for only BP5) and an interface condition
255
on the fault (
1
=
0
). At the free surface, all components of the traction vector are zeros, namely
256
3
¹
1
–푥
2
–
0
–푡
º
=
0
– 푗
=
1
–
2
–
3
•
(4)
257
Since the fault is always under compression in these benchmarks, there is no opening on the fault,
258
namely:
259
1
¹
0
̧
–푥
2
–푥
3
–푡
º
=
1
¹
0
–푥
2
–푥
3
–푡
º
•
(5)
260
We define the slip vector as the jump in horizontal and vertical displacements across the fault:
261
¹
2
–푥
3
–푡
º
=
¹
0
̧
–푥
2
–푥
3
–푡
º
¹
0
–푥
2
–푥
3
–푡
º
– 푗
=
2
–
3
–
(6)
262
with right-lateral motion yielding positive values of
2
. Positive values of
3
and
2
occur when the
263
̧
” or “
” side of fault moves in the positive or negative
3
and
2
directions, respectively.
264
We require that components of the traction vector be equal across the fault, which yields the
265
following conditions:
266
1
¹
0
̧
–푥
2
–푥
3
–푡
º
=
1
¹
0
–푥
2
–푥
3
–푡
º
– 푗
=
1
–
2
–
3
–
(7)
267
and denote the common values
11
,
21
, and
31
by
n
(positive in compression),
, and
,
268
respectively, i.e. one normal traction component and two shear traction components. Note that
269
positive values of
indicate stress that drives right-lateral faulting and positive values of
indicate
270
stress that tends to cause the “
̧
” side of the fault to move downward in the positive
3
direction and
271
the “
” side to move upward.
272
–8–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
We define the slip rate vector
in terms of its components,
=
¹
2
–푉
3
º
=
¹¤
2
–
¤
3
º
, where the
273
dot notation indicates the time derivative, and denote slip rate amplitude as the norm of the slip rate
274
vector,
=
jj
jj
. The shear stress vector is given by
=
¹
–휏
º
.
275
In both benchmark problems, we assign a frictional domain on the fault,
Ω
, with dimensions
276
of
¹
f
–푊
f
º
in the along-strike and along-dip directions, where fault slip is governed by a rate- and
277
state-dependent friction law (
Dieterich
, 1979;
Ruina
, 1983;
Marone
, 1998). The shear stress on the
278
frictional fault
is set to always equal the frictional strength
=
¹
2
– 퐹
3
º
, namely
279
=
¹
̄
n
–
–휃
º
–
(8)
280
where the effective normal stress is
̄
n
=
n
, with normal stress
n
and pore pressure
, and
281
is a state variable.
282
For quasi-dynamic problems (BP4-QD and BP5-QD),
=
0
̧
Δ
is the sum of the
283
prestress
0
, the shear stress change due to quasi-static deformation
Δ
, and the radiation damping
284
approximation of inertia
(
Rice
, 1993), where
=

2
s
is half the shear-wave impedance for
285
shear wave speed
s
=
√︁

, with the shear modulus
and density
. For fully dynamic problems,
286
=
0
̧
Δ
, where
Δ
includes all elastodynamic stress transfers due to prior slip on the fault.
287
The frictional resistance of the fault is the product of the effective normal stress,
̄
n
, and evolving
288
coefficient of friction,
, on the fault, namely
289
¹
̄
n
–
–휃
º
=
̄
n
¹
푉–휃
º

푉•
(9)
290
The effective normal stress is taken to be uniform in space and unvarying in time, which is valid
291
due to the symmetry across the planar fault and no fault opening. Since only the effective normal
292
stress, not the normal stress, matters in Eq. 9, we use
n
as a simpler notation for the effective
293
normal stress in the remainder of this paper. We adopt a regularized formulation for the rate-and-state
294
friction coefficient (
Lapusta et al.
, 2000)
295
¹
푉–휃
º
=

arcsinh

2

exp


̧
ln
¹


RS
º

–
(10)
296
where
RS
is the characteristic state evolution distance,

is the reference friction coefficient
297
determined at the reference slip rate

, and
and
are the parameters for the direct and evolution
298
effects, respectively. We couple Eq. 10 with the aging law for the evolution of the state variable
299
(
Dieterich
, 1979;
Ruina
, 1983):
300
푑휃
푑푡
=
1
푉휃
RS
–
(11)
301
The spatial distributions of parameters
and
are chosen to create a seismogenic zone with velocity-weakening
302
(VW;
푏 Ÿ
0
) frictional properties that is surrounded by regions with velocity-strengthening
303
(VS;
푏 ¡
0
) frictional properties, with a linear transition zone in-between. We use the same
304
–9–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
value for parameter
throughout the rate-and-state fault (denoted as
0
) and different values for
305
parameter
in the VW and VS regions (denoted as
0
and
max
, respectively).
306
Outside the frictional domain
Ω
, we impose a fixed long-term fault slip rate, which we refer
307
to as the plate loading rate
L
, giving rise to the interface conditions:
308
2
¹
2
–푥
3
–푡
º
=
L
–
(12a)
309
3
¹
2
–푥
3
–푡
º
=
0
–
(12b)
310
311
At an infinite distance from the fault (
j
1
j !1
), the far-field displacements should follow:
312

2
=

L
2
–
(13a)
313
1
=
3
=
0
–
(13b)
314
315
where the superscript “

" refers to the “+/–” sides of the fault, associated with positive and negative
316
displacement values, respectively. By imposing this boundary condition, we consider displacements
317
that are only caused by slip, excluding the deformation that produced the prestress
0
in the absence
318
of fault slip. As a result,
are essentially stress changes associated with the displacement field
319
relative to the prestress state. For the fully dynamic problem, Eq. 13 must be augmented with
320
radiation conditions that permit outgoing seismic waves (e.g.,
Bonnet
, 1999). We describe an infinitely
321
large domain in our benchmarks and leave choices of numerical implementation and approximation
322
to modelers (see Section 3.1).
323
2.2.2 Initial Conditions
324
We choose the initial values of the stress and state on the fault to enable a spatially uniform
325
distribution of initial fault slip rates, given by
326
=
¹
init
–푉
tiny
º
–
(14)
327
where we assign
init
=
L
for simplicity and
tiny
=
10
20
m/s to avoid infinity in logarithmic slip
328
rates. To achieve this, we prescribe the initial state over the entire fault with the steady-state value
329
at the slip rate
init
, namely
330
¹
2
–푥
3
–
0
º
=
RS

init
•
(15)
331
Accordingly, the initial stress vector takes the form
0
=
0

, where the scalar pre-stress
0
is
332
the steady-state stress:
333
0
=
푎휎
n

arcsinh

init
2

exp


̧
ln
¹


init
º

̧
휂푉
init
•
(16)
334
335
For quasi-dynamic problems, we need to specify an initial value for slip, which we take to be zero,
336
namely
337
¹
2
–푥
3
–
0
º
=
0
– 푗
=
2
–
3
•
(17)
338
–10–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
For fully dynamic problems, initial values for displacements and velocities in the medium need
339
to be specified. We spare the details here since our code comparisons below will be limited to
340
quasi-dynamic problems BP4-QD and BP5-QD.
341
To break the lateral symmetry of the fault and facilitate code comparisons, we add a square
342
zone within the VW region, with a width of
=
12
km and a center at (
22.5 km,
7.5 km) in BP4
343
and (
24 km,
10 km) in BP5, as a prescribed nucleation location for the first simulated earthquake.
344
To do that, we impose a higher initial slip rate,
i
, in the
2
direction within this square zone at
345
=
0
, while keeping the initial state variable
¹
2
–푥
3
–
0
º
unchanged. The resultant higher pre-stress
346
is calculated by replacing
init
with
i
in Eq. 16. This initial condition leads to an immediate initiation
347
of the first event. In BP5, we additionally use a smaller characteristic state evolution distance
RS
348
in this prescribed nucleation zone to promote the nucleation of subsequent earthquakes in the same
349
areas (see the next section). We note that future benchmarks can use a spatially smoother function
350
of the physical properties within the nucleation zone to minimize the influence of spatial discretizations
351
in numerical models (
Galis et al.
, 2015).
352
In simulations, the governing equations, Eqs. 1–3, are solved along with interface conditions,
353
Eq. 4 (for only BP5) and Eqs. 5–13, and initial conditions, Eqs. 14–17, over the period
0
6
6
f
,
354
where
f
is the maximum simulated time. Numerical methods that truncate model domain in the
355
fault-normal direction also need to explicitly incorporate the far-field boundary conditions on asymptotic
356
behavior of displacements at infinity (see Section 3.1). All model parameters in benchmarks BP4-QD
357
and BP5-QD are listed and compared in Table 1.
358
2.2.3 Critical Physical Length Scales
359
Numerical resolution is a critical issue for 3D benchmark problems, as we need to balance the
360
computational cost and adequate resolution to achieve acceptable model agreement. Two physical
361
length scales are generally important to consider in these problems. The first length scale, often
362
referred to as the process zone or cohesive zone,
Λ
, describes the spatial region near the rupture
363
front under which breakdown of fault resistance occurs, and shrinks as ruptures propagate faster
364
(
Freund
, 1990;
Palmer and Rice
, 1973). For faults governed by the rate-and-state friction, the quasi-static
365
process zone at a rupture speed of
0
̧
,
Λ
0
, can be estimated as follows (
Day et al.
, 2005;
Lapusta
366
and Liu
, 2009):
367
Λ
0
=
휇퐷
RS
푏휎
n
–
(18)
368
where
is a constant of order 1.
369
The second length scale that controls model behavior is the nucleation size

, which determines
370
the minimum size of the velocity-weakening region over which spontaneous nucleation may occur
371
(
Ampuero and Rubin
, 2008;
Rice and Ruina
, 1983;
Rubin and Ampuero
, 2005). For 3D problems,
372
–11–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.
manuscript submitted to
JGR: Solid Earth
the nucleation size can be estimated for the aging law for
0
•
5
Ÿ 푎

푏 Ÿ
1
as follows (
Chen and
373
Lapusta
, 2009):
374

=
2
휇푏퐷
RS
¹
º
2
n
•
(19)
375
Using Eqs. 18 and 19, we estimate that the nucleation size is
12
•
4
km and 12.5 km within the VW
376
region (outside the zone of frictional heterogeneity) in BP4 and BP5, respectively, whereas the
377
process zone is 2 and 6 km, respectively. This allows us to suggest
500
m and 1000 m for the grid
378
spacing,
Δ
, in low-order accurate methods for BP4 and BP5, respectively, which resolve
Λ
0
with at
379
least four cells in both benchmarks, following suggestions by
Day et al.
(2005).
380
The two benchmark problems are designed to produce a periodic sequence of spontaneous
381
earthquakes and slow slip, following the first event in which we impose higher local slip rates to
382
kickstart the earthquake rupture. BP5 is slightly different from BP4 in that the characteristic state
383
evolution distance
RS
is reduced within a square zone within the VW region, resulting in a smaller
384
nucleation size,

=
11
•
6
km. This form of persistent frictional heterogeneity is introduced to
385
favor (but not always determine) the initiation of subsequent earthquakes at the same location. We
386
choose the total simulated time to produce up to eight large earthquakes in the simulations, which
387
allows us to examine not only a few early events but also the seismic behavior of the fault in the
388
longer term.
389
2.3 Model Outputs
390
To assess model behavior over disparate spatial and temporal scales, we design several types
391
of simulation outputs for these benchmarks (Figure 3): (1) time series of local on-fault and off-fault
392
properties, (2) time series of global source properties, (3) a catalog of earthquake characteristics,
393
(4) profiles of slip accumulation and stress evolution, and (5) rupture times during the first event in
394
the sequence. The output formats for coseismic observables follow the practice in the code verification
395
of single-event dynamic rupture simulations (
Harris et al.
, 2009).
396
For local time series data, we are interested in resolving the time evolution of fault slip rates,
397
shear stress, and off-fault displacements throughout the coseismic, postseismic, and interseismic
398
periods. The global source properties refer to the evolving maximum slip rates and moment rates
399
over the entire seismogenic fault areas, which are useful for determining the precise time of initiation
400
and cessation of individual earthquakes. The catalog data contain key characteristics of simulated
401
earthquakes, including their initiation and termination times, coseismic slip, and static stress drop.
402
The beginning and end of the coseismic period are determined as the times at which any point
403
on the fault reaches above or all points drop below a threshold slip rate,
th
(chosen as
0
•
03
m/s),
404
respectively. We then estimate coseismic slip and stress drop as the change in the amplitude of fault
405
slip and shear stress (negative stress change corresponds to positive stress drop).
406
–12–
ESSOAr | https://doi.org/10.1002/essoar.10508582.2 | CC_BY_NC_ND_4.0 | First posted online: Thu, 3 Feb 2022 22:26:58 | This content has not been peer reviewed.