GEOPHYSICAL RESEARCH LETTERS
Supporting Information for ”Dynamic Emergence of
Plate Motions and Great Megathrust Earthquakes
Across Length and Time Scales”
Jiaqi Fang
1
, Michael Gurnis
1
, 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
Contents of this file
1. Texts S1 to S6
2. Tables S1
3. Figures S1 to S13
Introduction
This supplementary material contains texts about dislocation and diffusion creep laws
(Text S1), thermal structure used in the model (Text S2), additional models referred to
in the main text (Text S3–S4), and further explanations of 2.5D models (Text S5–S6).
Corresponding author: J. Fang, Seismological Laboratory, California Institute of Technology,
Pasadena, CA, USA. (jfang@caltech.edu)
September 26, 2024, 10:09pm
X - 2
:
Key parameters of models presented in the main text are provided in Table S1. This
supplement also contains figures for additional 2D models (Figures S1–S8), the stress
field of the 2D model in the main text (Figures S9–S11), and additional 2.5D models
(Figures S12–S13).
Text S1. Dislocation and diffusion creep laws
Viscosities in the model are governed by a composite creep law with both dislocation
and diffusion creep, which is parametrized as:
η
disl/dif f
=
A
−
1
n
̇
ε
1
−
n
n
II
exp
E
+
PV
nRT
,
(1)
where
η
disl/dif f
are viscosities determined by dislocation or diffusion creep law,
A
is a
prefactor,
n
is the stress exponent,
R
is the gas constant,
E
is the activation energy,
V
is the activation volumn, ̇
ε
II
is the second invariant of the strain rate tensor,
P
is the
pressure, and
T
is the temperature. The values of parameters are listed in Table S1.
The dislocation creep parameters are based on values for dry olivine-rich rocks for the
bulk mantle and lithosphere and wet quartzites for the shear zone (Ranalli, 1995). The
diffusion creep parameters are calculated based on values for dry olivine assuming a grain
size of 1.5 mm (Karato & Wu, 1993). The model viscosities are calculated as:
η
=
1
η
disl
+
1
η
dif f
!
−
1
,
(2)
Text S2. Thermal structure
Temperatures of the incoming plate from ridge to trench and the underlying mantle
are calculated using the half-space cooling model assuming a plate velocity of 5 cm/yr
(consistent with the actual plate velocity), a plate age of 50 Myrs at the trench, a thermal
September 26, 2024, 10:09pm
:
X - 3
diffusivity of 10
−
6
m
2
/s, and a mantle potential temperature of 1700 K. Temperatures for
each vertical profile on and below the overriding plate are calculated based on half-space
cooling at a cooling age of 80 Myrs. The slab interface is assumed to be straight at depths
between 100 and 400 km, with the slope equal to that at the depth of 100 km. For each
profile perpendicular to the slab interface, temperatures are calculated based on the 1D
thermal diffusion starting from the reference profile before subduction, i.e., the vertical
temperature profile of the incoming plate at the trench (50 Myrs old). The diffusivity
is 10
−
6
m
2
/s and the diffusion time is equal to the along-dip distance divided by a plate
velocity of 5 cm/yr. Linear smoothing is applied in the buffer zone surrounding the slab,
with a size that is twice the width of the zone where temperatures are below 1630 K. The
temperature data are available on CaltechDATA at doi: 10.22002/awg6a-x3t07.
Text S3. Additional 2D models
In order to demonstrate the ability to simultaneously model long-term plate motions
and rapid earthquake ruptures, we have computed several additional models with minor
modifications on the original model. The model temperatures or viscosities, as well as
their outcomes (akin to Figure 2 in the manuscript), are shown in Figures S1–S8.
1.
Model A1
(Figures S1–S2): The distance of the right boundary to the trench is
extended from the original 1500 km to 2885 km. In both Model A1 and the original model,
the surface horizontal velocity tapers to zero at distances
<
1500 km to the trench. The
plate velocities, interseismic period and coseismic slips all align closely, indicating that the
free-slip boundary condition on the right wall has minimal influence on model behavior.
September 26, 2024, 10:09pm
X - 4
:
2.
Model A2
(Figures S3–S4): A 40-km-wide buffer zone is introduced to detach the
overriding plate from the right boundary. The overriding plate is thus without direct
pinning at the boundary of the domain. Although the overriding plate can move slightly
leftward due to underlying mantle flow, all model behavior shows remarkable consistency
with the original model.
3.
Model A3
(Figures S5–S6): The model includes the lower mantle with a maximum
depth of 1500 km and a subducting slab that reaches the 660-km interface. The lower
mantle viscosities are governed by diffusion creep and the values of rheological parameters
are adopted from Behr, Holt, Becker, and Faccenna (2022), and an adiabatic temperature
gradient of 0.3 K/km is added following Behr et al. (2022). Under this configuration, the
convergence rate is reduced to
≈
3 cm/yr, which results in longer interseismic period,
but all first-order long-term and short-term behaviors are preserved, including the steady
interseismic plate motions, rapid coseismic slips and stress changes with reasonable values.
4.
Model A4
(Figures S7–S8): Wet olivine rheology in Ranalli (1995) is used in the
entire upper mantle, and the lithosphere yield stress is increased from 200 MPa to 300
MPa to maintain a reasonable plate velocity. In this extreme case, the asthenosphere
viscosity is reduced and the plate velocity accelerates substantially (
≈
10 cm/yr). The
evolving interseismic and coseismic patterns of megathrust slips, surface displacements
and megathrust stresses, remain clear, despite the shorter period due to faster loading.
Text S4. Additional 2.5D models
In the original model, the 2.5D models are expanded for 2 elements with a total width
of 5 km in the third dimension. To verify the robustness of 2.5D models, we conducted
September 26, 2024, 10:09pm
:
X - 5
additional computations: (1) 4 along-strike elements with a total width of 5 km; (2)
2 along-strike elements with a total width of 5 km. The results (akin to Figure 4 in
the manuscript) are shown in Figures S12 and S13, respectively. The curves of surface
displacements and megathrust slips versus the along-strike length scale
̃
L
are identical for
different numbers of along-strike elements and model widths, suggesting that these factors
do not affect their relationship.
Text S5. Evolving boundary condition in the 2.5D model
Since Robin boundary conditions are not implemented in Underworld2, we developed
an iterative framework as an alternative using an evolving Neumann boundary condition.
For one step, the traction is imposed continuously on the front wall based on the solved
velocity from the previous iteration:
τ
a
n
+1
=
φτ
e
n
+ (1
−
φ
)
τ
a
n
,
(3)
where
τ
a
i
is the applied traction in the
i
-th iteration,
τ
e
i
is the traction expected from the
velocity field solved in the
i
-th iteration, and
φ
is a weighting coefficient. The coefficient
φ
smooths the spatially heterogeneous change of traction between iterations and stabilizes
the update of this boundary condition, with
φ
= min(
L
L
0
,
0
.
5), where
L
0
= 10
4
km.
The iteration terminates when the misfit between the applied traction and the expected
traction calculated from the resulting velocity field is below 1%.
Text S6. Estimation of the rupture length scale
In elastic crack models where a crack in an elastic body is subjected to uniform stress
changes, the relative displacements for all the crack modes follow elliptical distributions
September 26, 2024, 10:09pm
X - 6
:
(Scholz, 2019). We use this analytical solution to approximate the slip distribution in the
y
direction along the strike of the plate boundary. The slip
U
as a function of
y
is given
by:
U
=
U
0
̃
L
(
̃
L
2
−
y
2
)
1
2
.
(4)
U
0
is the maximum slip occurring at the center of the rupture where
y
= 0, and
̃
L
: Half
the rupture length, from the center to the tip where
U
= 0.
The gradient of slip
U
with respect to
y
is:
d
U
d
y
=
−
U
0
y
̃
L
(
̃
L
2
−
y
2
)
1
2
.
(5)
Under the approximation that
y <<
̃
L
, which is reasonable near the rupture center
(
y
≈
0), Equation 5 simplifies to:
dU
dy
≈−
U
0
y
̃
L
2
.
(6)
In the model, a Robin boundary condition is imposed at
y
=
D
, where
D
is the domain
width and
D <<
̃
L
. It models the slip gradient as the difference between the slip at
y
=
D
and the far-field slip (negligible compared with
U
0
and assumed to be zero), divided by a
length scale parameter
L
:
dU
dy
y
=
D
=
−
U
0
D
̃
L
2
≈−
U
0
L
(7)
This leads to the scaling relation between
̃
L
and
L
:
̃
L
=
√
LD.
(8)
References
September 26, 2024, 10:09pm
:
X - 7
Behr, W. M., Holt, A. F., Becker, T. W., & Faccenna, C. (2022). The effects of plate
interface rheology on subduction kinematics and dynamics.
Geophysical Journal
International
,
230
(2), 796–812.
Karato, S.-i., & Wu, P. (1993). Rheology of the upper mantle: A synthesis.
Science
,
260
(5109), 771–778.
Ranalli, G. (1995).
Rheology of the earth
. Springer Science & Business Media.
Scholz, C. H. (2019).
The Mechanics of Earthquakes and Faulting
. Cambridge University
Press.
September 26, 2024, 10:09pm
X - 8
:
Table S1.
Model parameters.
Quantity
Symbol
Unit
Bulk value
Shear zone value
Dislocation creep
a
Exponent
n
—
3.5
2.3
Prefactor
A
Pa
−
n
·
s
−
1
2.83
×
10
−
16
3.13
×
10
−
17
Activation energy
E
J
·
mol
−
1
5.32
×
10
6
1.54
×
10
6
Activation volume
V
m
3
·
mol
−
1
1.10
×
10
−
5
0
Diffusion creep
b
Exponent
n
—
1.0
—
Prefactor
A
Pa
−
n
·
s
−
1
1.40
×
10
−
11
—
Activation energy
E
J
·
mol
−
1
3.00
×
10
6
—
Activation volume
V
m
3
·
mol
−
1
6.00
×
10
−
6
—
Elasticity
Shear modulus
G
GPa
30 – 60
30
Plastic yielding
Cohesion
C
MPa
100
10
Static friction
μ
s
—
0.6
0.6
Dynamic friction
μ
d
—
0.6
0.2 – 0.8
Reference slip velocity
V
ref
cm/year
2.1
2.1
Pore pressure factor
λ
—
0
0 – 0.95
Maximum yield stress
τ
max
MPa
200
200
a
Ranalli (1995).
b
Karato and Wu (1993).
September 26, 2024, 10:09pm