Earth and Planetary Science Letters 528 (2019) 115836

Contents lists available at ScienceDirect

Earth and Planetary Science Letters

www.elsevier.com/locate/epsl

Where does subduction initiate and cease? A global scale perspective

Martina M. Ulvrova a,∗ , Nicolas Coltice b , Simon Williams c,d , Paul J. Tackley a
a
Institute of Geophysics, Department of Earth Sciences, ETH Zürich, Zürich, Switzerland
b
Laboratoire de Géologie, École Normale Supérieure, CNRS-UMR 8538, PSL Research University, Paris, France
c
EarthByte Group, School of Geosciences, University of Sydney, Sydney, New South Wales, Australia
d
State Key Laboratory of Continental Dynamics, Department of Geology, Northwest University, Xi’an, China

a r t i c l e i n f o a b s t r a c t

Article history: The thermo-mechanical evolution of the Earth’s mantle is largely controlled by the dynamics of
Received 21 December 2018 subduction zones, which connect the surface tectonic plates with the interior. However, little is known
Received in revised form 18 July 2019 about the systematics of where subduction initiates and ceases within the framework of global plate
Accepted 9 September 2019
motions and evolving continental configurations. Here, we investigate where new subduction zones
Available online 2 October 2019
Editor: B. Buffett
preferentially form, and where they endure and cease using statistical analysis of large-scale simulations
of mantle convection that feature self-consistent plate–like lithospheric behaviour and continental drift
Keywords: in the spherical annulus geometry. We juxtapose the results of numerical modelling with subduction
subduction initiation histories retrieved from plate tectonic reconstruction models and from seismic tomography. Numerical
subduction zones models show that subduction initiation is largely controlled by the strength of the lithosphere and
mantle convection by the length of continental margins (for 2D models, the number of continental margins). Strong
plate tectonics lithosphere favours subduction inception in the vicinity of the continents while for weak lithosphere
tectonic reconstructions
the distribution of subduction initiation follows a random process distribution. Reconstructions suggest
numerical simulations
that subduction initiation and cessation on Earth is also not randomly distributed within the oceans, and
more subduction zones cease in the vicinity of continental margins compared to subduction initiation.
Our model results also suggest that intra-oceanic subduction initiation is more prevalent during times
of supercontinent assembly (e.g. Pangea) compared to more recent continental dispersal, consistent with
recent interpretations of relict slabs in seismic tomography.
© 2019 Elsevier B.V. All rights reserved.

1. Introduction respect to the position of the continental margins do subduction

zones cease?
Subduction of the rigid plates is a fundamental process in Earth Few examples of active subduction inception or cessation are
evolution, allowing chemical cycling between the surface and the available to study. Young subduction systems can be found at the
deep mantle (Kerrick, 2001). Indeed, the surface and interior of Mussau Trench (Hegarty et al., 1982) and Yap Trench (Lee, 2004)
the planet are interconnected within a self-organized system in in the western Pacific, but it is not clear if these will develop into
which subduction arises from an instability of the top boundary self-sustained subduction. Much of our knowledge on how subduc-
layer, while it also induces convective currents and pulls tectonic tion initiates and ceases is based on the geological record, includ-
plates (Lowman, 2011; Coltice et al., 2017). The evolution of the ing marine studies of forearcs (Reagan et al., 2010) and on-land
lithospheric plates including continents is then characterized by re- studies of ophiolites (Dilek and Furnes, 2009), and on numerical
peating Wilson cycles during which ocean basins periodically close modelling (e.g. Nikolaeva et al., 2010). These studies indicate that
and open while supercontinents assemble and disperse. However, subduction may initiate in a diverse range of tectonic settings; at
little is known about subduction inception in the setting of global passive margins (Nikolaeva et al., 2010), fracture zones such as for
tectonics with floating continental rafts. How far from the conti- example Aleutian subduction zone (Maffione et al., 2017), at ex-
nents do new subduction zones preferentially form? How do plate tinct spreading centres such as for example Puysegur (Lebrun et
motions influence subduction inception? At which locations with al., 2003), adjacent to fossil island arcs (Leng and Gurnis, 2015),
triggered by plumes (Gerya et al., 2015), or where mantle suction
flow occurs (Baes et al., 2018). Stern (2004) proposes two distinct
mechanisms for subduction initiation: spontaneous nucleation by
*
Corresponding author at: Department of Earth Sciences, Institute of Geophysics,
e.g. foundering at passive margins, or induced initiation involving
ETH Zürich, Sonneggstrasse 5, Zürich, 8092, Switzerland.
E-mail address: martina.ulvrova@erdw.ethz.ch (M.M. Ulvrova). forced convergence of existing plates.

https://doi.org/10.1016/j.epsl.2019.115836
0012-821X/© 2019 Elsevier B.V. All rights reserved.
2 M.M. Ulvrova et al. / Earth and Planetary Science Letters 528 (2019) 115836

Table 1
Dimensional and non-dimensional parameters of the convection model.

Variable Symbol Nondimensional Dimensional

value value
Gravitational acceleration g – 9.81 m s−2
Mantle thickness D 1 2890 km
Thermal expansivity α0 – 3 × 10−5 K−1
Thermal diffusivity κ 1 10−6 m2 s−1
Thermal conductivity k – 4 W m−1 K−1
Gas constant R – 8.314 J mol−1 K−1
Reference density ρ0 1 3300 kg m−3
Maximum viscosity cutoff ηmax 104 6 × 1026 Pa s
Internal heating rate H 20 5.44 × 10−12 W kg−1
Reference viscosity η0 1 6 × 1022 Pa s
Activation energy Ea 8 170 kJ mol−1
Activation volume Va 3 6.34 · 10−7 m3 mol−1
Surface temperature Ts 0.12 300 K
Superadiabatic temperature drop T 1 2500 K
Rayleigh number Ra 106 –
Surface yield stress in oceans σ0 103 to 8 × 103 7 MPa to 56 MPa
Yield stress depth derivative in oceans σY 3.3 × 105 810 Pa m−1

In common with subduction initiation, cessation of subduction 2. Method

has been attributed to a variety of mechanisms, including colli-
sion with continents or oceanic plateaus, interaction between the In order to investigate statistically the spatial relationships be-
subduction zone and spreading ridges and transforms, or within tween continents and the initiation, evolution and cessation of
the context of a broader-scale reorganisation of plate motions (e.g. subduction zones, we numerically calculate the solution of man-
Michaud et al., 2006). Some active subduction systems undergo a tle convection in a spherical annulus (Hernlund and Tackley, 2008)
so-called polarity reversal, when the overriding plate becomes the using the StagYY code (Tackley, 2008). The choice of geometry is
subducting plate and vice versa. In such a case, subduction initi- motivated by the necessity of having long temporal series of sev-
ation and termination are directly related. An example of reversal eral billions of years. Employment of the spherical annulus ensures
of an active convergent boundary is the New Hebrides with the re- similar scaling properties compared to the full 3D spherical shell.
versal of subduction of the Pacific beneath the Australian plate at The model features self-consistently generated plate-like surface
the Vitiaz trench (Auzende et al., 1988). tectonics and drifting continents.
A limiting factor in our current knowledge on subduction is
the reliance on geological evidence collected on land, for exam- 2.1. Physical and numerical model
ple where former intra-oceanic subduction products have been
accreted onto continents. Consequently, it is difficult to constrain We determine temperature, velocity, pressure and composition
where and when intra-oceanic arcs resided throughout their life within the mantle by solving the equations of conservation of
cycle, and tectonic reconstructions are often dominated by subduc- mass, momentum and energy and the advection of material com-
tion systems close to continents (e.g. Müller and Landgrebe, 2012). position considering an incompressible mantle under the Boussi-
Recently, studies mapping slab remnants imaged in seismic tomog- nesq approximation. Below, the equations are given in their di-
raphy (van der Meer et al., 2012; Domeier et al., 2017) point to mensionless form.
the existence of previously unrecognized intra-oceanic subduction
zones within the Pacific/Panthalassa domain that would have been ∇ · v = 0, (1)
active while Pangea was assembled or earlier in its dispersal, and   
much further from the continents than more recent examples. This ∇ · η ∇ v + (∇ v )T − ∇ p = Ra ( T + B C ) e r , (2)
raises the question of how important the presence of continents is 2
∂t T + v · ∇ T = ∇ T + H , (3)
to the life cycle of subduction systems, and whether this influence
varies between periods of supercontinent assembly and continen- ∂t C + v · ∇ C = 0 , (4)
tal dispersal.
with v the velocity, p the dynamic pressure, η the viscosity, T
In this paper, we investigate the pattern of subduction initia-
the temperature, C the composition, H the internal heating rate,
tion and cessation in time and space using numerical simulations
Ra the Rayleigh number, B the buoyancy ratio and e r the radial
of mantle convection. Numerical models are designed in a spher-
unit vector. ∂t is the partial time derivative. Length is scaled by
ical annulus geometry, and we vary the number of continents,
the thickness of the Earth’s mantle D, velocity v by κ / D, with κ
the strength of the lithosphere and its structure. We compare as-
the thermal diffusivity, time t by the thermal diffusion timescale
pects of the modelling results to reconstructed subduction histories
D 2 /κ , and viscosity η by the reference viscosity η0 .
based on both plate kinematics and analysis of slab remnants im-
Viscosity η follows the Arrhenius law and is strongly tempera-
aged by seismic tomography (van der Meer et al., 2010, 2012;
ture and pressure dependent
Müller et al., 2016). We focus on global scale models to infer how
 
the continental configuration and the plate layout change the dis- Ea + p V a
tribution of subduction initiation and cessation, and the lifespan η( T , p ) = η A exp , (5)
RT
of subduction zones. The models indicate that subduction initia-
tion is non-randomly distributed in the ocean, and cessation hap- where E a = 166 kJ mol−1 is the activation energy (kept constant
pens mostly in the vicinity of continents. Our calculations point to for all simulations), V a = 6.34 · 10−7 m3 mol−1 the activation vol-
different distributions of subduction inception between phases of ume (constant for all simulations) and R = 8.314 J mol−1 K−1 the
supercontinents and phases in which continents are dispersed. gas constant. We give all parameters in Table 1. η A is set such that
 M.M. Ulvrova et al. / Earth and Planetary Science Letters 528 (2019) 115836 3

η matches the reference viscosity η0 at zero pressure and at tem- Continental rafts are modelled using the tracer-ratio method
perature 1600 K, which is the expected temperature at the base of (Tackley and King, 2003). The detailed implementation is described
the lithosphere. We apply a viscosity cut off at 104 times η0 . The elsewhere (Tackley and King, 2003; Rolf and Tackley, 2011). We
viscosity varies over 6 orders of magnitude over the temperature consider continents with an interior that is 300 km thick sur-
variation  T , the superadiabatic temperature drop over the man- rounded by 140 km thick mobile belts in accordance with the
tle. Independently of temperature, viscosity increases exponentially thickness of the Archean cratons and Proterozoic belts. Continents
by an order of magnitude with depth. The lowest values of the vis- cover 30% of the model surface. In order to ensure the stability
cosity are in the asthenospheric mantle, while at the core mantle of the continents, two conditions must be fulfilled: positive buoy-
boundary viscosity is about 20 times lower than η0 . ancy and limited deformation within continents (e.g. Doin et al.,
To localize deformation in narrow zones and obtain realistic 1997; Lenardic and Moresi, 1999). Based on this, we choose a den-
plate boundaries at the surface, we use a pseudoplastic rheol- sity contrast between continental material and ambient mantle of
ogy (Moresi and Solomatov, 1998; Tackley, 2000). After reaching −100 kg m−3 , which gives the buoyancy ratio (ratio between the
a certain threshold value, the yield stress σY , the rocks undergo density contrast and the thermal density variation) of −0.4. Fur-
plastic yielding. σY is depth dependent and follows thermore, continents are 100× more viscous than the ambient
mantle and do not undergo pseudo-plastic deformation. Due to
σY = σ0 + d σY , (6) this high strength, continental erosion by mantle flow is negligi-
ble and the rafts are stable over billions of years.
where σ0 is the surface yield stress, d is the depth and σY is the The system is driven by heating from radioactive elements and
yield stress depth derivative. If the stress reaches σY , we calculate from heat conducted from the core. We keep the internal heating
the effective viscosity ηeff on the grid rate H constant at 5.44 × 10−12 W kg−1 throughout the simula-
tions. Core heat loss contributes around 20% to the total surface
ηeff = min [η( T , p ), ηY ] , (7) heat flow, falling into the 10%–40% estimate for the Earth (Jau-
part et al., 2015). Both the top and the bottom boundaries are
with ηY as isothermal and free-slip, except for the case with a free deformable
σY mantle surface, for which the top boundary of the computational
ηY = . (8) domain is no-slip above the air layer.
2 ˙II
The convective vigour of the system is measured by the
˙II is the second invariant of the strain rate tensor. We vary Rayleigh number Ra
the surface strength σ0 of the lithosphere between 7 MPa and
56 MPa while keeping the gradient constant for all simulations at ρ0 g α  T D 3
810 Pa m−1 . Ra = , (9)
κη0
Using this kind of rheology results in the self-consistent for-
mation of strong plate interiors moving with a uniform velocity where ρ0 is the reference density, g the gravitational acceleration,
delimited by narrow plate boundaries characterized by reduced α the thermal expansivity,  T is the superadiabatic temperature
viscosity and an abrupt velocity change (Moresi and Solomatov, drop across the mantle with the thickness D, κ is the thermal
1998; Tackley, 2000). Importantly, such rheology successfully re- diffusivity and η0 is the reference viscosity. For all experiments,
produces seafloor age distributions (Coltice et al., 2012) and is we keep Ra = 106 (calculated with η0 = 6 × 1022 Pa s). This is
predictive enough to investigate global surface tectonics (Coltice 10 − 100× lower than the Earth’s Rayleigh number, but lies at the
et al., 2017; Ulvrova et al., 2019). edge of the computational feasibility for the given rheology and
One of the important Earth-like features that is in general lack- the presence of sticky air.
ing in such convection models is single-sided subduction, i.e., sub- We use a resolution of 128 × 1024 cells in the radial and hor-
duction is often double-sided. This can be partially overcome by izontal directions, respectively, except for the case with the free
imposing a weak layer of oceanic crust at the surface, which gets surface when we use a grid with 256 × 2048 cells. Vertical grid re-
advected into the subduction channel and hence decouples to a finement close to the top and bottom limits is employed resulting
certain degree the sinking slab from the overriding plate, result- in 10 km and 15 km (5 km and 8 km for the case with the sticky
ing in strong subduction asymmetry (Gerya et al., 2008; Crameri air) thick cells at the surface and at the core-mantle boundary (cf.
and Tackley, 2014). In some of our simulations we include this more details in the supplementary material). To track composi-
weak crustal layer using tracers. The weak crustal layer is neu- tion, we use 4 × 107 tracers. This means that on average there
trally buoyant, and it follows the same viscosity law as ambient are around 300 tracers per cell except in simulations with a free
mantle but is characterized by a factor of 10 reduction in η A and it surface, which have higher resolution. In this case, the number of
is more easily deformable (surface yield stress 18 MPa, yield gra- tracers per cell is around 75.
dient 8.1 Pa m−1 ). The initial thickness of the weak crustal layer First, we run a model until it reaches a statistically steady state,
is 20 km and it is converted to a regular mantle after reaching at which the heat budget is balanced and characteristic properties
290 km depth. In one case, we further improve the model and ob- of the system such as mean temperature, mean velocity and av-
tain one-sided subduction zones by imposing a free surface that erage surface heat flux fluctuate around some constant values. In
can vertically deform. This is done using the so-called sticky-air this initial stage, we do not advect tracers and hence the compo-
method (Matsumoto and Tomoda, 1983) when vertical deforma- sition field remains the same throughout the initialization period,
tion of the mantle surface is allowed by prescribing an additional which means that continents do not move from their initial posi-
layer of “air” atop of the mantle. This layer is permanently forced tions. Once statistically steady state is achieved, we use a random
to be isothermal (300 K), has close to zero density and very low snapshot from the equilibrated evolution to start the calculation
viscosity, allowing it to be decoupled from the lithosphere. To ob- with continents that move freely. The choice of the particular
tain a valid solution, the air layer must be sufficiently thick while snapshot from the statistically steady state does not have any in-
having low viscosity (Crameri et al., 2012). We fix its initial thick- fluence on the statistics performed on the system. The model is
ness to 150 km and keep its viscosity at 10−3 times η0 . This results run until a sufficient number of subduction initiations N is col-
in a high viscosity contrast between the plates at the surface and lected. The shortest analyzed period is 1 Gy with 27 subduction
the air layer. The viscosity contrast is as high as 107 . initiation events. The longest analyzed period is 7 Gy with 288
4 M.M. Ulvrova et al. / Earth and Planetary Science Letters 528 (2019) 115836

detected subduction initiations. The length of the simulation does than 1%. When looking at the temperature structure of the sys-
not have any influence on the results as soon as N is large enough tem (Fig. 1b, e, h), one can note that there is a strong subadiabatic
and statistically representative of the system. The parameters of gradient (Fig. 1c, f, i). This is partly because the system is inter-
the calculations are listed in Table 1. We vary the strength of the nally heated but more importantly negative gradients arise due to
oceanic lithosphere and number of continents, and we test how pressure dependent viscosity (Christensen, 1985).
the presence of the weak crustal layer and presence of a free sur- Subduction zones (and in this study we refer to subduction
face alter the results. In total, we run 13 models having different zones as convective downwelling currents in the numerical sim-
parameterizations. ulations) next to continental margins are one-sided while a dis-
tinctive asymmetry is observed for intra-oceanic subduction zones
2.2. Model analysis (cf. Fig. 1 and animation S1 in the supporting material). The de-
gree of asymmetry in the latter case changes through time and
To detect subduction inception and follow its evolution, we an- differs among models. The dip angle with which plates sink into
alyze the divergence of the surface velocity. As soon as a negative the mantle is generally large (close to vertical, Fig. 1a) but more
peak of the surface velocity divergence appears (the threshold of realistic behaviour is observed (i.e. shallower dip angles) in mod-
detection being minus one tenth of the surface root mean square els with the weak crustal layer that allows partial decoupling be-
(rms) velocity), a new subduction zone is formed. The motion of tween the sinking and overriding plates (Fig. 1d), or in models
this peak is then tracked through time and we record the dis- with a free surface (Fig. 1g). However, these models are still far
tance to the closest continental margin until the peak disappears. from producing Earth-like subduction zones in their entirety, with
For each peak we specifically register this distance at subduction modelled sinking plates commonly experiencing phases of sym-
inception and at cessation. These values are further analyzed by metrical subduction during their lifespans, for example at their
calculating a cumulative distribution of the distance to the clos- inception or during polarity flips (Fig. 2). Crameri and Tackley
est continental margin. To construct these distributions we bin the (2014) have observed polarity reversals for intra-oceanic subduc-
range of all possible distances into 500 km wide intervals. We then tion in models with a similar parameterization but in 3D. Our
count the cumulative number of detected distances falling into a models also show polarity reversals at the continent-ocean bound-
specific bin. All subduction zones that initiate at a continental mar- ary (Fig. 2).
gin fall into the first bin. The long duration of each model allows
us to collect up to several hundreds of subduction initiations, typ- 3.1. Influence of lithospheric strength
ically several tens of them.
To investigate systematic patterns within the distributions re- To investigate the impact of the lithospheric strength on sub-
trieved from the models, we use a Monte Carlo method and cal- duction initiation we run several models with different yield stress
culate the synthetic distribution for subduction zones initiating at σ0 , which we vary between 7 and 56 MPa. The plate size hence
random positions within oceans. At each model time we gener- number of subduction zones is strongly influenced by how difficult
ate a set of randomly located points within the oceans, which it is for the lithosphere to localize deformation, which is directly
together have a spatial distribution equivalent to scenarios where related to the value of σ0 . Decreasing the surface yield stress re-
subduction initiation or cessation locations is randomly distributed. sults in weaker oceanic lithosphere and smaller plates (Moresi and
For models with one continent, the cumulative frequency of ran- Solomatov, 1998; Tackley, 2000; Langemeyer et al., 2018). For low
dom point locations as a function of distance to the continent is a yield stress (σ0 = 7 MPa), a large population of subduction zones
straight line since the distance between the two continental mar- exists at a given instant, fluctuating around 7 and 8 (cf. histogram
gins is constant through time and the probability that a subduction on Fig. 3a). The distribution of subduction initiation within the
zone initiates at a certain location is uniform. For more than one oceans can be characterized as a random process distribution as it
continent, the distribution is typically curved. This allows us to follows the distribution of subduction zones that are randomly ini-
show how the modelled subduction distributions deviate from the tiated in within oceans (Fig. 3a). There are small deviations of the
random one. However, it is necessary that a sufficiently large num- distribution from the mean synthetic random distribution, which
ber of random subduction initiations is generated at each time fall into the variance of the random distribution. This is due to fi-
level, in total at least several thousand. To estimate the variance of nite number of subduction zones collected for this model (N = 89).
the random distribution, we calculate distributions of N subduc- The fact that these deviations are small indicates that the number
tion initiation events that happen at random times and random of subduction zones collected is sufficient. The number of subduc-
positions. N is the number of initiations detected in the particular tion zones that are formed in the vicinity of a continental margin
model. (closer than 500 km) is small, around 5% (Fig. 3a). A stiffer litho-
sphere (σ0 = 35 MPa, Fig. 3b) promotes subduction initiation prox-
3. Results imal to continents, with around 30% of subduction zones being
formed close to continental margins. Within this case, the influ-
The organization of the system dictates its dynamic evolution: ence of the continents on subduction initiation is strongest close
sinking slabs drive plate motion while inducing the mantle flow to the continents, which is in accordance with previous studies
that in turn is at the origin of the plate motion (Fig. 1). Plate-like that showed stress focusing at the continental margins (Rolf and
behaviour is developed self-consistently using the yielding rheol- Tackley, 2011). Beyond a certain threshold distance from the conti-
ogy (cf. Section 2.2) with a surface velocity that is constant in plate nents the probability of subduction initiation is essentially random.
interiors while changing abruptly over plate boundaries. Since the For intermediate yield stress (σ0 = 35 MPa) this distance is around
system is heated mainly by internal heat sources (and we keep the 3000 km (Fig. 3b). For very stiff lithosphere (σ0 = 56 MPa), the
internal heating rate constant for all simulations), sinking slabs and distance within which continents influence the stress distribution
surface plates that compose the upper boundary layer dominate becomes greater still, with no clear threshold between random and
over plumes created at the core-mantle boundary. In particular, controlled behaviour (Fig. 3c). The frequency of subduction initia-
around 20% of the total surface heat flux is due to heating the tion in the immediate proximity of continental margins is further
mantle from the core. The core derived fraction of the total surface increased and is close to 40% for strong oceans with few (around
heat flux is very similar for all models, with differences smaller two) subduction zones on average (Fig. 3c and Fig. 4a).
 M.M. Ulvrova et al. / Earth and Planetary Science Letters 528 (2019) 115836 5

Fig. 1. A snapshot of the viscosity (left column) and temperature (middle column) fields at one particular time. (Right column) Azimuthally and temporally averaged temper-
ature profiles of minimum (blue), mean (black) and maximum (red) temperature. Model without (a)-(c) and with (d)-(f) weak crustal layer. (g)-(i) Model with free surface.
Air layer atop of the mantle is shown in light blue. Continents are emphasized in all panels in blue. (For interpretation of the colours in the figure(s), the reader is referred
to the web version of this article.)

Subduction zones that are formed next to a continental mar- subduction zones being able to migrate freely, but once they mi-
gin stay glued to the continent their whole existence (Fig. 5). In grate to a continental margin subduction zones stay there until
contrast, subduction zones that initiate as intra-oceanic reside in they cease (Fig. 3b and c).
the oceans where they either merge with another subduction zone
or migrate in the oceans before reaching a continental margin 3.2. Influence of number of continental margins
where they endure for some time before ceasing (Fig. 5). In some
cases, a subduction zone retreats toward the continent and once Repeated continental assembly and dispersal are observed on
it hits the margin, subduction continues with a reversed polarity Earth in the cycles that last for several 100 Myr (e.g. Rogers and
(Fig. 2). Rarely, a subduction zone is formed and ceases as intra- Santosh, 2004). The length of the continental margins thus varies
oceanic without colliding with another convergence zone (Fig. 5). according to the continental configuration, being minimal when
In the weak lithosphere case, termination is random just like initi- continents are aggregated and maximal when dispersed. To inves-
ation (Fig. 3a). Where the continents have an influence, subduction tigate the influence of the number of continental margins on the
zones are more likely to terminate adjacent to continents than to position of subduction initiation within our annulus models, we
initiate there, which is presumably a consequence of intra-oceanic perform a set of calculations with one, two and three continental
6 M.M. Ulvrova et al. / Earth and Planetary Science Letters 528 (2019) 115836

Fig. 2. (a)-(c) Subduction zone retreating toward the continent (emphasized in blue) and reaching the continental margin. (d)-(f) Subduction zone changes its polarity and
continues descending into the mantle next to the margin. The model has a weak crustal layer and intermediate yield stress σ0 = 35 MPa.

Fig. 3. Cumulative distribution of subduction initiation (blue) and cessation (black) as a function of distance from the nearest continent for three different yield stresses σ0
(increasing from left to right). Dark grey line represents the distribution of subduction initiation at random position for a large population of cases. Gray area designates
random distributions generated for N subduction initiations with N being the number of initiations detected for a particular model. Number of subduction zones detected
is (from left to right) 89, 93 and 78. Histograms in the bottom right corners show the distribution of the number of subduction zones in the respective models. The models
have a free slip top boundary but no weak crustal layer.

rafts (where each raft has two margins), while keeping the total strength; secondly, we systematically observe that more subduc-
cover of continents constant at 30% of the annulus. The number tion zones cease at a continental margin compared to the initia-
of continental margins in 2D models corresponds to the length of tion position (cf. Fig. 6). Both these relationships are weakest for
the continental margins in 3D. In these models, we keep the yield the case with a single continent (cf. Fig. 6a). The threshold dis-
stress fixed (σ0 = 35 MPa) and include a weak crustal layer at the tance below which the effect of continents is negligible and the
surface. The number of subduction zones fluctuates around four or distribution follows a random distribution of subduction initiation
five (cf. histograms on Fig. 6). With increasing number of conti- increases with increasing number of continental margins and is
nental margins, the number of intra-oceanic subduction initiations around 4000 km, 6000 km and 7000 km for cases with two, four
decreases, regardless of whether a weak crustal layer is present and six continental margins (Fig. 6).
(cf. Fig. 4b). In particular, subduction initiation at continental mar-
gins increases from 30% for the case with one continent to 67% for 3.3. Asymmetric subduction zones
three continents (cf. Fig. 4b and 6). Two fundamental results are
consistent across all model cases: firstly, all cases significantly dif- At the surface of the Earth, during during convergence of two
fer from initiation at random position for the chosen lithospheric oceanic plates, one of them subducts into the mantle while the
 M.M. Ulvrova et al. / Earth and Planetary Science Letters 528 (2019) 115836 7

Fig. 4. (a) Proportion of all subduction zones that initiate (blue) and cease (black) in the vicinity of the continent as a function of the lithospheric strength. One continental
raft is present (i.e., two margins) throughout the simulations. (b) Proportion of all subduction zones that initiate (blue) and cease (black) in the vicinity of the continent as
a function of the number of the continental margins. The yield stress is fixed at σ0 = 35 MPa. Solid lines on panel a) and b) are for models with compositionally uniform
oceanic lithosphere while dashed lines on panel b) are for runs with weak crustal layer.

Fig. 5. Position of the continents (blue) and subduction zones (coloured lines; one colour corresponds to individual subduction zone) through time together with the surface
heat flux (gray scale). The model has a weak crustal layer and intermediate yield stress σ0 = 35 MPa.

Fig. 6. The influence of the number of continental margins (increasing from left to right). (a) Two continental margins (N = 43), (b) four continental margins (N = 79), (c) six
continental margins (N = 288). N is the number of initiations detected for a particular model. Dark grey line represents the distribution of subduction initiation at random
position for a large population of cases. Gray area designates random distributions generated for N subduction initiations. Histograms in the bottom right corners show the
distribution of the number of subduction zones in the respective models. The yield stress is σ0 = 35 MPa and the models feature a weak crustal layer.
8 M.M. Ulvrova et al. / Earth and Planetary Science Letters 528 (2019) 115836

Fig. 7. Cumulative distribution of subduction initiation (blue) and cessation (black) for model with (a) no weak crustal layer (N = 132), (b) weak crustal layer (N = 79), and (c)
free surface (N = 27). N is the number of subduction initiations detected for a particular model. Dark grey line represents the distribution of subduction initiation at random
position for a large population of cases. Gray area designates random distributions generated for N subduction initiations. Histograms in the bottom right corners show the
distribution of the number of subduction zones in the respective models. The yield stress is fixed at σ0 = 35 MPa. There are two continents throughout the simulations.

overriding plate stays at the surface. Although numerical models evidence (van der Meer et al., 2010, 2012). Both reconstructions
still have limited ability to produce such behaviour, strong asym- attempt to reconcile geologic observations from arc remnants to
metry of sinking slabs is observed in our models as is described some extent, but nonetheless differ in many aspects – notably, slab
in the beginning of this section. The simplest model, with a free- remnants interpreted from seismic tomography suggest a much
slip surface and no weak crustal layer features the least realistic larger population of intra-oceanic subduction zones. Qualitatively,
subduction dipping angles and longer periods of vertical descent. subduction zones in the Müller et al. (2016) reconstruction are
In this case, about 45% of subduction inceptions are adjacent to typically closer to the continents compared to those in the in-
the continental margins (Fig. 7a). For more complex models with terpretation of van der Meer et al. (2010, 2012) (Fig. 8). Below
more realistic slab dip angles due to lubrication and partial de- we describe a first-order quantification of the proximity to con-
coupling at the interface between the two colliding plates (models tinents of subducting segments and their initiation and cessation
with the weak crustal layer), it is more likely that subduction initi- in both reconstructions, which offers some degree of comparison
ates at the continental margin; close to 60% of detected subduction with the results of the numerical simulations described in Sec-
zones start in the vicinity of the continent (Fig. 7b). This is be- tion 2.2. The supporting material contains the computer code and
cause the lateral density gradient between continental and oceanic data files used to perform this analysis and a detailed explanation
lithosphere produces additional compressive stresses that would of the steps we followed.
drive continents to spread below the free-slip boundary, if the vis- In computing distribution functions for reconstructions, uncer-
cosity was low enough. Hence compressive stresses focus at the tainty arises when we attempt to measure the distance from indi-
continent ocean boundary (Nikolaeva et al., 2010; Rolf and Tack- vidual segments along a subduction zone to the nearest continent
ley, 2011) and favour subduction initiation. In combination with – such measurements require a clear definition of what does and
a weak crustal layer having lower yield stress and hence localiz- does not constitute a continent within the reconstructions. The dis-
ing deformation more easily, more subduction zones initiate at the tinction is not clear-cut for the Earth, which contains a spectrum
continental margins. For the models with a free, deformable sur- of stretched, submerged continental fragments that have not con-
face that features more realistic slab dips and stress states, the ventionally been considered continents (e.g. Mortimer et al., 2017),
number of initiation events at the continental margins drops to and where the nature of the crust in the overriding plate can
around 30% (Fig. 7c). In this case, the free surface allows the con- change significantly along strike for individual arcs (for example
tinent to hamper continent spreading by generating a topography the present-day Aleutian arc). Using a set of reconstructable poly-
that can accommodate a fraction of the stresses at the continent gons defining major continental blocks (Fig. 8, see also supporting
ocean boundary and the lateral density difference between the information) we first derived the curves in Fig. 9 by generating
continents and the mantle through isostasy. sets of randomly distributed points within the oceans at different
reconstruction times. We then compute the distance of each point
3.4. Reconstructions of subduction, initiation, and cessation to the nearest continent boundary, and plot the cumulative distri-
bution functions of these distances subdivided into 25 Myr time
The timing and location of past subduction can be recon- windows. These results provide a visual baseline to show how the
structed from geological and geophysical constraints, though such distances to the continents of subduction segments compare to the
reconstructions are subject to large uncertainties over the time- pattern expected from a random process, similar to that used in
scales of supercontinent cycles (several hundred Myr). In partic- Figs. 3, 6 and 7 but on a spherical Earth rather than an annulus.
ular, the lengths and locations of plate boundaries within the For the Müller et al. (2016) reconstruction, we extract sub-
oceanic realm far from continents are poorly known, increasingly duction zone geometries at 1 Myr intervals and resample each
so further back in time, since the crust that comprised these re- line geometry to uniformly spaced half-degree segments. We com-
gions is scarcely preserved at the present day. We compare two pute the distance of segment mid-points to the nearest continent
reconstructions of subduction history (Fig. 8) that are derived by boundary, such that the distances that contribute to the cumula-
different methods, and use these as points of comparison with nu- tive distributions may vary significantly along each distinct line
merical model behaviour. The first reconstruction maps the extent feature. These distances are illustrated for each reconstruction time
of subduction zones within a globally self-consistent framework of in animation S2, and form the basis of the cumulative distribution
plate boundary configurations and plate kinematics (Müller et al., functions in Fig. 9b where the results are subdivided into 25 Myr
2016); the second reconstruction uses subducted slab signatures time bins. Subduction initiation and termination is not explicitly
mapped within seismic tomography models as the primary line of encoded into the kinematic reconstruction, and so we used simple
 M.M. Ulvrova et al. / Earth and Planetary Science Letters 528 (2019) 115836 9

Fig. 8. Position of the continents and subduction zones since the Triassic according to two alternative reconstructions (see text Section 3.4), subdivided into 4 distinct time
windows from Pangea times to recent. The detailed time-evolution of these reconstructions is illustrated in animations S2-S3. (a) Subduction zones and continent positions
for the M2016 model between 0 and 50 Ma, plotted at 10 Myr increments; locations of subduction zones are shown in colours corresponding to the colour legend, while the
continents are shown in gray with darker gray standing for younger positions within the 0-50 Myr period. (b) same as (a) for 50-100 Ma; (c) same as (a) for 100-150 Ma;
(d) same as (a) for 150-230 Ma; (e) V2012 model for times between 0-50 Ma; (f) same as (e) for 50-100 Ma; (g) same as (e) for 100-150 Ma; (h) same as (e) for 150-235
Ma.

criteria to detect initiating and ceasing segments. The main cri- are available and we take them as the timings of subduction termi-
terion is that initiating or ceasing segments will not lie within a nation and initiation respectively. The histories of additional slabs
threshold distance of any segment from the network of subduc- mapped in van der Meer et al. (2012) are not defined with the
tion zones from 1 Myr earlier or later. The threshold is required same level of detail. Consequently we do not include pre-Cenozoic
because trenches where subduction is ongoing will migrate by a subduction interpreted from slab remnants beneath the Pacific in
finite amount over 1 Myr (animation S2), and we use a threshold two of our distribution plots (Fig. 9e, f).
distance of 200 km to allow for the magnitude of trench migration With these limitations in mind, we ask the question as to
expected within 1 Myr (Williams et al., 2015). An exception to this whether the patterns of subduction initiation and cessation con-
criterion is segments that lie within this threshold distance, but tained within current reconstructions follow a similar distribution
where the pair of plates that bound the adjacent segments differs to that observed in our numerical models. The distribution of
between successive plate boundary snapshots 1 Myr apart (oc- subduction in relation to continents from the alternative subduc-
curring for example in cases of reconstructed subduction polarity tion histories is illustrated using cumulative distribution functions
reversal, which we consider as termination of existing subduction (Fig. 9), providing some analogy to the annulus model results.
and initiation of a new segment). Locations of subduction initiation For both subduction histories, the overall distribution of ongoing
and cessation defined in this way are highlighted in animation S2, subduction is typically closer to the continents than expected for
and the distances to the nearest continents at these points form randomly distributed points within the oceans (Fig. 9b, d, e). This
the basis of the cumulative distribution functions in Fig. 9c. pattern is particularly pronounced for the Müller et al. (2016) kine-
For subduction history interpreted from seismic tomography, matic reconstruction (Fig. 9b), which lacks many intra-oceanic sys-
we compute the distance to the nearest continent of the line tems interpreted by van der Meer et al. (2010, 2012) from seismic
geometries defined by these previous studies at discrete recon- tomography. When we isolate subduction initiation and cessation
struction times being 7 to 17 Myr apart. We make an important (Fig. 9c, f), a further trend emerges that is apparent in both kine-
distinction between slabs tabulated in van der Meer et al. (2010) matic and tomography-based subduction histories – subduction
(Fig. 9e) and the longer list of slabs considered in van der Meer et cessation is typically closer to the continents than subduction ini-
al. (2012) (Fig. 9d). For the former, the top and bottom slab ages tiation. This trend is only absent in the poorly resolved pre-100 Ma
10 M.M. Ulvrova et al. / Earth and Planetary Science Letters 528 (2019) 115836

Fig. 9. Cumulative distribution functions (CDFs) for distance between continents and points along subduction zones at different stages of their development for reconstructions
from the Triassic to present (see supporting text and animations S2-S3). Each coloured line represents the distribution for a specific time in panels (a), (b), (d) and (e). In
panels (c) and (f), the relatively small number of initiating and ceasing subduction segments are subdivided into broad time ranges encompassing the earlier and later stages
of Pangea breakup. See text for further explanation. (a) CDF for random points falling within reconstructed extent of ocean basins; the grey background shows the envelope
of these distributions based on random points, and is reproduced on the other panels for visual reference; (b) distance to continent for segments of active subduction zones
for the kinematic reconstruction of Müller et al. (2016); (c) distances to continent for initiating and ceasing subduction segments derived from Müller et al. (2016); (d)
distance to continent for remnants of past subduction mapped from seismic tomography (van der Meer et al., 2010, 2012); (e) As (d), but only including ‘primary’ subduction
according to the definition of van der Meer et al. (2010); (f) distance to continents for subduction zones in (e) at the beginning and end of their lifespans (assumed to
approximate initiation and cessation) for the slab remnant reconstruction.

section of the Müller et al. (2016) reconstruction and the results do the convective flow (e.g. Guillou and Jaupart, 1995) and influence
not show an obvious distinction between the periods before and heat loss out of the system as they act as thermal insulators (e.g.
after 100 Ma, broadly corresponding to the periods of initial and Lenardic and Moresi, 1999; Phillips and Coltice, 2010; Rolf et al.,
later stages of dispersal of Pangea. However, since the distributions 2012). Importantly, numerical simulations and laboratory experi-
for tomography-based initiation and termination do not include ments suggest that continents change the lithospheric stress dis-
additional slabs interpreted to have existed within the middle of tribution and facilitate subduction initiation (e.g. Nikolaeva et al.,
the Panthalassic ocean (Fig. 8), the proportion of cases far from 2010; Rolf and Tackley, 2011). However, systematic study of the lo-
continents while Pangea was assembled are likely to be underesti- cations of subduction initiation and their ensuing evolution taking
mated in our plots. into account global tectonic settings has received very little atten-
tion.
4. Discussion Comparison between the distribution of subduction in numeri-
cal simulations and those inferred from reconstructions offers in-
Previous studies into the effect of continents on mantle dy- sight into the most plausible model parameters. Models in which
namics have shown that continents increase the wavelength of the lithospheric strength is low, such that subduction is effectively
 M.M. Ulvrova et al. / Earth and Planetary Science Letters 528 (2019) 115836 11

randomly distributed within the oceans, are inconsistent with the of real-world plate tectonics. Nonetheless, even during superconti-
inferred patterns of subduction on Earth during the last >200 Myr nent assembly we would still expect a greater likelihood of intra-
(Fig. 9). Instead, the distributions from reconstructions are more oceanic subduction initiation proximal to continental margins in
compatible with scenarios in which the lithospheric strength is the geological record (e.g. Maffione et al., 2017).
relatively high, such that the continents generate a ‘strain shad- One must keep in mind that numerical simulations have certain
ow’ promoting subduction initiation closer to the continents. This limitations. Primarily, the study relies on 2D calculations, in which
region adjacent to margins is loaded by the density and topogra- only poloidal flow can exist. We do not consider any 3D charac-
phy contrast between continent and ocean lithosphere. Therefore, teristics of mantle flow such as toroidal motion around lateral slab
the lithosphere more readily yields upon experiencing additional edges that can slow down descending slabs (Li and Ribe, 2012).
stresses produced by convection, whereas the same convective The models operate at a Rayleigh number that is lower than the
stresses alone are less likely to reach the yield criterion further Earth’s and hence result in thicker boundary layers, less vigorous
from continents. convection and surface heat flow and plate velocities lower than
In models where sites of initiation of subduction are biased what is expected for our planet. Although basal heating in the
towards regions close to continents, the distribution of subduc- models (representing around 20% of the total surface heat flux)
tion throughout its duration and cessation is also naturally bi- falls into the estimated range of the Earth’s basal heat flux por-
ased towards these regions. Subduction that initiates at a con- tion (Jaupart et al., 2015), we fix the internal heating rate and we
tinental margin remains there until it ceases. Subduction sys- possibly underestimate the influence of plumes. Subduction can be
tems initiating in the oceans may migrate towards or away from initiated when plumes provide sufficient additional stresses caus-
the nearby continental margins; those that reach the continent ing local weakening of the lithosphere after hitting the subsurface
margin become continental arcs, and remain so until subduction (Gerya et al., 2008). However, this is not our case since the plumes
ceases. Hence, the control of continents on patterns of subduc- in the models are relatively weak and hence only rarely trigger
tion initiation influences the distribution of subduction cessation, the formation of new convergent boundaries. Importantly, the em-
such that subduction termination along continental margins occurs ployed rheology neglects any past deformation of the lithosphere
more frequently than subduction initiation, even in the absence of and reflects only the instantaneous stress distribution. However, on
continent-continent collision. Earth plates have memory of previous yielding and can be dam-
Comparison between models with different numbers of conti- aged or undergo healing (e.g Bercovici and Ricard, 2016).
nents (and therefore margins) offers an insight into the different
distributions of subduction that we may expect during different 5. Conclusions
phases of the supercontinent cycle. Cases with two or three dis-
tinct continents (analogous to periods of continental dispersal such We present an assessment of where subduction initiates and
as the last ∼100 Myr on Earth) are more favourable to subduc- ceases in global convection models with a plate-like surface and
tion initiation close to continental margins – in contrast, when the continental drift. We compare the results of numerical simulations
annulus includes only one continent (analogous to a superconti- with distributions of subduction initiation and cessation retrieved
nent such as Pangea), the area influenced by the continental strain from plate tectonics reconstructions and seismic tomography mod-
shadow effect is reduced and the proportion of subduction initi- els. We show that the location of subduction initiation and ces-
ations (and terminations) occurring within the oceans increases. sation is not randomly distributed within the oceans on Earth.
This result contrasts with the results for the Müller et al. (2016) Subduction zones that are formed at continental margins tend to
kinematic reconstruction, where very little subduction initiates, stay there, while subduction zones formed within the oceans mi-
evolves, or ceases far from continents prior to 100 Ma (Fig. 9b, grate and merge with other intra-oceanic subduction zones, or
c). However, distributions from kinematic reconstructions also con- reach continental margins where subduction usually continues,
trast with those for alternative subduction histories interpreted changing polarity before eventually ceasing. Hence, we systemat-
from seismic tomography – in particular when these take into ac- ically find that more subduction zones cease in the vicinity of
count Triassic-Jurassic, intra-Panthalassa subduction systems of van continental margins compared to subduction initiation. Numeri-
der Meer et al. (2012), illustrated in Fig. 8 and included in the cal models indicate that the critical parameters that influence the
distribution of Fig. 9d but not included in Fig. 9f. The ages and lo- position of subduction initiation are the lithospheric strength and
cations of subduction initiation for these systems is unknown, but the number of continental margins. Stronger lithosphere (which
it is reasonable to infer that some or all of these systems initiated implies larger plates and fewer subduction zones (Tackley, 2000))
far from Pangea. increases the probability of subduction initiation in the vicinity of
In addition to the global studies of van der Meer et al. (2010) continental margins. With our numerical simulations we also pre-
and van der Meer et al. (2012), regional studies have also inter- dict that intra-oceanic subduction initiation is more likely during
preted tomography to reveal previously unrecognised intra-oceanic the times of supercontinent assembly, while continental disper-
subduction zones existing within Panthalassa during the Creta- sal favours incipient subduction close to continents. These results
ceous and earlier (Domeier et al., 2017; Sigloch and Mihalynuk, favour interpretations of intra-oceanic subduction systems within
2013). Furthermore, numerical models of past mantle flow con- the Panthalassa Ocean during the time of Pangea based on seismic
strained by subduction histories similar to that of Müller et al. tomography (van der Meer et al., 2012, 2014), which are missing
(2016) produce present-day temperature fields that show good from earlier plate tectonic reconstructions.
agreement with deep mantle seismic velocity variations imaged
by tomography to first-order (Flament et al., 2017) but yield re- Acknowledgements
gional mismatches around parts of the Pacific which may be re-
solved when additional intra-oceanic subduction zones are incor- The authors thank the editor B. Buffet for handling the manu-
porated (e.g. Braz et al., 2018). The emerging view from observa- script and J. Lowman for constructive comments which enhanced
tions is that intra-oceanic subduction was more prevalent during the quality of this paper. The support for this research has been
Pangea’s existence and early breakup than previously recognised. provided by the European Union’s Horizon 2020 research and inno-
Our model results in which intra-oceanic subduction initiation and vation program under the ERC grant agreement no 617588 and the
evolution is favoured during supercontinental assembly are consis- Marie Skłodowska Curie grant agreement no 753755. Simulations
tent with this view, even though our models lack many features were performed on the AUGURY super-computer at P2CHPD Lyon.
12 M.M. Ulvrova et al. / Earth and Planetary Science Letters 528 (2019) 115836

The StagPy library was used in this study to process StagYY out- Leng, W., Gurnis, M., 2015. Subduction initiation at relic arcs. Geophys. Res. Lett. 42,
put data (https://github.com/StagPython/StagPy). The 7014–7021.
pygplates library was used to analyse plate tectonic reconstructions Li, Z.H., Ribe, N.M., 2012. Dynamics of free subduction from 3-d boundary el-
ement modeling. J. Geophys. Res., Solid Earth 117. https://doi.org/10.1029/
(https://www.gplates.org/docs/pygplates/).
2012JB009165.
Lowman, J.P., 2011. Mantle convection models featuring plate tectonic behavior: an
Appendix A. Supplementary material overview of methods and progress. Tectonophysics 510, 1–16.
Maffione, M., van Hinsbergen, D.J., de Gelder, G.I., van der Goes, F.C., Morris, A., 2017.
Supplementary material related to this article can be found on- Kinematics of Late Cretaceous subduction initiation in the neo-tethys ocean re-
line at https://doi.org/10.1016/j.epsl.2019.115836. constructed from ophiolites of Turkey, Cyprus, and Syria. J. Geophys. Res., Solid
Earth 122, 3953–3976.
Matsumoto, T., Tomoda, Y., 1983. Numerical simulation of the initiation of subduc-
References tion at the fracture zone. J. Phys. Earth 31, 183–194. https://doi.org/10.4294/
jpe1952.31.183.
Auzende, J.M., Lafoy, Y., Marsset, B., 1988. Recent geodynamic evolution of the north Michaud, F., Royer, J.Y., Bourgois, J., Dyment, J., Calmus, T., Bandy, W., Sosson, M.,
Fiji basin (southwest Pacific). Geology 16, 925–929. Mortera-Gutiérrez, C., Sichler, B., Rebolledo-Viera, M., et al., 2006. Oceanic-ridge
Baes, M., Sobolev, S.V., Quinteros, J., 2018. Subduction initiation in mid-ocean in-
subduction vs. slab break off: plate tectonic evolution along the Baja California
duced by mantle suction flow. Geophys. J. Int. 215, 1515–1522.
Sur continental margin since 15 Ma. Geology 34, 13–16.
Bercovici, D., Ricard, Y., 2016. Grain-damage hysteresis and plate tectonic states.
Moresi, L., Solomatov, V., 1998. Mantle convection with a brittle lithosphere:
Phys. Earth Planet. Inter. 253, 31–47. https://doi.org/10.1016/j.pepi.2016.01.005.
thoughts on the global tectonic styles of the Earth and Venus. Geophys. J.
Braz, C., Seton, M., Flament, N., Müller, R.D., 2018. Geodynamic reconstruction of
Int. 133, 669–682.
an accreted Cretaceous back-arc basin in the northern Andes. J. Geodyn. 121,
Mortimer, N., Campbell, H.J., Tulloch, A.J., King, P.R., Stagpoole, V.M., Wood, R.A.,
115–132.
Christensen, U.R., 1985. Heat transport by variable viscosity convection ii: pres- Rattenbury, M.S., Sutherland, R., Adams, C.J., Collot, J., et al., 2017. Zealandia:
sure influence, non-Newtonian rheology and decaying heat sources. Phys. Earth Earth’s hidden continent. GSA Today 27, 27–35.
Planet. Inter. 37, 183–205. https://doi.org/10.1016/0031-9201(85)90051-2. Müller, R., Landgrebe, T., 2012. The link between great earthquakes and the subduc-
Coltice, N., Gérault, M., Ulvrová, M., 2017. A mantle convection perspective on global tion of oceanic fracture zones. Solid Earth 3, 447–465.
tectonics. Earth-Sci. Rev. 165, 120–150. Müller, R.D., Seton, M., Zahirovic, S., Williams, S.E., Matthews, K.J., Wright, N.M.,
Coltice, N., Rolf, T., Tackley, P.J., Labrosse, S., 2012. Dynamic causes of the relation Shephard, G.E., Maloney, K.T., Barnett-Moore, N., Hosseinpour, M., Bower, D.J.,
between area and age of the ocean floor. Science 336, 335–338. https://doi.org/ Cannon, J., 2016. Ocean basin evolution and global-scale plate reorganization
10.1126/science.1219120. events since Pangea breakup. Annu. Rev. Earth Planet. Sci. 44, 107–138.
Crameri, F., Schmeling, H., Golabek, G.J., Duretz, T., Orendt, R., Buiter, S.J.H., May, Nikolaeva, K., Gerya, T.V., Marques, F.O., 2010. Subduction initiation at passive mar-
D.A., Kaus, B.J.P., Gerya, T.V., Tackley, P.J., 2012. A comparison of numerical sur- gins: numerical modeling. J. Geophys. Res., Solid Earth 115, B03406. https://
face topography calculations in geodynamic modelling: an evaluation of the doi.org/10.1029/2009JB006549.
‘sticky air’ method. Geophys. J. Int. 189, 38–54. Phillips, B.R., Coltice, N., 2010. Temperature beneath continents as a function of con-
Crameri, F., Tackley, P.J., 2014. Spontaneous development of arcuate single-sided tinental cover and convective wavelength. J. Geophys. Res. 115, B04408. https://
subduction in global 3-D mantle convection models with a free surface. J. Geo- doi.org/10.1029/2009JB006600.
phys. Res., Solid Earth 119, 5921–5942. https://doi.org/10.1002/2014JB010939. Reagan, M.K., Ishizuka, O., Stern, R.J., Kelley, K.A., Ohara, Y., Blichert-Toft, J., Bloomer,
Dilek, Y., Furnes, H., 2009. Structure and geochemistry of Tethyan ophiolites and S.H., Cash, J., Fryer, P., Hanan, B.B., et al., 2010. Fore-arc basalts and subduction
their petrogenesis in subduction rollback systems. Lithos 113, 1–20. initiation in the Izu-Bonin-Mariana system. Geochem. Geophys. Geosyst. 11.
Doin, M.P., Fleitout, L., Christensen, U., 1997. Mantle convection and stability of de- Rogers, J., Santosh, M., 2004. Continents and Supercontinents. Oxford University
pleted and undepleted continental lithosphere. J. Geophys. Res., Solid Earth 102, Press.
2771–2787. https://doi.org/10.1029/96JB03271.
Rolf, T., Coltice, N., Tackley, P., 2012. Linking continental drift, plate tectonics and
Domeier, M., Shephard, G.E., Jakob, J., Gaina, C., Doubrovine, P.V., Torsvik, T.H., 2017.
the thermal state of the Earth’s mantle. Earth Planet. Sci. Lett. 351, 134–146.
Intraoceanic subduction spanned the Pacific in the Late Cretaceous–Paleocene.
Rolf, T., Tackley, P., 2011. Focussing of stress by continents in 3D spherical mantle
Sci. Adv. 3, eaao2303.
convection with self-consistent plate tectonics. Geophys. Res. Lett. 38.
Flament, N., Williams, S., Müller, R., Gurnis, M., Bower, D.J., 2017. Origin and evo-
lution of the deep thermochemical structure beneath Eurasia. Nat. Commun. 8, Sigloch, K., Mihalynuk, M.G., 2013. Intra-oceanic subduction shaped the assembly of
14164. Cordilleran North America. Nature 496, 50.
Gerya, T., Stern, R., Baes, M., Sobolev, S., Whattam, S., 2015. Plate tectonics on the Stern, R.J., 2004. Subduction initiation: spontaneous and induced. Earth Planet. Sci.
earth triggered by plume-induced subduction initiation. Nature 527, 221–225. Lett. 226, 275–292.
Gerya, T.V., Connolly, J.A., Yuen, D.A., 2008. Why is terrestrial subduction one-sided? Tackley, P.J., 2000. Self-consistent generation of tectonic plates in time-dependent,
Geology 36, 43–46. https://doi.org/10.1130/G24060A.1. three-dimensional mantle convection simulations, 1: pseudoplastic yielding.
Guillou, L., Jaupart, C., 1995. On the effect of continents on mantle convection. J. Geochem. Geophys. Geosyst. 1.
Geophys. Res., Solid Earth 100, 24217–24238. Tackley, P.J., 2008. Modelling compressible mantle convection with large viscosity
Hegarty, K.A., Weissel, J.K., Hayes, D.E., 1982. Convergence at the Caroline- contrasts in a three-dimensional spherical shell using the Yin-Yang grid. Phys.
Pacific Plate Boundary: Collision and Subduction. American Geophysical Union, Earth Planet. Inter. 171, 7–18.
pp. 326–348. Tackley, P.J., King, S.D., 2003. Testing the tracer ratio method for modeling ac-
Hernlund, J.W., Tackley, P.J., 2008. Modeling mantle convection in the spherical an- tive compositional fields in mantle convection simulations. Geochem. Geophys.
nulus. Phys. Earth Planet. Inter. 171, 48–54. https://doi.org/10.1016/j.pepi.2008. Geosyst. 4, 8302. https://doi.org/10.1029/2001GC000214.
07.037. Ulvrova, M.M., Brune, S., Williams, S., 2019. Breakup without borders: how conti-
Jaupart, C., Labrosse, S., Lucazeau, F., Mareschal, J., 2015. 7.06 – Temperatures, heat nents speed up and slow down during rifting. Geophys. Res. Lett. 46, 1338–1347.
and energy in the mantle of the earth. In: Schubert, G. (Ed.), Treatise on Geo- https://doi.org/10.1029/2018GL080387.
physics, vol. 7, second edition. Elsevier, Oxford, pp. 223–270.
van der Meer, D.G., Spakman, W., Van Hinsbergen, D.J., Amaru, M.L., Torsvik, T.H.,
Kerrick, D.M., 2001. Present and past nonanthropogenic CO2 degassing from the
2010. Towards absolute plate motions constrained by lower-mantle slab rem-
solid Earth. Rev. Geophys. 39, 565–585. https://doi.org/10.1029/2001RG000105.
nants. Nat. Geosci. 3, 36–40.
Langemeyer, S.M., Lowman, J.P., Tackley, P.J., 2018. The sensitivity of core heat flux
van der Meer, D.G., Torsvik, T.H., Spakman, W., Van Hinsbergen, D.J.J., Amaru, M.L.,
to the modeling of plate-like surface motion. Geochem. Geophys. Geosyst. 19,
2012. Intra-Panthalassa ocean subduction zones revealed by fossil arcs and
1282–1308.
mantle structure. Nat. Geosci. 5, 215–219.
Lebrun, J.F., Lamarche, G., Collot, J.Y., 2003. Subduction initiation at a strike-slip
plate boundary: the Cenozoic Pacific-Australian plate boundary, South of New van der Meer, D.G., Zeebe, R.E., van Hinsbergen, D.J., Sluijs, A., Spakman, W., Torsvik,
Zealand. J. Geophys. Res., Solid Earth 108. T.H., 2014. Plate tectonic controls on atmospheric CO2 levels since the Triassic.
Lee, S.M., 2004. Deformation from the convergence of oceanic lithosphere into Yap Proc. Natl. Acad. Sci. USA 111, 4380–4385.
trench and its implications for early-stage subduction. J. Geodyn. 37, 83–102. Williams, S., Flament, N., Müller, R.D., Butterworth, N., 2015. Absolute plate mo-
Lenardic, A., Moresi, L.N., 1999. Some thoughts on the stability of cratonic litho- tions since 130 ma constrained by subduction zone kinematics. Earth Planet.
sphere: effects of buoyancy and viscosity. J. Geophys. Res., Solid Earth 104, 12. Sci. Lett. 418, 66–77.
https://doi.org/10.1029/1999JB900035.