of 10
Dynamic
Emergence
of
Plate
Motions
and
Great
Megathrust
Earthquakes
Across
Length
and
Time
Scales
Jiaqi
Fang
1
,
Michael
Gurnis
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
Abstract
The
slow
motion
of
tectonic
plates
over
thousands
of
kilometers
is
intermittently
interrupted
by
great
earthquakes
with
sudden
slips
localized
near
convergent
plate
boundaries.
We
developed
a
subduction
model
that
self‐consistently
integrates
buoyancy
forces,
diffusion
and
dislocation
creep,
and
inter‐plate
friction.
From
the
nonlinear
dynamics
emerge
long‐term
plate
motions
that
achieve
velocities
of
5
cm/year,
effective
viscosities
of
10
19
Pa⋅s
below
plates,
and
sudden
slips
up
to
10
m
repeating
every
several
hundred
years.
Along‐strike
resistance
arising
from
long‐wavelength
variation
of
coseismic
slip
is
naturally
incorporated
with
a
rupture
length
scale,
̃
L
.
Computations
with
̃
L
10
3
km
generate
events
with
M
w
9.
When
̃
L
decreases,
there
is
a
commensurate
decrease
in
the
effective
moment
of
rupture
events.
Predicted
long‐term
plate
velocities,
mantle
viscosities,
cycles
of
stress
loading
and
release,
and
rupture
event
size
and
magnitude
all
show
good
agreement
with
observations.
Plain
Language
Summary
Subduction
zones
are
regions
where
one
tectonic
plate
slides
under
another
at
a
velocity
of
several
centimeters
per
year,
but
the
plates
can
slip
suddenly
during
great
earthquakes
in
the
convergent
plate
boundary.
We
developed
a
model
that
self‐consistently
integrates
tectonic
forces
and
material
properties.
This
model
replicates
steady
plate
motions
of
several
centimeters
per
year
and
rupture
events
with
slips
up
to
10
m
and
a
period
of
several
hundred
years.
Mantle
viscosities
and
space–time
patterns
of
stress
accumulation
and
release
also
align
with
observations.
During
great
earthquakes,
the
rupturing
segment
advances
rapidly
but
is
held
by
neighboring
non‐rupturing
segments.
Resisting
forces
arising
from
this
along‐
strike
variation
are
naturally
incorporated
and
tested
with
different
rupture
length
scales,
and
Magnitude
7–9
earthquakes
emerge.
The
models
highlight
the
significance
of
along‐strike
slip
variations
in
modulating
earthquake
processes.
1.
Introduction
Megathrust
faults
at
convergent
plate
boundaries
host
the
largest
earthquakes
with
slips
of
meters
that
interrupt
the
regular
several
cm
/year
plate
motions.
Historical
records
of
megathrust
earthquakes
show
a
complicated
variability
in
space,
time
and
source
properties
(Bilek
&
Lay,
2018
;
Sykes
&
Quittmeyer,
1981
;
Thatcher,
1990
;
Ye
et
al.,
2018
).
This
variability
is
thought
to
be
intimately
linked
to
physical
properties
of
the
megathrust,
for
example,
the
interface
roughness,
material
composition,
sediment
thickness,
fluid
content
and
thermal
structure,
which
directly
influence
rupture
mechanics
(Heuret
et
al.,
2012
;
Hyndman
et
al.,
1997
;
J.
Li
et
al.,
2018
;
K.
Wang
&
Bilek,
2014
;
Wirth
et
al.,
2022
).
Larger‐scale
geodynamic
controls,
such
as
subduction
zone
geometry
and
mantle
rheology,
also
play
a
crucial
role
in
shaping
the
complexity
of
megathrust
earthquake
occurrence,
as
they
modulate
the
tectonic
loading,
material
transport,
and
ultimately
stresses
delivered
to
the
megathrust
(Dal
Zilio
et
al.,
2018
;
Nishikawa
&
Ide,
2014
;
Ruff
&
Kanamori,
1980
;
K.
Wang
et
al.,
2012
;
Wirth
et
al.,
2022
).
Meanwhile,
background
plate
motions,
resulting
from
a
balance
between
the
driving
slab
pull
and
resisting
forces,
including
the
viscous
shear
in
the
asthenosphere,
plate
bending
and
plate
boundary
friction,
are
regulated
by
the
same
factors
inside
and
outside
the
megathrust
(Alisic
et
al.,
2012
;
Behr
et
al.,
2022
;
Conrad
&
Hager,
1999
;
King,
2016
;
Y.
Li
&
Gurnis,
2023
;
Zhong
et
al.,
1998
).
Cycles
of
steady
interseismic
motions
and
sudden
coseismic
slip
reflect
an
intricate
interplay
between
tectonic
structures
and
material
properties.
As
plate
motions
and
earthquake
ruptures
span
a
large
range
of
length
and
time
scales
and
are
thought
to
be
governed
by
different
constitutive
relationships,
it
has
been
challenging
to
integrate
them
into
a
unified
computational
model.
Models
of
long‐term
plate
deformation
(Stadler
et
al.,
2010
;
Tackley,
2000
;
Zhong
&
Gurnis,
1995
)
and
seismic
cycles
(Barbot
et
al.,
2012
;
Lapusta
&
Liu,
2009
;
D.
Li
&
Liu,
2021
)
are
usually
RESEARCH
LETTER
10.1029/2024GL110821
Key
Points:
Development
of
a
self‐consistent
sub-
duction
model
that
integrates
the
dy-
namics
of
plate
tectonics
and
great
megathrust
earthquakes
With
realistic
model
parameters,
predicted
long‐term
plate
motions
and
seismic
cycle
behaviors
align
with
ob-
servations
in
subduction
zones
The
along‐strike
resistance
arising
from
slip
variations
plays
a
critical
role
in
influencing
the
magnitude
of
a
rupture
event
Supporting
Information:
Supporting
Information
may
be
found
in
the
online
version
of
this
article.
Correspondence
to:
J.
Fang,
jfang@caltech.edu
Citation:
Fang,
J.,
Gurnis,
M.,
&
Lapusta,
N.
(2024).
Dynamic
emergence
of
plate
motions
and
great
megathrust
earthquakes
across
length
and
time
scales.
Geophysical
Research
Letters
,
51
,
e2024GL110821.
https://doi.
org/10.1029/2024GL110821
Received
3
JUL
2024
Accepted
9
NOV
2024
Author
Contributions:
Conceptualization:
Jiaqi
Fang,
Michael
Gurnis,
Nadia
Lapusta
Data
curation:
Jiaqi
Fang
Formal
analysis:
Jiaqi
Fang
Funding
acquisition:
Michael
Gurnis,
Nadia
Lapusta
Methodology:
Jiaqi
Fang,
Michael
Gurnis,
Nadia
Lapusta
Project
administration:
Michael
Gurnis,
Nadia
Lapusta
Resources:
Jiaqi
Fang,
Michael
Gurnis
Software:
Jiaqi
Fang
Supervision:
Michael
Gurnis,
Nadia
Lapusta
Validation:
Jiaqi
Fang,
Michael
Gurnis,
Nadia
Lapusta
Visualization:
Jiaqi
Fang
Writing
original
draft:
Jiaqi
Fang
©
2024.
The
Author(s).
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.
FANG
ET
AL.
1
of
10
implemented
separately.
For
subduction
zones,
there
have
been
a
few
studies
using
2D
cross‐sections
with
visco‐
elasto‐plasticity
that
generate
earthquake‐like
ruptures
and
seismic
cycles
consistent
with
observations
(Dal
Zilio
et
al.,
2022
;
Herrendörfer
et
al.,
2015
;
Oryan
&
Buck,
2020
;
Sobolev
&
Muldashev,
2017
;
van
Dinther
et
al.,
2013
,
2014
).
They
provide
valuable
insight
into
the
interplay
between
tectonic
loading,
rapid
rupture
processes,
lithosphere
stress
state
and
aftershock
seismicity.
Critically,
the
imposed
loading
velocity
boundary
conditions
are
not
self‐consistent
with
the
rheology
in
the
system
and
ignore
the
impact
of
great
earthquakes
on
long‐term
plate
motions
and
stress
states.
Being
2D,
existing
approaches
miss
along‐strike
slip
variations.
A
great
earthquake
only
ruptures
a
finite
length
along
the
plate
boundary
and
displacements
of
the
rupturing
segment
will
be
resisted
by
forces
imposed
by
adjacent
non‐rupturing
segments.
Typically,
the
size
in
this
third
dimension
governs
how
large
a
megathrust
earthquake
becomes
(Kanamori,
1978
;
Scholz,
1982
).
We
formulated
a
visco‐elasto‐plastic
subduction
model
that
treats
the
dynamics
of
plate
tectonics
self‐
consistently
with
great
megathrust
earthquakes.
The
model
is
driven
by
negative
buoyancy
with
a
realistic
thermal
structure
instead
of
imposed
kinematics
while
employing
experimentally‐derived
rheological
parameters.
The
model
is
capable
of
simultaneously
predicting
long‐term
plate
motion
and
mantle
convection
as
well
as
short‐
term
slip
and
relaxation
behaviors.
Given
the
computational
difficulties,
instead
of
implementing
a
full
3D
model,
we
efficiently
incorporate
the
along‐strike
resistance
by
expanding
the
model
into
a
small
thickness
and
applying
tractions
proportional
to
the
incremental
displacements
that
occur
during
a
megathrust
rupture.
The
results
highlight
how
the
length
scale
of
this
along‐strike
resistance
from
long‐wavelength
slip
variations
modulates
rupture
behaviors.
2.
Methods
2.1.
The
2D
Model
The
models
are
computed
using
Underworld2,
open‐source
finite‐element
software
that
solves
the
conservation
of
mass
and
momentum
for
an
incompressible
medium
(Mansour
et
al.,
2020
;
Moresi
et
al.,
2007
).
Originally
tailored
for
long‐time‐scale
mantle
convection,
Underworld2
has
been
adapted
to
shorter‐time‐scale
problems
with
Maxwell
viscoelastic
and
visco‐elasto‐plastic
rheology
(Farrington
et
al.,
2014
;
Moresi
et
al.,
2002
).
This
computational
method
has
expanded
its
applicability
to
earthquake‐related
problems
(Yang
et
al.,
2023
),
which
has
also
been
demonstrated
with
several
other
numerical
implementations
(Dal
Zilio
et
al.,
2022
;
Sobolev
&
Muldashev,
2017
;
van
Dinther
et
al.,
2013
).
We
set
up
a
2D
model
of
a
4040
×
660
km
2
trench‐normal
section
with
the
whole
incoming
plate
(Figure
1a
).
The
megathrust
is
approximated
by
a
6‐km‐wide
shear
zone
between
the
incoming
and
overriding
plates
(Figure
1b
).
The
domain
boundaries
are
free‐slip,
and
the
deformation
is
driven
entirely
by
internally
generated
forces.
The
system
has
initial
conditions
based
on
results
of
large‐time‐step
computations,
which
reflect
the
long‐term
dy-
namics
of
this
system,
and
progresses
through
multiple
cycles
of
coseismic,
postseismic
and
interseismic
phases
with
a
time
step
size
of
5
years.
Displacements
and
stress
changes
are
averaged
over
each
5‐year
step
and
the
inertia
term
is
omitted.
The
mesh
has
1024
×
256
linear
elements
and
is
refined
in
the
area
near
the
megathrust,
with
a
minimum
element
size
of
1.2
km.
The
rheology
of
the
lithosphere
and
mantle
is
a
composite
of
creeping,
elastic
and
plastic
processes,
with
key
parameters
listed
in
Table
S1
in
Supporting
Information
S1
.
The
viscosity
(Figure
1d
)
is
determined
by
diffusion
and
dislocation
creep
(Karato
&
Wu,
1993
;
Ranalli,
1995
).
The
shear
modulus
is
30
GPa
for
the
shear
zone,
top
6
km
of
the
incoming
plate
and
20
km
of
the
overriding
plate,
and
60
GPa
for
other
parts
of
the
bulk
lithosphere
and
mantle
throughout.
Plastic
yielding
follows
the
Drucker‐Prager
criterion
in
which
the
second
invariant
of
the
deviatoric
stress,
τ
I
I
,
is
limited
by
a
yield
stress
defined
as:
τ
y
=
C
+
μP
(
1
λ
)
,
(
1
)
where
C
is
the
cohesion,
μ
is
the
friction
coefficient,
P
is
the
pressure,
and
λ
is
a
pore
pressure
factor.
Fluids
in
the
megathrust
are
approximated
with
a
fixed
pore
pressure
factor
λ
=
0
.
95
in
the
shear
zone
shallower
than
60
km
and
a
linear
transition
from
0.95
to
0
from
60
to
80
km.
Rate‐dependent
friction
coefficients
from
van
Dinther
et
al.
(
2013
)
are
applied
in
the
shear
zone:
Writing
review
&
editing:
Jiaqi
Fang,
Michael
Gurnis,
Nadia
Lapusta
Geophysical
Research
Letters
10.1029/2024GL110821
FANG
ET
AL.
2
of
10
19448007, 2024, 22, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2024GL110821 by California Inst of Technology, Wiley Online Library on [27/11/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
μ
=
μ
d
yn
+
μ
st
μ
d
yn
1
+
V
V
ref
,
(
2
)
where
μ
st
and
μ
dyn
are
static
and
dynamic
friction
coefficients,
V
ref
is
a
reference
slip
velocity,
and
V
is
the
local
slip
velocity,
defined
as
twice
the
shear
zone
width
multiplying
the
second
invariant
of
the
strain
rate
(van
Dinther
et
al.,
2013
;
Yang
et
al.,
2023
).
The
temperature
dependence
of
frictional
parameters
in
the
shear
zone
(Figure
1c
)
is
guided
by
previous
studies
(Den
Hartog
et
al.,
2012
;
Liu
&
Rice,
2005
).
The
top
and
bottom
sections
are
rate‐
strengthening
μ
d
yn
>
μ
st
)
,
favoring
stable
aseismic
slip.
The
strong
rate‐weakening
friction
μ
dyn
<
μ
st
)
allows
spontaneous
slip
instability
in
the
seismogenic
zone
at
an
intermediate
depth.
The
modeling
technique
developed
in
this
work
will
allow
us
to
consider
other
laboratory‐derived
constitutive
formulations
for
faults
that
have
been
shown
to
reproduce
a
number
of
earthquake
source
phenomena,
such
as
rate‐and‐state
friction
and
additional
effects
due
to
frictional
shear
heating,
evolution
of
pore
fluid
pressure,
and
coupling
between
on‐fault
and
off‐
fault
deformation
(Allison
&
Dunham,
2021
;
Ampuero
&
Ben‐Zion,
2008
;
Dieterich,
1979
;
Di
Toro
et
al.,
2011
;
Faulkner
et
al.,
2011
;
Lambert
et
al.,
2021
;
Marone,
1998
;
Rice,
2006
;
Segall
et
al.,
2010
).
2.2.
The
2.5D
Model
With
Along‐Strike
Resistance
Rapid
slips
in
the
rupturing
segment
will
encounter
resistance
from
neighboring
non‐rupturing
segments
during
a
great
megathrust
earthquake,
an
effect
missing
in
2D.
While
a
full
3D
model
is
computationally
prohibitive,
we
compensate
for
the
lack
of
along‐strike
resistance
by
expanding
the
domain
in
the
x
z
plane
to
two
elements
with
a
thickness
of
5
km
in
the
along‐strike
y
direction
and
initializing
with
state
variables
from
the
step
before
the
first
rupture
event.
The
thin
2.5D
model
is
located
at
the
center
of
the
rupturing
zone
and
a
traction
is
imposed
on
the
front
wall,
which
is
caused
by
the
strain
between
the
rapidly
advancing
segment
in
the
model
domain
and
distal
segments
with
an
interseismic
velocity.
The
expected
shear
traction
is
proportional
to
a
long‐wavelength
displacement
difference
and
inversely
proportional
to
an
along‐strike
length
scale,
L
:
τ
e
=
G
u
u
far
L
=
G
Δ
t
v
v
far
L
,
(
3
)
Figure
1.
(a)
Thermal
structure
of
the
model
(see
details
in
Text
S2
in
Supporting
Information
S1
).
x
,
z
and
T
denote
the
horizontal
length,
depth
and
temperature.
(b)
Enlargement
of
the
thermal
structure
near
the
megathrust
with
shear
zone
(width
is
6
km,
5
×
local
element
size)
surrounded
by
silver
dashed
lines.
(c)
The
temperature‐dependent
frictional
parameters,
μ
st
and
μ
dyn
,
in
the
shear
zone.
The
top
secondary
axis
shows
the
depth
z
corresponding
to
the
temperature
T
.
(d)
Effective
viscosity
η
initially
for
the
model,
from
a
viscously‐dominated
computation
with
a
large
time
step
size.
Geophysical
Research
Letters
10.1029/2024GL110821
FANG
ET
AL.
3
of
10
19448007, 2024, 22, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2024GL110821 by California Inst of Technology, Wiley Online Library on [27/11/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
where
τ
e
is
the
traction
on
the
front
wall,
u
u
far
and
v
v
far
denotes
the
displacement
and
velocity
difference
from
a
distal
reference,
τ
e
xy
depends
on
v
x
v
far
x
,
τ
e
zy
depends
on
v
z
v
far
z
,
and
τ
e
yy
=
0.
Equation
3
is
essentially
a
Robin
boundary
condition
that
specifies
a
linear
combination
of
velocity
and
stress
on
the
front
wall.
With
Robin
boundary
conditions
not
available
in
Underworld2,
we
developed
an
iterative
approach
using
evolving
Neumann
boundary
conditions
within
each
step.
The
traction
applied
on
the
front
wall
is
constantly
updated
in
each
iteration
based
on
Equation
3
,
until
the
misfit
between
the
applied
traction
and
expected
traction
calculated
from
the
velocity
solution
is
below
1%.
The
parameter
L
reflects
the
gradient
of
slips
at
the
front
wall
and
describes
a
length
scale
of
this
along‐strike
resistance
due
to
the
long‐wavelength
displacement
difference
between
the
rupture
center
and
non‐rupturing
distal
segments.
A
small
L
indicates
rupture
confined
to
a
narrow
along‐strike
width,
and
conversely
for
a
large
L
.
Considering
the
nonlinear
decay
of
slip
in
the
along‐strike
direction,
L
scales
with,
but
is
usually
larger
than
the
actual
along‐strike
size
of
an
earthquake.
If
we
take
an
analytical
elastic
crack
model—an
elliptical
slip
distribution
along
the
strike
with
a
maximum
slip
of
U
0
at
the
rupture
center
and
extending
for
a
length
of
2
̃
L
(Scholz,
2019
),
our
boundary
condition
translates
into
a
slip
gradient
of
U
0
L
at
a
distance
of
the
domain
width
(
D
=
5
km)
from
the
rupture
center.
By
equating
the
gradient
in
this
analytical
solution
(
U
0
D
̃
L
̅̅̅̅̅̅̅̅̅̅
̅
̃
L
2
D
2
U
0
D
̃
L
2
)
with
the
boundary
condition,
we
obtain:
̃
L
̅̅̅̅̅̅
̅
LD
.
(
4
)
Despite
uncertainties
from
geometric
complexity
and
material
heterogeneity,
this
relation
converts
the
parameter
L
in
Equation
3
to
a
rupture
length
scale
̃
L
.
In
realistic
scenarios,
the
rupture
length
can
be
controlled
by
the
lithosphere
thickness
and
the
along‐strike
variability
in
the
plate
geometry,
and
investigation
of
these
compli-
cations
will
require
3D
models
that
incorporate
such
along‐strike
changes.
3.
Results
With
a
realistic
thermal
structure
derived
from
the
cooling
of
plates
and
slabs
(Figure
1a
)
and
a
set
of
experimentally‐constrained
creeping
and
frictional
parameters
(Table
S1
in
Supporting
Information
S1
),
the
slow
motion
of
plates
emerges
but
is
disrupted
by
abrupt
slips
of
10
rupture
events
in
2,000
years,
with
a
period
of
200
years
(Figure
2
).
For
the
entire
incoming
plate,
the
average
surface
displacements
over
these
seismic
cycles
has
the
incoming
plate
moving
at
a
steady
velocity
of
5
cm/year,
consistent
with
the
initial
thermal
structure
as
well
as
the
long‐term
plate
velocity
from
viscously
dominated
computations
with
a
large
time
step
size.
The
surface
velocity
shows
different
patterns
co‐
and
inter‐seismically
(Figure
2a
):
Prior
to
a
rupture
event,
the
surface
velocity
on
both
incoming
and
overriding
plates
show
in‐plane
compression,
which
quickly
flips
to
tension
during
the
event.
On
the
incoming
plate,
the
pace
at
locations
near
the
trench
is
up
to
1
cm
/year
slower
than
the
average
plate
velocity
between
rupture
events,
but
accelerates
substantially
in
areas
near
the
trench
during
a
coseismic
phase,
leading
to
a
leap
in
meters
of
displacement
(Figure
2c
).
The
plate
maintains
slow
velocities
(equal
to
the
long‐term
average)
at
a
distance
larger
than
1,500
km
from
the
trench.
The
overriding
plate
is
compressed
by
the
advancing
incoming
plate
and
areas
near
the
trench
move
slightly
forward
at
2–3
cm
/year.
During
the
rupture
event,
the
overriding
plate
becomes
decoupled
and
experiences
rebound
of
more
than
1
m
in
the
opposite
di-
rection.
Both
the
interseismic
movement
and
coseismic
rebound
diminish
with
distance
from
the
trench.
The
results
generally
align
with
the
elastic
rebound
theory
where
the
elastic
strain
is
gradually
built
up
by
interseismic
displacements
and
abruptly
released
by
coseismic
slips,
despite
non‐negligible
permanent
viscous
deformation
on
both
incoming
and
overriding
plates.
Over
the
2000
years
of
the
model,
the
incoming
plate
shows
an
extension
of
2
m
(1
km
/Myrs)
and
the
overriding
plate
is
compressed
by
0.5
m
(0.25
km
/Myrs)
at
locations
100
km
away
from
the
trench,
compared
with
the
average
plate
motion
(Figure
2c
).
Changes
in
far‐field
boundary
conditions,
rheological
parameters
and
plate
geometry
do
not
influence
first‐order
model
outcomes,
although
the
interseismic
period
and
coseismic
slip
may
vary
when
convergence
rates
are
different
(Figures
S1–S8
in
Supporting
Information
S1
).
The
space–time
variations
of
slip
(Figure
2b
)
and
stresses
(Figure
2d
)
show
cyclical
patterns.
The
coseismic
stress
changes
near
the
megathrust
are
further
shown
in
Figures
S9–S11
in
Supporting
Information
S1
.
In
the
inter-
seismic
phase,
the
megathrust
remains
locked
in
the
section
shallower
than
30
km
and
the
portion
deeper
than
Geophysical
Research
Letters
10.1029/2024GL110821
FANG
ET
AL.
4
of
10
19448007, 2024, 22, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2024GL110821 by California Inst of Technology, Wiley Online Library on [27/11/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
40
km,
where
the
temperature
is
higher,
steadily
creeps
at
a
rate
close
to
the
plate
velocity.
The
main
thrust
event
is
characterized
by
extensive
slip
triggered
along
the
megathrust,
with
a
maximum
of
13
m
at
a
depth
of
10–
30
km.
The
section
at
a
depth
of
8–40
km
can
be
defined
as
the
seismogenic
zone,
where
stresses
accumulate
during
the
interseismic
period
and
suddenly
released
by
rupture.
With
high
pore
fluid
pressure
imposed,
shear
stresses
in
the
seismogenic
zone
are
10–30
MPa
and
the
maximum
coseismic
drop
is
8
MPa
(Figure
2d
).
Conversely,
there
is
a
prominent
increase
in
deviatoric
stress
in
the
updip
and
downdip
portions
of
the
shear
zone.
Due
to
a
rate‐strengthening
friction
below
380
K
(depth
<
7
km),
stresses
in
the
stable‐sliding
updip
section
exhibits
an
increase
of
0.6
MPa
during
rupture
but
drop
immediately
afterward
due
to
notable
afterslip.
The
coseismic
stress
increase
is
7
MPa
at
the
downdip
rupture
tip
and
decays
exponentially
when
the
temperature
increases
with
depth,
although
the
rate‐weakening
parameter
setting
continues
to
a
depth
of
50
km
(700
K).
In
this
section,
the
loaded
stresses
gradually
relax
after
the
rupture.
This
depth
variation
underscores
a
brittle‐to‐ductile
Figure
2.
(a)
The
surface
horizontal
velocity,
v
x
,
along
a
profile
perpendicular
to
the
trench
at
an
interseismic
time
step
(blue
solid),
at
a
coseismic
time
step
(red
solid),
and
averaged
over
the
time
span
of
2,000
years
(black
dashed).
As
the
time
step
size
is
5
years,
coseismic
velocities
represent
integrated
coseismic
displacements
over
a
step
of
5
years.
(b)
The
slip
velocity
V
at
an
interseismic
time
step
(blue)
and
the
slip
U
at
a
coseismic
time
step
(red)
in
the
shear
zone
versus
the
depth
z
.
(c)
The
evolution
of
the
mean
slip
velocity
in
the
shear
zone,
V
mean
,
is
depicted
by
the
teal
solid
line.
Each
peak
represents
a
rupture
event.
The
cumulative
displacement
relative
to
the
average
plate
motion,
Δ
u
x
,
is
plotted
with
magenta
solid
(at
x
=
100
km,
on
the
incoming
plate)
and
magenta
dashed
(
x
=
+
100
km,
on
the
overriding
plate)
lines.
(d)
Space–time
pattern
of
the
relative
stress
change,
Δ
τ
I
I
,
in
the
shear
zone.
Δ
τ
I
I
is
calculated
by
subtracting
the
2000‐year
mean
of
τ
I
I
at
each
depth
from
the
second
invariant
of
the
deviatoric
stress
τ
I
I
.
Geophysical
Research
Letters
10.1029/2024GL110821
FANG
ET
AL.
5
of
10
19448007, 2024, 22, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2024GL110821 by California Inst of Technology, Wiley Online Library on [27/11/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License
transition:
While
plastic
strength
limits
the
stresses
in
the
cold
section,
viscous
strength
regulates
the
shear
zone
deformation
at
temperatures
above
620
K.
When
a
traction
(proportional
to
the
long‐wavelength
displacement
difference
between
the
rupture
center
and
distal
segments)
is
applied
on
the
front
wall,
the
rupture
length
scale,
̃
L
,
modulates
the
size
and
magnitude
of
rupture
(Figures
3
and
4
).
These
results
are
generally
insensitive
to
the
number
of
elements
and
the
domain
width
in
the
along‐strike
direction
(see
Figures
S12–S13
in
Supporting
Information
S1
).
When
̃
L
1000
km
(
L
2
×
10
5
km),
the
rupture
behavior
aligns
with
results
in
2D
(effectively
̃
L
=
+
),
because
the
along‐strike
resistance
associated
with
such
a
large
scale
is
negligible.
Under
these
conditions,
plates
can
experience
accel-
eration
across
spans
greater
than
1,000
km
from
the
trench,
with
coseismic
slips
reaching
up
to
13
m.
As
̃
L
de-
creases,
resisting
forces
become
stronger,
reducing
slips
along
the
entire
megathrust,
with
the
maximum
slip
dropping
from
13
to
4
m
and
near‐surface
displacements
vanishing
from
6
m
to
almost
zero
(Figure
3
).
The
extent
of
plate
acceleration
narrows
from
a
range
of
>
1000
km
to
a
confined
zone
close
to
the
trench.
On
the
incoming
plate,
areas
near
the
trench
exhibit
a
coseismic
displacement
of
1.8
m
for
̃
L
=
1000
km
(
L
=
2
×
10
5
km),
which
converges
to
0.2
m
for
̃
L
=
32
km
(
L
=
2
×
10
2
km),
which
is
the
interseismic
movement
over
each
compu-
tational
step
of
5
years.
The
coseismic
rebound
on
the
overriding
plate
reaches
7
m
for
̃
L
=
1000
km
but
gradually
decreases
for
smaller
along‐strike
length
scales.
When
̃
L
32
km,
these
displacements
toward
the
trench
are
eventually
impeded
at
the
surface,
and
ruptures
are
restricted
near
the
shear
zone
at
a
depth
of
8–40
km
(380–
620
K)
with
limited
slip
(Figure
3c
).
The
reduction
in
both
the
incoming
plate
displacement
and
the
maximum
slip
is
most
evident
when
̃
L
decreases
in
the
critical
range
between
700
and
20
km,
which
is
mapped
from
Equation
4
for
L
decreasing
from
10
5
to
10
2
km
(Figure
4
).
This
range
represents
a
critical
threshold
for
a
pronounced
in-
fluence
of
̃
L
on
the
rupture
behavior.
The
displacements
at
the
surface
tend
to
converge
to
pre‐rupture
values
for
Figure
3.
(a)
The
coseismic
surface
displacement,
u
x
,
along
a
profile
perpendicular
to
the
trench
for
cases
with
different
rupture
length
scales
̃
L
.
(b)
The
coseismic
slip
U
in
the
shear
zone
versus
the
depth
z
for
the
same
cases.
(c)
Coseismic
displacement
near
the
trench
in
representative
events
with
different
values
of
̃
L
.
Silver
dashed
lines
show
the
location
of
the
shear
zone.
Geophysical
Research
Letters
10.1029/2024GL110821
FANG
ET
AL.
6
of
10
19448007, 2024, 22, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2024GL110821 by California Inst of Technology, Wiley Online Library on [27/11/2024]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License