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
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 “