of 10
A
Giant
Impact
Origin
for
the
First
Subduction
on
Earth
Qian
Yuan
1,2
,
Michael
Gurnis
2
,
Paul
D.
Asimow
1
,
and
Yida
Li
2
1
Division
of
Geological
and
Planetary
Sciences,
California
Institute
of
Technology,
Pasadena,
CA,
USA,
2
Seismological
Laboratory,
California
Institute
of
Technology,
Pasadena,
CA,
USA
Abstract
Hadean
zircons
provide
a
potential
record
of
Earth's
earliest
subduction
4.3
billion
years
ago.
It
remains
enigmatic
how
subduction
could
be
initiated
so
soon
after
the
presumably
Moon‐forming
giant
impact
(MGI).
Earlier
studies
found
an
increase
in
Earth's
core‐mantle
boundary
(CMB)
temperature
due
to
the
accumulation
of
the
impactor's
core,
and
our
recent
work
shows
Earth's
lower
mantle
remains
largely
solid,
with
some
of
the
impactor's
mantle
potentially
surviving
as
the
large
low‐shear
velocity
provinces
(LLSVPs).
Here,
we
show
that
a
hot
post‐impact
CMB
drives
the
initiation
of
strong
mantle
plumes
that
can
induce
subduction
initiation
200
Myr
after
the
MGI.
2D
and
3D
thermomechanical
computations
show
that
a
high
CMB
temperature
is
the
primary
factor
triggering
early
subduction,
with
enrichment
of
heat‐producing
elements
in
LLSVPs
as
another
potential
factor.
The
models
link
the
earliest
subduction
to
the
MGI
with
implications
for
understanding
the
diverse
tectonic
regimes
of
rocky
planets.
Plain
Language
Summary
Plate
tectonics
remains
unique
to
Earth,
but
when
and
how
it
started
is
debated.
Earth's
oldest
minerals
indicate
a
clement
surface
by
4.3
Ga,
resembling
the
modern
Earth
with
its
granitic
crust
and
oceans.
Granite
is
most
easily
explained
as
originating
from
subduction.
However,
the
mechanisms
for
subduction
initiation,
especially
so
soon
after
the
Moon‐forming
giant
impact,
remain
elusive.
Earlier
studies
indicate
that
the
core‐mantle
boundary
(CMB)
temperature
is
increased
due
to
accumulation
of
the
impactor's
core
during
the
impact.
Our
recent
work
further
shows
that
the
lower
half
of
Earth's
mantle
remains
mostly
solid
after
this
impact
and
that
parts
of
the
impactor's
mantle
might
have
survived
as
the
two
seismically‐observed
large
low‐shear
velocity
provinces
(LLSVPs).
In
this
study,
we
perform
whole‐mantle
convection
models
to
illustrate
that
strong
mantle
plumes
can
arise,
weaken
the
lithosphere,
and
eventually
initiate
subduction
200
Myr
after
the
giant
impact.
Our
systematic
computations
reveal
that
the
hot
CMB
temperature
after
the
impact
is
the
primary
factor
determining
whether
there
is
early
subduction
initiation,
with
enrichment
of
heat‐producing
elements
in
LLSVPs
as
another
potential
contributor.
Our
model
ties
the
Moon's
formation
to
incipient
subduction,
providing
insights
for
understanding
the
diverse
tectonic
regimes
of
rocky
planets.
1.
Introduction
As
planetary
exploration
continues,
Earth
remains
unique
as
the
only
planet
that
is
known
to
operate
in
the
plate
tectonic
regime.
In
plate
tectonics,
subduction
is
the
major
process
affecting
the
geodynamic
and
geochemical
evolution
of
our
planet,
but
when
and
how
the
first
subduction
occurred
remains
contentious.
Although
their
interpretation
is
debated
(Kemp
et
al.,
2010
),
analyses
of
Hadean
detrital
zircons
show
geochemical
signals
consistent
with
subduction
as
early
as
4.3
Ga
(Valley
et
al.,
2002
;
Watson
&
Harrison,
2005
).
Thus,
the
earliest
subduction
may
have
followed
less
than
200
million
years
after
the
Moon‐forming
giant
impact
(MGI)
(Canup
&
Asphaug,
2001
).
The
collision
of
a
protoplanet
called
Theia
with
the
proto‐Earth
at
4.51
Ga
is
expected
to
largely
reset
the
initial
conditions
of
Earth's
evolution,
creating
a
need
for
systematic
and
quantitative
exploration
of
the
planetary
state
established
by
MGI
and
its
lasting
influence
on
Earth's
tectonic
evolution.
Here
we
propose
that
the
thermal
and
chemical
structure
established
by
a
canonical
MGI
creates
favorable
conditions
for
the
initiation
of
subduction
in
Earth's
earliest
period,
exemplifying
the
pivotal
influence
of
initial
conditions
set
by
giant
impact
processes
for
the
tectonic
evolution
of
terrestrial
planets.
A
number
of
mechanisms
have
been
proposed
to
initiate
early
subduction
(Stern
&
Gerya,
2018
),
but
all
except
three
of
these
models
require
a
pre‐existing
weakness
in
the
lithosphere.
The
exceptions
are
grainsize
evolution
(Foley
et
al.,
2014
),
impact‐driven
(Borgeat
&
Tackley,
2022
;
O’Neill
et
al.,
2017
)
and
plume‐induced
(Gerya
et
al.,
2015
;
Ueda
et
al.,
2008
)
subduction
initiation.
However,
large
(
>
700
km)
and
energetic
impacting
bolides
capable
of
initiating
subduction
are
rare
and
only
occurred
very
early
in
Earth's
history
(Marchi
et
al.,
2014
;
RESEARCH
LETTER
10.1029/2023GL106723
Key
Points:
The
mantle
thermochemical
structure
left
by
the
Moon‐forming
impact
trig-
gers
strong
plumes
that
may
have
initiated
the
first
subduction
The
strong
mantle
plumes
may
rise
from
a
heated
core
or
from
large
low‐
shear
velocity
provinces
enriched
in
heat‐producing
elements
Early
giant
impact
events
may
have
a
profound
influence
on
the
diverse
tectonic
regimes
of
rocky
planets
Supporting
Information:
Supporting
Information
may
be
found
in
the
online
version
of
this
article.
Correspondence
to:
Q.
Yuan,
qyuan@caltech.edu
Citation:
Yuan,
Q.,
Gurnis,
M.,
Asimow,
P.
D.,
&
Li,
Y.
(2024).
A
giant
impact
origin
for
the
first
subduction
on
Earth.
Geophysical
Research
Letters
,
51
,
e2023GL106723.
https://doi.org/10.1029/2023GL106723
Received
10
OCT
2023
Accepted
14
APR
2024
Author
Contributions:
Conceptualization:
Qian
Yuan
Formal
analysis:
Qian
Yuan
Funding
acquisition:
Qian
Yuan,
Michael
Gurnis
Investigation:
Qian
Yuan,
Yida
Li
Methodology:
Qian
Yuan
Software:
Qian
Yuan
Supervision:
Michael
Gurnis,
Paul
D.
Asimow
Validation:
Qian
Yuan,
Michael
Gurnis,
Paul
D.
Asimow
Visualization:
Qian
Yuan
Writing
original
draft:
Qian
Yuan
Writing
review
&
editing:
Qian
Yuan,
Michael
Gurnis,
Paul
D.
Asimow,
Yida
Li
©
2024.
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.
YUAN
ET
AL.
1
of
10
O’Neill
et
al.,
2017
),
which
may
coincide
with
the
period
of
magma
ocean
solidification.
Recently,
plume‐
induced
subduction
has
been
increasingly
proposed
as
a
viable
mechanism
for
subduction
initiation
in
the
early
Earth
(Baes
et
al.,
2021
).
Nevertheless,
it
is
unclear
how
to
generate
strong
mantle
plumes
in
a
vigorous
convecting
Hadean
mantle
(Korenaga,
2008
).
The
MGI
is
typically
thought
violent
enough
to
melt
a
substantial
portion
or
even
the
entire
mantle,
especially
in
scenarios
involving
high
energy
and
high
angular
momentum
(Lock
&
Stewart,
2017
;
Nakajima
&
Steven-
son,
2015
).
Using
two
different
giant
impact
computational
methods
with
improved
equations
of
state
(Stewart
et
al.,
2020
)
and
at
high
(Deng
et
al.,
2019
)
and
ultra‐high
(Kegerreis
et
al.,
2022
)
resolution,
our
recent
work
(Yuan
et
al.,
2023
)
demonstrates
that
the
lower
half
of
Earth's
mantle
would
have
remained
mostly
solid
after
a
canonical
MGI.
It
also
suggests
that
large
intact
domains
of
Theia's
likely
Fe‐rich
mantle
are
candidates
of
the
two
seismically‐observed
large
low‐shear
velocity
provinces
(LLSVPs).
Furthermore,
the
concomitant
core
formation
process
during
the
MGI
may
substantially
heat
the
core
to
raise
the
temperatures
at
the
core‐mantle
boundary
(CMB)
(Canup,
2004
;
Deng
et
al.,
2019
).
This
unique
thermochemical
structure,
featuring
a
solid
lower
mantle
and
a
hotter
CMB
resulting
from
the
MGI,
provides
ideal
conditions
for
the
formation
of
strong
mantle
plumes,
a
scenario
less
likely
when
the
mantle
and
core
are
in
equilibrium.
Accretion
scenarios
built
from
only
smaller‐scale
impacts
dissipate
more
of
their
gravitational
potential
energy
in
the
mantle
and
yield
a
more
stable
initial
thermal
structure.
Here,
we
test
the
hypothesis
that
such
strong
mantle
plumes
can
be
developed
in
the
early
Hadean,
weakening
the
lithosphere
and
eventually
causing
subduction
initiation
(Figure
1
).
We
model
this
scenario
using
2‐D
and
3‐D
whole‐mantle
thermo‐chemical‐mechanical
models
with
a
visco‐elasto‐plastic
rock
rheology.
2.
Geodynamic
Models
We
compute
whole‐mantle
convection
models
in
2D
and
3D
Cartesian
geometries
using
Underworld2
,
an
open‐
source,
particle‐in‐cell,
finite
element
code
(Mansour
et
al.,
2020
;
Moresi
et
al.,
2007
).
2.1.
Governing
Equations
and
Material
Description
The
computations
solve
the
conservation
equations
of
mass,
momentum,
and
energy:
·
u
=
0
(
1
)
·
σ
P
+
ρg
ˆ
z
=
0
(
2
)
ρc
p
(
T
t
+
u
·
T
)
·
k
T
H
=
0
.
(
3
)
Figure
1.
A
schematic
illustration
of
the
path
from
the
Moon‐forming
giant
impact
(MGI)
to
plume‐induced
subduction
in
the
early
Hadean.
A
canonical
MGI
(a)
leads
to
a
two‐layered
mantle
structure
with
a
mostly
solid
lower
mantle
and
an
elevated
core‐mantle
boundary
temperature
(b)
(Deng
et
al.,
2019
;
Yuan
et
al.,
2023
).
After
solidification
of
the
upper
magma
ocean
within
several
million
years
(Hamano
et
al.,
2013
;
Nikolaou
et
al.,
2019
),
we
assume
a
proto‐lithosphere
gradually
formed
over
50–100
million
years
(c),
which
was
then
destabilized
and
segmented
by
a
mantle
plume
from
an
large
low‐shear
velocity
province
made
of
the
Theia
mantle
remnant
(Yuan
et
al.,
2023
)
(d).
Geophysical
Research
Letters
10.1029/2023GL106723
YUAN
ET
AL.
2
of
10
19448007, 2024, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023GL106723 by California Inst of Technology, Wiley Online Library on [31/05/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
u
is
velocity,
σ
is
deviatoric
stress,
P
is
pressure,
T
is
temperature,
ρ
is
density,
g
is
gravitational
accel-
eration,
ˆ
z
is
the
vertical
unit
vector,
c
p
is
heat
capacity,
k
is
thermal
conductivity
and
H
is
heat
production.
An
incompressible
Maxwell‐rheological
model
is
adopted
to
reflect
viscoelastic
behavior:
̇
ε
ij
=
1
2
K
e
̇
σ
ij
+
1
2
η
σ
ij
(
4
)
where
̇
ε
ij
is
the
strain
rate
tensor,
̇
σ
ij
is
the
stress
rate
tensor,
K
e
is
the
shear
modulus
for
elasticity,
and
η
is
the
viscosity.
Non‐Newtonian
and
temperature‐dependent
rheology
is
specified
by:
η
(
T
,
̇
ε
)
=
η
0
e
(
E
nRT
E
nRT
0
)
(
̇
ε
̇
ε
0
)
1
n
n
(
5
)
where
E
is
the
activation
energy,
R
is
the
gas
constant,
and
n
is
the
exponent
for
shear
thinning
in
dislocation
creep.
Material
yielding
is
included
by
imposing
an
upper
limit
on
the
stress,
thus
the
effective
viscosity
becomes
η
ef
f
=
min
(
η
(
T
,
̇
ε
)
,
τ
y
̇
ε
I
I
)
(
6
)
where
τ
y
is
the
yield
stress
and
̇
ε
I
I
is
the
second
invariant
of
the
strain
rate
tensor.
The
effective
viscosity
is
restricted
to
the
range
between
10
18
and
10
24
Pa
s.
The
yield
stress
follows
the
Drucker–Prager
yielding
criteria
with
an
upper
cut‐off:
τ
y
=
min
C
+
μ
y
p
,
τ
y
0
)
(
7
)
where
τ
y
is
the
yield
stress,
C
and
μ
y
are
cohesion
and
friction
coefficient,
p
is
pressure,
and
τ
y
0
is
the
maximum
yield
stress.
Similar
to
previous
studies
(Li
&
Gurnis,
2022
),
weakening
processes
are
approximated
by
reducing
the
yield
stress
with
plastic
strain,
following
a
two‐stage
process:
prior
to
a
strain
saturation,
C
and
μ
y
linearly
decrease
as
follows:
C
=
(
τ
yf
C
0
)
min
(
ε
P
ε
P
0
,
1
)
+
C
0
(
8
)
μ
y
=
μ
y
0
μ
y
0
min
(
ε
P
ε
P
0
,
1
)
(
9
)
where
the
variables
are
the
accumulated
plastic
strain
ε
P
,
the
reference
plastic
strain
ε
P
0
that
determines
the
rate
of
weakening
(thus,
the
inverse
of
ε
P
0
is
the
weakening
rate),
minimum
yield
stress
τ
yf
,
initial
cohesion
C
0
,
and
initial
friction
coefficient
μ
y
0
.
Three
materials
are
present
within
the
model
domain:
a
30‐km‐thick
mafic
crust
on
top
that
can
metamorphose
to
eclogite
at
depth
(Hacker
et
al.,
2003
),
LLSVP
material,
and
background
mantle.
2.2.
Model
Setup
in
2D
Most
2D
models
have
a
domain
that
is
5,780
×
2,890
km
with
256
×
128
linear,
quadrilateral
elements,
with
30
Lagrangian
particles
in
each
element.
We
have
verified
the
validity
of
2D
results
with
a
higher‐resolution
grid
of
512
×
256
elements,
and
in
longer
aspect
ratio
of
14,450
×
2,890
km
with
640
×
128
elements.
A
free‐slip
boundary
condition
is
applied
to
all
model
domain
boundaries.
The
temperature
at
the
upper
boundary
is
maintained
at
273
K,
and
the
bottom
boundary
is
fixed
at
a
varying
temperature
(5273
K
down
to
3273
K).
The
lithosphere's
maximum
yield
strength
in
2D
is
set
between
80
and
180
MPa.
As
seismic
observations
require
that
Geophysical
Research
Letters
10.1029/2023GL106723
YUAN
ET
AL.
3
of
10
19448007, 2024, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023GL106723 by California Inst of Technology, Wiley Online Library on [31/05/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
LLSVPs
have
a
distinct
bulk
modulus
(Tan
&
Gurnis,
2005
),
we
use
a
depth‐dependent
density
profile
within
the
LLSVP.
The
values
of
physical
parameters
are
listed
in
Table
S1
of
Supporting
Information
S1
.
2.3.
Model
Setup
in
3D
Three‐dimensional
plume‐induced
subduction
computations
are
used
to
verify
the
2D
results
and
consist
of
a
Cartesian
domain
2,890
km
deep,
2,000
×
2,000
km
horizontally.
We
use
a
resolution
of
128
×
80
×
80
elements
with
a
uniform
grid
spacing
in
each
coordinate
direction.
Each
element
initially
has
30
material
points
randomly
located
in
each
element.
Free
slip
conditions
are
applied
to
all
the
boundaries.
Maximum
yield
stresses
of
lith-
osphere
between
60
and
100
MPa
were
tested
in
3D.
3.
Results
Following
the
MGI,
a
surface
magma
ocean
extending
to
a
certain
depth
is
expected.
Tracking
magma
ocean
solidification
is
beyond
the
scope
of
this
study,
but
it
is
generally
accepted
that
the
solidification
process
would
be
mostly
complete
within
several
million
years
(Hamano
et
al.,
2013
;
Nikolaou
et
al.,
2019
).
For
simplicity,
our
models
begin
50
Myr
after
the
magma
ocean
solidification;
the
initial
state
is
derived
from
a
2D
thermo-
mechanical
model
with
a
proto‐lithosphere
defined
by
a
half‐space
cooling
model
with
a
prescribed
age
of
50
Ma.
Most
model
runs
include
an
intrinsically
dense
layer
of
Theia
mantle
materials
above
the
CMB.
Current
estimates
of
the
present‐day
CMB
temperature
that
account
for
a
lower
boundary
layer
range
from
3473
to
4573
K
(Deschamps
&
Cobden,
2022
),
but
because
of
Earth's
secular
cooling,
the
Hadean
CMB
temperature
should
have
been
higher.
Moreover,
a
hot
CMB
is
anticipated
after
the
MGI
due
to
the
deposition
of
Theia's
iron
core,
as
shown
in
previous
MGI
simulations
(Canup,
2004
).
Consequently,
for
most
cases
we
used
a
higher
CMB
tem-
perature
(5273
K)
and
a
boundary‐layer
contrast
of
1273
K
between
the
CMB
and
the
initial
LLSVP
temperature.
Such
high
CMB
temperature
may
lead
to
partial
or
complete
melting
in
the
boundary
layer,
but
temperature
should
drop
quickly
below
the
solidus
of
likely
lower
mantle
compositions
above
the
CMB.
A
thin
partial
melting
zone
within
the
mantle
does
not
substantially
alter
the
dynamics,
since
the
fully
molten
core
sets
the
lower
free
slip
boundary
condition
in
all
cases.
Furthermore,
we
note
that
similar
or
higher
CMB
temperatures
were
used
in
previous
convection
simulations
(Nakagawa
&
Tackley,
2010
),
and
we
confirm
our
results
with
lower
CMB
temperatures
(4273
and
3773
K).
We
mostly
used
a
80
MPa
yield
strength
for
the
proto‐lithosphere
(Rudi
et
al.,
2022
).
Model
runs
with
this
initial
configuration
consistently
show
that
a
mantle
plume
nucleates
atop
the
LLSVP‐like
thermochemical
pile,
ascends
through
the
mantle,
ruptures
the
lithosphere,
and
initiates
subduction
after
120
Myr
(Figure
2
).
In
detail,
typical
model
evolution
can
be
divided
into
several
stages.
First,
the
initial
layer
of
intrinsically
dense
material
(Figure
2a
)
at
the
CMB
concentrates
into
discrete
piles
due
to
a
combination
of
strong
CMB
heating
and
its
depth‐dependent
density
(Tan
&
Gurnis,
2005
).
Meanwhile,
the
lithosphere
grows
thicker
due
to
inefficient
cooling
while
in
a
stagnant‐lid
convection,
although
the
mode
of
magmatism
may
affect
a
planet's
tectonic
regime
(Lourenço
et
al.,
2018
).
At
99
Myr,
a
strong
mantle
plume
develops
from
the
top
of
the
central
pile
(Figure
2b
).
This
plume
quickly
rises
to
the
upper
mantle,
becoming
narrower
due
to
the
viscosity
decrease
in
the
transition
zone.
As
the
plume
approaches
the
surface,
the
lithosphere
is
locally
weakened
through
yielding
above
the
plume
and,
at
109
Myr,
the
persistent
localized
thinning
leads
to
plume
penetration
through
the
lithosphere,
forming
a
high‐temperature,
low‐viscosity
wedge
structure
(Figure
2c
).
Meanwhile,
the
proximal
lithosphere
thickens
around
the
plume
head,
resulting
in
two
lithospheric‐scale
shear
bands
on
each
side
(Figure
2c
).
The
plume
then
starts
to
spread
near
the
surface,
overcoming
the
yield
strength
of
the
surrounding
lithosphere
segments
and
forming
shear
zones
along
the
subducting
plate
boundary
that
promote
self‐sustaining
and
retreating
subduction
at
123
Myr
(Figure
2d
).
Our
reference
case
above
shows
that
strong
plumes
can
arise
from
the
LLSVPs
and
initiate
subduction
when
the
lithosphere's
strength
is
less
than
100
MPa,
a
threshold
yield
stress
consistent
with
previous
geodynamic
studies
(Moresi
&
Solomatov,
1998
;
Van
Heck
&
Tackley,
2008
).
Here,
we
aim
to
explore
the
maximum
lithosphere
strength
that
the
MGI‐caused
mantle
plumes
can
breach
to
initiate
subduction
and
assess
the
influence
of
LLSVPs
on
this
capability
(Figure
3
).
Initially,
we
increase
the
lithosphere
yield
strength
from
80
to
100
MPa
and
find
that
it
becomes
more
challenging
to
induce
subduction;
subduction
only
occurs
on
one
side
of
the
plume
head
(Figure
3a
).
When
the
lithosphere
yield
stress
is
further
increased
to
150
MPa,
plumes
from
LLSVPs
do
not
penetrate
through
the
lithosphere
to
induce
subduction
but
instead
spread
at
the
base
of
the
lithosphere
Geophysical
Research
Letters
10.1029/2023GL106723
YUAN
ET
AL.
4
of
10
19448007, 2024, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023GL106723 by California Inst of Technology, Wiley Online Library on [31/05/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
(Figure
3b
).
Next,
for
the
same
calculation
after
removing
LLSVPs
materials,
we
find
that
strong
plumes
rising
from
the
CMB
can
initiate
subduction,
as
shown
in
Figure
3c
for
a
case
with
a
yield
stress
of
150
MPa,
and
in
another
case
with
a
yield
stress
of
180
MPa
(Figure
3d
).
As
studied
in
previous
research
on
plume‐induced
subduction,
temperature,
size,
and
buoyancy
of
plumes
play
a
major
role
in
subduction
initiation.
Therefore,
we
systematically
explore
the
influence
of
CMB
temperature,
which
significantly
affects
all
these
factors
in
models
where
plumes
are
self‐consistently
generated.
We
lower
the
CMB
temperature
to
4773
K
(Figure
3e
)
and
3773
K
(Figure
3f
),
and
find
an
increase
in
the
time
delay
for
plume‐
induced
subduction
initiation
(Figures
3i
and
3j
)
due
to
the
lower
buoyancy
of
the
plumes.
Nonetheless,
sub-
duction
still
occurs
as
long
as
CMB
temperature
is
3773
K.
When
CMB
temperature
is
3273
K,
no
subduction
is
initiated
(Figures
3g
and
3k
).
However,
if
we
assume
that
LLSVP
materials
are
enriched
in
heat‐producing
el-
ements
for
the
same
calculation
(Citron
et
al.,
2020
),
the
results
indicate
that
plumes
rising
from
hotter
LLSVPs
can
still
generate
subduction
(Figures
3h
and
3l
).
Figure
2.
Temporal
evolution
for
the
temperature
(left
column)
and
viscosity
(right
column)
fields
of
the
reference
case
showing
a
LLSVP‐sourced
plume
inducing
subduction
initiation.
(a)
Initial
set
up
of
the
model;
(b)
Plume
rising
from
the
LLSVP;
(c)
Plume
penetrates
the
lithosphere;
(d)
Slab
retreat
and
melting
of
crust.
The
yellow
contours
in
temperature
field
marks
the
boundary
of
different
materials
including
crust,
lithosphere,
molten
materials,
and
large
low‐shear
velocity
province
materials.
Geophysical
Research
Letters
10.1029/2023GL106723
YUAN
ET
AL.
5
of
10
19448007, 2024, 9, Downloaded from https://agupubs.onlinelibrary.wiley.com/doi/10.1029/2023GL106723 by California Inst of Technology, Wiley Online Library on [31/05/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