Confidential manuscript submitted to
GRL
Supporting Information for
“Subduction initiation with vertical lithospheric heterogeneities and new
fault formation”
Xiaolin Mao
1
, Michael Gurnis
1
, Dave A. May
2
1
Seismological Laboratory, California Institute of Technology, Pasadena, CA, 91125, USA
2
Department of Earth Sciences, University of Oxford, Oxford, United Kingdom
Contents
1. Method details
2. Tables S1 to S2
3. Figures S1 to S8
Method details
Surface process
To incorporate the influence of surface processes on subduction zone evolution, we
consider a simple 2D law of topographic changes that could simulate erosion and sedimenta-
tion at the scale of several tens of kilometers [
Avouac and Burov
, 1996]. Assuming the rate
of downslope transport of debris,
q
e
, is proportional to the local slope, we have
q
e
=
−
k
∂
h
∂
x
,
(1)
where
k
is the mass diffusivity coefficient. From the mass conservation law,
h
obeys
∂
h
∂
t
=
−
∂
q
e
∂
x
.
(2)
With constant
k
, equations (1) and (2) lead to the linear diffusion equation
∂
h
∂
t
=
k
∂
2
h
∂
x
2
.
(3)
Equation (3) is applied to the free surface between time steps and after the topography
updating related to subsurface dynamics. A layer of sediment tracers are added on the free
surface before applying Eq. (3). These sediment tracers become out of the model domain
Corresponding author: Xiaolin Mao,
xlmao@gps.caltech.edu
–1–
Confidential manuscript submitted to
GRL
and are deleted in the erosion region after the topographic changes, while they survive at the
sedimentation region. The mass diffusivity coefficient is set to be
2
×
10
4
m
2
/yr in our model
[Mao et. al., manuscript in preparation].
Phase change
Discrete 4D (temperature, pressure, rock type and total water content) phase maps for
density and free water content are calculated using Perplex [
Connolly
, 2005] and stored for
lookup during the computation. Four rock types are considered: pyrolite, harzburgite (only
in case SI8), basalt and sediment. The ranges of temperature and pressure cover the whole
upper mantle. The total water content varies from 0 to 10-20 percent. During the compu-
tation, rock type and total water content are recorded on tracers, and temperature and pres-
sure are obtained from solution of the numerical model. After each Stokes solve, density and
free water content on a tracer are updated through refer to stored phase maps for the corre-
sponding rock type. Linear interpolation is used to obtain the value between adjacent nodes
of temperature, pressure and total water content [Mao et. al., manuscript in preparation].
Elastic Viscous Stress Splitting (EVSS)
In our models, we employ the Elastic Viscous Stress Split (EVSS) visco-elastic model
[
Keunings
, 2000]. The EVSS model is built upon the following deviatoric stress decomposi-
tion
τ
=
2
η
v
¤
ε
(
u
)
+
τ
e
,
(4)
where
η
v
is the viscosity of the purely viscous component,
¤
ε
(
u
)
is the strain rate tensor, and
τ
e
is the visco-elastic contribution to the stress tensor. The visco-elastic stress
τ
e
evolves
according to the following
τ
e
+
η
e
G
(
D
τ
e
Dt
+
J
(
u
,
τ
e
))
=
2
η
e
¤
ε
(
u
)
,
(5)
where
η
e
and
G
are the viscosity and shear modulus of the elastic component of the stress
respectively, and
D
τ
e
Dt
+
J
(
u
,
τ
e
)
is the Jaumann corotational stress rate, where
J
(
u
,
τ
e
)
=
τ
e
W
−
W
τ
e
,
(6)
and
W
i j
=
1
2
(
∂
u
i
∂
x
j
−
∂
u
j
∂
x
i
)
.
(7)
–2–
Confidential manuscript submitted to
GRL
The partition of the total viscosity is given by
η
=
(
1
−
φ
)
η
v
+
φη
e
, φ
∈[
0
,
1
]
,
(8)
where
η
is the conventional viscosity, and
φ
is the partition coefficient.
φ
decreases linearly
from close 1 in the lithosphere to 0 in the asthenosphere with increasing temperature, which
results visco-elastic and viscous behaviors in the two regions respectively. When
φ
is set to
1, the EVSS formulation reverts back to the visco-elastic model commonly used in geody-
namics [e.g.,
Moresi et al.
, 2003]. The benefit of using EVSS model is that when time step
continuously decreases,
τ
converges to the stress from the purely viscous part, which helps to
stabilize numerical simulation [Mao et. al., manuscript in preparation].
References
Avouac, J.-P., and E. Burov (1996), Erosion as a driving mechanism of intracontinental
mountain growth,
Journal of Geophysical Research: Solid Earth
,
101
(B8), 17,747–
17,769, doi:10.1029/96JB01344.
Connolly, J. A. (2005), Computation of phase equilibria by linear programming: a tool for
geodynamic modeling and its application to subduction zone decarbonation,
Earth and
Planetary Science Letters
,
236
(1), 524–541, doi:10.1016/j.epsl.2005.04.033.
Keunings, R. (2000), Advances in the computer modeling of the flow of polymetric liquids,
Computational Fluid Dynamics Journal
,
9
(1), 449–458.
Moresi, L., F. Dufour, and H.-B. Mühlhaus (2003), A lagrangian integration point finite el-
ement method for large deformation modeling of viscoelastic geomaterials,
Journal of
Computational Physics
,
184
(2), 476–497.
–3–
Confidential manuscript submitted to
GRL
Tables
Table S1.
Thermal and rheological parameters
Parameter
Value
surface temperature
0
◦
C
bottom temperature
1400
◦
C
thermal diffusivity (
κ
)
10
−
6
m
2
s
−
1
maximum viscosity cutoff (
η
max
)
10
24
Pa s
minimum viscosity cutoff (
η
min
)
10
18
Pa s
shear modulus (
G
)
30 GPa
maximum yield stress (
τ
yield
)
150 MPa
initial water content in mantle
100 ppm
preexponential factor in mantle (
A
m
) 1.6 x 10
−
15
Pa
−
n
s
−
1
exponent in mantle (
n
m
)
3.2
activation energy in mantle (
E
m
)
540 kJ mol
−
1
initial water content in crust
2.68
%
preexponential factor in crust (
A
c
)
2 x 10
−
23
Pa
−
n
s
−
1
exponent in crust (
n
c
)
3.2
activation energy in crust (
E
c
)
238 kJ mol
−
1
initial water content in sediment
7.29
%
preexponential factor in sediment (
A
s
) 5 x 10
−
21
Pa
−
n
s
−
1
exponent in sediment (
n
s
)
2.3
activation energy in sediment (
E
s
)
154 kJ mol
−
1
–4–
Confidential manuscript submitted to
GRL
Table S2.
Table of models. WZ is short for weak zone.
Model WZ dip angle WZ width Harzburgite layer
SI1
90
°
10 km
No
SI2
30
°
10 km
No
SI3
45
°
10 km
No
SI4
60
°
10 km
No
SI5
75
°
10 km
No
SI6
90
°
15 km
No
SI7
90
°
5 km
No
SI8
90
°
10 km
15 km
–5–
Confidential manuscript submitted to
GRL
Figures
Figure S1.
Model results (case SI2). (a) Effective viscosity evolution. (b) Accumulation of plastic strain.
The black box shows the corresponding region for Figure 4b. (c) Density evolution. Black lines show direc-
tion and magnitude of maximum principal stress. Rock types and free water contents are shown in the insets
with different colors and white contours. (d) Topography changes. Blue and red lines show initial and final
topography for each time interval.
–6–
Confidential manuscript submitted to
GRL
Figure S2.
Model results (case SI3). (a) Effective viscosity evolution. (b) Accumulation of plastic strain.
(c) Density evolution. Black lines show direction and magnitude of maximum principal stress. Rock types
and free water contents are shown in the insets with different colors and white contours. (d) Topography
changes. Blue and red lines show initial and final topography for each time interval.
–7–
Confidential manuscript submitted to
GRL
Figure S3.
Model results (case SI4). (a) Effective viscosity evolution. (b) Accumulation of plastic strain.
(c) Density evolution. Black lines show direction and magnitude of maximum principal stress. Rock types
and free water contents are shown in the insets with different colors and white contours. (d) Topography
changes. Blue and red lines show initial and final topography for each time interval.
–8–
Confidential manuscript submitted to
GRL
Figure S4.
Model results (case SI5). (a) Effective viscosity evolution. (b) Accumulation of plastic strain.
(c) Density evolution. Black lines show direction and magnitude of maximum principal stress. Rock types
and free water contents are shown in the insets with different colors and white contours. (d) Topography
changes. Blue and red lines show initial and final topography for each time interval.
–9–
Confidential manuscript submitted to
GRL
Figure S5.
Model results (case SI6). (a) Effective viscosity evolution. (b) Accumulation of plastic strain.
(c) Density evolution. Black lines show direction and magnitude of maximum principal stress. Rock types
and free water contents are shown in the insets with different colors and white contours. (d) Topography
changes. Blue and red lines show initial and final topography for each time interval.
–10–
Confidential manuscript submitted to
GRL
Figure S6.
Model results (case SI7). (a) Effective viscosity evolution. (b) Accumulation of plastic strain.
(c) Density evolution. Black lines show direction and magnitude of maximum principal stress. (d) Topogra-
phy changes. Blue and red lines show initial and final topography for each time interval.
–11–
Confidential manuscript submitted to
GRL
Figure S7.
Along depth variation of initial density structure. hz: harzburgite.
–12–