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