The effect of metastable pyroxene on the slab dynamics

Seismic studies show that some subducting slabs penetrate straight into the lower mantle, whereas others seem to flatten near the base of the mantle transition zone. Slab stagnation is often attributed to an increase in viscosity and phase transformations in the olivine system. However, recent mineral physics studies showed that due to extremely low transformational diffusion rates, low‐density metastable pyroxene may persist into the transition zone in cool slabs. Here we use a dynamically fully self‐consistent subduction model to investigate the influence of metastable pyroxene on the dynamics of subducting oceanic lithosphere. Our results show that metastable pyroxene affects slab buoyancy at least as much as olivine metastability. However, unlike metastable olivine, which can inhibit slab penetration in the lower mantle only for cold, old, and fast slabs, metastable pyroxene is likely to also affect sinking of relatively young and slow slabs.


Introduction
Seismic studies reveal a large range of subducting slab morphologies. Most notable are the variations near the base of the mantle transition zone, between 400 and 800 km depth, where some subducted slabs seem to flatten to form stagnant slabs, for example, beneath the Izu-Bonin region and Tonga, while others sink further into the lower mantle, for example, beneath the Mariana Arc and Central America [van der Hilst et al., 1997;Fukao et al., 2001;Grand, 2002;Li et al., 2008;Kawakatsu and Yoshioka, 2011;Fukao and Obayashi, 2013]. The volumes and duration of slab stagnation may have important implications for the efficiency of the transport of heat and chemistry through the mantle [van Keken et al., 2002]. Hence, numerous studies have investigated what governs the ability of the slab to penetrate or to lie down in the transition zone. Christensen [1996] showed that high trench migration rates promote slab flattening in the transition zone. Trench motion is, in turn, encouraged by high downgoing plate age and strength and may be impeded by coupling with the overriding plate [Bellahsen et al., 2005;Tagawa et al., 2007;Garel et al., 2014;Sharples et al., 2014]. Another relevant factor is slab dip angle, where a low dip tends to produce slab stagnation [Torii and Yoshioka, 2007]. A low dip angle may be the result of fast trench retreat but can also result from interaction with the upper plate [Tagawa et al., 2007]. Furthermore, slabs can encounter resistance to further sinking due to density and viscosity increases associated with phase changes of the major mantle minerals which encourage slab stagnation [Christensen and Yuen, 1984;Gurnis and Hager, 1988;Torii and Yoshioka, 2007;Yoshida, 2013;Garel et al., 2014].
Phase transformations result in additional local buoyancy effects. At equilibrium conditions, the exothermic olivine-wadsleyite (ol-wd) phase transition (around 400 km depth) occurs at shallower depths in the cool slab than in the surrounding mantle, adding extra negative buoyancy to the downgoing plate and encouraging sinking, while the endothermic ringwoodite to perovskite and magnesiowüstite (rw-pv + mw) transition (near 700 km depth) leads to an increment of positive buoyancy hampering slab penetration into the lower mantle [Christensen and Yuen, 1984;Billen, 2010]. However, if the ol-wd transition is kinetically slow at slab conditions, this transition also provides positive slab buoyancy and may retard or prevent slabs from penetrating into the lower mantle, although olivine metastability is not likely to affect younger slabs [Marton et al., 1999;Schmeling et al., 1999;Tetzlaff and Schmeling, 2000;Bina et al., 2001;Mosenfelder et al., 2001;Bina and Kawakatsu, 2010].

Citation:
Agrusta, R., J. van Hunen, and S. Goes (2014) , 1977]. At cold slab temperatures, this dissolution of pyroxene is kinetically inhibited, and hence, pyroxene may also survive as a metastable phase [Hogrefe et al., 1994;Nishi et al., 2008Nishi et al., , 2013. The dissolution rate is controlled by how fast the majorite component can diffuse through the garnet structure, which has recently been found to be one of the slowest diffusion processes in the Earth's mantle [van Mierlo et al., 2013]. The slow diffusion can promote the presence of low-density pyroxene even at relatively high slab temperatures and can thus contribute to the deceleration and potential stagnation of slabs in the transition zone [Bina, 2013a[Bina, , 2013bvan Mierlo et al., 2013].
No dynamic models have yet investigated what the influence of metastable pyroxene on the dynamics of subducting oceanic lithosphere is and what the relative importance of the olivine phase transitions and a kinetically hindered pyroxene transformation. This is the aim of our study.

Governing Physics
The slab dynamics are studied using the finite-element code Citcom [Moresi and Gurnis, 1996;Zhong et al., 2000;van Hunen et al., 2005]. The code solves the incompressible extended Boussinesq approximation of the governing equations of thermally driven Stokes' convection [Christensen and Yuen, 1985], without radiogenic heating. The mantle phase transitions are included in the model using a harmonic phase function [Čížková et al., 2007]. The effective temperature, pressure, and stress-dependent viscosity is defined as a composite rheological model in the upper mantle [Karato et al., 1995], combining diffusion and dislocation creep [Hirth and Kohlstedt, 2003] and a brittle deformation mechanism [Di Giuseppe et al., 2008] and pure diffusion creep in the lower mantle (see Text S1 and Table S1 in the supporting information for details).
The potentially kinetically slow phase transitions have been modeled using a growth rate model for the ol-wd transformation [Rubie and Ross, 1994] and a diffusion model for pyroxene-garnet (px-gt) [Van Mierlo et al., 2013].
For the growth rate model, the volume fraction ξ transformed at a given temperature T in a time t is where ẋ is the temperature-dependent growth rate of wadsleyite at the interface between olivine grains with diameter d g [Rubie and Ross, 1994], ΔH ol-wd is activation enthalpy for growth for the ol-wd transformation, ΔG is the free-energy change of reaction, k is a constant, and R is the gas constant. The volume fraction of pyroxene that has transformed into majorite-garnet is with D the diffusion rate of pyroxene into garnet, t diff the time necessary to diffuse the majorite component at a distance of a half grain size d g [Van Mierlo et al., 2013], and ΔH px-gt is activation enthalpy for this transformation (see Text S1 and Table S1 for details and parameter values).

Model Setup
The model domain is 5000 km wide and 2000 km high, and the box is discretized into 1088 × 219 elements. The element sizes range from 3 km to 9 km. Vertically, the grid is refined between 0 and 210 km depth (lithosphere and transition to asthenosphere) and around the mantle phase transitions between 354 and 432 km and between 648 and 720 km depth. The mesh is also refined horizontally around the boundary between the subducting and the overriding plate between x = 810 and 3204 km.
The mechanical boundary conditions are free slip along all boundaries, so mantle and slabs are driven purely internally by buoyancy forces only. The thermal boundary conditions are the following: constant temperature at the surface (273 K) and the bottom (potential temperature of 1573 K) and zero heat flux on the vertical boundaries. The subducting plate extends from a ridge located in the upper left corner to the trench initially situated at x = 2500 km. The plate's initial geotherm follows a half-space cooling model to produce plate ages at the trench between 50 Myr and 150 Myr. 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 50 Myr. To avoid Geophysical Research Letters

10.1002/2014GL062159
complications associated with subduction initiation, a proto slab long enough to allow self-sustained subduction is prescribed as an initial condition. The initial slab has a curvature radius of 500 km and an initial slab tip depth of 180 km. On the top of the subducting plate, a 15 km wide low-viscosity layer (10 20 Pa s) is present down to 200 km depth to facilitate the decoupling of the two plates (see Text S1 and Figure S1 for details).
Our focus is on the effect of the phase transformations, but we test the sensitivity to rheology by using three different realistic viscosity profiles. In all cases, the rheological parameters have been chosen to obtain mantle viscosity at the base of the upper mantle of 10 21 Pa s [Karato and Wu, 1993;Hirth and Kohlstedt, 2003;Mitrovica and Forte, 2004]. The average viscosity of the lower mantle has been suggested to range from 2 to 10 times of the upper mantle viscosity [Mitrovica and Forte, 2004;Quinteros et al., 2010;Čížková et al., 2012]. We test two profiles for the lower mantle viscosity, one in which there is no jump across the upper-lower mantle transition (profiles i and ii, Table 1) and another where the viscosity increases by a factor of 2 (profile iii, Table 1). A viscosity jump of 5 has been tested, but this does not allow the slabs to penetrate into the lower mantle, even for a model where thermodynamic equilibrium is assumed for the phase transformations. For the upper mantle, we test profiles with (ii) and without (i and iii) dislocation creep (Table 1).

Mantle Phase Transitions
To test the effect of pyroxene metastability, we used three compositions: a reference case of 100% olivine, a mixture of 10% pyroxene and 90% olivine (PX1), and one of 20% pyroxene and 80% olivine (PX2). These amounts of pyroxene span the range of pyroxene volume fractions expected in a peridotitic/harzburgitic mantle composition [Xu et al., 2008].
The models include the ol-wd and rw-pv + mw phase transitions at equilibrium depths of 410 km and 670 km, respectively. The density contrast associated with ol-wd and rw-pv + mw transformations is 250 kg m À3 and 350 kg m À3 [Xu et al., 2008]. The Clapeyron slope of the exothermic ol-wd transition has been estimated to be between 2 MPa K À1 and 4 MPa K À1 [Bina and Helffrich, 1994;Katsura et al., 2004], whereas for the endothermic rw-pv + mw, it ranges from À2 MPa K À1 to À0.5 MPa K À1 [Katsura et al., 2003;Fei et al., 2004;Litasov et al., 2005].
In this study, we use a slope of 3 MPa K À1 for the ol-wd transition and slopes of À2 MPa K À1 and À1 MPa K À1 for the rw-pv + mw transformation (Table 1). The density increase related to the wadsleyte-ringwoodite transition at~550 km is relatively small 2% [Marton et al., 1999;Bina et al., 2001]. Therefore, we approximate both the exothermic transitions ol-wd and wd-rw together by one single transition [e.g., Tetzlaff and Schmeling, 2009].
The px-gt solid solution at equilibrium conditions starts at~200 km and completes at~450-480 km depth [Akaogi and Akimoto, 1977;Irifune and Isshiki, 1998;Bina and Helffrich, 2014]. We model the equilibrium px-gt transition to start at 200 km and be completed at 400 km depth. We use a shallower depth to complete the transformation for convenience so that we can keep the density contrast below the ol-wd transition the same for pure olivine models and models with olivine and pyroxene. Extending the transition a further~80 km The reference viscosity profile i has no jump across the upper-lower mantle transition and no dislocation creep in the upper mantle; profile ii is as i but with upper mantle dislocation creep; and profile iii as i but with a factor 2 jump in viscosity at the upper-lower mantle boundary.

Geophysical Research Letters
10.1002/2014GL062159 deeper will not have a significant effect on the results, because the additional volume fraction of px produced in this interval would be smaller than 1%. The px-gt Clapeyron slope is positive [Akaogi and Akimoto, 1977;Akaogi et al., 1987] but is set to a zero in the models. Given the significant metastability of this transition in all models, this simplification has a negligible effect. At depths greater than~660 km garnet transforms into perovskite [Frost, 2008], but this transition is spread out over a larger depth interval than rw-pv + mw, and the contributions to local density anomalies and hence buoyancy forces should be considerably more diffuse [Christensen and Yuen, 1985;Bina et al., 2001]. Hence, we do not include the gt-pv transition in the models.
The density increase for pure px during the phase transition of px-gt is assumed to be 300 kg m À3 [Xu et al., 2008], and its effect on slab buoyancy has been scaled depending on the amount of pyroxene considered, to 30 kg m À3 or 60 kg m À3 for 10% or 20% of pyroxene, respectively. At pressures above 20 GPa, and low-temperature conditions, pyroxene transforms into its higher-pressure polymorph akimotoite [Hogrefe et al., 1994], resulting in a large-density increase of~16%. We model this transformation simply by transforming any remaining pyroxene to akimotoite below 670 km depth. We compensated for the addition of pyroxene into the mantle composition by reducing the density contrast of the ol-wd transition to 220 kg m À3 and 190 kg m À3 , for models with 10% and 20% of pyroxene, respectively.
The coefficients used to define the growth rateẋ of the ol-wd transition (equation (1)) have been taken from laboratory experiments [Rubie and Ross, 1994] and from previous numerical studies [Tetzlaff and Schmeling, 2009] (see Table S1). The activation enthalpy (ΔH px-gt ) and preexponential coefficient (D 0 ) used to define the diffusion rate D of the px-gt transition (equation (2)) are from a recent laboratory experiment [Van Mierlo et al., 2013]. We chose their preferred values, ΔH px-gt = 291 ± 51 kJ mol À1 and D 0 = 2.3 × 10 À11 m 2 s À1 , for the models. Uncertainties in ΔH px-gt ± 51 kJ mol À1 and the confidence interval for D 0 lie between the two limiting values 9.4 × 10 À14 and 5.4 × 10 À10 m 2 s À1 , leading to an uncertainty of 5 orders of magnitude in the diffusion rates [Van Mierlo et al., 2013]. We will discuss how that may affect our results in section 4.

Results
A total of 48 simulations is used to study the effects of metastable phases in the descending slab (Table 1). There are 12 groups of models for different viscosity profiles, rw-pv + mw Clapeyron slope and age of the subducting plate. For each of these we perform the following four simulations: 100% olivine without kinetically slow phase transition (EQ), a kinetically slow ol-wd model (OL), and two with kinetically slow px-gt and ol-wd transitions (PX1 with 10% pyroxene and PX2 with 20% pyroxene).
We illustrate the resulting subduction dynamics for the reference Subd1b set of models, i.e., for a 100 Myr old subducting plate, Newtonian viscosity throughout the mantle, without a viscosity jump at the upper-lower mantle transition, and with a lower bound of the rw-pv + mw Clapeyron slope. Figures 1a and 1b show the time evolution for the thermodynamic equilibrium case EQ and case PX2, respectively, while Figure 1c compares slab sinking through time, evolution of buoyancy forces, and convergence velocity as a function of max slab depth. Negative thermal buoyancy force drives the slab sinking in all cases. In the EQ case, the convergence accelerates after the slab reaches the ol-wd phase transition (Figure 1a, left), due to additional negative buoyancy from the upward deflection of the phase boundary. As this slab arrives at the base of the upper mantle (Figure 1a, middle), the downward equilibrium deflection of the rw-pv + mw transition reduces the forces acting to drive subduction, and the convergence slows down (Figure 1c). In this model, a viscosity jump at the upper-lower mantle transition is absent, and the resisting force from the endothermic phase transition cannot prevent the slab from penetrating deeply into the lower mantle (Figure 1a, right). The slowdown of convergence as the slab sinks deeper is due to increased resistance as viscosity gradually increases with depth.
If a kinetically slow ol-wd phase transition is assumed (equation (1)) (OL, Figure 1c), the slab behavior is similar to the equilibrium case, due to the fairly insignificant reduction of buoyancy forces from~60 × 10 12 N m À1 tõ 55 × 10 12 N m À1 , which reduces the convergence velocity from~13 cm yr À1 to~11 cm yr À1 (Figure 1c).
The persistence of metastable pyroxene has a much stronger effect (Figures 1b and 1c). Inside the PX1 and PX2 slabs, metastable pyroxene extends down to 670 km depth, before transforming into its higher-pressure polymorph. By contrast, olivine metastability is limited to the cold inner core of the slab with temperatures below 650°C (Figure 1b). Depending on the amount of metastable pyroxene present, the slab either penetrates into the lower mantle (PX1) or the slab dip becomes smaller (PX2, Figure 1b, middle), ultimately promoting stagnation in the transition zone (Figure 1b, right). For PX2, the resulting buoyancy forces are only F tot ≈ À30 × 10 12 N m À1 when the slab tip is at the rw-pv + mw phase boundary, compared with F tot ≈ À60 × 10 12 N m À1 for the EQ case (Figure 1c, middle). This "parachute effect" of metastable pyroxene is highlighted by the strong reduction of the convergence velocity when the slab reaches 670 km depth, from 13 cm yr À1 to~7 cm yr À1 (Figure 1c, right).
Other initial slab ages, deformation mechanisms in the upper mantle and lower mantle properties, can lead to slab dynamics and subduction styles different from those shown in Figure 1 (as previously shown by Christensen [1996]). We show some examples in Figure 2. Metastable pyroxene can promote slab stagnation not only for old but also for younger slabs, as illustrated in Figure 2a, where a slab with an initial age at the trench of 50 Myr descends through a composite rheology upper mantle (Subd2a). As illustrated by case Subd3b, shown in Figure 2b, for a set of 100 Myr old slabs, a more negative rw-pv + mw Clapeyron slope helps the stagnation also for PX1, whereas for the cases without pyroxene metastability (EQ and OL), the slab still penetrates into the lower mantle, albeit with a flattening in the transition zone. A viscosity increase in the lower mantle also promotes slab stagnation for PX1 models, except for the oldest (150 Myr) slabs, for which stagnation and absence of penetration deep into the lower mantle are only observed for PX2 (Figure 2c, Subd4c).
An overview of how metastable phases can change slab dynamics is given in Figure 3, which shows sinking slab velocities (averaged between 500 km and the maximum depth of the slab, Figure 3a) and the maximum slab depth (Figure 3b) for all simulations.
All models with 20% pyroxene (PX2) have a sinking velocity close to zero and a maximum depth that does not exceed 800 km, clear evidence of slab stagnation. The only exception is the Subd2b where the slab reaches the depth of~900 km. Although the difference between the models with dislocation creep (Subd2) and without it (Subd1) is not large, the slabs sink somewhat faster when upper mantle dislocation creep is included. As a result the Subd2b case retains more negative buoyancy than case Subd1b. The same models with younger (50 Myr, Subd2a) and older (150 Myr, Subd2c) slabs both result in a flat slab without penetration in the lower mantle. The slab stagnation for the younger slab is attributed to positive buoyancy associated with the metastable pyroxene, whereas for the oldest plate, positive buoyancy is enhanced by metastable olivine that due to the high convergence velocity affects a larger part of the slab's cold core than for younger slabs. Increasing the resistance to sinking into the lower mantle through a more negative rw-pv + mw Clapeyron slope allows the slab to lie down more easily at the transition even with less metastable pyroxene. Despite a larger Clapeyron slope, the effect of only considering metastable olivine does not show a significant contribution to stagnation for 150 Myr lithosphere (subd3c OL), as was reported in previous studies [Tetzlaff and Schmeling, 2000;Mosenfelder et al., 2001]. This discrepancy is due to our slower slab sinking velocity than in the mentioned previous studies, where velocities faster than 15 cm yr À1 occurred. Instead, increasing the lower mantle viscosity only favors slab stagnation for the PX1 case for younger subducting slabs (Subd4a).

Discussion and Conclusion
In this study, the effect of the metastable persistence of pyroxene on slab stagnation in the mantle transition zone has been analyzed using fully dynamically self-consistent subduction models. Pyroxene can survive as a metastable phase in cold slabs and due to its low density favors slab deceleration and stagnation in the  transition zone, as was proposed by Van Mierlo et al. [2013]. Unlike metastable olivine, which can inhibit slab penetration in the lower mantle only for cold, old, and fast slabs [Tetzlaff and Schmeling, 2000;Bina et al., 2001], metastable pyroxene is likely to also affect sinking of relatively young and slow slabs.
Old slabs are normally associated with more negative buoyancy than young slabs and from a dynamical perspective are therefore expected to sink faster and deeper into the mantle. However, observations suggest that many slabs, including old ones such as Izu-Bonin, Tonga, and Japan, stagnate in the transition zone [Fukao and Obayashi, 2013]. Our model results show that the combined effects of the pyroxene and olivine metastability can lead to the stagnation of both younger and the oldest slabs ( Figure 3).
The reaction kinetics of the ol-wd and px-gt phase transitions are likely to be dependent on several parameters that are not considered in this study. Large grain sizes can enhance metastability significantly and increase the positive buoyancy forces acting on the slab. On the other hand, slabs may transport water deep in the mantle and into the transition zone [Ferot and Bolfan-Casanova, 2012;Bercovici and Karato, 2003], which could speed up the kinetics by perhaps several orders of magnitude [Rubie, 1986]. We tested those effects for model Subd3b for the PX2 case by changing the preexponential factor D 0 . These tests show that faster reaction kinetics confine the metastable pyroxene wedge more to the colder inner part of the slab. A 2 order-of-magnitude increase in diffusion rates (2.3 × 10 À9 m 2 s À1 ) still leaves enough metastable pyroxene in the slab to promote slab stagnation, but if diffusion rates are increased by another order of magnitude, then pyroxene metastability becomes incapable of producing slab stagnation. This illustrates that the results related to slab stagnation by pyroxene metastability in this study are robust over a significant range of phase change kinetics. Therefore, metastable pyroxene in subducting slabs is likely to play an important role in their stagnation in the mantle transition zone.  To compare the different simulations, we measured these quantities for each model at a time that is 3 times the time that the slab takes to reach the 670 km depth. The symbols distinguish groups with different upper mantle deformation mechanism, Clapeyron slope of the rw-pv + mw transition, and lower mantle viscosity. For each group, cases a, b, and c correspond to subducting plates of 50, 100, and 150 Myr, respectively. See Table 1 for details.