Geophys. J. Int.
(2022)
229,
1098–1121
https://doi.org/10.1093/gji/ggab521
Advance Access publication 2021 December 27
GJI Seismology
Subduction earthquake sequences in a non-linear visco-elasto-plastic
megathrust
Luca Dal Zilio
,
1
,
2
,
3
Nadia Lapusta
,
1
,
2
Jean-Philippe Avouac
1
,
2
and Taras Gerya
3
1
Seismological Laboratory, California Institute of Technology (Caltech),
Pasadena, CA 91125
, United States. E-mail:
luca.dalzilio@erdw.ethz.ch
2
Department of Mechanical and Civil Engineering, California Institute of Technology (Caltech),
Pasadena, CA 91125
, United States
3
Institute of Geophysics, Department of Geosciences, Swiss Federal Institute of Technology (ETH Z
̈
urich), CH-
8093
Z
̈
urich, Switzerland
Accepted 2021 December 22. Received 2021 December 10; in original form 2021 August 22
SUMMARY
We present a 2-D thermomechanical computational framework for simulating earthquake
sequences in a non-linear visco-elasto-plastic compressible medium. The method is developed
for a plane-strain problem and incorporates an invariant formulation of the classical rate-
and state-dependent friction law and an adaptive time-stepping, which allows the time step
to vary by many orders of magnitude during a simulation. Long-term tectonic convergence is
imposed by displacing a boundary at a constant rate, whereas temperature-dependent viscosity
is solved using a rapidly converging Newton–Raphson scheme. The 2-D volume is discretized
using finite differences on a fully staggered grid and marker-in-cell techniques. An adaptive
free-surface approximation is used to modulate the air viscosity with the time step, which
allows stresses to vanish on the free surface during the propagation of fast slipping events. We
present a set of increasingly complex models in which we investigate how inertia, radiation
damping, thermally activated non-linear rheology and off-megathrust splay-fault events affect
sequences of
seismic
and
aseismic
slip on a simplified subduction megathrust. The new method
provides a unique computational framework to analyse earthquake sequences and to connect
forearc deformation with the dynamic properties of the megathrust, thus providing a physical
link between observations spanning from slow interseismic strain accumulation to localized
coseismic slip of individual earthquakes and post-seismic viscoelastic relaxation.
Key words:
Earthquake dynamics; Numerical modelling; Seismicity and tectonics; Subduc-
tion zone processes; Rheology and friction of fault zones.
1 INTRODUCTION
Subduction interfaces—often called megathrusts—produce the largest earthquakes on Earth (e.g. Wang
et al.
2012
;Lay
2015
). Despite
decades of observational, laboratory and theoretical studies, the processes leading to large megathrust earthquakes are still a matter of debate
(e.g. Kawamura
et al.
2012
; Kato & Ben-Zion
2020
). In particular, it remains unclear what controls earthquake complexity; persistent fault
features, failure physics and dynamics of the failure process can all play a role (e.g. Ye
et al.
2018
). Continuous geophysical monitoring has
shown that large megathrust earthquakes occur in the context of irregular interseismic strain accumulation, with static and dynamic stress
changes causing multiple effects on the fault zone, such as post-seismic slip (often called afterslip), viscoelastic relaxation and eventually
relocking of the fault (Wang
et al.
2012
; Bilek & Lay
2018
). Furthermore, decades of geophysical observations have demonstrated that
subduction megathrusts are the locus of a wide spectrum of slip modes, ranging from the dynamic propagation of seismic rupture to slow-slip
events of various durations and sizes, and to stable (aseismic) slip (e.g. B
̈
urgmann
2018
; Jolivet & Frank
2020
).
Seismological and geodetic data from different subduction zones, including rupture velocity, radiated energy, and the distribution
of
seismic
and
aseismic
slip (Lay
et al.
2012
; Bilek & Lay
2018
), provides clear evidence of depth-dependent behaviour of subduction
megathrusts (e.g. Lay
et al.
2012
;Ye
et al.
2018
). This depth dependence has been attributed to changes in fault properties (Huang
et al.
2012
; Noda & Lapusta
2013
), fluid pressure (Tobin & Saffer
2009
), the presence of weak sediments (Heuret
et al.
2012
), fault roughness and
geometry (Wang & Bilek
2011
) and rigidity of the upper plate (Sallar
`
es & Ranero
2019
). In particular, temperature is thought to constrain the
downdip extent of the seismogenic zone. The latter is particularly important when the
∼
350
◦
C isotherm is reached before the intersection
of the megathrust with the forearc Moho (Oleskevich
et al.
1999
), consistent with the common view that temperature is the dominant factor
controlling friction and viscoelastic behaviour of subduction megathrusts (e.g. Sibson
1982
; Blanpied
et al.
1991
; Hyndman & Wang
1993
;
Avouac
2015
). Furthermore, the transition zone from seismogenic to aseismic deformation often involves a mix of episodic slow-slip events,
1098
C
The Author(s) 2021. Published by Oxford University Press on behalf of The Royal Astronomical Society.
Downloaded from https://academic.oup.com/gji/article/229/2/1098/6484796 by California Institute of Technology user on 24 February 2022
Subduction earthquake cycles
1099
low frequency earthquakes and seismic tremors (e.g. Peng & Gomberg
2010
;Beroza&Ide
2011
; Houston
et al.
2011
; Avouac
2015
;Obara
& Kato
2016
;B
̈
urgmann
2018
; Jolivet & Frank
2020
).
Both geodetic and seismological observations have increased dramatically over the past two decades, yet their short instrumental record
makes it difficult to understand the maximum expected earthquake magnitude (McCaffrey
2008
). Records of the relative sea level change
extracted from corals may help to constrain the sequence of large events in the last centuries (e.g. Sieh
et al.
2008
; Philibosian
et al.
2017
), but defining the recurrence intervals and the long-term seismic pattern for a megathrust fault have proven to be problematic. In the
meantime, progress in computing technology has led to the development of numerical models of sequences of seismic and aseismic slip that
incorporate the physics governing earthquakes and slow-slip events (e.g. Tse & Rice
1986
; Ben-Zion & Rice
1995
; Lapusta
et al.
2000
;
Lapusta & Liu
2009
). Further advance of this class of models have incorporated variable material properties (Kaneko
et al.
2011
;Erickson
& Dunham
2014
), poroelastic effects (Heimisson
et al.
2019
), pore pressure evolution (Zhu
et al.
2020
), off-fault plasticity (Erickson
et al.
2017
; Thomas & Bhat
2018
), dilatant strengthening (Segall
et al.
2010
), semi-brittle rheologies (Goswami & Barbot
2018
; Lavier
et al.
2021
)
and viscoelasticity (Lambert & Barbot
2016
; Barbot
2018
; Allison & Dunham
2018
; Viesca & Dublanchet
2018
; Miyake & Noda
2019
).
The increasing complexity in resolving multiscale earthquake source processes has also called for extensive efforts to verify codes and their
reproducibility (Erickson
et al.
2020
;Jiang
et al.
2021
).
More recently, a new class of models has been developed to combine long-term geodynamic processes and short-term earthquake cycles
(van Dinther
et al.
2013b
,
a
; Sobolev & Muldashev
2017
). In particular, this approach allows for simulating seismic behaviour in relation to the
evolving tectonic forces, including elastic bending of the lithosphere (van Dinther
et al.
2014
), slab retreat and buoyancy forces (D’Acquisto
et al.
2020
; Dal Zilio
et al.
2020a
) in realistic subduction settings. However, due to simplifying assumptions, these simulations yield slip
transients that are many orders of magnitude slower than earthquakes in nature. In a recent study, Herrend
̈
orfer
et al.
(
2018
) significantly
improved this modelling approach by incorporating (i) an invariant formulation of the classical rate- and state-dependent friction equations
and (ii) an adaptive time-stepping (Lapusta
et al.
2000
) to capture the evolution of dynamic ruptures. However, some other simplifications
adopted in these simulations, for example linear rheology and the absence of temperature and gravity, have limited the possibility of comparing
the model predictions with observations of the earthquake sequences in nature.
In this study, we present VELO2CYCLEs (Visco-ELastO 2-D Cycles of Earthquakes), a continuum mechanics modelling approach to
simulate sequences of earthquakes and aseismic slip (SEAS) in a visco-elasto-plastic medium. The computational method is developed for a
plane-strain problem of a thrust fault governed by rate-and-state friction. The fault and bulk rheology is controlled by temperature-dependent
viscosity, which is given by a combination of diffusion and dislocation creep. Within the context of a time-stepping method, we solve the
Navier–Stokes equations by incorporating inertia effects. We propose a finite difference method, as it is easy to program, efficient, and supports
the marker-in-cell technique that can be applied to simulate large strains, as required to model long-term tectonic deformation (Gerya & Yuen
2007
;Gerya
2019
).
Here, we present our newly developed methodology and demonstrate its versatility to tackle a broad range of challenging problems in
computational earthquake physics. To test and illustrate the novel capabilities of this method, herein we present a set of increasingly complex
models in which we investigate how inertia, radiation damping, thermally activated non-linear rheology and off-megathrust splay-fault events
affect sequences of
seismic
and
aseismic
slip on a simplified subduction megathrust. In particular, the relatively simple model setup allows
us to compare the modelling results against other computational methods in the literature and simple analytical solutions. In Section 2,
we present the conceptual and mathematical model and state the continuum problem solved in this work. In Section 3, we present results
and compare solutions for quasi-dynamic (QD) and fully dynamic (FD) simulations (Section 3.1), non-linear rheology (Section 3.2) and
off-megathrust splay-fault events (Section 3.3). Finally, Section 4 analyses and interprets the results in the context of observations of large
megathrust earthquakes and discusses future avenues for research through the use of this methodology.
2 METHODOLOGY
2.1 Visco-elasto-plastic rheology and governing equations
We assume a continuum, isotropic, Maxwell-type model. The mathematical description of a deforming continuum in a stationary gravity field
is given by the following set of conservation equations for mass—eq. (
A1
), linear and angular momentum—eq. (
A4
), and energy—eq. (
A5
)
under the plane-strain assumption:
D
ρ
D
t
+
ρ
∇·
v
=
0
,
(1)
ρ
D
v
D
t
−∇·
σ
−
ρ
g
=
0
,
(2)
ρ
C
p
D
T
D
t
−∇·
(
k
∇
T
)
=
0
,
(3)
where
ρ
is density,
v
the velocity vector,
σ
is the Cauchy stress tensor,
g
is the gravity vector,
C
p
the isobaric heat capacity,
T
is the temperature
and
k
is the thermal conductivity.
Downloaded from https://academic.oup.com/gji/article/229/2/1098/6484796 by California Institute of Technology user on 24 February 2022
1100
Luca Dal Zilio
et al
.
The stress is decomposed into a deviatoric and volumetric component:
σ
=
τ
−
p
I
(4)
where
τ
is the deviatoric stress tensor,
p
=−
1
3
Tr
σ
is the pressure or the mean compressive stress, and
I
is the identity tensor.
The volumetric deformation is assumed to be purely elastic, with the bulk modulus
K
defined as
1
K
=
1
ρ
D
ρ
D
p
,
(5)
such that the time derivative of density in eq. (
A1
) satisfies the following:
D
ρ
D
t
=
ρ
K
D
p
D
t
.
(6)
The deviatoric deformation occurs as a combination of elastic, viscous, and plastic deformation. Components of the deviatoric strain-rate
tensor are defined via spatial gradients of the velocity under the assumption of
infinitesimal
strain as follows:
̇
ε
ij
=
1
2
∂v
i
∂
x
j
+
∂v
j
∂
x
i
−
1
3
Tr ̇
ε
ij
,
(7)
where
i
and
j
are coordinate indexes and
x
i
and
x
j
are spatial coordinates such that in 2-D we can define four tensor components. We assume
that strain rate can be additively decomposed into its elastic, viscous and plastic components as
̇
ε
ij
=
[ ̇
ε
ij
]
viscous
+
[ ̇
ε
ij
]
elastic
+
[ ̇
ε
ij
]
plastic
.
(8)
The elastic part is defined as
[ ̇
ε
ij
]
elastic
=
1
2
μ
̃
D
̃
D
t
(
τ
ij
)
,
(9)
where
μ
is the shear modulus, while the differential operator
̃
D
/
̃
D
t
denotes the corotational time derivative—that is, the local rates of change
in the deviatoric stress with respect to translations, rotations, and/or deformation.
Viscous creep obeys either Newton’s law of viscous friction or is computed in terms of deformation invariants and depends on strain
rate, temperature and pressure (Ranalli
1995
). The viscous behaviour of a rock is defined here via a single effective viscosity according to
η
eff
=
τ
II
2 ̇
ε
II(v)
,
(10)
where the
τ
II
and ̇
ε
II(v)
denote the square root of the second invariant of the deviatoric stress and the viscous strain-rate tensors, respectively.
The unified effective viscosity (
η
eff
) is calculated as a combination of diffusion (
η
diff
) and dislocation creep (
η
disl
):
η
diff
=
1
2
A
d
τ
1
−
n
tr
exp
E
+
pV
RT
,
(11)
η
disl
=
A
d
τ
n
−
1
II
exp
−
E
+
pV
RT
.
(12)
where
A
d
is the pre-exponential factor,
n
is the stress exponent,
E
is the activation energy,
V
is the activation volume,
R
is the gas constant and
τ
tr
defines the stress transition between diffusion creep and dislocation creep, which is assumed to occur at 30 kPa (Turcotte & Schubert
2002
).
Dislocation creep is dominant for high-stress levels, while diffusion creep—a thermally activated mechanism—is dominant at low-stress
conditions. When both diffusion and dislocation creep mechanisms are simultaneously active, the two viscous mechanisms are combined
(Appendix A) while a Newton–Raphson scheme with quadratic convergence is used to obtain an effective viscosity (
η
eff
) (Appendix B). The
effective viscosity is then used to define the viscous contribution to the strain-rate tensor:
[ ̇
ε
ij
]
viscous
=
τ
ij
2
η
eff
.
(13)
Plastic (irreversible) deformation and frictional sliding is governed by rate- and state-dependent friction (RSF), conclusively documented
in laboratory experiments and used to reproduce, both qualitatively and quantitatively, a number of earthquake-source observations (Dieterich
1979
; Dieterich
et al.
1981
; Ruina
1983
; Marone
1998
; Blanpied
et al.
1991
; Tullis
1996
; Dieterich
2007
, and references therein). Such laws
express the dependence of frictional shear strength
τ
y
on the confining pressure
p
and on steady-state friction
f
∗
at the reference slip rate
V
0
,
which depends on sliding velocity
V
(
t
) and state variable
θ
(
t
), with the latter obeying a state evolution equation. Notably, in our continuum
approach physical quantities are
invariant
of the coordinate system, which allows the solution to adapt to spontaneous bulk evolution. As
such, our approach contains an
invariant
reformulation of the
regularized
version of RSF (Nakatani
2001
; Herrend
̈
orfer
et al.
2018
)
τ
y
=
τ
II
=
ap
arcsinh
⎡
⎣
V
(
t
)
2
V
0
exp
⎛
⎝
f
∗
+
b
ln
θ
(
t
)
V
0
D
RS
a
⎞
⎠
⎤
⎦
+
μ
2
c
s
V
(
t
)
,
(14)
where the
invariant
shear resistance (
τ
y
) depends on the slip rate
V
(
t
), the effective confining pressure
p
, the characteristic slip distance
D
RS
,
the rate-and-state parameters
a
,
b
, and the state variable
θ
(
t
). Note that eq. (
14
) imposes that
τ
y
=
τ
II
; hence, the yielding condition is always
Downloaded from https://academic.oup.com/gji/article/229/2/1098/6484796 by California Institute of Technology user on 24 February 2022
Subduction earthquake cycles
1101
fulfilled and some plastic deformation always occurs. This way some slip always occurs if shear stress is larger than zero. Compared to the
classical formulation of rate-and-state friction (e.g. Lapusta
et al.
2000
), our
invariant
formulation replaces the effective normal stress
σ
n
with the effective confining pressure
p
, and the shear traction
τ
with the shear strength
τ
y
.
The slip rate
V
(
t
) is computed as the magnitude of the second invariant of deviatoric plastic strain rate ̇
ε
II(p)
over the thickness of the
fault zone
W
z
V
(
t
,
x
)
=
2 ̇
ε
II(p)
W
z
.
(15)
The state variable
θ
describes the evolutionary effects in response to changing slip rates
d
θ
d
t
=
1
−
V
θ
D
RS
.
(16)
This state evolution law, commonly referred to as the Dieterich–Ruina or
aging law
, has experimental support (Beeler
et al.
1996
)and
partial theoretical justification (Rice
et al.
2001
). Other formulations for the evolution of the state variable exist, such as the
slip law
(Ruina
1983
) as well as various composite laws, and the formulation that best describes various laboratory experiments remains a topic of ongoing
research (Bhattacharya
et al.
2017
). The last term in eq. (
14
), (
μ/
2
c
s
)
·
V
(
t
), represents the
radiation damping
term (Rice
1993
; Cochard &
Madariaga
1996
), which depends on the shear modulus (
μ
) and the shear wave speed (
c
s
) of the medium where the fault is embedded. The
radiation damping
term is exclusively used in eq. (
14
) only when inertia is not solved in eq. (
A4
) (i.e. QD formulation). From a physical point
of view, the radiation damping approximation is used to mimic the energy lost due to seismic-wave propagation (e.g. Crupi & Bizzarri
2013
).
For a constant slip velocity
V
(
t
), the state variable
θ
evolves to its steady-state value
θ
ss
=
L
/
V
(
t
), transforming the shear strength
τ
y
into
its steady-state form
τ
ss
=
p
f
∗
+
(
a
−
b
)
ln
V
(
t
)
V
0
.
(17)
Eq. (
17
) highlights the importance of the non-dimensional parameter (
a
−
b
) for the rate dependence of friction; stable slip results for
velocity-strengthening (VS) friction, (
a
−
b
)
>
0, and conditionally unstable slip for velocity-weakening (VW) friction, (
a
−
b
)
<
0(Rice&
Ruina
1983
; Scholz
1998
).
In this study, we assume the dilation angle to be zero, which implies that the volume does not change during plastic yielding. Accordingly,
the deviatoric plastic strain is defined as
[ ̇
ε
ij
]
plastic
=
χ
∂τ
ij
2
∂τ
II
,
(18)
with
χ
being the
plastic multiplier
, which serves as
variable scaling coefficient
to ensure that the deviatoric stress always satisfies the yield
criteria (
τ
II
=
τ
y
), and it is equal to the double of the square root of the second invariant of the deviatoric plastic strain rate tensor (
χ
=
2 ̇
ε
II(p)
)
(Gerya
2019
).
The deviatoric components (
τ
ij
) of the stress tensor
σ
in eq. (
A4
) are obtained from the visco-elasto-plastic constitutive relationships
(eq.
8
) by using an implicit first-order finite-difference scheme in time in order to represent objective time derivatives of viscoelastic stresses
(e.g. Moresi
et al.
2003
)
τ
ij
=
2
η
vp
Z
̇
ε
ij
+
τ
0
ij
·
(1
−
Z
)
,
(19)
where
Z
is the viscoelasticity factor (Schmalholz
et al.
2001
)
Z
=
μ
t
μ
t
+
η
vp
,
(20)
and
η
vp
is the effective viscoplastic viscosity that characterizes the combined viscoplastic deformation:
η
vp
=
η
eff
τ
II
2
η
eff
̇
ε
II(p)
+
τ
II
,
(21)
where ̇
ε
II(p)
is the square root of the second invariant of the deviatoric plastic strain rate tensor.
2.2 Adaptive time-stepping
To resolve different timescales—ranging from long-term (slow) deformation to periodic (fast) slip events—a variable time stepping is required
during computation. We use the approach of Lapusta
et al.
(
2000
), with a variable time step (
t
) chosen as
t
=
max [
t
min
,
t
dyn
]
,
(22)
where
t
min
is the minimum time step to resolve dynamic rupture propagation, whereas
t
dyn
is the dynamic time step, which varies with
slip velocity. The minimum time step is given by
t
min
=
γ
x
/
c
s
,
(23)
Downloaded from https://academic.oup.com/gji/article/229/2/1098/6484796 by California Institute of Technology user on 24 February 2022
1102
Luca Dal Zilio
et al
.
where
γ
=
1/4 is a constant suggested by dynamic rupture simulations (Day
et al.
2005
),
x
is the minimum grid size. The dynamic time
step
t
dyn
is set to be inversely proportional to the maximum slip velocity (
V
max
):
t
dyn
=
θ
max
D
RS
V
max
,
(24)
where the coefficient
θ
max
is determined in Lapusta
et al.
(
2000
) by first computing a prescribed parameter that depends on the RSF
parameters and the grid size
ξ
=
1
4
KD
RS
aP
−
b
−
a
a
2
−
KD
RS
ap
,
(25)
with
K
=
2
π
μ
∗
x
defining the stiffness (assuming that slip is homogeneous throughout a single cell), and with the shear modulus
μ
∗
=
μ
/(1
−
ν
), where
ν
is the Poisson’s ratio. Accordingly
θ
max
=
⎧
⎨
⎩
min
ap
KD
RS
−
(
b
−
a
)
p
,
0
.
2
→
ξ>
0
min
1
−
(
b
−
a
)
p
KD
RS
,
0
.
2
→
ξ<
0
.
(26)
Following Lapusta & Liu (
2009
),
θ
max
is not allowed to be larger than 0.2. As the result of this adaptive time-stepping, large time-steps
(years) are used in the interseismic period, while small time-steps (fraction of a second) are used to capture the evolution of each dynamic
rupture. Note that the stability of this adaptive time-stepping procedure relies on the presence of the positive direct effect (
a
>
0) in the RSF
formulation (eq.
14
).
2.3 Model of a subduction megathrust
To study the behaviour of a generic subduction megathrust (Fig.
1
a), we design a simple and computationally efficient setup that mimics
a subduction zone in which a slab subducts beneath a continental plate (Fig.
1
a). We consider a 2-D, plane-strain model, with a 1-D fault
embedded in a 2-D isotropic, visco-elasto-plastic medium (Fig.
1
b). This model setup assumes a megathrust interface represented as a
finite-width fault zone. In order to promote depth dependency of the physical parameters, gravity acceleration (
g
) is applied with a slab dip
angle of 10
◦
(Fig.
1
b).
A key element in models of sequences of earthquakes and aseismic slip is the spatial distribution of the rate-and-state dependency of
friction (
a
−
b
)(Fig.
1
c). To allow for seismic slip, we consider a segment with steady-state velocity-weakening properties surrounded by
an updip and downdip VS regions. To avoid the formation of splay-faults and allow the comparison of our modelling results with outputs
from other models, we assume velocity-neutral frictional behaviour within the wedge. We then relax this assumption to assess the impact
of splay-faults when they are allowed to form. The seismogenic zone incorporates gradual variations of rate and state parameters near the
trench (12
<
x
<
8 km), and in the downdip region of the seismogenic zone (62
<
x
<
58 km). The length of the seismogenic segment (
x
2
−
x
1
, Fig.
1
c) is 50 km. In the downdip region (
x
>
62 km), the VS properties of the fault allow the megathrust to slip at a rate similar or
equal to plate rate. While the depth distributions of rate and state parameters assumes constant
a
and
L
,
b
increases with depth resulting in a
transition to stable steady-state sliding. Pressure (
p
), which plays the role of the normal traction on the fault interface (
σ
n
), is set to 30 MPa.
The constant value of
p
close to the free surface is chosen for numerical tractability (Lapusta & Liu
2009
), and to investigate several issues
unrelated to the free surface, such as interaction of rupture with stress heterogeneities over several earthquake sequences. Model parameters
are given in Tables
1
and
2
.
2.4 Initial and boundary conditions
The system of equations has a unique solution given appropriate boundary conditions for the velocity and pressure (Fig.
1
b). All initial
unknowns in eq. (
A1
)andeq.(
A4
) are set to zero and they evolve as a result of the boundary conditions. The applied velocity boundary
conditions resemble those of an analog sandbox model. For the specifics of the model setup presented in this study, the conditions on the left
(
x
=
0), right (
x
=
L
x
) and upper boundaries (
y
=
0) are defined
free slip
as follow:
∂v
y
∂
x
=
0
x
=
0
,v
x
=
0;
(27a)
∂v
y
∂
x
=
0
x
=
L
x
,v
x
=
0;
(27b)
∂v
x
∂
y
=
0
y
=
0
,v
y
=
0
.
(27c)
A
free slip
condition requires that the normal velocity component on the boundary is zero and the other component do not change
across the boundary. On the other hand, the fault and the wedge are loaded from the lower boundary (
y
=
L
y
) by applying Dirichlet boundary
conditions for horizontal velocity
v
x
and zero velocity on the vertical component
v
y
(Fig.
1
b):
v
x
=
V
pl
|
y
=
L
y
,v
y
=
0
|
y
=
L
y
,
(28)
Downloaded from https://academic.oup.com/gji/article/229/2/1098/6484796 by California Institute of Technology user on 24 February 2022