Interaction of subducted slabs with the mantle transition-zone: A regime diagram from 2-D thermo-mechanical models with a mobile trench and an overriding plate

Transition zone slab deformation inﬂuences Earth’s thermal, chemical, and tectonic evolution. However, the mechanisms responsible for the wide range of imaged slab morphologies remain debated. Here we use 2-D thermo-mechanical models with a mobile trench, an overriding plate, a temperature and stress-dependent rheology, and a 10, 30, or 100-fold increase in lower mantle viscosity, to investigate the effect of initial subducting and overriding-plate ages on slab-transition zone interaction. Four subduction styles emerge: (i) a ‘‘vertical folding’’ mode, with a quasi-stationary trench, near-vertical subduction, and buckling/folding at depth (VF); (ii) slabs that induce mild trench retreat, which are ﬂattened/‘‘horizontally deﬂected’’ and stagnate at the upper-lower mantle interface (HD); (iii) inclined slabs, which result from rapid sinking and strong trench retreat (ISR); (iv) a two-stage mode, displaying backward-bent and subsequently inclined slabs, with late trench retreat (BIR). Transitions from regime (i) to (iii) occur with increasing subduct-ing plate age (i.e., buoyancy and strength). Regime (iv) develops for old (strong) subducting and overriding plates. We ﬁnd that the interplay between trench motion and slab deformation at depth dictates the sub-duction style, both being controlled by slab strength, which is consistent with predictions from previous compositional subduction models. However, due to feedbacks between deformation, sinking rate, temperature, and slab strength, the subducting plate buoyancy, overriding plate strength, and upper-lower mantle viscosity jump are also important controls in thermo-mechanical subduction. For intermediate upper-lower mantle viscosity jumps ( 3 30), our regimes reproduce the diverse range of seismically imaged slab morphologies.


Introduction
Much of the cold material at Earth's surface slides back into the underlying mantle at subduction zones, where one tectonic plate (''slab'') subducts beneath a second (''overriding'') plate, providing the main driving force for plate tectonics [e.g., Forsyth and Uyeda, 1975].By imaging seismic structure, it has been inferred that, in some cases, subducted material extends below 700 km into Earth's lower mantle, while in other regions, slabs bend, buckle, or flatten out in the transition zone, in response to phase and viscosity changes in this depth range [e.g., Isacks and Molnar, 1971;van der Hilst et al., 1991;Fukao et al., 1992;Gudmundsson and Sambridge, 1998;Li et al., 2008;Hayes et al., 2012].This variability in slab deformation is accompanied by significant variations in earthquake potential and basin formation or mountain building [Uyeda and Kanamori, 1979].It has also been linked to the record of plate kinematics throughout the geological past [e.g., van der Hilst and Seno, 1993;Goes et al., 2008] and has likely played a key role in Earth's thermal and chemical evolution [Davies, 1999].
The end-member slab morphologies imaged seismically [e.g., Bijwaard et al., 1998;Li et al., 2008;Fukao et al., 2009] range from: (i) an almost constant dip from upper to lower mantle with possibly some broadening in the transition zone (e.g., Central America, Aegean); (ii) near-vertical slabs that appear to thicken upon interaction with the base of the upper mantle (Marianas, Kermadec); (iii) slabs that have a constant small dip in the upper mantle and flatten at the base of the transition zone (e.g., Japan); to (iv) steep slabs through the upper mantle, which flatten in the transition zone (e.g., Izu-Bonin, Tonga).

PUBLICATIONS
slab flattening in the transition zone [e.g., Kincaid and Olson, 1987;van der Hilst and Seno, 1993;Griffiths et al., 1995;Christensen, 1996].However, only recently have dynamic models started to shed light on the parameters that control the velocity and direction of trench motion.The most systematic modeling undertaken, thus far, has been with a compositional/multimaterial, athermal, free-subduction approach, where subduction is driven only by the downgoing plate's negative buoyancy, while it is resisted by the viscous mantle, without any influence of an overriding plate.Such models have illustrated that high slab density, high slab viscosity, and small slab width all increase the trench's ability to move and that trench motion exerts an important control on slab morphology [e.g., Becker et al., 1999;Funiciello et al., 2003;Bellahsen et al., 2005;Enns et al., 2005;Stegman et al., 2006;Capitanio et al., 2007;Schellart et al., 2007;Di Giuseppe et al., 2008;Ribe, 2010;Stegman et al., 2010a].More recent models demonstrate that the overriding plate also affects trench velocity, possibly introducing a time-dependent trench evolution [e.g., Clark et al., 2008;Capitanio et al., 2010;Leng and Gurnis, 2011].
In addition to the important role of trench motion, it has been shown that a slab's ability to penetrate into the lower mantle is affected by: (i) slab strength and density, relative to the viscosity and density contrast between upper and lower mantle; (ii) the angle at which the slab reaches the transition zone; and (iii) phase transitions [e.g., Gurnis and Hager, 1988;Kincaid and Olson, 1987;Guillou-Frottier et al., 1995;Torii and Yoshioka, 2007;B ehounkov a and C ı zkov a, 2008;Alisic et al., 2012;C ı zkov a and Bina, 2013].Slab density, viscosity, thickness, and the Clapeyron slope of phase transitions are all temperature-dependent.Accordingly, to study slab-transition zone interaction self-consistently, models must account for the role of temperature [e.g., Gurnis and Hager, 1988;Zhong and Gurnis, 1995a;Schmeling et al., 1999;Ci zkov aetal., 2002].With the exception of some recent studies [e.g., Tagawa et al., 2007;Magni et al., 2012;Nakakuki and Mura, 2013], most thermo-mechanical models have considered a fixed trench [e.g., B ehounkov a and C ı zkov a, 2008;Billen, 2010;Lee and King, 2011;Rodr ıguez-Gonz alez et al., 2012], thus constraining the subduction force-balance and providing only limited insight into the dominant controls on the diversity of slab morphologies.
The aim of our study is to obtain an improved understanding of slab-transition zone interaction modes, in a selfconsistent thermo-mechanical system, with a freely moving trench and an overriding plate.We use a 2-D model setup that captures most of the key physics, excluding phase transitions and 3-D effects.Our models include both subducting and overriding plates, and the upper-lower mantle transition is represented by a viscosity jump.A weak decoupling layer, at the interface between subducting and overriding plates, and a free-surface facilitate trench motion [e.g., Schmeling et al., 2008;De Franco et al., 2006;Quinquis et al., 2011;Crameri et al., 2012].Plate and mantle rheology depend on temperature and stress.Thermal structure controls slab thickness, density, and viscosity, and, hence, slab strength (resistance to bending and/or stretching deformation dependent on both effective viscosity and thickness) and buoyancy (the integral of density over slab thickness) evolve self-consistently with the underlying thermal state.Our numerical approach, which uses an adaptive-mesh methodology, allows us to consider very large domains, thereby reducing the control of boundary conditions on the resulting dynamics.We perform a systematic thermo-mechanical study, in line with recent compositional/ multimaterial models, and derive a set of diagnostic outputs to quantitatively distinguish different modes of subduction, allowing us to produce a regime diagram as a function of initial subducting and overriding-plate ages.
Models where the upper-lower mantle transition is represented by a viscosity jump, like those examined herein, are a necessary first step to gain an understanding of the system before adding additional complexities that will induce further (nonlinear) feedbacks.We examine cases where subducting and overridingplate ages are varied, spanning the range observed on Earth.With this setup, we assess (i) how thermomechanical subduction, in which density and strength evolve with temperature, differs from compositionally defined subduction; (ii) the role of the overriding plate, which is a natural part of any thermally defined subduction system; and (iii) how subduction dynamics and the resulting slab morphologies are affected by the viscosity increase from upper to lower mantle.Future work will incorporate the complexities associated with phase transitions, in addition to the third dimension.

Model Description
In this section, we describe our model setup.We summarize the governing equations, geometry, and imposed boundary and initial conditions in section 2.1, the choice of temperature, pressure, and strain ratedependent composite rheology in section 2.2, and the numerical solution strategies in section 2.3.

Model Setup and Governing Equations
We solve the standard equations describing the conservation of mass, momentum, and energy, for an incompressible Stokes fluid, under the Boussinesq approximation @ i u i 50; (1) (2) @T @t 1u i @ i T5j@ 2 i T; (3) where u and g denote velocity and gravity vectors, respectively, r ij the stress tensor, T the temperature, j the thermal diffusivity, and Dq52aq s T2T s ðÞ the density difference due to temperature, with a the coefficient of thermal expansion and q s the nominal density at the surface temperature T s (Table 1).
The full stress tensor r ij can be decomposed into deviatoric and lithostatic components where s ij is the deviatoric stress tensor, p is dynamic pressure, and d ij is the Kronecker delta function.The deviatoric stress tensor and the strain rate tensor _ e ij are related via where l represents the viscosity.
Our model setup is presented in Figure 1.We limit the influence of side and bottom boundary conditions [Enns et al., 2005;Chertova et al., 2012] by running our simulations in a very wide computational domain (10,000 km), which spans the full 2900 km mantle depth.Thermal boundary conditions are isothermal at the domain's top (T 5 T s ) and bottom surfaces (T 5 T m ), with insulating sidewalls (zero heat flux).Velocity boundary conditions are free-slip everywhere, except for the free-surface top boundary.
Our initial condition consists of two plates: the subducting plate (SP, also referred to as downgoing plate) and the overriding plate (OP).A half-space cooling model [Turcotte and Schubert, 2002] yields the initial plate temperature via with z the depth, x the horizontal coordinate, and t the time.The initial plate age, Age 0 (x), evolves linearly from 0 Myr at the left (right) corner to Age 0 SP ðAge 0 OP Þ at the trench (at x 0 trench 55000 km), for the subducting (overriding) plate (Figure 1).We investigate a wide range of initial plate ages, from young 20 Myr old plates to older plates with a thermal age of 100 Myr.We explore a wide range of Age 0 SP and Age 0 OP combinations, to capture the diverse subduction scenarios found on Earth [e.g., Lallemand et al., 2005;M€ uller et al., 2008;Seton et al., 2012].
To circumvent the problem of subduction initiation, a slab long enough to allow for self-sustaining subduction is prescribed as an initial condition.The initial slab has a bending radius of 250 km and an initial slap tip depth of 194 km (see Figure 1), which is sufficient to initiate subduction without external forcing.
A 5 km thick low-viscosity decoupling layer is positioned between the two plates and extends along the top of the subducting plate.This layer is similar in thickness to oceanic crust.It ensures decoupling between the two plates and ensures asymmetric subduction.The interface between the two plates (the ''trench'') is free to move in response to the underlying dynamics.

Rheology
Laboratory experiments have shown that olivine (and its polymorphs), the dominant upper-mantle mineral(s), can deform through different processes under upper-mantle conditions.At high temperatures, either diffusion creep (low stress) or dislocation creep (high stress) dominates [Ranalli, 1995].In addition, some form of stress limiter is expected to occur at low temperature (i.e., in the cold core of the slab) [Karato, 2008].At very shallow depths, yielding (either brittle or plastic) limits the shear stress.At greater depths, it is expected that the plastic yield strength depends on temperature, and although we may not know deeper mantle (or even lithospheric) plastic deformation mechanisms well, the Peierls mechanism, a low-temperature dislocation glide, gives a sensible formulation of a temperature-dependent strength limit at low temperatures [e.g., Evans and Goetze, 1979].
We implement a composite viscosity based upon these four deformation mechanisms.In our simulations, the entire computational domain is governed by identical rheological laws: there is no compositional boundary between slab and mantle material.In this way, plate age determines both buoyancy (through a temperature-dependent density and thickness) and strength (through a temperature-dependent viscosity and thickness).This differs from the compositional/multimaterial approach, where buoyancy and strength are only coupled through thickness and can also be varied independently.Note that we will use the term We use a generic prefactor with no explicit dependency on water content and grain size, assumed constant in the model.c In the lower mantle, we set minute prefactors of 10 242 and 10 2300 for dislocation and Peierls creeps, respectively, so that they do not occur below 660 km.d A friction coefficient of 0.2 is intermediate between lower values of previous subduction models [Di Giuseppe et al., 2008;Crameri et al., 2012] and the actual friction coefficient of the Byerlee law [Byerlee, 1978].
e A very high value is taken for cases with Peierls mechanism, for which the yield strength mechanism dominates only at shallow depth.We discuss in section 5.2 lower values of the yield strength used when decoupling of strength and buoyancy.
Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 strength to refer to a plate's resistance to deformation, which is proportional to effective viscosity times thickness in the case of uniaxial stretching or compression and to effective viscosity times thickness cubed in the case of bending [Ribe, 2001].In our models, strength varies spatially and temporally within plates and the surrounding mantle.
Diffusion, dislocation, and Peierls creep viscosity are calculated via a generic relationship between stress and strain rate for each mechanism with A a prefactor, n the stress exponent, E, V the activation energy and volume, respectively, P5q s gz the lithostatic pressure, R the gas constant, and _ e II the second invariant of the strain rate tensor.T r is the temperature obtained by adding to the Boussinesq solution an adiabatic gradient of 0.5 K/km in the upper mantle and 0.3 K/km in the lower mantle [Fowler, 2005].
Yielding at low pressure is implemented through a brittle-failure type yield-stress law with l y the yielding viscosity and s y the yield strength, given by s y 5min s 0 1f c P; s y;max ÀÁ ; with s 0 the surface yield strength, f c the friction coefficient, P the lithostatic pressure, and s y;max the maximum yield strength.
For Peierls creep we follow the simplification in Ci zkov aetal.[2002] by using equation ( 7) with a high exponent n to approximate the exponential strain rate and the temperature dependency of the Peierls mechanism.With the parameters given in Table 1, this implementation yields a similar variation of effective viscosity as a function of stress and temperature as published Peierls laws [e.g., Goetze, 1978;Rubie, 1984;Kameyama et al., 1999;Karato et al., 2001].
The composite viscosity of the material is calculated via Viscosity is subsequently capped with both upper and lower limits (see Table 1).In the lower mantle, the same formulations are used, but we set negligible prefactors for dislocation and Peierls creep such that The material in the whole domain has the same rheology: the distinction between plates and mantle is not imposed, but arises self-consistently from temperature structure.The heat flux q is null at the side boundaries.Free-slip (FS) conditions are imposed on bottom and sides, with a free-surface at the top.Age 0 SP and Age 0 OP are the initial ages of the subducting (SP) and overriding plates (OP) at the trench, respectively.Dl is the jump between diffusion creep upper (UM) and lower mantle (LM) viscosities at 660 km.Dl is 30 in most cases and is varied in section 4.5.The coordinate origin is on the top left of the domain.The initial hook geometry of the subducting plate is prescribed using a bending radius of 250 km (including the weak layer, shown in light gray) and an angle b of 77 .The star indicates the rightmost location of the slab tip at x 5 5246 km and z 5 194 km.
The activation parameters and stress-dependent exponent used (Table 1) are consistent with those derived from experimental data on olivine [e.g., Karato and Wu, 1993;Hirth and Kohlstedt, 1995;Ranalli, 1995;Hirth and Kohlstedt, 2003;Korenaga and Karato, 2008] and are also similar to those used by previous studies [e.g., Schmeling et al., 1999;Ci zkov aetal., 2002;Billen and Hirth, 2007;Buffett and Becker, 2012;Nakakuki and Mura, 2013].The prefactors A for the different deformation mechanisms are tailored based on the arguments hereafter, which allows us to obtain a plausible range of subduction dynamics.
The prefactor A for diffusion creep is set to yield a mantle viscosity profile that is compatible with geoid and postglacial rebound observations [e.g., Mitrovica and Forte, 2004].Differences of A between upper (UM) and lower mantle (LM) for diffusion creep are set to yield a viscosity increase Dl of a factor 10, 30, or 100 (Table 1).Most models are run with Dl 5 30, but we investigate the influence of viscosity jumps of Dl 5 10 and Dl 5 100 in section 4.5.The relative parameters for dislocation and diffusion creep are chosen to allow dislocation creep to occur down to 100-200 km below the plate and around the slab, consistent with observations of seismic anisotropy [e.g., Karato et al., 2008;Debayle and Ricard, 2013].Seismic evidence of slab deformation throughout the upper mantle [e.g., Isacks and Molnar, 1971;Bijwaard et al., 1998] allows us to set the relative strength of dislocation creep and stress-limiting rheology, as dislocation creep would lead to exceedingly high viscosities at the low temperatures in the slab's centre.Our parameters are set such that for effective yield strengths of around 100-500 MPa, yield-stress or Peierls deformation mechanism is activated in the slab core.
In the weak decoupling layer, the viscosity follows the same rheological laws (equation ( 10)).However, a friction coefficient f c that is 10 times lower than the remainder of the domain is prescribed, while viscosities are truncated above 10 20 Pa s to ensure decoupling between downgoing and overriding plates.The prescribed friction coefficient is in agreement with studies that quantify the degree of coupling at plate boundaries [e.g., Iaffaldano, 2012;Arcay, 2012].In contrast to a weak layer with a constant viscosity, the composite rheology used allows strain rate weakening, such that strain naturally localizes.The unique rheological properties of the weak layer are graded out below 200 km depth.

Numerical Solution Strategy
Numerical simulation of dynamic subduction is challenging.The thermo-mechanical approach utilized herein demands (i) extreme local resolution, particularly at interfaces between the slab and background mantle/overriding plate; (ii) solvers that can handle the sharp, order-of-magnitude viscosity contrasts that arise under this framework; (iii) a free-surface, to allow the slab to decouple from the surface.The recently developed computational modeling framework, Fluidity [Davies et al., 2011;Kramer et al., 2012], a finiteelement, control-volume code that is built upon adaptive, unstructured discretizations [e.g., Davies et al., 2007Davies et al., , 2008]], is ideal for such simulations.Fluidity has been extensively validated for geodynamical problems against a range of analytical and benchmark solutions [Davies et al., 2011;Kramer et al., 2012;Le Voci et al., 2014].The adaptive mesh strategies underlying Fluidity are outlined in detail by Davies et al. [2011].These allow us to run multiresolution simulations with element size varying between 400 m at the plate interface to 200 km in less active regions of the domain, such as the lowermost mantle.A series of tests at different resolutions demonstrate that such a mesh spacing yields well-resolved solutions.
The solution strategies employed to solve the discretized form of equations ( 1) and ( 2) are identical to those outlined in Davies et al. [2011], with the addition of a modification to incorporate an implicit treatment of the free-surface [see Kramer et al., 2012, for further details].Equation (3) is discretized using finite volumes, constructed around the vertices of the triangular finite element mesh.The temperature is evaluated at the boundaries of the volumes using the finite element basis functions.A flux limiter [Sweby, 1984] is then applied to avoid spurious oscillations.The location of the weak decoupling layer, at the interface between the subducting and overriding plates, is tracked using a volume fraction, the evolution of which is described by a linear advection equation.In order to avoid excessive numerical diffusion of the weak layer into neighboring regions, this is discretized on the control volume mesh using the minimally diffusive HyperC face-value scheme (see Leonard [1991] and Wilson [2009] for further details).

Simulation Diagnostics: Inputs and Outputs
We use several diagnostic quantities to analyze model output, beyond the visual observations of slab morphology.The quantities chosen are, in part, motivated by previous insights from compositional models of free-subduction [e.g., Bellahsen et al., 2005;Capitanio et al., 2007;Di Giuseppe et al., 2008;Ribe, 2010;Stegman et al., 2010a].These are introduced below: Initial lithosphere thickness.The lithosphere is defined as material colder than 1300 K: following equation ( 6), the initial thermal thickness h 0 SP of the subducting plate at the trench is 45, 55, 63, 80, and 100 km for Age 0 SP of 20, 30, 40, 65, and 100 Myr, respectively.
Sinking times.t 660 is defined as the time when the slab's 1300 K isotherm reaches 660 km depth.Similarly, we define t 800 as the time when the slab reaches 800 km depth.These provide estimates of average sinking rates during transit through the upper mantle and penetration into the lower mantle, respectively.
Diffusion efficiency.The efficiency of thermal diffusion from slab to mantle over t 660 is characterized by the ratio of diffusion length L diff to h 0 SP .L diff equals 2 ffiffiffiffiffiffiffiffiffiffi jt 660 p .
Plate velocities.Horizontal plate velocities V SP and V OP of the subducting and overriding plates are measured at the surface at x 5 2000 and 8000 km, respectively, within the rigid part of the plates (Figure 2b).Note that V OP has a positive value when the overriding plate velocity is moving toward the trench (leftward), as shown in Figure 2b.The convergence velocity V conv is the sum of V SP and V OP .Subducting-plate velocity V SP reaches a peak value during upper-mantle sinking, at a time termed t5t max : V tmax SP and V tmax OP are measured at this time.Stokes sinking velocity.In compositional free-subduction, retreating slabs were found to sink with their Stokes velocity [Guillou-Frottier et al., 1995;Capitanio et al., 2007;Ribe, 2010].Furthermore, subducting-plate velocities tended toward slab sinking velocities for steeply dipping slabs [Capitanio et al., 2007], consistent with observations of Pacific subduction velocities [Faccenna et al., 2007;Goes et al., 2011].We calculate a Stokes velocity, V Stokes , according to the expression given by Capitanio et al. [2009]: ffiffi ffi 2 p l UM Þ, with upper-mantle background viscosity l UM 54310 20 Pa s, slab length L 0 of 1000 km, and excess density jDqj, calculated from the equation of state, of 57 kg/m 3 .The Stokes velocity is then compared to an estimate of the sinking velocity V sink to identify different styles of subduction (Table 2).V sink is calculated at t max as V conv sin ðdÞ, where the dip d is determined between 150 and 400 km depth.
Trench motion.Trench location x trench is tracked as the boundary between the weak layer and the overriding plate at the surface.In particular, we evaluate Dx 02t660 trench and Dx t6602t800 trench , which denote the amount of trench retreat at t 660 and t 800 relative to the initial trench location x 0 trench or to x t660 trench , respectively.These provide estimates of the mean trench velocity during the upper-mantle subduction phase, V 02t660 trench , and during SP 5100 Myr and Age 0 OP 520 or 100 Myr, at times 3.2 and 17.4 Myr, respectively (trench location of 4884 and 4917 km, respectively).Note that a transition in velocities direction occurs in the trench region.The back-arc region (right of the trench) displays different velocities than the bulk overriding plate for the younger overriding plate age.The different quantities calculated to quantify morphology, and subduction dynamics are defined in section 3.
Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 lower-mantle sinking, V t6602t800 trench , respectively.Note that Dx trench is positive when trench retreats (i.e., the overriding plate move toward the left of the domain (see Figure 2a)).
Overriding-plate deformation.Figure 2b presents the variation of the horizontal surface velocity along the x axis, from two cases that are later examined.In one of the cases, velocity is constant in both subducting and overriding plates away from the trench.However, in the other, a higher velocity is observed 100-200 km away from the trench in the overriding plate.The velocity is termed V back-arc , and the difference V tmax back-arc 2V tmax OP roughly quantifies the extensional deformation within the overriding plate.
Slab geometry.As shown in Figure 2a, the coordinates of the rightmost location of the slab in the mantle are used to calculate slab advance Dx max .The increase in the slab's horizontal extent from its initial state, Dx slab , is the sum of trench retreat and slab advance (i.e., Dx slab 5Dx trench 1Dx max ).
Role of trench retreat in slab morphology.Subduction in our thermo-mechanical approach is achieved by a combination of advance of the subducting plate toward the trench and motion of the trench over the downgoing plate (i.e., positive V OP values).The ratio V OP =V conv quantifies the relative contribution of subducting plate advance and trench retreat to the length of slab subducted.Similarly, the ratio Dx trench =Dx slab indicates the importance of trench retreat (rather than forward push) in the total horizontal slab extent.This ratio is calculated over time interval [0-t 660 ] for the slab's upper-mantle evolution.

Results
4.1.Self-Consistent Subduction Dynamics in a Thermo-Mechanical Model: An Example We illustrate the generic dynamics of a simulation using a case where the initial ages of subducting and overriding plates at the trench are 100 and 20 Myr, respectively.Figure 3 shows the simultaneous temporal evolution of temperature, viscosity, the dominant deformation mechanism, and the underlying computational mesh.Figure 4 illustrates the associated subducting plate velocity V SP, trench motion Dx trench , and V OP =V conv ratio, as a function of time.
In this example, the slab sinks rapidly through the upper mantle, with the 1300 K isotherm reaching 660 km depth after 4 Myr.The initial subducting-plate acceleration(phase1)showedinFigure4ais a See Figure 2a and section 3 for a definition of the different parameters.
Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 predominantly driven by increasing slab pull as the slab lengthens.However, it is enhanced by strain rate weakening in the dislocation creep regime that prevails in the mantle around the slab (Figure 3c).Similarly, strain rate weakening naturally forms an asthenosphere below the unsubducted part of the plate.As the slab starts to interact with the high-viscosity lower mantle (phase 2), V SP reduces to 2cm/ yr (around 10% of the maximum velocity attained in the upper mantle see Figure 4a).The slab then deforms and flattens along the upper-lower mantle interface and subsequently penetrates into the lower mantle (Figure 3b).
Convergence is the result of both motion of the subducting plate into the trench (V SP in Figure 4a) and motion of the trench toward the subducting plate (Dx trench in Figure 4b).During phase 1, subduction occurs predominantly through downgoing-plate advance, as illustrated by the minima of around 20% in V OP =V conv (Figure 4c).In phase 2, V SP decreases and trench retreat becomes the dominant contribution to convergence (V OP =V conv increases to >50% between 3 and 8 Myr).As the slab lengthens in the upper mantle, the increased pull leads to a decrease of V OP =V conv and, subsequently, penetration into the lower mantle.The relative contribution of trench retreat to convergence continues oscillating between phases of slab flattening and lower-mantle penetration, on a 10 Myr period from 45%-75%.
With our thermo-mechanical approach, plate viscosities and thicknesses evolve with temperature.Additional time dependence in plate strength occurs due to yielding (at shallow depths) or through Peierls creep (Figure 3c).This leads to low viscosity in high strain rate regions: (i) in the subducting plate bending region before the trench; (ii) in the slab unbending region at depth; and (iii) in the deformed part of the slab adjacent to the mantle viscosity jump (Figure 3b).Other thermo-mechanical subduction models that utilize a composite rheology have found similarly complex and time-dependent viscosity patterns [e.g., Billen and  3b and 3c and red lines in Figure 3d mark the location of the 1300 K isotherm.The images display only a zoom-in of the total simulation domain (see Figure 1).Note that we have saturated the lower bound of the viscosity (log) scale at 10 20 Pa s (the viscosity can be as low as 10 18 Pa s).There is no compositional difference between the slab and the mantle, only a density difference due to temperature.Note the weakening (low viscosity) of the subducting plate in the bending region left of the trench.The computational mesh is adapted at fixed intervals during the simulation to provide high resolution in regions of high solution curvatures.Metric advection [Wilson, 2009] is used to ensure that sufficient resolution is available between those intervals.
Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 Hirth, 2007;C ı zkov aetal., 2007].Overriding-plate extension during upper-mantle sinking V tmax back-arc > V tmax OP ÀÁ is enhanced by strain rate weakening, facilitating trench retreat.Nonetheless, the plates behave rigidly away from these deformation zones, and a strong slab core at maximum viscosity is maintained throughout the upper mantle, acting as a down-dip stress guide.Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257

Influence of Subducting-Plate Age
Here we illustrate how the dynamics of subduction change as a function of subducting-plate age, using models with an initially young (20 Myr) overriding plate (Table 2 and Figures 4 and 5a-5c).We identify three different subduction modes, termed ''inclined-strong retreat'' (ISR), ''horizontally deflected'' (HD), and ''vertical folding'' (VF).The first style (ISR, Figure 5a, Age 0 SP 5100 Myr ), described in the previous section, is characterized by rapid sinking and trench retreat, leading to an inclined slab that partly flattens upon interaction with the upper-lower mantle interface (phase 2).In the second mode (HD, Figure 5b, Age 0 SP 530 Myr ), the slab sinks slowly, at first, and near-vertically through the upper mantle (phase 1).It is then deflected at the upperlower mantle interface, where it flattens (phase 2) and stagnates for a short period.In the third style (VF, Figure 5c, Age 0 SP 520 Myr ), the slab subducts vertically through the upper mantle (phase 1) before buckling and folding upon interaction with the viscosity jump (phase 2) and sinking into the lower mantle.In this case, the trench is quasi-stationary throughout.
We find that the older downgoing plates sink, subduct, and retreat faster (Table 2 and Figure 4).This is as expected from previous free-subduction models in which both velocities increase with increasing negative buoyancy.In addition, those models showed that high slab resistance to bending together with high density (which both increase with Age 0 SP ) leads to stronger retreat [Bellahsen et al., 2005;Schellart, 2008;Capitanio et al., 2007;Ribe, 2010;Stegman et al., 2010a].Sinking rates as a function of age vary more strongly than would be expected purely based on variations in slab buoyancy, as can be seen from the increase in V tmax sink = Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 V Stokes with increasing Age 0 SP (Table 2).The highest ratio, for Age 0 SP 5100 Myr , is 1.9 and 3 times higher than the lowest value observed, for Age 0 SP 520 Myr .This discrepancy is beyond the effect of variable slab length in the upper mantle (around a factor of 1=sin ð60 Þ1:2 between dipping and vertical slabs).We attribute it to larger mantle lubrication around old slabs compared to younger slabs Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 (compare, e.g., Figure 5a at 2.4 Myr and Figure 5c at 9.5 Myr).The lowest ratio is below 1.0, which may indicate some other factors play a role, e.g., thermal equilibration of the slab and or loss of energy in bending (discussed in sections 4.3 and 5.3, respectively).However, it may simply reflect the uncertainty in calculating a reference V Stokes .
The three styles exhibit distinct temporal evolutions of the contribution of trench retreat to total subduction rates (Figure 4c).The older, ISR, cases (100 and 65 Myr) have high retreat rates, but also low V OP =V conv ratios in phase 1 and oscillating ratios in phase 2 that generally exceed 40-50%.The ratios for intermediate age slabs (Age 0 SP 530 and 40 Myr), in HD mode, display a less distinct phase 1 and lower-amplitude phase-2 oscillations with maxima generally below 40%, as for the first 30-40 Myr the slab remains confined to the upper mantle.For the young VF slab (20 Myr), trench retreat, initially already small, rapidly diminishes, and subduction subsequently occurs only through downgoing-plate advance into the trench.
Slab morphology in the transition zone is the result of trench motion in phase 1 and 2 plus deep slab deformation in phase 2. For the strongly retreating ISR slabs, most of the horizontal extent of the slab in the upper mantle is due to trench retreat, evidenced by a value close to 1 for the ratio Dx trench =Dx slab at t 660 (Table 2).For younger HD slabs, low values of Dx trench =Dx slab indicate that less than half of the increased horizontal slab extent is due to trench retreat: the remainder arises due to advance of the slab tip during deflection along the upper-lower mantle interface.This behavior is intermediate between that of freesubduction models without an overriding plate, where flattening was accommodated by trench retreat [e.g., Kincaid and Olson, 1987;Griffiths et al., 1995;Capitanio et al., 2007], and models with a fixed trench [Billen, 2010], in which slow slab sinking in the lower mantle resulted in flattening exclusively by forward propagation of the subducted slab tip.
Deep slab deformation in our thermo-mechanical models differs from that in purely compositional subduction simulations, due to the thermal weakening of the slab that occurs during sinking.The ratio L diff =h 0 SP shows that slab heating is a significant effect, with values that approach or exceed unity for the sinking velocities of the younger slabs.As a result, all our slabs deform significantly within the transition zone.

Influence of Overriding-Plate Age
By increasing the initial age of the overriding plate, we evaluate how the overriding plate influences the style of subduction (Figures 5d, 5e, and 6 and bottom half of Table 2).Older overriding plates hamper trench retreat (Dx trench in Figure 6b and V 02t660 trench in Table 2).Restricted trench retreat forces plates subducting below older overriding plates to bend more strongly at the trench [Capitanio et al., 2007].This leads to steeper and, for Age 0 SP 5100 Myr and the oldest (100 Myr) overriding plate, rolled-over slab shapes during upper-mantle sinking (Figure 5d 5e), and a new style develops for old slabs, once Age 0 SP > 40 Myr (Figure 5d).This latest style is termed ''bent-inclined retreat'' (BIR) and is characterized by low sinking and retreat rate in phase 1 while the slab bends strongly, encountering the viscosity jump with a rolled-over shape.Interaction with the viscosity jump promotes unbending through late trench retreat (phase 2), eventually resulting in an inclined slab.
The role of the overriding plate is particularly important during the first phase of subduction.Due to feedbacks between deformation and strength in the upper plate, the decrease in phase-1 trench motion with Age 0 OP can be almost as large as the decrease in trench motion with Age 0 SP (Table 2).Stress transmitted across the plate contact and by slab-induced flow in the mantle wedge couple overriding plate deformation or translation to slab pull at depth [e.g., De Franco et al., 2006;Arcay et al., 2007;Capitanio et al., 2010;Stegman et al., 2010b;Husson, 2012].Older and, hence, stronger overriding plates offer more resistance to the stretching promoted by slab sinking and hamper trench retreat [Yamato et al., 2009;Capitanio et al., 2010;Butterworth et al., 2012].In several of our cases, additional overriding-plate weakening occurs through yielding, thus enhancing the variation of trench motion with Age Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 (between 1.2 and 1.5 cm/yr; Table 2), indicating that long-term trench retreat is dominated by lower-mantle slab sinking rather than strength of the overriding plate.Smaller trench retreat rates lead to lower convergence rates, which triggers several thermo-mechanical feedbacks (Table 2): (i) less mantle weakening around the slab and thus lower sinking rates (higher t 660 ), (ii) larger thermal diffusion L diff =h 0 SP at t 660 ÀÁ , which may further slow sinking.Although diffusion, by itself, should not affect the slab's net buoyancy [e.g., Davies, 1999;Turcotte and Schubert, 2002], less of the negative buoyancy may be effective in driving flow, due to the low viscosity of the slab's edges.For old slabs Age 0 SP 5100 Myr ÀÁ , slab heating does not vary much with Age 0 OP ðL diff =h 0 SP 0:5 in Table 2), but the consequent slab weakening during sinking does facilitate phase-2 unbending of the BIR slab.For young subducting plates (e.g., Age 0 SP 530 Myr and Age 0 OP 565 Myr in Figure 5e), lower sinking rates modify slab strength sufficiently to change the subduction style.

Summary of the Controls: Toward a Regime Diagram
Before analyzing the influence of the upper-lower mantle viscosity jump, we summarize subducting and overriding plate effects.Figure 7 presents some of the diagnostic outputs, introduced in section 3, for all simulations examined with an upper-lower mantle viscosity jump, Dl 5 30.These are combined to discriminate between the different subduction regimes defined in sections 4.2 and 4.3.
The ''vertical folding'' (VF) morphology is easily discriminated by (i) a quasi-stationary trench throughout its evolution (Figure 7e); (ii) a vertical slab morphology (Figures 5c and 5e at t 40 Myr ); (iii) slow uppermantle sinking that results in very high L diff =h 0 SP ratios ( 0.9; Figure 7a), leading to a significant loss of strength during descent; and (iv) low sinking velocities relative to V Stokes ( 0.6; Figure 7b).''Horizontally deflected'' (HD) slabs are identified by steep initial descent followed by flattening at 660 km depth (Figure 5b).A typical characteristic is their large Dx 02t660 max ( 60 km; Figure 7c).They exhibit mild trench retreat between 0.5 and 1.5 cm/yr (Figures 7e and 7f), with a smaller loss of strength and buoyancy through diffusion than VF slabs (L diff =h 0 SP ratios between 0.6 and 1.2 (Figure 7a)).Flattening is, in roughly equal proportions, the result of trench retreat and deflection at the viscosity jump resulting in a forward push, with Dx trench =Dx slab ratios, at t 660 ,of<0.7 (Figure 7d).Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 ''Inclined-strong retreat'' (ISR) and ''bent then inclined and retreat'' (BIR) modes display larger Dx trench =Dx slab at t 660 than the HD regime (Figure 7d): trench retreat is the main factor increasing the slab's horizontal extent in the two styles.The discrimination between ISR and BIR is made through the ratio V tmax sink =V Stokes (Figure 7b), the upper-mantle trench retreat velocity V 02t660 trench (Figure 7e), and from comparisons between V 02t660 trench and V t6602t800 trench (Figure 7f): acceleration (BIR) or slow-down (ISR) of trench retreat following interaction with the lower mantle (Figure 6b).We observe overriding plate extension only for ISR cases (Table 2).
Guided by the dynamic, kinematic, and geometric characteristics of Figure 7, we build a regime diagram (Figure 8), as a function of initial subducting and overriding-plate ages.Note that plate ages evolve during the simulations as a result of the competition between thermal diffusion from the surface and convergence velocity.Transitions between different regimes are gradual rather than sharp and we find some cases that are intermediate between the end-member regimes identified.As Figure 8 illustrates, subducting-plate age clearly exerts the dominant control on subduction style.However, a strong overriding plate can modulate subduction mode.

Effects of the Viscosity Jump
All results presented, thus far, are for cases with an intermediate upper-lower mantle viscosity jump, Dl 5 30. Figure 9 illustrates the influence of lower (Dl 5 10) or higher (Dl 5 100) viscosity contrasts, for cases with initial downgoing and overriding-plate ages of 40 and 20 Myr, respectively, while Figure 10 illustrates the temporal evolution of V SP ; Dx trench , and V OP =V conv , for the same cases.These figures demonstrate that slab morphologies can change substantially by varying Dl.While Dl 5 30 resulted in a horizontally deflected (HD) slab, for Dl 5 10, the evolution resembles that of an ISR case, while Dl 5 100 starts off like a HD case, but eventually forms a folded vertical slab, of similar morphology to the VF mode.The latter mode of deformation is new: an addition to the four end-member styles defined previously.In this regime, the slab is deflected at depth (Figure 9 at 16 Myr), with subduction subsequently evolving toward a stagnant trench regime, as shown in Figures 10b and 10c (regime ''HD-VF'').Similarly, for older slabs, we find styles that initially experience retreat before the trench becomes stationary (''R-VF mode,'' during which the slab can also exhibit some flattening at the bottom of the upper mantle).
One of the main reasons for modified subduction regimes with different viscosity jumps is a change in upper-mantle sinking time (t 660 ).Although initially the evolution of the three cases in Figure 9 is similar, as soon as the slab starts to interact with the higher-viscosity lower mantle, slab sinking slows down to a different degree, depending on Dl: upper-mantle sinking time, t 660 , is 7, 11, and 27 Myr for jumps of 10, 30, and Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 100, respectively.Slower sinking rates in cases with a higher Dl lead to weaker (lower-viscosity) slabs encountering the viscosity jump.Lower slab strength due to increased Dl leads to lower trench retreat rates (e.g., upper/lower mantle rates for Age 0 SP of 40 Myr and Age 0 OP of 20 Myr, V 02t660 trench and V t6602t800 trench , range from 1.3/1.1 cm/yr at Dl 5 10 to 1.0/0.7 cm/yr at Dl 5 30 and to 0.5/0 cm/yr at Dl 5 100).These effects, together with the increased resistance to penetration associated with a more viscous lower mantle, produce Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 increasingly deformed and stalling slabs at depth, as Dl increases (Figure 9).This behavior was also observed by Loiselet et al. [2010] and C ı zkov a and Bina [2013].
Results from cases with different viscosity jumps, for a suite of downgoing and overriding-plate ages, are presented in the supporting information.This information is summarized in Figure 11, where regime diagrams are presented for simulations with Dl 5 10 and Dl 5 100.In general, we find an extension of the strong-retreat cases at the expense of the HD and VF domain when Dl is reduced from 30 to 10.In contrast, for Dl 5 100, slow slab sinking results in weak-plate modes for most slabs, with a cessation of trench retreat when the slab is delayed by the viscosity jump (HD-VF and R-VF modes).This results in deflected and vertical slab morphologies at the expense of ISR modes.The BIR mode, for strong slabs and thick overriding plates, is found for all values of Dl, as is the VF mode for younger slabs and older overriding plates.

Model Assumptions and Limitations
Our models are advanced in terms of (i) the inclusion of both subducting and overriding plates; (ii) plate, trench, and mantle motions being fully dynamic; (iii) the inclusion of a free-surface; and (iv) their ability to self-consistently capture multiscale feedbacks between temperature, strength, and buoyancy.However, they have been simplified in other respects, which are next discussed.
The assumption of two-dimensionality.Our results reiterate the key role of free trench motion in controlling the dynamics of subduction.However, several studies have demonstrated that trench motions are facilitated by mantle flow around 3-D slabs [e.g., Funiciello et al., 2003;Bellahsen et al., 2005;Piromallo et al., 2006;Morra and Regenauer-Lieb, 2006;Stegman et al., 2006;Schellart et al., 2007] and, hence, absolute values of trench motion from our 2-D models must be treated with caution.Nonetheless, the subduction regimes predicted here, in many ways, agree with those from 3-D compositional models (see section 5.3), while the mechanisms controlling slab behavior in response to a viscosity jump predicted herein should remain valid in 3-D.The 2-D models examined here likely behave similarly to the centre of wide 3-D slabs.As a consequence, for narrower slabs we would expect boundaries between the subduction modes in Figure 8 to shift, enlarging the regimes with significant trench retreat to younger Age 0 SP and older Age 0 OP .Boundary conditions.All the results presented are from cases with free-slip sidewalls and the same initial trench location (5000 km from the sidewalls; see Figure 1) and initial plate lengths.However, additional Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 simulations, which are not presented herein, demonstrate that these parameters can modify mantle flow and the intensity of trench retreat, consistent with the predictions of previous studies [e.g., Gurnis and Hager, 1988;Enns et al., 2005;Quinquis et al., 2011;Chertova et al., 2012].For example, in a computational domain of identical dimensions, cases with shorter subducting and overriding plates exhibit increased trench motion.Furthermore, we found that box depth and basal boundary conditions also play an Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 important role: plate velocities are lower in cases with a 2900 km deep, free-slip bottom (simulations presented here) when compared to a simulation with 1000 km deep, free-stress bottom.
Interface properties.Properties of the weak decoupling layer have also been shown to influence subducting dynamics [e.g., De Franco et al., 2006, 2008;Arcay, 2012;C ı zkov a and Bina, 2013].For a thicker or less viscous weak layer, we would expect more trench motion and, accordingly, more HD or ISR modes at the expense of VF and BIR modes in Figure 8, respectively.Subduction initiation.In our models, subduction initiates from an undeformed, stationary system with prescribed slab shape and length.Self-consistent subduction should initiate from an active, convecting mantle.However, this likely requires additional external forcing [e.g., McKenzie, 1977;Davies, 1988;Zhong and Gurnis, 1995b;Moresi et al., 2000], which is beyond the scope of our study.Ribe [2010] has shown that slab shape upon interaction with the transition zone (which is strongly modulated by the choice of initial condition) plays an important role in controlling the passage of material into the lower mantle, thus affecting model evolution.To minimize this artificial constraint, our initial slab is prescribed to be as short as possible.
Absence of phase transitions.We do not model the mineralogical phase transitions encountered by the slab during its descent through the mantle.These will modify slab buoyancy and rheology (e.g., through grain size reduction) and have been shown to affect slab morphology [e.g., Karato et al., 2001;Tagawa et al., 2007;B ehounkov aand C ı zkov a, 2008].For example, slab buoyancy and, hence, sinking velocity is increased at the olivine-wadsleyite transition near 410 km depth.This will decrease thermal reequilibration of the slab.On the other hand, the addition of an endothermic phase transition at 660 km provides an additional resistance to penetration and helps slab flattening and deflection along the upper-lower mantle boundary [e.g., Christensen and Yuen, 1984;Ci zkov aetal., 2002;Torii and Yoshioka, 2007;Fukao et al., 2009;C ı zkov a and Bina, 2013].

Implementation of Yielding Rheology
There are substantial uncertainties in how the rheology of lithospheric and mantle material evolves as a function of pressure, temperature, strain rate, and grain size.The parameterization of slab viscosity used herein, through a temperature-dependent Peierls creep mechanism, yields weaker slabs than in some other models and leads to predominantly retreating subduction styles.Furthermore, many previous studies have assumed a temperature-independent stress-limiting mechanism.To quantify the effect of our viscosity parameterization, we have examined a few simulations with a different yielding rheology: the Peierls mechanism is removed and the temperature-independent maximum yield strength, s y;max , is varied between 100 and 1000 MPa.Note that the yielding rheology (equation ( 8)) only dominates over other creep mechanisms in cold regions (i.e., in plates and slab).
Figure 12 illustrates coeval slab morphologies for five cases with different yielding formulation, and otherwise identical simulation parameters.When the Peierls creep mechanism is removed, we observe, for the case with a high maximum yield strength of 1000 MPa (Figure 12e), that the slab is very strong (broad strong core at maximum viscosity), conserves an ''umbrella handle,'' rolled-over shape in the upper mantle and experiences little trench motion.This is due to both more difficult unbending of the strong slab [e.g., Bellahsen et al., 2005;Ribe, 2010;Billen, 2010] and a stronger overriding plate, which resists trench motion [e.g., Butterworth et al., 2012].Similar morphologies are predicted at s y;max 5500MPa (Figure 12d), with a thinner high-viscosity core (the slab is somewhat weakened by the yield-strength mechanism when encountering the viscosity jump).The case with s y;max 5100MPa (Figure 12b) weakens more in bending and unbending, which allows the trench to retreat more easily, with the slab also weakening substantially in the vicinity of the viscosity jump (strong folding at an obtuse angle).For the case of s y;max 5300MPa (Figure 12c), slab morphology is similar to the slab with an additional Peierls rheology (Figure 12a), although the Peierls slab is slightly weaker (lower viscosity) upon interaction with the lower mantle.
Most mantle deformation mechanisms are thermally activated, so thermal weakening of the slab during sinking, as occurs in our Peierls formulation, is plausible.The maximum yield strength value is uncertain.However, we find that yield strengths above 500 MPa confine stress-limiting mechanisms to the shallow parts of the slab only and yield a very strong slab core, which is dominated by diffusion creep and maximum viscosity.
Ci zkov aetal.[2002] and B ehounkov a and C ı zkov a [2008] also found that yield strengths of Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 1 GPa make slabs essentially undeformable under the mantle stresses expected.Finally, we emphasize that a temperature-dependent stress-limiting mechanism yields various morphologies as a function of plate ages rather than as a consequence of independently chosen yield-stress values.

Comparison With Previous Subduction Models
A comparison of our results with previous studies yields further insight into (i) the mechanisms controlling which subduction regimes arise and (ii) the role of temperature in these mechanisms.
Such compositional models have illustrated how partitioning of energy dissipation between sinking and bending at the trench can control the dynamics of subduction [e.g., Becker et al., 1999;Bellahsen et al., 2005;Ribe, 2010].Conrad and Hager [1999] first highlighted that bending of a strong downgoing plate at the trench can require a substantial part of a slab's potential energy and even stall subduction.Capitanio et al. [2007] found that, in free-subduction, slabs adjust their dip and partitioning of subduction by plate advance and trench motion to try to minimize the energy dissipated in the bending process and allow the slab to sink at Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 its Stokes velocity.In some of the modes (e.g., advance), however, bending consumes so much of the potential energy that slabs sink at lower velocities.In our models, we also found that sinking velocities may be less than Stokes velocity (in VF mode) or may be enhanced by mantle lubrication through dislocation creep.
In a scaling analysis, Ribe [2010] showed that this energy partitioning makes slab stiffness (resistance to bending) the main control on the subduction mode.As both sinking and bending are driven by slab pull (i.e., slab buoyancy), stiffness is the parameter that modulates the relative importance of the two.If there are feedbacks between slab strength and sinking velocity (as occurs in strain rate-dependent rheologies like yielding, or visco-elasticity), then buoyancy also affects the subduction mode [Ribe, 2010;Stegman et al., 2010a;Leng and Gurnis, 2011;Capitanio and Morra, 2012].In our thermo-mechanical approach, the effect of buoyancy is further amplified as sinking time also modulates slab strength.
For low Age 0 OP , our subduction modes change from VF to HD to ISR (Figure 8) with increasing Age 0 SP .This is fully consistent with the compositionally established mode controls given that Age 0 SP increases buoyancy and strength in tandem and, hence, favors trench retreat and high sinking velocities over bending.None of our models exhibit significant trench advance, because it is discouraged by slab weakening during sinking and by the fact that strong slabs also have high density, which drives trench retreat.Even when our slabs retain their strength at the base of the upper mantle (Figure 12e), trench advance is hampered by the overriding plate.Due to these factors, our BIR mode only occurs with the additional forcing from a thick overriding plate, and, unlike compositional models, does not solely depend on downgoing plate properties.C ı zkov a and Bina [2013], in thermo-mechanical models where they considered only old subducting and overriding plates (70-150 Myr), also obtained subduction in the BIR mode.
Previous thermo-mechanical models proposed a range of controls on the subduction style: upper-lower mantle viscosity contrast, slab strength, relative density contrast of the slab in upper and lower mantle, trench retreat rate, dip of the slab through the transition zone [e.g., Guillou-Frottier et al., 1995;Karato et al., 2001;Ci zkov aetal., 2002;Billen and Hirth, 2007;Tagawa et al., 2007;Torii and Yoshioka, 2007;B ehounkov a and C ı zkov a, 2008;Nakakuki and Mura, 2013;Androvic ˇov aetal., 2013].These results can be understood by their influence on energy partitioning between slab deformation and sinking.
We sketch this framework in Figure 13.Our phase-1 trench motions and slab geometry are controlled by slab buoyancy B SP , slab resistance to bending S SP , plus overriding plate resistance to stretching S OP ,asinthecompositional models with an overriding plate [Yamato et al., 2009;Capitanio et al., 2010;Butterworth et al., 2012].In phase 2, when the slab interacts with the upper-lower mantle transition, sinking is still driven by overall slab buoyancy, while slab deformation is now driven by the difference in upper and lower-mantle sinking velocities.Hence, also in phase 2, slab strength S SP (resistance to bending and thickening) is expected to be the dominant parameter determining subduction style and will modulate slab dip, slab folding, and trench motion.S SP is influenced through Peierls temperature-dependent viscosity by sinking rate during phase 1, hence by slab buoyancy and Dl.Additionally, S OP has a continued effect on trench motion in phase 2, although the similar behavior of BIR and ISR slabs in phase 2 illustrate that this is a secondary effect.Note that the models yield two styles of reasonably efficient lower-mantle penetration (high sinking rate; see t 800 2t 660 in Table 2): (i) the VF style in which the thickening sufficiently increases slab buoyancy to facilitate lower-mantle sinking, and (ii) the strong-slab ISR/BSR styles where slabs enter the lower mantle driven by their high density.

Potential Relevance for Earth
A detailed comparison between model predictions and observations is beyond the scope of this study, and caution should be exercised in comparison with Earth because of model limitations (section 5.1).However, a few points are worth noting: Our subduction modes comprise the range of seismically imaged slab morphologies [e.g., Bijwaard et al., 1998;Van der Voo et al., 1999;Li et al., 2008;Fukao et al., 2009].It should be noted, however, that slabs that are Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 deformed and inclined, either flattening in the transition zone or penetrating into the lower mantle, are most common in present-day subduction zones.Only two vertical slabs (Mariana and Kermadec) and one with a bentover shape (Himalaya) are imaged.Only a few imaged slabs are in phase 1 (i.e. , not yet interacting with the transition zone) and that HD, BIR, and ISR slabs display similar inclined morphologies over longer evolution times.
Models with Dl 5 100 do not reproduce the full range of imaged morphologies, tending to produce vertically folded slabs with stationary trenches.For Dl 5 10, the majority of slabs penetrate through the transition zone without significant deformation.Tomographic subduction surveys [e.g., Bijwaard et al., 1998;Li et al., 2008;Fukao et al., 2009] are thus most compatible with our intermediate Dl 5 30 models.
We find that the initial age of the subducting plate, through its influence on the evolution of slab strength and buoyancy, is the major control on subduction modes (Figure 8).However, as subduction systems evolve on Earth, so does the age of the subducting plate and,h e n c e ,p r e s e n t -d a y ages do not correlate with slab morphology [e.g., King, 2001;Lallemand et al., 2005;Billen, 2010].Pacific subduction exhibits morphological contrasts between eastern (ISR shapes such as under Central America) and western subduction zones (HD-like slabs in Izu-Bonin, VF-like slabs in Mariana).These features do not correlate with present-day ages at the trench (e.g., 80-100 Myr for Tonga-Kermadec; 10-30 Myr in Northern Andean subduction), but may correlate better with early Cenozoic ages of young (35-50 Myr) subduction zones in the west and older (50-60 Myr) subduction zones in the east [Sdrolias and M€ uller, 2006;Seton et al., 2012].
In our models, the HD mode only occurs across a narrow parameter space, while HD slabs only fatten over a few 100 km.By contrast, stagnating morphologies are relatively common on Earth with slabs flattening over much larger dista n c e s( e .g .,I z u -B o n i n ,T o n g a ,C a l a b r i a [ Bijwaard et al., 1998;Li et al., 2008;Fukao et al., 2009].The addition of an endothermic phase transition may increase the occurrence of the HD regime and lead to further flattened morphologies for the ISR mode [e.g., C ı zkov a and Bina, 2013;Zhong and Gurnis, 1995a].

Conclusion
We have presented a new 2-D thermo-mechanical subduction model that includes an overriding plate, a mobile trench, and evolving temperature, pressure, and strain rate-dependent viscosities.
Together with the renewal of lithospheric material through thermal diffusion at the surface, this setup ensures the self-consistency of dynamics and deformation arising from feedbacks between rheology, flow, temperature, and density.The ability of the trench to move in response to the underlying flow field is an essential ingredient to understand the evolution of subduction motions and morphology, in response to interaction with the boundary between upper and lower mantle, here represented by a viscosity jump.
We developed a set of quantitative diagnostics to distinguish four modes of subduction, characterized by different plate-motion histories and resulting slab morphologies: (i) vertical folding (VF); (ii) horizontally deflected (HD); (iii) inclined-strong retreat (ISR); and (iv) bent then inclined and retreat (BIR).Slab deformation and dynamics are controlled by the initial ages of subducting and overriding plates at the trench: the high strength and buoyancy of old subducting plates encourage trench retreat and slab sinking (ISR and BIR modes), while thick overriding plates decrease trench motion and slab velocities (BIR versus ISR, VF versus HD).Slabs that have retained sufficient strength upon reaching the base of the upper mantle will continue to retreat, or will induce retreat, resulting in flattened (young, weak slabs and large viscosity jump-HD regime) or inclined morphologies (old slabs and mild viscosity jump-ISR and BIR).Weaker slabs, which are unable to induce significant trench retreat throughout their evolution, fold upon interaction with the higher viscosity lower mantle (VF regime).
We have built a new regime diagram for thermo-mechanical models with interdependent strength and buoyancy.Although results agree well with those from previous compositional models, evolution of slab strength, and buoyancy during upper-mantle sinking have important effects on interaction with, and penetration into, the lower mantle.Variations in the initial ages of subducting and overriding plates allow us to reproduce the wide range of observed slab morphologies in the transition zone, through their control on Geochemistry, Geophysics, Geosystems 10.1002/2014GC005257 (i) trench motion and (ii) the slab's ability to resist deformation at depth.On Earth, these processes are likely further modified through a number of additional controlling parameters.

Figure 1 .
Figure1.Setup and initial geometry of the subduction simulations.The material in the whole domain has the same rheology: the distinction between plates and mantle is not imposed, but arises self-consistently from temperature structure.The heat flux q is null at the side boundaries.Free-slip (FS) conditions are imposed on bottom and sides, with a free-surface at the top.Age 0 SP and Age 0 OP are the initial ages of the subducting (SP) and overriding plates (OP) at the trench, respectively.Dl is the jump between diffusion creep upper (UM) and lower mantle (LM) viscosities at 660 km.Dl is 30 in most cases and is varied in section 4.5.The coordinate origin is on the top left of the domain.The initial hook geometry of the subducting plate is prescribed using a bending radius of 250 km (including the weak layer, shown in light gray) and an angle b of 77 .The star indicates the rightmost location of the slab tip at x 5 5246 km and z 5 194 km.

Figure 2 .
Figure 2. (a) Sketch of a slab (contour of 1300 K isotherm) when hitting the boundary between upper (UM) and lower mantle (LM).The diagnostic outputs are defined in section 3. The trench retreat Dx t660 trench is equal to x 0 trench 2x t660 trench .Slab advance Dx 02t660 max is defined as the difference between the rightmost x coordinate of the slab and the initial horizontal location of the slab tip.The slab horizontal extent Dx slab is Dx trench 1Dx max .(b) Horizontal velocity along the x axis at the surface of the domain for cases with Age 0SP 5100 Myr and Age 0 OP 520 or 100 Myr, at times 3.2 and 17.4 Myr, respectively (trench location of 4884 and 4917 km, respectively).Note that a transition in velocities direction occurs in the trench region.The back-arc region (right of the trench) displays different velocities than the bulk overriding plate for the younger overriding plate age.The different quantities calculated to quantify morphology, and subduction dynamics are defined in section 3.

Figure 3 .
Figure 3. Snapshots of (a) temperature, (b) viscosity, (c) dominant deformation mechanism, and (d) underlying computational mesh at various times during a simulation with Age 0 SP and Age 0 OP of 100 and 20 Myr, respectively (initial ages).Black squares indicate initial trench location.White lines in Figure 3a indicate isotherms from 600 to 1400 K with a 200 K interval.Black lines in Figures3b and 3cand red lines in Figure3dmark the location of the 1300 K isotherm.The images display only a zoom-in of the total simulation domain (see Figure1).Note that we have saturated the lower bound of the viscosity (log) scale at 10 20 Pa s (the viscosity can be as low as 10 18 Pa s).There is no compositional difference between the slab and the mantle, only a density difference due to temperature.Note the weakening (low viscosity) of the subducting plate in the bending region left of the trench.The computational mesh is adapted at fixed intervals during the simulation to provide high resolution in regions of high solution curvatures.Metric advection[Wilson, 2009] is used to ensure that sufficient resolution is available between those intervals.

Figure 4 .
Figure 4. Influence of Age 0 SP on the surface dynamics: evolution of (a) subducting plate velocity V SP , (b) trench retreat Dx trench , and (c) V OP =V conv as a function of time for simulations with different initial subducting plate ages, and Age 0 OP 520 Myr .The dotted gray lines indicate t 660 for the different simulations.See Figure 2a and section 3 for a definition of the different parameters, and Figure 5 and section 4.2 for a description of the subduction modes ISR (inclined-strong retreat), HD (horizontally deflected), and VF (vertical folding).

Figure 6 .
Figure 6.Influence of Age 0 OP on the surface dynamics: evolution of (a) subducting plate velocity V SP , (b) trench retreat Dx trench , and (c) V OP =V conv as a function of time for simulations with different initial overriding plate ages, and Age 0 SP 5100 Myr .The dotted gray lines indicate t 660 for the different simulations.See Figure 2a and section 3 for a definition of the different parameters, and Figure 5 and section 4.2 for a description of the subduction modes ISR (inclined-strong retreat) and BIR (bent then inclined and retreat).

Figure 8 .
Figure 8. (left) Regime diagram for subduction dynamics as a function of initial subducting and overriding plate ages, with the modes underlined in Figure 5 for a viscosity jump Dl of 30 between upper and lower mantle.More discussion on the discrimination between the different modes is made in section 4.4.Gray symbols indicate a slab morphology intermediate between two subduction modes.Increasing Dl yields more HD cases, as explained in section 4.5.(right) Slab morphologies at t 660 (black contour) and at t 800 (gray fill).Black squares indicate initial trench location, length scales correspond to 200 km, and the dashed line indicates the 660 km viscosity jump.

Figure 9 .
Figure 9. Influence of the viscosity jump Dl on slab morphology: snapshots of viscosity fields during simulations with a 10, 30, or 100-fold viscosity increase between upper and lower mantle.Black squares indicate trench initial location.ISR 5 inclined-strong retreat, HD 5 horizontally deflected, HD-VF 5 horizontally deflected then vertical folding.

Figure 10 .
Figure 10.Influence of the viscosity jump Dl on subduction dynamics: evolution of (a) subducting plate velocity V SP , (b) trench retreat Dx trench , and (c) V OP /V conv as a function of time for simulations with the initial subducting and overriding plate ages of 40 Myr and 20 Myr, respectively, and Dl 5 10, 300, or 100.The dotted lines indicate t 660 for the different simulations.See Figure 2a and section 3 for a definition of the different parameters.

Figure 11 .
Figure 11.Regime diagram of subduction styles for (a) Dl 5 10 and (b) Dl 5 100, to be compared with Figure 8 for Dl 5 30.VF 5 vertical folding, HD 5 horizontally deflected, ISR 5 inclined-strong retreat, BIR 5 bent then inclined and retreat.New modes include HD-VF and R-VF (retreat then vertical folding) (see section 4.5 for more details).

Figure 12 .
Figure 12.Snapshots of (left) viscosity field and (right) dominant deformation mechanism, 16 Myr after the beginning of simulations with Age 0 SP and Age 0 OP of 65 and 20 Myr, respectively.The panels display simulation outputs with either both (a) a Peierls and a yield strength deformation mechanism or only (b-e) a yield strength mechanism with a maximum yield strength s y;max of 100, 300, 500, or 1000 MPa, respectively.Black squares indicate initial trench location.The black contour marks the 1300 K isotherm.

Figure 13 .
Figure13.Schematic summary of the controls on slab morphology during phase 1 and phase 2. Slab kinematics and morphology are controlled by slab strength S SP , slab buoyancy B SP , overriding plate resistance to stretching S OP , and the viscosity increase Dl between upper and lower mantle.Other control parameters may additionally influence final slab morphology (section 5.3).

Table 1 .
Physical Parameters Used in the Simulations a b

Table 2 .
Quantitative Outputs of Subduction Simulations Used to Discriminate Between Subduction Regimes a at 24 Myr).
Once the 100 Myr slab has penetrated into the lower mantle, the mean trench retreat velocity between t 660 and t 800 is almost the same for the different Age 0 OP cases Dx trench =Dx slab at t 660 ,(e)V 02t660 trench ,and(f)V 02t660 trench =V t6602t800 trench .Initial overriding plate ages and subduction regimes are labeled by symbols and colors, respectively.VF cases (with quasi-stationary trenches) have been removed for Figures7d and 7f.
Figure 7. Evolution of diagnostic outputs as a function of initial subducting plate age Age 0 SP for all simulations with Dl530: (a) L diff =h 0 SP at t 660 ,(b)V tmax sink =V Stokes ,(c)Dx 02t660 max ,(d)