This manuscript is a preprint
and has been submitted for publication in
Seismological Research Letters
. Please note that the manuscript has yet to be
peer-reviewed and has not been formally accepted for publication. Subsequent versions of
this manuscript may have slightly different content. If accepted, the final version of this
manuscript will be available via the
’Peer-reviewed Publication DOI’
link on the right-hand
side of this webpage. Please feel free to contact any of the authors; we welcome feedback.
1
The Community Code Verification Exercise for
1
Simulating Sequences of Earthquakes and Aseismic Slip
2
(SEAS)
3
Brittany A. Erickson
1
*
, Junle Jiang
2
*
, Michael Barall
3
, Nadia Lapusta
4
,
Eric M. Dunham
5
, Ruth Harris
6
, Lauren S. Abrahams
5
, Kali L. Allison
7
,
Jean-Paul Ampuero
8
, Sylvain Barbot
9
, Camilla Cattania
5
, Ahmed Elbanna
10
,
Yuri Fialko
11
, Benjamin Idini
4
, Jeremy E. Kozdon
12
, Valère Lambert
4
,
Yajing Liu
13
, Yingdi Luo
4
, Xiao Ma
10
, Maricela Best McKay
14
,
Paul Segall
5
, Pengcheng Shi
15
, Martijn van den Ende
8
, Meng Wei
15
1
University of Oregon, 1477 E. 13th Ave., Eugene , OR 97403-1202
2
Cornell University
3
Invisible Software
4
California Institute of Technology
5
Stanford University
6
United States Geological Survey
7
University of Maryland
8
Université Côte d’Azur, IRD, CNRS,
Observatoire de la Côte d’Azur, Géoazur, France
9
University of Southern California
10
University of Illinois Urbana-Champaign
11
University of California San Diego
12
Naval Postgraduate School
13
McGill University
14
Portland State University
15
University of Rhode Island
4
September 4, 2019
5
*
Correspondence to B.A.E (bae@uoregon.edu) and J.J. (jjiang@cornell.edu)
2
Any use of trade, firm, or product names is for descriptive purposes only and does not
6
imply endorsement by the U.S. Government.
7
3
Abstract
8
Numerical simulations of Sequences of Earthquakes and Aseismic Slip (SEAS) have
9
made great progress over the past decades to address important questions in earthquake
10
physics and fault mechanics. However, significant challenges in SEAS modeling remain
11
in resolving multiscale interactions between aseismic fault slip, earthquake nucleation,
12
and dynamic rupture; and understanding physical factors controlling observables such
13
as seismicity and ground deformation. The increasing capability and complexity of
14
SEAS modeling calls for extensive efforts to verify codes and advance these simulations
15
with rigor, reproducibility, and broadened impact. In 2018, we initiated a community
16
code-verification exercise for SEAS simulations, supported by the Southern California
17
Earthquake Center (SCEC). Here we report the findings from our first two benchmark
18
problems (BP1 and BP2), designed to test the capabilities of different computational
19
methods in correctly solving a mathematically well-defined, basic problem in crustal
20
faulting. These benchmarks are for a 2D antiplane problem, with a 1D planar vertical
21
strike-slip fault obeying rate-and-state friction, embedded in a 2D homogeneous, linear
22
elastic half-space. Sequences of quasi-dynamic earthquakes with periodic occurrences
23
(BP1) or bimodal sizes (BP2) and their interactions with aseismic slip are simulated.
24
The comparison of >70 simulation results from 11 groups using different numerical
25
methods, uploaded to our online platform, show excellent agreements in long-term and
26
coseismic evolution of fault properties. In BP1, we found that the truncated domain
27
boundaries influence interseismic fault stressing, earthquake recurrence, and coseismic
28
rupture process, and that agreement between models is only achieved with sufficiently
29
large domain sizes. In BP2, we found that complexity of long-term fault behavior
30
depends on how well important physical length scales related to spontaneous nucleation
31
and rupture propagation are resolved. Poor numerical resolution can result in the
32
generation of artificial complexity, impacting simulation results that are of potential
33
interest for characterizing seismic hazard, such as earthquake size distributions, moment
34
release, and earthquake recurrence times. These results inform the development of more
35
advanced SEAS models, contributing to our further understanding of earthquake system
36
4
dynamics.
37
5
Introduction and Motivation
38
When we develop models of physical systems, credible and reproducible results are
39
essential to scientific progress. Robust predictive models of earthquake source pro-
40
cesses have become important means for studying fundamental questions in earthquake
41
science. Models of single earthquakes (known as
dynamic rupture simulations
) have
42
emerged as powerful tools for understanding the influence of fault geometry, friction
43
and prestress on rupture propagation, and for explaining observations of high-frequency
44
ground motions and damage zones (
Day
, 1982;
Olsen et al.
, 1997;
Nielsen et al.
, 2000;
45
Duan and Oglesby
, 2006;
Ripperger et al.
, 2007;
Bhat et al.
, 2007;
Dunham et al.
,
46
2011a,b;
Lozos et al.
, 2011;
Gabriel et al.
, 2012;
Shi and Day
, 2013;
Kozdon and Dun-
47
ham
, 2013;
Xu et al.
, 2015;
Wollherr et al.
, 2018;
Ma and Elbanna
, 2019). Many of
48
the codes used for these studies incorporate advanced features such as 3D domains and
49
complex fault geometries, leading to very large problems for which rigorous convergence
50
tests can be too computationally expensive. An alternative means for verifying model
51
results are code comparisons made across the different modeling groups, using cell sizes
52
at the limit of computational feasibility. Over the past decade, the SCEC/USGS Spon-
53
taneous Rupture Code Verification Project has made significant progress in using code
54
comparison studies to provide confidence in model outcomes (
Harris et al.
, 2009;
Barall
55
and Harris
, 2015;
Harris et al.
, 2018).
56
Although these dynamic rupture simulations have contributed greatly to our un-
57
derstanding of the physical factors that govern ground motion, they are limited to
58
single-event scenarios with imposed artificial prestress conditions and
ad hoc
nucle-
59
ation procedures. In order to understand earthquake source processes and how fault
60
slip history influences subsequent events, it has been widely recognized that we need
61
models that simulate behavior over multiple seismic events and the intervening periods
62
of aseismic deformation. To address this need, models of Sequences of Earthquakes
63
and Aseismic Slip (SEAS) have emerged that consider all phases of earthquake fault-
64
ing, from slow tectonic loading to earthquake nucleation (under self-consistent prestress
65
6
conditions), propagation and termination. However, so far codes for SEAS simulations
66
remain untested. Inspired by the success of the SCEC/USGS Spontaneous Rupture
67
Code Verification Project, this paper describes the efforts of the SEAS initiative – a
68
SCEC (Southern California Earthquake Center) funded working group who has initi-
69
ated the first code-verification study for earthquake sequence simulations. In this pa-
70
per, we present the initial benchmark problems and results from the code comparisons
71
submitted to our online platform (
http://scecdata.usc.edu/cvws/seas/
). Through
72
these exercises, we aim to provide confidence in SEAS model outcomes, determine best
73
practices for improvement of accuracy and efficiency of SEAS simulations, and provide
74
other scientists strategies for verification during code development.
75
In SEAS models the goal is to capture the interplay of interseismic periods and
76
the associated aseismic fault slip that ultimately lead to earthquake nucleation and
77
earthquakes (dynamic rupture events) themselves, in an effort to understand which
78
physical factors control the full range of observables such as aseismic deformation,
79
nucleation locations of earthquakes, ground shaking during dynamic rupture, recurrence
80
times and magnitudes of major earthquakes, see Figure 1. These features distinguish
81
SEAS models from both dynamic rupture models which only consider single events,
82
and the so-called earthquake simulators (
Tullis et al.
, 2012). Earthquake simulators are
83
capable of simulating seismicity patterns over millennium time scales in complex fault
84
network systems (
Richards-Dinger and Dieterich
, 2012) but are missing key physical
85
features that could potentially dominate earthquake and fault interaction, such as stress
86
transfer generated by dynamic waves, aseismic slip within fault segments, and inelastic
87
responses.
88
SEAS modeling is not without significant challenges, due to the varying tempo-
89
ral and spatial scales that characterize earthquake source behavior. For computational
90
efficiency the vast majority of SEAS models do not consider full dynamics during earth-
91
quake rupture, but rather take a "quasi-dynamic" approach, where inertia is only ap-
92
proximated (see section for further details). Computations are further complicated
93
when material heterogeneities, bulk inelastic responses and fault nonplanarity are in-
94
7
cluded. However, accounting for such complexity is widely recognized as crucial for
95
understanding the real Earth and predicting seismic hazards. Significant developments
96
in SEAS models over the past decade have incorporated some of these complexities and
97
connected model outcomes to geophysical observations. For example, seismological and
98
geodetic observations have been combined with modeling of coseismic and quasi-static
99
(aseismic) deformation to infer the spatial distribution of fault frictional properties
100
(
Johnson et al.
, 2006;
Barbot et al.
, 2009;
Mitsui and Hirahara
, 2011;
Dublanchet et al.
,
101
2013;
Floyd et al.
, 2016;
Jiang and Fialko
, 2016), the decay rate of aftershocks (
Per-
102
fettini and Avouac
, 2004, 2007), the role of tremor and slow slip (
Mele Veedu and
103
Barbot
, 2016;
Dublanchet
, 2017;
Luo and Ampuero
, 2017), and long-term models have
104
been used to reproduce characteristics of multiple and/or repeating events (
Chen and
105
Lapusta
, 2009;
Barbot et al.
, 2012). The framework of earthquake cycle modeling is
106
also adopted to explain geodetic and geologic data (
Meade et al.
, 2013;
Kaneko et al.
,
107
2011;
Wei et al.
, 2013, 2018), study subduction zones (
Hori et al.
, 2004;
van Dinther
108
et al.
, 2013;
Noda and Lapusta
, 2013;
Liu and Rice
, 2005, 2007;
Li and Liu
, 2016,
109
2017), collision zones (
Qiu et al.
, 2016;
Michel et al.
, 2017), and explore induced seis-
110
micity phenomena (
McClure and Horne
, 2011;
Dieterich et al.
, 2015), among many
111
applications.
112
While SEAS models are being used to explain, reproduce, and predict earthquake
113
behavior and other geophysical phenomena, a critical step must be to ensure that these
114
methodologies are accurate. The SEAS initiative is also taking the step to improve
115
and promote a new generation of verified numerical SEAS models that can simulate
116
much longer periods of earthquake activity than single-event dynamic rupture simula-
117
tions but with the same level of computational rigor, while incorporating qualitatively
118
different features such as (a) pre-, inter-, and post-seismic slip and the resulting stress
119
redistribution, (b) spontaneous earthquake nucleation, and (c) physical processes rele-
120
vant to long-term slip such as interseismic healing of the fault zone, viscoelasticity, and
121
fluid flow. Such SEAS models can provide physics-based approximations for larger-scale
122
and longer-term earthquake simulators. In addition they can inform the initial condi-
123
8
tions and nucleation procedures for dynamic rupture simulations, however our vision
124
for SEAS models is to develop them all to include full dynamic ruptures, capturing the
125
range of processes and heterogeneities known to be essential for realistic ground motion
126
modeling.
127
SEAS Modeling Challenges and Initial Benchmark
128
Problems
129
Although the ultimate SEAS modeling framework would naturally include dynamic rup-
130
ture modeling, current methods for simulating SEAS problems require computational
131
codes that are fundamentally different from those used in single-event dynamic rupture
132
simulations. The use of variable time stepping and possible switching between different
133
computational schemes is required in order to resolve sub-seconds to year-long changes.
134
The interaction between the highly nonlinear nature of the problems and round-off er-
135
rors can lead to model divergence. The need to distinguish between legitimate solution
136
differences due and improper choices of algorithm and modeling procedures necessitates
137
new and more suitable comparison metrics.
138
SEAS models are unique in that they cover a wide range of numerical methodologies
139
and applications in earthquake science. Methods based on spectral boundary integral
140
formulations (BIEM) are efficient in solving for earthquake ruptures with quasi-dynamic
141
or full inertial effects (
Lapusta and Rice
, 2003;
Lapusta and Liu
, 2009;
Jiang and La-
142
pusta
, 2016). Methods based on the finite difference method (FDM) or a hybrid finite
143
element/spectral BIEM have been used to simulate quasi-dynamic ruptures on faults
144
with more complex bulk rheologies (
Erickson and Dunham
, 2014;
Erickson et al.
, 2017;
145
Allison and Dunham
, 2018;
Mckay et al.
, 2019;
Abdelmeguid et al.
, 2019). Other SEAS
146
modeling approaches include boundary element methods (BEM) for simulating slow slip
147
and tremor (e.g.,
Tse and Rice
, 1986;
Rice and Tse
, 1986;
Ong et al.
, 2019;
Goswami
148
and Barbot
, 2018;
Luo and Ampuero
, 2011;
Nakata et al.
, 2012;
Liu
, 2013;
Wei et al.
,
149
9