1.
Introduction
Earthquakes occur in the context of fault networks and many large earthquakes span several fault segments.
This reality brings about the issue of fault interaction and highlights the need for simulating earthquake
source processes over several fault segments and regional-scale fault networks. How dynamic ruptures
navigate fault segmentation has strong implications for seismic hazard analysis (Biasi & Wesnousky,
2021
;
Field,
2019
). Earthquakes are capable of jumping fault segments. For example, the 1992 Landers earthquake
succeeded in rupturing across at least four fault segments, amounting to a Mw 7.3 event (Sieh et al.,
1993
).
Abstract
Physics-based numerical modeling of earthquake source processes strives to predict
quantities of interest for seismic hazard, such as the probability of an earthquake rupture jumping
between fault segments. How to assess the predictive power of numerical models remains a topic of
ongoing debate. Here, we investigate how sensitive the outcomes of numerical simulations of sequences of
earthquakes and aseismic slip are to choices in numerical discretization and treatment of inertial effects,
using a simplified 2-D crustal fault model with two co-planar segments separated by a creeping barrier.
Our simulations demonstrate that simplifying inertial effects and using oversized cells significantly affects
the resulting earthquake sequences, including the rate of two-segment ruptures. We find that fault models
with different properties and modeling assumptions can produce similar frequency-magnitude statistics
and static stress drops but have different rates of two-segment ruptures. For sufficiently long faults, we
find that long-term sequences of events can substantially differ even among simulations that are well
resolved by standard considerations. In such simulations, some outcomes, such as static stress drops, are
similar among adequately resolved simulations, whereas others, such as the rate of two-segment ruptures,
can be highly sensitive to numerical procedures and modeling assumptions. While it is possible that the
response of models with additional ingredients -Realistic fault geometry, fluid effects, etc. -Would be less
sensitive to numerical procedures, our results emphasize the need to examine the potential dependence
of simulation outcomes on the modeling procedures and resolution, particularly when assessing their
predictive value for seismic hazard assessment.
Plain Language Summary
There is growing interest in using computer simulations of
long-term earthquake sequences to determine quantities of interest for seismic hazard, such as the
probability of an earthquake rupture jumping from one fault segment to another. This is because large
earthquakes are rare, hence the need to assess potential future earthquake scenarios through numerical
modeling based on all available field observations and knowledge of fault physics. However, the outcomes
of numerical simulations can depend on choices in modeling approximations and numerical procedures.
Here, we numerically simulate earthquake sequences in a model with two fault segments separated by a
barrier and study how the resulting earthquake sequences depend on common modeling choices. We find
that different treatment of the inertial effects and small changes in physical and numerical parameters
can result in different simulated long-term sequences, including significant changes in the rate of multi-
segment ruptures. This is true even when certain properties of the earthquake sequences are similar,
such as the earthquake frequency-magnitude statistics and the average stress drop. Our results emphasize
the need to examine how simulation outcomes may depend on modeling choices when assessing their
predictive value and explore the sensitivity of hazard parameters to model uncertainty.
LAMBERT AND LAPUSTA
© 2021. The Authors.
This is an open access article under
the terms of the
Creative Commons
Attribution
License, which permits use,
distribution and reproduction in any
medium, provided the original work is
properly cited.
Resolving Simulated Sequences of Earthquakes and Fault
Interactions: Implications for Physics-Based Seismic
Hazard Assessment
Valère Lambert
1
and Nadia Lapusta
1,2
1
Seismological Laboratory, California Institute of Technology, Pasadena, CA, USA,
2
Department of Mechanical and
Civil Engineering, California Institute of Technology, Pasadena, CA, USA
Key Points:
•
Long-term interactions of fault
segments qualitatively differ
among fully versus quasi-dynamic
simulations, and when using
oversized cells
•
Simulations with similar stress drops
and power-law frequency-magnitude
statistics can have different rates of
multisegment ruptures
•
Simulated earthquake sequences can
differ due to compounded effects
of numerical errors, even when
individual ruptures are well resolved
Supporting Information:
Supporting Information may be found
in the online version of this article.
Correspondence to:
V. Lambert,
valerelambert@gmail.com
Citation:
Lambert, V., & Lapusta, N. (2021).
Resolving simulated sequences of
earthquakes and fault interactions:
Implications for physics-based
seismic hazard assessment.
Journal
of Geophysical Research: Solid Earth
,
126
, e2021JB022193.
https://doi.
org/10.1029/2021JB022193
Received 7 APR 2021
Accepted 10 OCT 2021
Author Contributions:
Conceptualization:
Valère Lambert,
Nadia Lapusta
Formal analysis:
Valère Lambert,
Nadia Lapusta
Writing – original draft:
Valère
Lambert, Nadia Lapusta
10.1029/2021JB022193
RESEARCH ARTICLE
1 of 33
Journal of Geophysical Research: Solid Earth
LAMBERT AND LAPUSTA
10.1029/2021JB022193
2 of 33
The 2016 Mw 7.8 Kaikoura earthquake ruptured at least 21 segments of the Marlborough fault system
(Ulrich et al.,
2019
). Increasingly, seismological observations show that it is not uncommon to see ruptures
navigating and triggering subsequent ruptures within fault networks, including the recent 2019 Mw 6.4 and
7.1 Ridgecrest earthquakes (Ross et al.,
2019
), and the 2012 Mw 8.6 and 8.2 Indian Ocean earthquakes (Wei
et al.,
2013
), the largest and second-largest recorded strike-slip earthquakes to date. Yet, in any given seis
-
mogenic region, the record of past large events is not long enough to forecast the behavior of ruptures with
respect to the existing fault segments, specifically how likely would the rupture be to jump between nearby
segments, prompting the discussion on whether and how physics-based models may inform this and other
questions important for seismic hazard assessment (Field,
2019
).
Determining what conditions allow a dynamic rupture to propagate or arrest are key to understanding
the maximum potential magnitude of an earthquake. Previous modeling of single fully dynamic ruptures
have shown great success in investigating earthquake propagation in nonplanar and multi-segment fault
models, including step-overs and branched geometries (Ando & Kaneko,
2018
; Douilly et al.,
2015
; Duan
& Oglesby,
2006
; Dunham et al.,
2011b
; Galvez et al.,
2014
; Harris & Day,
1993
,
1999
; Harris et al.,
1991
;
Hu et al.,
2016
; Kame et al.,
2003
; Lozos et al.,
2015
; Ulrich et al.,
2019
; Withers et al.,
2018
; Wollherr
et al.,
2019
). In particular, such modeling has shown that the ability of a rupture to propagate across seg
-
ments depends on the stresses before the rupture and shear resistance assumptions, as well as the geometry
of the fault system. However, single-rupture simulations need to select initial conditions and need addition
-
al assumptions to incorporate the effect of previous seismic and aseismic slip.
Fault processes involve both sequences of dynamic events and complex patterns of quasi-static slip. Simu
-
lating this behavior in its entirety is a fascinating scientific problem. However, even for the more pragmatic
goal of physics-based predictive modeling of destructive large dynamic events, it is still important to consid
-
er sequences of earthquakes and aseismic slip (SEAS), since prior slip events, including aseismic slip, may
determine where earthquakes would nucleate as well as modify stress and other initial conditions before
dynamic rupture. Furthermore, such simulations provide a framework for determining physical properties
consistent with a range of observations including geodetically recorded surface motions, microseismicity,
past (including paleoseismic) events, and thermal constraints, and hence may inform us about the current
state of a fault segment or system and potential future rupture scenarios (e.g., Allison & Dunham,
2018
;
Barbot et al.,
2012
; Ben-Zion & Rice,
1997
; Cattania,
2019
; Chen & Lapusta,
2009
; Erickson & Day,
2016
;
Erickson & Dunham,
2014
; Jiang & Lapusta,
2016
; Kaneko et al.,
2010
; Lambert & Barbot,
2016
; Lambert,
Lapusta, Perry,
2021
; Lapusta & Rice,
2003
; Lapusta et al.,
2000
; Lin & Lapusta,
2018
; Liu & Rice,
2005
; Noda
& Lapusta,
2013
; Perry et al.,
2020
; Segall et al.,
2010
). However, simulating long-term slip histories is quite
challenging because of the variety of temporal and spatial scales involved.
Recently, several earthquake simulators have been developed with the goal of simulating millions of earth
-
quake ruptures over regional fault networks for tens of thousands of years (Richards-Dinger & Dieter
-
ich,
2012
; Shaw et al.,
2018
; Tullis et al.,
2012
). The term “simulators” typically refers to approaches that em
-
ploy significant simplifications, compared to most SEAS simulations, in solution procedures and physical
processes, in order to simulate earthquake sequences on complex, regional scale 3-D fault networks for long
periods of time. For example, earthquake simulators typically account only for the quasi-static stress trans
-
fer due to earthquake events, ignoring wave-mediated stress changes, aseismic slip/deformation, and fluid
effects; employ approximate rule-based update schemes for earthquake progression instead of solutions of
the governing continuum mechanics equations; and use oversized numerical cells. Such simplifications
are currently necessary to permit simulations of hundreds of thousands of events over hundreds of fault
segments that comprise the regional networks (Shaw et al.,
2018
). Earthquake simulators have matched a
number of regional-scale statistical relations, including the Gutenberg-Richter frequency-magnitude scal
-
ing (Shaw et al.,
2018
), and highlighted the importance of large-scale fault and rupture interactions.
There is growing interest in using earthquake simulators to directly determine quantities of interest for
seismic hazard, such as the probability of an earthquake rupture jumping from one fault segment to anoth
-
er (Field,
2019
; Shaw et al.,
2018
). However, assessing the predictive power of numerical models remains
a topic of active research. Determining how sensitive simulated outcomes may be to modeling choices and
how reliably they can be determined from a given numerical model are topics of great importance for phys
-
ics-based hazard assessment. Here, we examine the sensitivity of the long-term interaction of fault segments
Journal of Geophysical Research: Solid Earth
LAMBERT AND LAPUSTA
10.1029/2021JB022193
3 of 33
to choices in numerical discretization and representations of inertial ef
-
fects in simulated sequences of earthquakes and aseismic slip, using a
relatively simple 2-D model of two co-planar strike-slip fault segments
separated by a velocity-strengthening (VS) barrier. We explore how con
-
siderations for adequate numerical resolution and convergence depend
on the physical assumptions and complexity of earthquake sequences as
well as on the modeling outcome of interest. We especially focus on the
rate of earthquake ruptures jumping across the VS barrier and examine
whether reproducing comparable earthquake frequency-magnitude sta
-
tistics and static stress drops provides sufficient predictive power for the
jump rate, a quantity of interest to seismic hazard studies (Field,
2019
).
2.
Model Setup and Numerical Resolution
Our simulations are conducted following the methodological devel
-
opments of Lapusta et al. (
2000
), Noda and Lapusta (
2010
), and Lam
-
bert, Lapusta, Perry (
2021
). We consider a one-dimensional (1-D) fault
embedded into a 2-D uniform, isotropic, elastic medium (Figure
1
). The
2-D model approximates a faulted crustal plate coupled to a moving sub
-
strate using the idea of a constrained continuum (Johnson,
1992
; Lehner
et al.,
1981
). Fault slip may vary spatially along-strike but it is depth-av
-
eraged through a prescribed seismogenic thickness
퐴퐴퐴퐴
푆푆
=15
km, beneath
which the elastic domain is coupled to a substrate moving at the pre
-
scribed loading rate (
퐴퐴퐴퐴
pl
=10
−9
m/s). The elastodynamic equation for the
depth-averaged displacement along-strike
퐴퐴퐴퐴퐴
(
푥푥푥푥푥푥푥푥
)
is given by (Kaneko &
Lapusta,
2008
; Lehner et al.,
1981
):
푍푍
휕휕
2
̄푢푢
휕휕휕휕
휕휕
2
̄푢푢
휕휕휕휕
휆휆
eff
(
1
휕휕
)
푉푉
pl
푡푡
−
̄푢푢
)
=
1
푐푐
푠푠
휕휕
2
̄푢푢
휕휕푡푡
2
,
(1)
where
퐴퐴퐴퐴
eff
=(
휋휋
∕4)
퐴퐴
S
and
퐴퐴퐴퐴
=1∕(1−
휈휈
)
, with
퐴퐴퐴퐴
being the Poisson's ratio. The effective wave speed along-
strike for the crustal plane model is
퐴퐴퐴퐴
퐿퐿
=
푍푍퐴퐴
푠푠
, where
퐴퐴퐴퐴
푠푠
is the shear wave speed. The along-strike slip is then
given by
퐴퐴퐴퐴
(
푥푥푥푥푥
)=
̄푢푢
(
푥푥푥푥푥
=0
+
푥푥푥
)−
̄푢푢
(
푥푥푥푥푥
−
푥푥푥
)
.
Our simulations resolve sequences of earthquakes and aseismic slip (SEAS) in their entirety, including the
gradual development of frictional instability and spontaneous nucleation, dynamic rupture propagation,
post-seismic slip that follows the event, and the interseismic period between events (Figure
2
). In all mod
-
els, frictional resistance along the fault interface is governed by the standard laboratory-derived rate-and-
state friction law with the state evolution described by the aging law (Dieterich,
1979
; Ruina,
1983
):
휏휏
̄휎휎휎휎
=(
휎휎
−
푝푝
)
[
휎휎
∗
+
푎푎
ln
푉푉
푉푉
∗
+
푏푏
ln
푉푉
∗
휃휃
퐷퐷
RS
]
,
(2)
푑푑푑푑
푑푑푑푑
푉푉푑푑
퐷퐷
RS
,
(3)
where
퐴퐴퐴퐴퐴
=(
퐴퐴
−
푝푝
)
is the effective normal stress,
퐴퐴퐴퐴
is the normal stress,
퐴퐴퐴퐴
is the pore pressure,
퐴퐴퐴퐴
is the shear
stress,
퐴퐴퐴퐴
is the friction coefficient,
퐴퐴퐴퐴
is the slip velocity,
퐴퐴퐴퐴
is the state variable,
퐴퐴퐴퐴
RS
is the characteristic slip
for the evolution of the state variable,
퐴퐴퐴퐴
∗
is the reference steady-state friction coefficient corresponding to
a reference slip rate
퐴퐴퐴퐴
∗
, and
퐴퐴퐴퐴
and
퐴퐴퐴퐴
are the direct and evolution effect constitutive parameters, respectively.
At steady-state (constant slip velocity), the shear stress and state variable evolve to their steady-state values
퐴퐴퐴퐴
푠푠푠푠
and
퐴퐴퐴퐴
푠푠푠푠
given by:
휏휏
푠푠푠푠
(
푉푉
)=(
휎휎
−
푝푝
)
[
푓푓
∗
+(
푎푎
−
푏푏
)ln
푉푉
푉푉
∗
]
,
(4)
휃휃
푠푠푠푠
(
푉푉
)=
퐷퐷
RS
푉푉
.
(5)
Figure 1.
Schematic of a strike-slip fault with two co-planar velocity-
weakening fault segments separated by a velocity-strengthening barrier.
In our simulations, we use a 2D approximation of the problem with a
1D along-strike depth-averaged fault, in which the fault is assumed to
be creeping at the loading plate rate
퐴퐴퐴퐴
pl
=10
−9
m/s below the depth of
퐴퐴퐴퐴
푆푆
=
15 km.
Journal of Geophysical Research: Solid Earth
LAMBERT AND LAPUSTA
10.1029/2021JB022193
4 of 33
The combination of frictional properties such that
퐴퐴
(
푎푎
−
푏푏
)
>
0
results in steady-state velocity-strengthening
(VS) behavior, where the shear resistance increases with an increase in slip velocity and where stable slip is
expected. If
퐴퐴
(
푎푎
−
푏푏
)
<
0
then the fault exhibits velocity-weakening (VW) behavior, in which case an increase
in slip velocity leads to a decrease in shear resistance, making these regions of the fault potentially seismo
-
genic if their size exceeds a critical nucleation size.
Figure 2.
Interaction of two co-planar fault segments in well-resolved simulations of model M1 demonstrating convergence of simulated earthquake
sequences. (a and b) History of cumulative slip over 4,000 years in well-resolved fully dynamic simulations of fault model M1 with initial conditions S1 using
(a) 12.5-m and (b) 25-m cell size. Contours for seismic slip are plotted every 0.5 s, with ruptures that jump across the VS barrier colored blue. The simulated
fault behavior is virtually indistinguishable between the two resolutions. (c) Frequency-magnitude histograms of events, on top of each other for the two
resolutions. The well-resolved simulations produce the same relatively simple and quasi-periodic behavior. (d and e) The evolution of local shear stress and
slip velocity at a point (
퐴퐴퐴퐴
=−20
.
5
km, shown by star in (a and b), practically indistinguishable even after over 3,800 years of simulated time. (f–h) Spatial
distribution of shear stress at the rupture front for three locations (
퐴퐴퐴퐴
=−20
km, 5 and 20 km) throughout the first rupture in (a and b). While the quasi-static
estimate of the cohesive zone
퐴퐴
Λ
0
is about 1.1 km, the actual size of the cohesive zone varies with the local rupture speed throughout the rupture. In these well-
resolved simulations, the cohesive zone is always resolved by at least 10 cells.
Journal of Geophysical Research: Solid Earth
LAMBERT AND LAPUSTA
10.1029/2021JB022193
5 of 33
Two theoretical estimates of the nucleation size in mode II are (Rice & Ruina,
1983
; Rubin & Ampuero,
2005
):
ℎ
푅푅푅푅
=
휋휋
4
휇휇휇휇
(1−
휈휈
)(
푏푏
−
푎푎
)(
휎휎
−
푝푝
)
;
ℎ
∗
푅푅푅푅
=
2
휋휋
휇휇휇휇푏푏
(1−
휈휈
)(
푏푏
−
푎푎
)
2
(
휎휎
−
푝푝
)
,
(6)
where
퐴퐴퐴퐴
is the shear modulus. The estimate
퐴퐴퐴
∗
푅푅푅푅
was derived from the linear stability analysis of steady fric
-
tional sliding by Rice and Ruina (
1983
). It also represents the critical cell size for steady-state quasi-static
sliding such that larger cells can become unstable on their own. Thus
퐴퐴퐴
∗
푅푅푅푅
represents a key length scale
to resolve for slow interseismic processes and earthquake nucleation (Lapusta et al.,
2000
; Rice & Rui
-
na,
1983
). The estimate
퐴퐴퐴
∗
푅푅퐴퐴
was determined in the parameter regime
퐴퐴퐴퐴
∕
푏푏푏
0
.
5
using the energy balance of
a quasi-statically expanding crack (Rubin & Ampuero,
2005
), and provides an estimate of the minimum size
for a slipping region that releases enough stored energy to result in the radiation of waves.
We aim to explore the impact of numerical resolution on the long-term simulated slip behavior of sequences
of earthquakes and aseismic slip. The nucleation size,
퐴퐴퐴
∗
, estimated by either
퐴퐴퐴
∗
푅푅푅푅
or
퐴퐴퐴
∗
푅푅퐴퐴
from Equation
6
,
is one length-scale that clearly needs to be well resolved. Early resolution studies for sequences of events
showed that resolution of the nucleation scale
퐴퐴퐴
∗
푅푅푅푅
by 20–40 cells is required for stable numerical results
(Lapusta et al.,
2000
). Later, the need to resolve the nucleation size by at least 20 cells was shown to be due to
the more stringent criterion of resolving the region where shear resistance breaks down at the rupture front,
often referred to as the cohesive zone. The cohesive zone can be an order of magnitude smaller than the nu
-
cleation size, depending on the constitutive description (Day et al.,
2005
; Lapusta & Liu,
2009
). The size of the
cohesive zone depends on the weakening rate
퐴퐴퐴퐴
of shear stress with slip associated with the constitutive law.
The quasi-static estimate
퐴퐴
Λ
0
of the cohesive zone size at near-zero rupture speed and constant
퐴퐴퐴퐴
is given by:
Λ
0
=
퐶퐶
1
휇휇
′
푊푊
,
(7)
where
퐴퐴퐴퐴
1
is a constant,
퐴퐴퐴퐴
′
=
퐴퐴
for mode III, and
퐴퐴퐴퐴
′
=
퐴퐴
∕(1−
휈휈
)
for mode II (Rice,
1980
). For standard
rate-and-state friction with the aging form of the state variable evolution, the weakening rate is given by
퐴퐴퐴퐴
=
퐷퐷
RS
∕(
푏푏푏푏푏
)
(Lapusta & Liu,
2009
) and:
Λ
0
=
퐶퐶
1
휇휇
′
퐷퐷
RS
푏푏푏푏푏
.
(8)
If one assumes that the traction distribution within the cohesive zone is linear, then the constant
퐴퐴퐴퐴
1
can be
approximated as
퐴퐴퐴퐴
1
=9
휋휋
∕32
(Rice,
1980
).
For fully dynamic rupture simulations, continuously resolving the breakdown process at the rupture front
becomes even more challenging as the cohesive zone size
퐴퐴
Λ
exhibits a contraction with increasing rupture
speed
퐴퐴퐴퐴
푅푅
(e.g., Rice,
1980
):
Λ=Λ
0
퐴퐴
−1
(
푣푣
푅푅
);
퐴퐴
−1
퐼퐼퐼퐼
=
(1−
휈휈
)
푐푐
2
푠푠
푣푣
2
푅푅
(1−
푣푣
2
푅푅
∕
푐푐
2
푠푠
)
1∕2
;
퐴퐴
−1
퐼퐼퐼퐼퐼퐼
=(1−
푣푣
2
푅푅
∕
푐푐
2
푠푠
)
1∕2
,
(9)
where
퐴퐴
=4(1−
푣푣
2
푅푅
∕
푐푐
2
푠푠
)
1∕2
(1−
푣푣
2
푅푅
∕
푐푐
2
푝푝
)
1∕2
−(2−
푣푣
2
푅푅
∕
푐푐
2
푠푠
)
1∕2
with
퐴퐴퐴퐴
푝푝
=
√
2(1−
휈휈
)(1−2
휈휈
)
퐴퐴
푠푠
. Note that
퐴퐴퐴퐴
−1
(0
+
)=1
, giving the quasi-static cohesive zone estimate
퐴퐴
Λ
0
when
퐴퐴퐴퐴
푅푅
=0
+
. As the rupture speed approach
-
es the limiting wave speed,
퐴퐴퐴퐴
푅푅
→
푐푐
푅푅
(Rayleigh wave speed) for mode II and
퐴퐴퐴퐴
푅푅
→
푐푐
푠푠
(shear wave speed) for
mode III, one has
퐴퐴퐴퐴
−1
(
푣푣
푅푅
)
→
0
and the width of the breakdown region approaches zero. Hence, it becomes
increasingly more challenging to resolve the rupture front during fully dynamic simulations if the rupture
accelerates toward the limiting speeds. Such acceleration typically occurs during long enough propagation
of dynamic rupture over favorable prestress, unless impeded by additional factors such as unfavorable pre
-
stress or situations with increasing effective breakdown energy, for example, due to off-fault inelasticity
or navigating fault roughness (Andrews,
2005
; Dunham et al.,
2011b
; Lambert & Lapusta,
2020
; Okubo
et al.,
2019
; Perry et al.,
2020
; Poliakov et al.,
2002
; Rice,
2006
). Simulations of faults with rate-and-state
friction and the aging form of the state variable evolution embedded in elastic bulk result in ruptures with
near-constant breakdown energy (Perry et al.,
2020
) and this holds for most cases considered in this study.
In Section
7
, we show that adding an approximation of off-fault inelasticity to our simulations that reduces
the rupture speeds does not alter our conclusions.
In our model, the fault contains a frictional domain consisting of two VW regions of length
퐴퐴퐴퐴
푉푉푉푉
= 32 km
that are separated by a 2-km-long VS region that impedes rupture propagation. We select large enough