Earth and Planetary Science Letters Subducting-slab transition-zone interaction: Stagnation, penetration and mode switches

Seismic tomography shows that subducting slabs can either sink straight into the lower mantle, or lie down in the mantle transition zone. Moreover, some slabs seem to have changed mode from stagnation to penetration or vice versa. We investigate the dynamic controls on these modes and particularly the transition between them using 2D self-consistent thermo-mechanical subduction models. Our models conﬁrm that the ability of the trench to move is key for slab ﬂattening in the transition zone. Over a wide range of plausible Clapeyron slopes and viscosity jumps at the base of the transition zone, hot young slabs (25 Myr in our models) are most likely to penetrate, while cold old slabs (150 Myr) drive more trench motion and tend to stagnate. Several mechanisms are able to induce penetrating slabs to stagnate: ageing of the subducting plate, decreasing upper plate forcing, and increasing Clapeyron slope (e.g. due to the arrival of a more hydrated slab). Getting stagnating slabs to penetrate is more diﬃcult. It can be accomplished by an instantaneous change in the forcing of the upper plate from free to motionless, or a sudden decrease in the Clapeyron slope. A rapid change in plate age at the trench from old to young cannot easily induce penetration. On Earth, ageing of the subducting plate (with accompanying upper plate rifting) may be the most common mechanism for causing slab stagnation, while strong changes in upper plate forcing appear required for triggering slab penetration. scenar-ios, because in a global system, plates are hardly ever fully free (except possibly in the case of an actively opening back arc) or fully ﬁxed. The models show that such strong changes in upper plate forcing are very eﬃcient in inducing a mode change, within few Myr, both from penetrating to stagnating and from stagnating to penetrating.


Introduction
Subducted slabs in today's mantle display a wide variation in morphology, where some slabs seem to stagnate in the transition zone, e.g., beneath the Izu-Bonin region, South-Kurile and Japan, whereas others seem to penetrate into the lower mantle, for example beneath Peru, the Marianas and Central America (van der Hilst et al., 1997;Fukao and Obayashi, 2013). Furthermore, it seems that some slabs that are stagnant in the transition zone today, switched from penetrating to stagnating, for example the Tonga slab . Comparisons of plate reconstructions with tomography indicate that most slabs that subducted in the last ∼200 Myr are now in the lower mantle (e.g., Ricard et al., 1993;van der Meer et al., 2010). This suggests that to avoid longterm accumulation of slabs in the transition zone, slabs that previously stagnated, as many slabs do today, must have changed from stagnation to penetration. Indeed, this has been suggested to have occurred in subduction of the Tethys ocean, where lower-mantle tomography has been interpreted as evidence of slab fragments deposited under significant trench motion (Hafkenscheid et al., 2006).
The mechanisms controlling the behaviour of the slab in the transition zone have been extensively studied. Most models investigated these mechanisms in isolation, and prescribed part of the system, e.g. trench position, trench motion, slab dip. In this way, it has been shown that slab stagnation is often accompanied by trench retreat (van der Hilst and Seno, 1993) and that high trench migration rates promote slab flattening in the transition zone (Christensen, 1996;Čížková et al., 2002;Griffiths et al., 1995), and low dip angles hamper penetration (Tagawa et al., 2007;Torii and Yoshioka, 2007). However, some models indicate that high trench migration rates are only effective at stagnating slabs if they sufficiently weaken in the transition zone (Čížková et al., 2002). It is clear also that a strong Clapeyron slope for the postspinel transition as well as increasing viscosity at the upperlower mantle boundary may inhibit slab penetration (Billen, 2010;Torii and Yoshioka, 2007). Similarly, the persistence inside the slab of relatively low-density metastable olivine and pyroxene can hamper penetration (Tetzlaff and Schmeling, 2000;Agrusta et al., 2014). In a dynamic system, trench motion and slab dip will evolve as slabs of different density and strength interact with the increasing resistance to sinking at the base of the upper mantle. Although intuitively one might expect old slabs to sink more readily into the lower mantle because of their higher negative buoyancy and larger resistance to bending, free subduction models find that young (weak and relatively buoyant) plates drive less trench retreat than older (denser and stronger) slabs (Capitanio et al., 2007;Garel et al., 2014;Goes et al., 2008;Ribe, 2010;Stegman et al., 2010). Karato et al. (2001) suggested that both young and very old, cold slabs are weak and stagnate, the latter weakened by a metastable olivine region leading to small grain sizes, while intermediate-aged slabs are expected to be strong and hence penetrate. Č ížková et al. (2002), however, showed that stagnation by weakening due to grain-size reduction is only effective if subduction is accompanied by fast trench retreat (>4 cmyr −1 ). Moreover, Nakakuki et al. (2010) showed that slab weakening can lead to penetration in the long-term. Weak slabs have also been suggested to accumulate in the transition zone and cause so-called flushing events, which represent the sinking of the accumulated slab material into the lower mantle (Machetel and Weber, 1991;Tackley et al., 1993;Zhong and Gurnis, 1995). So a range of behaviour may be possible as a function of slab strength and age.
Forcing by the upper plate has also been proposed as a possible trigger or control on subduction-transition-zone interaction. For example, several studies Čížková et al., 2007;Clark et al., 2008;Tagawa et al., 2007) illustrate the importance of upper plate forcing (free vs. fixed) on trench motion. Holt et al. (2015) show that a thicker and more buoyant upper plate leads to decreased trench motions, and Iaffaldano et al. (2006) proposed that building of the Andean mountain belt put a break on Nazca subduction motions.  proposed that the start of back-arc opening, i.e., upper-plate weakening, may have initiated a mode change for Tonga subduction from penetration to stagnation. Clark et al. (2008) demonstrated that slab-transition zone interaction can induce episodic back-arc opening.
In this paper we aim to address the following questions, using a dynamically self-consistent numerical subduction model: (i) Which parameters most efficiently allow slabs to penetrate into or to stagnate above the lower mantle? (ii) What are the plausible mechanisms to change the slab morphology from one mode to the other? and (iii) Which mechanisms are most likely on Earth?

Model set-up
Numerical simulations of the subducting slab are performed using the finite-element code CITCOM (Moresi and Gurnis, 1996;Wang et al., 2015). The system of equations, conservation of mass, momentum, and energy -for an incompressible fluid, infinite Prandtl number, under the extended Boussinesq approximation (Christensen and Yuen, 1985), without internal heating -describing the flow are: where u, P , T , T S , t, z, σ , ε and μ are the dimensionless velocity, pressure, temperature, surface temperature, time, depth, stress, strain rate, and viscosity, respectively. The dimensionless parameters are Ra, the thermal Rayleigh number, Rb i , the boundary Rayleigh associated with the ith phase transition, and Di, the dissipation number. They are defined as follows: with α 0 the surface thermal expansion coefficient, T the temperature drop over the box of height H , g the gravity acceleration, c p the specific heat, κ the thermal diffusivity, μ 0 the reference viscosity, ρ 0 the surface reference density, and ρ i the density anomaly related to the ith mantle phase transition. The coefficient of thermal expansion is considered to be depth dependent, ranging from 3 × 10 −5 K −1 at the surface to 1 × 10 −5 K −1 at the core-mantle boundary (e.g. Tosi et al., 2013) as: The relative fraction of the heavier phase is described by the phase function Γ , varying from 0 and 1 as a function of pressure and temperature, as: where d i is the width of the transformation in depth, γ i is the Clapeyron slope, and z i and T i are the depth and temperature of the ith mantle phase transition at equilibrium conditions, respectively. Used values are listed in Table S1 in the supplementary material. The rheological model is assumed to be a combination of linear diffusion creep and a pseudo-brittle Byerlee rheology. The effective viscosity μ eff is calculated from the viscosities of the individual mechanisms as: . The factor μ lower/upper defines the viscosity jump at 670 km depth, and is 1 in the upper mantle. P in this case corresponds to hydrostatic pressure. The meaning and values of the other parameters are listed in Table S1 (supplementary material). A lower and an upper viscosity cutoff of 10 18 Pa s and 10 24 Pa s, respectively, are imposed for numerical stability.
The model domain is 6000 km wide (with x ranging from −3000 to 3000 km) and 3000 km high, and the box is discretized into 1392 × 390 elements. The element sizes range from 3 km to 9 km. The grid is refined vertically between 0 and 240 km depth, and horizontally between x = −2100 km and x = 1050 km.
The mechanical boundary conditions are free slip along all boundaries, so mantle and slabs are driven purely internally by buoyancy forces. The thermal boundary conditions are: constant temperature, of 273 K, at the surface, a constant temperature at the bottom, corresponding to that at the base of an adiabat with a potential surface temperature of 1573 K, and a mid-ocean ridge temperature profile on the vertical boundaries, to keep mid-ocean ridges at the two model corners. The subducting plate extends from the ridge located in the upper-left corner to the trench, initially situated at x = 0 km, and its initial temperature follows the half-space cooling model. To avoid complications associated with subduction initiation, at the trench the plate follows a path with radius of curvature 500 km into the mantle to a depth of 200 km, which allows self-sustained subduction from the start. The overriding plate extends from a ridge at the upper-right corner to the trench, with a plate age at the boundary with the subducting plate of 100 Myr. On the top of the subducting plate, an 8 km wide low-viscosity crustal layer (μ crust = 10 20 Pa s) is present down to 200 km depth to facilitate the decoupling of the two plates (Fig. 1a).

Model cases
Our first aim is to investigate how slabs with different thermal buoyancy and strength behave during their sinking and how they interact with the transition zone and lower mantle. Because buoyancy and strength are associated with slab age, where young slabs are more buoyant and less viscous than old ones, we performed a total of 32 simulations (Table 1) considering two subducting plate ages at the trench, 25 Myr and 150 Myr.
All transition-zone and lower-mantle parameters investigated are listed in Table 1. In the second part of the paper, we will con-sider cases where one of the parameters, e.g., Clapeyron slope or subducting plate age, changes during the simulation.

Role of slab age
Previous models of free subduction have shown that older slabs (dense and strong) tend to retreat more than younger slabs (Capitanio et al., 2007;Capitanio, 2013;Garel et al., 2014;Ribe, 2010). In light of this, we initially focus on the difference of the sinking behaviour between slabs of different ages. The evolution of all models, up to the time that the slab first reaches the base of the upper mantle (hereafter defined as τ 670 ), is comparable in terms of their general behaviour, regardless of transition-zone and lower-mantle parameters. Indeed, all models begin with the slab tip sinking almost vertically and followed by slow trench retreat for old plates and almost no trench motion for young slabs. The time it takes the slabs to reach 670 km depth, τ 670 , is ∼10 Myr for the plate of 150 Myr old and ∼20 Myr for the 25 Myr old plate (Table 1).
After this time, the slabs behave differently, and we will initially discuss the cases where young and old slabs show the greatest difference in their interaction with the transition zone and lower mantle. Fig. 2 shows the evolution of Simulation 7, which has a transition zone characterised by Clapeyron slopes of 4 MPa K −1 for the ol-wd transition and −2 MPaK −1 for the post-spinel transition, and a viscosity jump μ lower/upper = 10. When the slab reaches the base of the upper mantle, it encounters resistance from the post-spinel transition and the increase in viscosity. The young slab initially folds inside the transition zone and afterwards continues its sinking into the lower mantle ( Fig. 2a), whereas the old one initially penetrates into the lower mantle up to ∼1000 km, and afterwards bends and lies down above the lower mantle, stagnating in the transition zone (Fig. 2b).
Fig. 2c-f show some parameters that can be related to the interaction of the sinking slab with the transition zone and lower mantle. The subducting plate velocity (Fig. 2c) shows, for both models, a maximum that is due to the additional negative buoy-   ancy from the upward deflection of the ol-wd transformation, and a smaller peak linked with the slab folding and/or bending upon reaching the base of the upper mantle (Čížková and Bina, 2013). Afterwards, the plates move at 1-2 cm yr −1 until the end of the simulations. Both slabs are accompanied by motion of the overriding plate and the related trench motion, which is faster for the old plate, 2-3 cm yr −1 , and relatively slow for the young slab, <1 cm yr −1 (Fig. 2d, e). The old slab attains a smaller subduction dip angle of ∼50 • (estimated between 100 km and 400 km), compared to the ∼60 • of the young slab. These observations agree with the previous studies that showed that old slabs drive faster trench retreat with shallower slab dip angle than younger slabs (e.g. Capitanio et al., 2007;Capitanio, 2013).
The ability of the slab to sink into the lower mantle or accumulate at the base of the upper mantle can be visualised by showing the evolution through time of the maximum depth of the slab, and by the amount of the slab material that is in the wadsleyiteringwoodite stability field, i.e. the slab material in the transition zone (Fig. 2g, h). Indeed, after ∼20 Myr (≈2τ 670 ), the old slab does not continue its sinking below ∼1000 km and the amount of slab material in the transition zone increases roughly linearly with time. By contrast, the young slab continues to increase its maximum depth and the amount of slab material trapped in the transition zone is almost constant during the simulation time. The relatively low sinking velocity of the young slab leads to an ageing of the subducting plate during the model evolution. This ageing induces the slab to sink under more trench retreat and with a shal- lower dip angle (after ≈3.5τ 670 in Fig. 2d, f), but despite this, the young slab continues to penetrate (Fig. 2g).

Modes as a function of transition-zone properties
Next, we analyse whether the slab behaviour observed in section 3.1 is true for the range of transition-zone resistance properties. We vary the viscosity jump at the base of the transition zone and both Clapeyron slopes, where the ol-wd and post-spinel Clapeyron slopes are increased or decreased together, considering that the slope depends on water content, with stronger values for a hydrated compositions (2 wt%) and lower slopes for a dry mineralogy (Litasov and Ohtani, 2007).
To discriminate the interaction mode for each simulation, we evaluate the change of the maximum slab depth and the amount of slab material in transition zone through time (Fig. 2g, h). A slab that stagnates in the transition zone must show a relatively constant value of the maximum slab depth and an increase of slab material in the transition zone (Fig. 2g, h black line), whereas a slab sinking into the lower mantle must show an increase of the maximum slab depth and a relatively constant amount of slab material in the transition zone (Fig. 2g, h, red line). We choose the time of 2τ 670 (twice the time it takes the slab to reach 670 km depth) as the initial time of the quasisteady state, from which the mode of the slabs does not change significantly through time. We use the value of the maximum slab depth, Slab max-depth (2τ 670 ), and the area of the slab material, Slab T Z (2τ 670 ), at this time, to normalise the corresponding averaged values over the quasi-steady state time window, which last until 7τ 670 , which corresponds to ∼50 Myr for the simulation with the old slab and ∼100 Myr in the case with the young plate. For a theoretical slab case where thermal diffusion is negligible and the slab is not deformable, Slab max-depth /Slab max-depth (2τ 670 ) = 1 and Slab T Z /Slab T Z (2τ 670 ) > 1 would represent stagnation into the transition zone, whereas Slab max-depth /Slab max-depth (2τ 670 ) > 1 and Slab T Z /Slab T Z (2τ 670 ) = 1 would indicate penetration. Fig. 3a shows the two ratios for all simulations performed. The simulations lie in two clusters that represent the two modes of transition-zone interaction. Young slabs (open symbols) and old slabs (filled symbols) are found in both clusters, indicating that young slabs do not always penetrate and that old slabs do not always stagnate. At intermediate Clapeyron slopes (green and red symbols), both modes are present, i.e. penetrating for young slabs and stagnating for old, as in Simulation 7 above. Instead, at high Clapeyron slopes (blue symbols), slabs of all ages are prone to stagnation, irrespective of the viscosity jump, whereas at low Clapeyron slopes (purple symbols), slabs of all ages penetrate. Increasing the viscosity of the lower mantle does not hamper slab penetration for these low Clapeyron slopes, but slab morphology varies with young slabs sinking more vertically, driving less trench retreat than old slabs, that flatten even as they keep sinking into the lower mantle.
A regime diagram summarising how the relative behaviour of young and old slabs changes with phase-transition strength and viscosity jump at the base of the upper mantle is shown in Fig. 3b. There is quite a wide range of plausible conditions over which penetrating young slabs and stagnating olds slabs are expected to co-exist.

Role of trench retreat and sinking velocity
It has long been clear that trench motion exerts a control on the mode of slab-transition zone interaction (Christensen, 1996;Cížková et al., 2002;Griffiths et al., 1995;van der Hilst and Seno, 1993). However, if we analyse our models, there is not an immediately apparent correlation between absolute trench motion and lower mantle slab sinking rates (Fig. 4a), e.g. there is no simple limiting upper-plate (trench) velocity above which slabs stagnate. Strong trench retreat occurs both in cases that stagnate and in cases that penetrate. Nonetheless, it is clear that for stagnating cases, upper plate motion is large relative to lower mantle sinking velocity.
We investigate this further by performing a scaling exercise. One might expect that velocities would scale with mantle viscosity, as Stokes sinking velocities do, and that trench velocities too would scale with sinking velocity (Capitanio et al., 2007;Ribe, 2010). In Fig. 4b, we therefore have scaled (i) trench (upper-plate) velocity with lower-mantle sinking velocity, (ii) lower-mantle sinking velocity with an estimate of Stokes sinking velocity. The first scaling (ratio of trench retreat over sinking velocity) directly links to the average dip angle of the transition-zone slab, which strongly links to slab stagnation. The second scaling should account for the effect of the changing lower-mantle viscosity between models. Because Stokes velocity cannot easily be estimated for the deforming slabs in the models, we take the cases with the lowest viscosity jump and lowest Clapeyron slopes (i.e. Simulation 1 in Table 1) as our reference, and assume they are sinking with something close to their Stokes velocity into the lower mantle. . Stokes velocities depend on slab density, shape and mantle viscosity (Turcotte and Schubert, 1982). We obtain approximate Stokes velocities for the other cases by correcting the reference velocity for the viscos- . Some of the scatter in the resulting scaled sinking velocities (including values >1) is due to differences in slab geometry between the scaled and reference cases. Nonetheless, scaling the velocities in this way separates out nicely the stagnating and penetrating cases (Fig. 4b).
Our models show that the slabs that accumulate in the transition zone, are actually not fully stalled, but rather sink slowly into the lower mantle compared to the rate at which slab material is delivered to the transition zone. These slabs have a flattened morphology, and this is how they have been identified in seismic tomography. They are stagnant in the sense that they sink inefficiently into the lower mantle and hence accumulate in the transition zone. For slab stagnation, there appears to be a critical lower mantle sinking velocity of about 0.5 times the sinking velocity expected if the slabs could sink freely (as a Stokes sinker) into the lower mantle. The main factor hampering free sinking is the negative Clapeyron slope (see Fig. S1 in supplementary material).
If there is a threshold retreat velocity above which slabs stagnate, then it is higher for old (about 3.5 times V LM sink ) than for young slabs (about 1.5 V LM sink ). This difference in threshold retreat velocity may be due to differences in slab strength between old and young slabs, making it easier to flatten weaker (young) than stronger (old) slabs, i.e., young slabs require less relative retreat (and hence have a higher threshold dip below which slabs stagnate) than old slabs.

Mechanisms to induce mode changes
Few previous studies have looked at how subduction mode can switch from a slab initially trapped into the transition zone to sinking into the lower mantle, or vice versa. In earlier studies, slabs accumulated in large piles above the endothermic phase transition leading to flushing events (Machetel and Weber, 1991;Tackley et al., 1993;Zhong and Gurnis, 1995). However, such slabs had relatively low strength, and more recent subduction models in which slabs behave more plate-like do not easily generate such large-scale piles (Garel et al., 2014;Stegman et al., 2010). A recent study showed that penetration could be triggered if the viscosity of a stagnant slab would reduce by heating over time (Nakakuki et al., 2010). Here we investigate several other mechanisms that may allow a slab to switch its interaction with the transition zone and lower mantle from one mode to another.

Change in Clapeyron slope
Clapeyron slope has been shown to be one of the dominant controls on whether slabs stagnate in the transition zone. If slabs carry sufficient water, and comprehensive slab dehydration occurs while it is in the transition zone (Iwamori and Nakakuki, 2013), then the strength of the regional Clapeyron slope could reduce significantly over time.
Motivated by the regime diagram obtained in the Section 3.2, we test if such reduction of the strength of the Clapeyron slope can induce a change from stagnation to penetration. For this test, we choose the Simulation 4 as starting point, with a subducting plate of 150 Myr old and initial Clapeyron slopes of 5 MPa K −1 for the ol-wd transition and −3 MPa K −1 for the post-spinel transition (i.e., for ∼2 wt% water according to Litasov and Ohtani, 2007).
This case leads to slab stagnation (Figs. 3 and 5a). During the calculation, we change the values of the Clapeyron slopes to conform to those of Simulation 1, with 2.5 MPa K −1 for the ol-wd transition and −0.5 MPa K −1 for the post-spinel transition (dry composition, Litasov and Ohtani, 2007), where the slab sinks into the lower mantle (Figs. 3 and 5b). For simplicity, we assume that slab dehydration takes place instantaneously at ∼48 Myr from the beginning of the simulation (∼6τ 670 ).
This change induces a reduction of the driving force due to the ol-wd transition and of resisting force due to the post-spinel transition. The reduction of resisting forces of the post-spinel transition is enough to induce the sinking of an initially stagnant slab (Fig. 5c). Within approximately 5 Myr after the change, the slab starts to sink into the lower mantle as observed by the slab tip which begins to increase its depth and the decrease in the amount of material trapped in the mantle transition zone (black line in Fig. 5g, h). The mode switch does not seem to strongly change the trench retreat, but reduces the slab dip angle (black line in Fig. 5e, f). It takes about 2τ 670 , i.e. about 20 Myr, for the slab dip and the rate of material accumulation in the transition zone to become similar to that of the continuously penetrating Simulation 1.
An opposite change in Clapeyron slope (from low dry values to higher wet values) is not a very likely scenario, but it may induce a change from penetration to stagnation. We perform the same exercise, suddenly changing the Clapeyron slopes from the Simulation 1 to the Simulation 4 at ∼43 Myr (∼5τ 670 ). The slab begins to lie down into the transition zone (Fig. 5d), although the leading part of the slab that is already in the lower mantle continues its sinking (green line in Fig. 5g, h). Also in this case, there is no strong change in the roll back, but the slab dip increases (green line in Fig. 5e, f). The response in slab dip and accumulation rate in the transition zone occurs more rapidly, within a few Myr, than for the transition from stagnant to penetrating.

Change in subducting-plate age
Since the penetration mode depends on plate age (i.e. the subducting-plate buoyancy and strength) at the trench, and this is a parameter that varies considerably through time at most trenches (e.g. Sdrolias and Müller, 2006), we next test the possibility to switch from one mode to another by changing the age of the subducting plate. We build on Simulation 7, which results in different interaction modes as a function of plate age, with a stagnant old slab and a penetrating young slab (Figs. 2 and 3b). We impose for this case an initial subducting plate structure with a change in age at x = −1000 km. Fig. 6a shows the evolution of the subducting plate that begins to subduct with a plate of 25 Myr old. While the subducting plate is young, the slab penetrates into the lower mantle and the evolution of plate velocities, trench retreat and dip angle (Fig. 6b-e) largely follows the same trend as the 25 Myr old subducting plate shown in Fig. 2. When the old part of the subducting plate reaches the trench and begins to subduct at ∼70 Myr (∼3τ 670 ), an increase of subducting plate velocity, trench retreat, and a reduction of dip angle is observed (Fig. 6b-e). As the older plate continues its sinking, and as soon as the oldest part reaches the base of the upper mantle (∼10 Myr) the slab lies down and stagnates in the transition zone (Fig. 6g).
The opposite change, where old-plate subduction is followed by subduction of a young plate is shown in Fig. 7. Like for the previous case, the evolution follows the same trend as the 150 Myr old plate shown in Fig. 2 (Fig. 7b-g), until the first part of the young slab starts sinking. At ∼25 Myr (∼2.5τ 670 ), when the young part of the plate begins to subduct, a reduction of trench retreat rate accompanied by an increase of the slab dip angle is observed and the accumulation of the amount of slab in the transition zone almost stalls (Fig. 7c-e). However, as the young slab subducts, plate motion slows down and the age of the subducting plate at the trench increases. This leads to a resumption of trench retreat and slab accumulation in the transition zone ( Fig. 7d-g). Hence, in this case, the change to a more buoyant and weaker sinking plate does not change the mode fast enough to overcome the effect of the previous history of subduction, and the slab remains trapped in the transition zone.
We further test whether a stagnant slab can be made to penetrate by making the subducting plate younger, using another variation on the simulation in Fig. 7, a case where the mid-ocean ridge is free to move. This allows for a more protracted time period of young-plate subduction (over ∼100 Myr). However, once again the initially stagnant slab does not start to penetrate. Only, as soon as the mid-ocean ridge enters into the trench, the slab breaks off and the remaining trapped slab, slowly starts to sink in the lower mantle (Fig. 8).

Change in upper-plate forcing
Changing the Clapeyron slope by (de)hydration may not be a common mechanism for a mode change, and changing slab age may not be effective enough. Since most of the stagnant slabs are associated with significant trench retreat (Čížková and Bina, 2013;Tagawa et al., 2007;Torii and Yoshioka, 2007) (Section 3.3), we next test how a change in the trench mobility forced by the upper (rather than subducting) plate will affect slab-transition zone interaction.
To change the trench mobility, we change the surface velocity boundary condition of the overriding plate (from x = 0 to 3000 km), from free-slip to no-slip at ∼40 Myr, starting from simulation 8 with a subducting plate of 150 Myr. Fig. 9 shows the difference between the reference model, where the trench is free to move throughout the simulation (Fig. 9a), and the case where the mobility of the overriding plate is fixed after ∼40 Myr (Fig. 9b). As soon as the trench motion is inhibited, the slab be-gins to sink under a steep angle and starts to penetrate into the lower mantle ( Fig. 9e-f) vertically below the trench, leaving the previous stagnant section trapped in transition zone, from where it will eventually be dragged down by the sinking slab portion. The opposite change is also possible. If we start from Simulation 7 with a plate that is 150 Myr old but with a fixed upper plate, this induces the slab to sink almost vertically into the lower (Fig. 9c, h). As soon as the upper plate is left free to move, the slab starts to lie down and stagnates in the transition zone (Fig. 9d, l).
These changes in upper plate mobility are end-member scenarios, because in a global system, plates are hardly ever fully free (except possibly in the case of an actively opening back arc) or fully fixed. The models show that such strong changes in upper plate forcing are very efficient in inducing a mode change, within few Myr, both from penetrating to stagnating and from stagnating to penetrating.

Possible role of rheology
In the models presented for evaluating the interaction modes of a subducting slab with the mantle transition zone, we assumed that the upper mantle is deforming only by diffusion creep. However, there is evidence that the upper mantle, at the least in the upper 250 km, may deform preferably by dislocation creep (Hirth and Kohlstedt, 2003;Karato and Wu, 1993). The strain-rate dependence of dislocation creep would lead to a weaker asthenosphere which will induce a stronger decoupling of the subducting plate from the mantle which would enhance the tendency of the slabs to sink, driving less trench retreat (e.g., Capitanio et al., 2007). Previous numerical models show that the mobility of the trench may also be reduced by increasing the viscosity of the weak crustal layer (Čížková and Bina, 2013) and by the use of a free-surface boundary instead of a free-slip condition (Garel, pers. comm. 2016). Such rheological changes might enhance the tendency of the slabs to penetrate, and thus somewhat shift the boundaries in Fig. 3b, but would not change the main results of this paper.
In Section 4.2 we showed that younging of the subducting plate may not be effective enough to induce penetration of an initially stagnant slab. In the presented models we do however not account for mechanisms that can induce weakening of the plate at greater depth. For example, if the plastic yield strength depends on temperature (as in Peierls mechanism, Karato et al., 2001), then the slab could weaken as it spends time in the transition zone. This mechanism has been suggested to facilitate folding and buckling of the slab in the transition zone and lower mantle (Garel et al., 2014), a process that may promote slab penetration (Garel et al., 2014;Goes et al., 2008). To test the potential effect of a temperature-dependent yield stress, we redid the simulation shown in Fig. 7, but now using a slab that weakened in the transition zone, by reducing the maximum cutoff viscosity below the ol-wd transition depth. Despite the fact that in this model with drastic slab weakening more deformation is observed when the youngest slab reaches the transition zone, the initially stagnant slab still does not easily penetrate into the lower mantle (see Video S7 in the supplementary material). As Č ížková et al. (2002) previously found, slab strength in the transition zone appears to play only a secondary role in controlling the mode of subductiontransition zone interaction.

Mechanisms for switching slab-transition zone interaction on Earth
Can these model results help explain any observations of subduction evolution? The tomographic images of the Tonga and Japan slabs look like they switched from penetrating to stagnation with trench retreat (e.g., Li et al., 2008). For Tonga,  linked this to back-arc opening (i.e., a changed overriding plate forcing). In addition, back-arc opening has been previously linked to the age of the downgoing plate (Molnar and Atwater, 1978), i.e. to a change in the subducting plate forcing. For Japan, the slab also aged and this too may have induced back-arc extension (Goes et al., 2008). So these appear to be cases where the dominant forcing is by subducting plate buoyancy, although upper-plate weakening in response to the slab pull (e.g., Magni et al., 2014) aids and possibly forms the trigger of the switch in modes.  Goes et al. (2008) proposed that the subduction of very young lithosphere caused lower mantle penetration for the early Cenozoic phase of subduction below Japan and for much of the Cenozoic Farallon subduction below Central and northern South America. Our models indicate that indeed one might expect young plates to cause penetration, although this is most efficient if the subducting plate is young at the start and less so once a stagnant slab has been formed.
The Tethys slab has been inferred to go from a scenario of significant back-arc basin formation (to account for current slab volume and flattened geometry in lower-mantle tomography) to current lower-mantle sinking (Hafkenscheid et al., 2006). The material in the shallow lower-mantle is estimated to have subducted around 50 Myr, which more or less coincides with the time of India-Asia collision, i.e. a change in subducting plate buoyancy more extreme than what occurs by changing plate age may have triggered this switch.
In the above discussed examples, evidence points to an important role of subducting plate buoyancy. In our models changing upper plate forcing appears to be the most efficient mechanism to induce a mode switch. However, note that the Eurasian plate is very large, i.e., difficult to shift by a single slab, and appears to not have moved very much through the Cenozoic. Yet, this has not impeded old western Pacific slabs from flattening in the transition zone. So upper plate strength may limit how efficiently the upper plate can impede trench motion. Nonetheless, it is possible that large-scale plate reorganisations, such as those associated with amalgamation and break-up of Pangaea, led to phases of strong upper-plate forcing that may have triggered mode switches.

Conclusions
In this study, using a series of 2D self-consistent thermomechanical subduction models, we investigate which parameters are most efficient in allowing slabs to penetrate into or stagnate above the lower mantle and if there are plausible mechanisms to change the slab-transition zone interaction from one mode to the another.
A systematic study shows that within a plausible range of olivine-system Clapeyron slopes and upper-lower-mantle viscosity jumps, there are many cases where more buoyant and weaker slabs (younger plates) are more prone to penetration than dense and strong (old) slabs. The coexistence of the both modes disappears for the strongest Clapeyron slopes (−3 MPa/K for postspinel), which favour slab stagnation, whereas the lowest Clapeyron slope (−0.5 MPa/K for post-spinel) allows plates of all ages to sink into the lower mantle.
A hampering of sinking into the lower mantle appears to be the primary factor controlling how slabs interact with the transition zone. However, the increase of lower mantle viscosity over the range considered (factor 5 to 30) is by itself not sufficient to induce slab stagnation. Trench migration clearly also plays a role in slab stagnation in the transition zone, however there is no simple threshold trench migration rate above which (or dip below which) slab flattening occurs.
An efficient and quick mechanism for switching the subducting mode (few Myr in response to an instantaneous change) is a strong change in upper plate forcing, e.g. from free to fixed which can induce an initially stagnant slab to sink. Whether actual upper plate force changes would be strong enough to induce such rapid changes is unclear from observations, but possibly large-scale plate reorganisations might be able to trigger mode switches.
A decrease in the Clapeyron slope of the post-spinel transition due to slab dehydration also seems to be a possible mechanism to switch from stagnation to penetration. Inducing a change from penetration to stagnation is easier than vice versa. A change in subducting plate age is effective if the change is from young to old, which allows an initially penetrating slab to stagnate into the transition zone. But a change from old to young results only in a slow change of slab geometry, and in most of our models is not enough to overcome the effect of the previous history of subduction.
Changes in subducting plate age at the trench occur very commonly throughout the history of subduction zones, and in several cases increases in subducting plate age or buoyancy appear to correlate with changes in subduction mode. Hence, this may be the most common mechanism to trigger a change in mode from penetrating to the preferred free-subduction mode by trench retreat that leads to slab flattening in the transition zone.