## Summary

In this study, we calculate timescales for the growth of gravitational instabilities forming in the sediment layer on the downgoing slab at subduction zones. Subducted metasediments are buoyant with respect to the overlying mantle and may form diapirs that detach from the slab and rise upwards into the mantle wedge. We use a particle-in-cell, finite-difference method to calculate growth rates for instabilities forming within a buoyant, wet-quartz metasediment layer underlying a dense mantle half-space composed of wet olivine. These growth rates are used to determine where sediment diapirs initiate and detach from the slab over a range of subduction zone thermal structures. We find that, given a sufficient layer thickness (200–800 m, depending on slab-surface and mantle-wedge temperatures), sediment diapirs begin to grow rapidly at depths of ∼80 km and detach from the slab within 1–3 Myr at temperatures ≤900 °C and at depths roughly corresponding to the location of the slab beneath the arc. Diapir growth is most sensitive to absolute slab temperature, however it is also affected by the viscosity ratio between the sediment layer and the mantle wedge and the length-scale over which viscosity decays above the slab. These secondary affects are most pronounced in colder subduction systems with old slabs and faster subduction rates. For a broad range of subduction zone thermal conditions, we find that diapirs can efficiently transport sediments into the mantle wedge, where they would melt and be incorporated into arc magmas. Thus, we conclude that sediment diapirism is a common feature of many subduction zones, providing a potential explanation for the ‘sediment signature’ in the chemistry of arc magmas.

## Introduction

It is well established that subducted sediments are incorporated into magmatic systems beneath arcs and that these sediment melts influence the production and chemical composition of new continental crust. Despite an abundance of isotopic and trace-element evidence for a sedimentary component in arc magmas (e.g. Armstrong 1968, 1991; Sun 1980; Karig *et al.* 1981; White *et al.* 1985; Tera *et al.* 1986; Plank & Langmuir 1993, 1998), as well as geochemical (e.g. Elliott *et al.* 1997; Hawkesworth *et al.* 1997; Behn *et al.* 2011) and geophysical (e.g. Vogt 1974; Marsh 1979; de Bremond d’Ars *et al.* 1995) constraints on subduction zone melting processes, the mechanism for sediment melting and the transport of these melts to the surface remains uncertain. Two primary classes of models have been proposed: (1) slab-surface melting of sediments with subsequent porous or diapiric flow of the melts to the surface (e.g. Spiegelman & McKenzie 1987; Hall & Kincaid 2001) and (2) ‘cold' diapirism in which the sediment layer detaches from the downgoing slab and melts as it ascends through the hot mantle wedge (e.g. Kelemen *et al.* 2003). In the latter case, diapirism can be driven by hydration and partial melting of the mantle wedge at the slab surface (e.g. Gerya & Yuen 2003; Gerya *et al.* 2006; Castro & Gerya 2008) or by the intrinsic buoyancy of the sediment layer itself (Currie *et al.* 2007; Behn *et al.* 2011).

Modern subduction zone thermal models that include temperature- and stress-dependent viscosities (van Keken *et al.* 2002; Kelemen *et al.* 2003; Conder 2005; Arcay *et al.* 2007; Wada & Wang 2009; England & Katz 2010; Syracuse *et al.* 2010) and recent geothermometry (Plank *et al.* 2009) predict slab-surface temperatures in excess of the fluid-saturated solidus for metasediments at slab depths corresponding to the location of the volcanic arc (∼650 °C at 2 GPa; Nichols *et al.* 1994; Schmidt *et al.* 2004), implying that sediments may melt at the slab surface. However, trace-element-depletion trends in metasedimentary rocks that endured subduction to pressures in excess of 2 GPa suggest that key trace elements associated with the observed ‘sediment-melt signature' are not released until temperatures exceed ∼1050 °C (Behn *et al.* 2011)-significantly hotter than slab-surface temperatures at similar pressures in subduction zone thermal models (e.g. van Keken *et al.* 2002; Wada *et al.* 2008). Behn *et al.* (2011) interpreted this discrepancy to imply that, although melting may commence at lower temperatures, the sediment-melt signature is retained in accessory mineral phases (e.g. phengite, monazite, allanite) to higher temperatures, and proposed that sediments detach from the slab and rise diapirically into the overlying mantle wedge. These diapirs would undergo melting as they ascend into the hot mantle wedge (e.g. Gerya & Yuen 2003), rapidly cycling sediment-derived melts and volatiles into arc magma systems.

The growth of diapiric, Rayleigh-Taylor-type instabilities in a buoyant sediment layer underlying a dense mantle half-space is controlled by the layer buoyancy and the absolute and relative viscosities of both materials (e.g. Rayleigh 1883; Taylor 1950; Chandrasekhar 1961; Whitehead & Luther 1975). Increasing density contrasts and/or layer thicknesses increase buoyancy, promoting faster instability growth. Diapirs form more readily at lower absolute viscosities, and the dependence of viscosity on strain rate has a significant effect on instability growth. Specifically, diapirs in linear (i.e. Newtonian) rheologies grow exponentially whereas materials with strain-rate-dependent, nonlinear (i.e. non-Newtonian) rheologies weaken as they deform, giving rise to superexponential growth (e.g. Conrad & Molnar 1997; Houseman & Molnar 1997; Jull & Kelemen 2001). The relative viscosities of the sediment layer and mantle half-space also play a key role in controlling the growth of diapirs, with a larger layer-to-half-space viscosity ratio promoting faster growth (e.g. Conrad & Molnar 1997; Houseman & Molnar 1997; Jull & Kelemen 2001).

Previous scaling analyses (Behn *et al.* 2011) and geodynamic models (Currie *et al.* 2007), suggest that subducted metasediments are buoyant in the mantle and may form diapiric instabilities. Assuming a fixed ratio of sediment and mantle viscosities, background strain rate and density contrast, Behn *et al.* (2011) showed that, for typical subduction zone thermal structures (Wada & Wang 2009), instabilities would form over a broad range of sediment layer thicknesses. However, the densities and viscosities of the sediment layer and overlying mantle wedge are expected to evolve independently as the slab warms during subduction, causing the key controls on diapir growth (i.e. density contrast, viscosity ratio and half-space viscosity structure) to vary with depth along the slab interface. Currie *et al.* (2007) included these depth-dependent effects in subduction-zone-scale geodynamic models that characterized the fate of subducted sediments. Although this work explored the sensitivity of diapir growth to sediment and mantle material parameters, the effect of far-field parameters (e.g. slab age/temperature, convergence rate and depth of mechanical decoupling between the slab and mantle wedge) was not investigated. Further, these models lacked the resolution required to accurately quantify the growth of small-scale instabilities in sediment layers <2 km thick.

In this study, we used a finite-difference method (Gerya 2010) to determine timescales for the growth of individual instabilities initiating within a wet-quartz metasediment layer underlying an olivine mantle wedge, with independent temperature- and strain-rate-dependent viscosities in each layer. To predict the location of diapir formation, we integrated these timescales over time- and depth-varying thermal structures and background strain rates from a range of modern subduction zone thermal models (Wada & Wang 2009; Syracuse *et al.* 2010). We find that, while diapir growth depends strongly on the temperature of the subducting slab and mantle wedge, sediments can form diapirs that detach from the slab at locations corresponding to the location of arc volcanism over a wide range of slab thermal structures. In ‘hot' subduction zones (young slabs and/or slow convergence rates, for example, Cascadia), effective viscosities are relatively small, enabling instabilities to form over a broad range of sediment thicknesses and mantle-wedge thermal and strain rate structures. In colder subduction zones (old slabs and/or fast convergence rates, for example, Izu-Bonin), effective viscosities are greater, slowing instability growth, and secondary controls such as the ratio of sediment to mantle viscosities and the length scale of viscous decay in the mantle play a more significant role in the development of sediment diapirs.

## Methods

We modelled the growth of metasedimentary diapirs in subduction zones as Rayleigh-Taylor-type instabilities forming in a buoyant, water-rich, quartz layer that underlies a dense, hydrous, olivine half-space (Fig. 1). We first used a 2-D, particle-in-cell, finite- difference method (Gerya 2010) to calculate diapir growth rates for Newtonian and non-Newtonian viscosities at a range of temperatures and strain rates. From these growth rates, we calculated timescales for diapir formation, which we quantified through the ‘instability time' (defined below for Newtonian and non-Newtonian rheologies). Next, we exploited previously derived scaling relationships found by linear analysis and numerical experiments (e.g. Rayleigh 1883; Taylor 1950; Chandrasekhar 1961; Whitehead & Luther 1975; Conrad & Molnar 1997; Houseman & Molnar 1997; Molnar *et al.* 1998) to scale instability times for the range of temperature and strain rate conditions in subduction-zone thermal models (Wada & Wang 2009; Syracuse *et al.* 2010). Finally, we integrated these instability times over the history of subduction to estimate the timescales for diapir formation in individual thermal models. Comparing calculated timescales from a range of subduction zone settings reveals the sensitivity of sediment-diapir growth to key parameters such as sediment-layer thickness, slab temperature (slab age and subduction rate) and the strength of the slab/mantle interface (e.g. Wada & Wang 2009; Syracuse *et al.* 2010). The symbols and scaling terms defined below and used throughout this paper are listed in Table 1.

### Rheologic flow law

The absolute and relative viscosities of the sediments and overlying mantle determine how easily a buoyant sediment layer can flow upward into the dense mantle. In most rock-forming minerals, viscosity decreases with increasing strain rate and temperature (e.g. Hirth & Kohlstedt 2003). To model deformation with this behaviour, we use the flow law (e.g. Kirby 1983)

where*A*is an empirical constant, is the differential stress,

*n*is the stress exponent,

*Q*is the activation energy,

*R*is the gas constant (= 8.31 J K

^{−1}mole

^{−1}) and

*T*is temperature. Assuming that deviatoric stress is related to strain rate by the constitutive relationship we define an effective viscosity as where is the second invariant of the strain-rate tensor . The rheologic parameter

*B*describes the sensitivity of viscosity to changes in strain rate and temperature and is defined (Molnar

*et al.*1998): For

*n*= 1 (Newtonian), , viscosity varies with temperature alone, and strain rate is a linear function of stress. For

*n*> 1 (non-Newtonian), viscosity depends on strain rate and temperature, causing strain rate to vary nonlinearly with stress.

### Theoretical model: scaling relationships and non-dimensionalization

Solutions found by linear analysis of the Stokes equation (i.e. the Navier-Stokes equations with inertial terms ignored) can describe the growth of infinitesimal-amplitude instabilities in a buoyant layer underlying a half-space with Newtonian or non-Newtonian viscosities in each. For instabilities to form, the layer boundary must be initially perturbed. In prior analytic solutions and in our models, the height of this initial boundary is defined by

where*Y*

_{0}is the perturbation amplitude and is the wavelength of the perturbation.

For Newtonian materials, Chandrasekhar 1961) showed that maximum instability displacement *Y* increases exponentially from this initial deflection according to

*t*is time and

*q*is a growth rate. Here, maximum vertical velocity

*w*is given by where

*w*

_{0}=

*qY*

_{0}is the initial velocity.

To quantify the timescale for instability formation, we used the ‘instability time' (e.g. Houseman & Molnar 1997), which we refer to by the symbol *t _{b}*. For the Newtonian (i.e. exponential growth) case, we define

*t*as the time at which diapir displacement reaches a height equal to the initial layer thickness

_{b}*h*, and setting

*Y*=

*h*in eq. (6) yields

Non-dimensionalization reveals the sensitivity of diapir growth to changes in layer buoyancy, initial perturbation wavelength and viscosities. By choosing the initial layer thickness (*h*), the density contrast between the layer and half-space (), and the rheologic parameter and stress exponent in the sediment layer (*B _{L}* and

*n*, respectively) as characteristic scales, we can define a characteristic timescale (Houseman & Molnar 1997)

_{L}*r*=

*B*/

_{L}*B*is the ratio of the layer to half-space (

_{H}*B*) rheologic parameters.

_{H}For cases with non-Newtonian rheologies in which *B _{L}* and

*B*are constant throughout each layer, Houseman & Molnar (1997) showed that the maximum vertical displacement of the layer and the velo-city change with time according to the superexponential functions:

_{H}*t*is defined as the time at which maximum vertical velocity accelerates to infinity (Houseman & Molnar 1997), given by As the velocity accelerates to infinity,

_{b}*Y*rapidly approaches

*h*, and the instability time given by eq. (14) can thus be compared directly to the instability time defined by eq. (8).

At depth in subduction zones, temperature increases with distance from the slab surface (e.g. Wada & Wang 2009), causing viscosities to decay exponentially (eq. 3) into the mantle wedge. This viscous decay can be quantified in terms of the *e*-folding length scale *L* (Fletcher & Hallet 1983; Conrad & Molnar 1997), defined by

*T*

_{0}and describe a linear, mantle-wedge geotherm, . The ratio of

*L*to the layer thickness

*h*also controls instability growth (Conrad & Molnar 1997; Houseman & Molnar 1997). If

*L*/

*h*is small, the half-space weakens rapidly with distance from the interface and the layer ‘senses' the lower viscosity within the half-space, shortening instability times. For such cases, Molnar

*et al.*(1998) showed that the relationships in eqs (6-14) still apply, although

*L*replaces

*h*as the controlling length scale.

In settings with strong background strain rates, such as those associated with corner flow in the mantle wedge, Molnar *et al.* (1998) showed that diapirs first grow exponentially according to eq. (6), but with viscosities given by eq. (3) with defined by the background strain rate. For this initial growth stage, Molnar *et al.* (1998) defined a second timescale

### Model parameters for subduction zones

We calculated growth rates and instability times for the range of thermal and background strain rate conditions in the subduction-zone modelling results of Wada & Wang (2009) and Syracuse *et al.* (2010). To calculate viscosities, we assume that the subducting sediments and overlying mantle are both hydrated (Hacker 2008). We apply values for rheologic parameters (*A, n* and *Q*) reported for hydrous quartz (Hirth *et al.* 2001) to the metasediment layer and those for hydrous olivine (Hirth & Kohlstedt 2003) to the overlying mantle wedge (Table 2). While slab-top conditions can exceed the fluid-saturated solidus for metasediments (Nichols *et al.* 1994; Schmidt *et al.* 2004), we assume that any melting has a negligible effect on viscosity. If melt migration is rapid compared to the timescale of the matrix flow, only a small amount (3 per cent) of melt is likely to be retained in the sediment layer and/or overlying mantle (McKenzie 1984), reducing viscosities by at most an order of magnitude (Zimmerman & Kohlstedt 2004).

Subducted metasediments are likely similar in composition to ultrahigh pressure (UHP) metapelites (Behn *et al.* 2011). These rocks have between 47 and 76 wt% SiO_{2} and, at pressures less than ∼2.5 GPa, would be composed of ∼35 per cent quartz by volume, with smaller volumes of phengite and plagioclase (Behn *et al.* 2011). Of these minerals, quartz is the weakest phase (Fig. 2a). Multiphase rocks have been shown to deform according to the rheology of the weakest phase when even small amounts (<15 per cent) of that mineral are present (Dell’Angello & Tullis 1996), and we thus model deformation in the sediment layer with a quartz rheology. At pressures greater than ∼2.5 GPa, quartz undergoes a phase change to coesite (Fig. 2b), which is stronger than quartz (Fig. 2a). As discussed below in Section 4, this viscosity change will have only a small effect on instability times for diapirs forming in subduction zones, and thus we ignore the quartz-to-coesite transition in our calculations of layer viscosity.

Hirth *et al.* (2001) used microstructures and thermochronology in naturally deformed quartzite samples, along with laboratory experiments, to constrain the rheologic parameters of quartz at strain rates of ∼10^{−14} s^{−1} and over a temperature range of ∼600-1200 °C. For strain rates typical of active deformation (10^{−15} to 10^{−12} s^{−1}), this flow law has been shown to fit stress profiles inferred from naturally deformed, metamorphosed rocks from the mid-crust (Behr & Platt 2011). These conditions are similar to the strain rates and temperatures for diapir formation in our models, suggesting that this flow law can accurately describe deformation in subducted metasediments. For olivine, Hirth & Kohlstedt (2003) used experimental results to determine rheologic parameters for temperatures of ∼1100-1300 °C and strain rates of ∼10^{−4}-10^{−6}. Viscosity profiles extrapolated from these experimental conditions to upper mantle conditions agree well with constraints on viscosities from geophysical observations (Hirth & Kohlstedt 2003), indicating that these parameters accurately describe the rheology of the bulk upper mantle over the range of conditions in the mantle wedge at subduction zones.

For our calculations of instability times for exponential growth, we imposed a constant background strain rate in the mantle wedge () based on average values in the Wada & Wang (2009) thermal models. Strain rates in the sediments depend on assumptions about the nature of mantle-wedge flow and strain localization at the slab surface. Both the Wada & Wang (2009) and Syracuse *et al.* (2010) models include a shallow (<50-90 km) region in which the subducting slab descends faster than the overall flow in the mantle wedge, requiring that strain be accommodated within a weak zone at the slab surface. This shallow ‘decoupling' creates a stagnant, cold, mantle-wedge nose beneath the forearc, a feature needed to explain surface heat flow and seismic attenuation measurements (e.g. Kneller *et al.* 2005, 2007; Abers *et al.* 2006; Wada & Wang 2009). At shallow depths, we assume that strain is accommodated within the sediment layer, and we set the background strain rate in the sediments based on the layer thickness and the slab/arc convergence velocity. At greater depths, we set the background strain rate in the sediment layer to be equal to the strain rate in the mantle wedge . That is,

*z*is depth in the subduction zone model,

*V*

_{slab}is the slab velocity, and

*z*

_{DC}is the depth of the decoupling/coupling transition in the thermal models. In this way, the strain rate in the sediment layer can differ from the strain rate governing the non-Newtonian viscosity of the overlying mantle wedge. is a value less than one that determines how much of the relative motion between the slab and mantle wedge is accommodated within sediments on the slab surface. For consistency with the parametrization of the thermal models, we set to 1.0 (full decoupling) in the Wada & Wang (2009) models and 0.98 (partial decoupling) in the Syracuse

*et al.*(2010) models.

The absolute and relative viscosities of hydrous quartz and olivine vary by several orders of magnitude over the range of subduction zone temperatures (∼200-1000 °C) and background strain rates (∼10^{−16}-10^{−12} s^{−1}) predicted by thermal models (Fig. 2a), suggesting that growth rates for sediment diapirs may vary significantly. While quartz and olivine are brittle at low temperatures (<280-350 °C for quartz, Stöckert *et al.* 1999; <600 °C for olivine, Boettcher *et al.* 2007), we are interested in the formation of sediment diapirs at depths greater than ∼50 km in subduction zones, where diapirs would be incorporated into arc magma systems. At these depths, temperatures are near or above the temperature of the brittle/ductile transition for quartz and olivine, and we thus only consider deformation by viscous flow. At ∼600 °C, olivine viscosities are high (10^{20}-10^{23} Pa s for strain rates of 10^{−12}-10^{−16} s^{−1}) and the viscosity ratio of quartz to olivine (i.e. *r = B*_{qtz}/*B*_{ol}) at a given strain rate is <10^{−3}. As temperature increases to ∼1000 °C, olivine viscosities decrease to ∼10^{18}-10^{21} Pa s for strain rates of 10^{−12}-10^{−16} s^{−1}, approaching the viscosity of quartz and causing *r* to increase to 0.1. This decrease in absolute viscosities with increasing temperature promotes diapir growth, shortening instability times. Further, the relative weakening of olivine with respect to quartz (i.e. larger *r*) with temperature allows quartz diapirs to more readily ascend into the overlying olivine half-space, decreasing instability times (e.g. Houseman & Molnar 1997; Conrad & Molnar 1997).

The effect of the viscous decay length on the growth of sediment diapirs depends on the ratio of *L* to the layer thickness *h*. In subduction zone thermal models, *L* is typically minimal at depths of ∼80 km (Fig. 3e) and acts to weaken the mantle with respect to the sediment layer, promoting diapir growth. This decrease in *L* is a consequence of steepening, inverted geotherms near the slab/mantle decoupling/coupling transition (e.g. Wada & Wang 2009). The abruptness of the transition determines both the shape of the mantle geotherms and the depth and magnitude of the minimum in *L* (Syracuse *et al.* 2010), which in turn influences the time it takes sediment instabilities to form.

Diapirism is driven by the buoyancy of the sediment layer, which depends on the density contrast between the sediments and the overlying mantle and the thickness of the sediment layer. We assign a density contrast between the sediments and mantle wedge based on Behn 's (2011) calculations of temperature- and pressure-dependent densities for metasediments and mantle-wedge harzburgite along a slab-top geotherm (Fig. 2b). These densities assume 2 wt% H_{2}O in both compositions, which is representative of the amount of H_{2}O retained by subducting sediments at 700 °C and 3 GPa (Hacker 2008). The resulting density contrast varies from −340 to −80 kg m^{−3} at depths <100 km, but is nearly constant at ∼−200 kg m^{−3} at depths >100 km (Fig. 3c).

The geometry and thickness of the subducting sediment layer is largely unconstrained at depth, and we therefore calculate timescales for diapirism over a range of layer thicknesses and initial perturbation sizes. The average thickness of the sediment layer on the incoming plate at trenches varies from ∼100 to ∼4000 m globally (Plank & Langmuir 1998), although much of this of sediment is scraped off the slab at the trench, leaving at most ∼1000 m that is transported to depths beyond the forearc (von Huene & Scholl 1991; Clift & Vannucchi 2004). For completeness, we consider instabilities forming in 10-2000-m-thick layers, given that sediment layers could be >1000 m locally.

Deposition over rough topography and deformation during subduction likely cause sediments to vary in thickness across the slab, with undulations that are potentially similar in amplitude and wavelength to abyssal hill topography on the subducting plate. These variations in thickness would provide the initial perturbations in the layer surface (i.e. ) that are required for diapirism. While large portions of the sediment layer may be removed at the trench during subduction, sediments are deposited in lows on the abyssal hill topography and erosion from the top of the layer would leave significant variations in sediment thickness. Folding and faulting during shallow subduction would also create variations in sediment thickness. Abyssal hills can be 100s of meters tall (e.g. Goff & Jordan 1988) and shallowly deformed slab-top sediments have fold amplitudes of up to ∼500 m (Chauhan *et al.* 2009), leading us to conservatively (instability time decreases with increasing ) consider instability growth from initial perturbations of 10 and 30 per cent of the layer thickness ().

Diapirs grow fastest at a preferred wavelength that depends on the viscosities of the layer and half-space (e.g. Whitehead & Luther 1975), and this wavelength of fastest growth may vary significantly from the initial perturbation wavelength (Schmeling 1987). In our models and scaling for the growth of sediment diapirs in subduction zones, we find that a majority of instability growth occurs at temperatures greater than ∼800 °C and that, at these temperatures, diapirs grow fastest at wavelengths near or . For the sediment thicknesses that we consider (10-2000 m), wavelengths corresponding to are within the spectra of typical, largely stochastic, abyssal hill topography (Goff & Jordan 1988), suggesting that diapirs forming in subducted, abyssal sediments may initiate at this preferred spacing. In predicting timescales and detachment depths for instabilities forming in subduction zones, we therefore focus on diapir growth at .

### Numerical model

To calculate instability times and (eqs (18) and (14)), we first calculated the growth rates and . Since no analytic growth rate solution exists for the case of non-Newtonian rheologies and an exponentially decaying half-space viscosity (i.e. finite *L*), we calculated these growth rates using a particle-in-cell, finite-difference method (Gerya 2010). The code assumes plane strain and solves equations for two-dimensional Stokes flow and continuity on a fixed (Eulerian) grid with regularly spaced, advecting (Lagrangian) tracer points (particles). Material properties (i.e. density and rheologic parameters) are carried by the particles and transferred to the grid by bilinear interpolation. The particles are advected over small time steps by first-order Runge-Kutta integration of the gridded velocities.

We setup the initial material geometry by defining a buoyant layer beneath a dense half-space (Fig. 1). The layer and half-space were distinguished by assigning particle densities and , rheologic parameters B _{L} and B _{H} and stress exponents n _{L} and n _{H}, respectively. The material boundary was defined by eq. (5). We assigned a temperature *T*_{0} to particles within the layer, and imposed a temperature profile in the half-space. Our models do not explicitly include the effects of thermal diffusion, a simplification that does not significantly alter the accuracy of the model results. The subduction zone models of Wada & Wang (2009) and Syracuse *et al.* (2010) include heating of the subducting slab by thermal diffusion, and, by assigning values for *T*_{0} and from these subduction zone models, we implicitly include this effect as a starting condition in our diapir models.

We set the width of the model domain to , ignoring possible interactions between adjacent perturbation peaks. We numerically simulated a half-space by setting the model height to be four times this width-the minimum height found to prevent interaction of the layer with the top boundary during initial instability growth in our models. For mechanical boundary conditions, we assumed that the sediment layer is rigidly coupled to deeper slab rocks (i.e. u =*w*= 0 at *y*= 0) and imposed free-slip side and top boundaries (i.e. at *x*= [0, *x*_{2}] and *y*=*y*_{2}).

In all of our models, we used a 64 × 256-node grid (resolution of in the horizontal and vertical) with 25 particles in each grid cell, resolutions that we found to produce accurate results as compared to predictions of growth rates from analytic solutions. We determined growth rates by tracking the maximum vertical velocity of the layer particles as a function of time. For a layer and half-space with homogenous, Newtonian viscosities, plotting as a function of produces a linear relationship with a slope of , verifying that the numerical instabilities grow according to eq. (7) (Fig. 4a). Similarly, plotting versus for models with non-Newtonian rheologies yields a linear relationship with slope , indicating that the code can model the superexponential growth predicted by eq. (13) (Fig. 4b).

The effective viscosities of metasediments and the mantle vary by several orders of magnitude over the range of temperatures and strain rates in subduction zone thermal models (Fig. 2a), requiring that the models be able to handle significant viscosity contrasts (i.e. *r*). A comparison of predicted by analytic solutions (eq. (33) in Conrad & Molnar 1997) to maximum values found numerically (Fig. 5) shows that we can calculate growth rates for *r* < 0.1 (i.e. the range for quartz and olivine under subduction conditions) to within ∼0.03 of the analytic solutions. In the scaling for subduction zone conditions that follows, this *O*(10^{−2}) uncertainty in non-dimensional growth rate translates to a <1 Myr uncertainty in instability time or a <25 km uncertainty in the depth for diapir detachment for typical subduction rates.

### Calculating instability times for sediment diapirs in subduction zones

To evaluate the likelihood of sediment diapir growth in subduction zones, we first calculated non-dimensional exponential () and superexponential () growth rates for quartz diapirs in an olivine half-space over a range of temperatures (*T*_{0}= 200-1200 °C), thermal gradients (*L*/*h*= 1-1000) and perturbation wavenumbers ( 0.5-5). Examples of these growth rates are shown in Fig. 6 at 600 and 1000 °C and all values are reported in Table 3.

We then used these growth rates to calculate dimensional instability times and using eqs (18) and (14), respectively, for a range of thermal conditions (*T*_{0} and ) and, for , a range of background strain rates ( and ). In the calculations of exponential growth (eq. 18), enters only through the timescale (eq. 16). For a constant value of , increasing reduces the effective viscosity of the mantle wedge with respect to the viscosity of the layer, increasing the viscosity ratio *r*=*B _{L}*/

*B*, and we thus included via the scaling for

_{H}*r*. We found that and differ by several orders of magnitude over the range of temperatures predicted along the slab top, implying that, at low temperature, instability times would be controlled by exponential growth, but, at higher temperatures, would be dominated by superexponential growth (Fig. 7). We thus approximate instability times for diapirs forming in the presence of a background strain rate by writing:

We then scaled these instability times for different layer thicknesses, initial perturbation amplitudes, and density contrasts (Fig. 8). From eq. (14), it can be shown that differences in *h*, and are related to a difference in instability time by

*n*= 1 (exponential growth) or

*n*=

*n*(superexponential growth) and the subscripts 1 and 2 refer to the original and scaled values, respectively.

_{L}Finally, to determine the point at which diapirs will detach from the downgoing slab for specific subduction zone models, we integrated the scaled instability times calculated for slab-surface conditions as a function of depth. If instability time as a function of time elapsed during subduction (i.e. ‘subduction time', *t _{s}*) is

*t*(

_{b}*t*), then the time at which diapir displacement approaches the initial layer thickness is found by solving for when the sum of fractional displacement

_{s}*n*is the total number of time (depth) intervals along the slab surface and

*i*is the index of a discrete subduction time. After reaching this doubling height, diapirs ascend vertically at a nearly constant velocity, effectively detaching from the downgoing slab and the depth of detachment can thus be inferred from

*t*.

_{s}## Sensitivity of Diapir Growth To Temperature and Strain Rate

### Effect of temperature and half-space thermal structure

Calculated dimensional instability times show the strong dependence of diapir growth on temperature (Fig. 8). At low temperatures (800 °C), diapir growth is slow, viscosities are dominated by background strain rates, and superexponetial growth does not occur until after displacement reaches the doubling height (i.e. *Y*=*h*). At temperatures greater than ∼800 °C, lower viscosities enable faster instability growth, causing local strain rates (from the diapiric flow) to exceed background strain rates and enabling superexponential growth to initiate before diapir displacement reaches the doubling height, thereby reducing instability times.

Non-dimensionalizing growth rates by shows that diapir growth is also sensitive to secondary controls such as the ratio of the sediment-to-mantle viscosity, *r*, which is a function of temperature and the length scale of viscous decay in the mantle wedge, *L* (Fig. 6). At low temperatures (∼600 °C), olivine is much more viscous than quartz (*r*= 5.2 × 10^{−4}) and the ascent of weak quartz diapirs into the strong, overlying olivine half-space is slow . At these low temperatures, diapirs grow fastest at small wavenumbers () (Fig. 6). As temperature increases, olivine weakens with respect to quartz (*r*= 4.6 × 10^{−2} at 1000 °C), enabling faster diapir growth . By 800 °C, the preferred wavenumber for diapir growth begins to converge on (Table 3).

The thermal gradient in the mantle wedge also influences the relative viscosities of the sediment layer and mantle wedge and has the strongest effect on diapir growth at lower temperatures and longer wavelengths. For example, decreasing *L* from *L*≫*h* to *L*=*h* increases by a factor of ∼3 at 600 °C and (Fig. 6); this same change in the decay length has only a minimal effect at hotter temperatures and/or short wavelengths (). In terms of dimensional instability time, this same change in *L*/*h* results in a shortening of by two orders of magnitude at 600 °C, but only one order of magnitude at 1000 °C (Fig. 8).

### Effect of background strain rate

Finite background strain rates reduce effective viscosities, shortening diapir instability times and effectively reducing the temperature and/or sediment layer thickness required to produce diapirs. The scaled instability times in Fig. 8 show that increasing the background strain rate from 10^{−16} to 10^{−12} s^{−1} can reduce the time required to form an instability by up to three orders of magnitude (∼10^{9} to 10^{6} yr) at temperatures of 600-800 °C. Similarly, for and a temperature of 900 °C, a background strain rate of 10^{−12} s^{−1} allows a 100-m-thick sediment layer to reach a doubling height in approximately the same time it takes a 500-m-thick layer to grow by the same amount with 10^{−16} s^{−1}. The relative value of the layer and mantle background strain rates also has an effect on instability times, with a factor of 100 increase in slowing instability times by up to half an order of magnitude (Fig. 8). As shown above in instability times for variable *L*/*h* and *T*_{0}, these results demonstrate the importance of the relative viscosities of the sediment layer and the mantle half-space in diapir formation.

## Timescales for Sediment Diapirs in Subduction Zones

### Instability times for subduction zone thermal models

The sensitivity of sediment diapir growth to the thermal and strain rate structure of the slab and mantle wedge is revealed by scaling instability times to the conditions found in a range of subduction zone thermal models (Fig. 9). Diapir growth is principally controlled by the viscosity-reducing effect of temperature, with instabilities initiating when slab-surface temperatures reach ∼500 °C (η_{eff}≃ 10^{18}-10^{24} Pa s) and detaching (i.e. growing to a height equal to the initial layer thickness) within ∼1–3 Myr. The depth of initiation depends on the temperature of the incoming slab, which is a function of plate age and, at depth in subduction zones, subduction rate (Molnar & England 1990; Molnar & England 1995). In arcs with young slabs subducting at moderate rates, such as Cascadia, the slab reaches ∼500 °C at ∼30 km and diapirs begin to grow. Here, instabilities forming in a 1500-m-thick sediment layer, for example, would detach at a depth of ∼90 km. By contrast, sediments on older, faster subducting slabs, such as Izu, do not reach ∼500 °C until a depth of ∼70 km. At this depth, diapirs forming in the same 1500-m-thick layer would begin to grow rapidly and detach from the slab at a depth of ∼140-170 km.

The maximum depth of full or partial mechanical decoupling between the slab and mantle wedge significantly influences the thermal structure of subduction zone models, and thus affects viscosities and the growth of sediment diapirs in our calculations of instability times. The decoupling/coupling transition marks an abrupt onset of mantle-wedge flow that brings heat from the back arc to the slab (e.g. van Keken *et al.* 2002; Wada *et al.* 2008). The shallower the depth of this transition, the sooner subducting sediments heat up, leading to lower absolute viscosities and faster diapir growth. This effect can be seen in instability times scaled for subduction zone thermal models that make different assumptions about what controls the depth of the decoupling/coupling transition (e.g. D80: prescribed 80 km depth, X25: distance from arc, T550: depth of the brittle/ductile transition, W1300: subarc mantle temperatures; see Syracuse *et al.* 2010 for full discussion; Fig. 9). The change in temperature across the transition, and therefore the effect of the transition depth on diapir growth, is most pronounced in subduction zones with colder ambient temperatures. In models of diapirs forming on a cold slab such as that subducting at Izu, increasing the depth of the decoupling/coupling transition from 80 to 120 km delays the onset of hotter slab-surface temperatures and produces a similar difference in the predicted depth of diapir detachment. By contrast, a 30 km difference in the depth of the decoupling/coupling transition in models of the hotter Cascadia arc produces only small differences in slab temperatures, and thus diapir detachment depths are similar in these models. Changes in the degree of slab/mantle coupling also have strong effects on the thermal structure of the mantle wedge, with partial coupling in the D80 Syracuse *et al.* (2010) models enabling mantle flow to reach further into the mantle-wedge nose than in the Wada & Wang (2009) models, which assume full slab/mantle decoupling at depths less than 80 km. This enhanced mantle flow produces hotter slab-top and mantle temperatures, promoting diapirism (Fig. 10).

A shortening of the viscous decay length in the mantle wedge near the depth of the decoupling/coupling transition promotes the development of instabilities, enabling diapir growth in even the coldest subduction systems. In the mantle wedge, the decoupling/coupling transition is marked by a zone of steep, inverted geotherms (i.e. short *L*) at depths of ∼50-180 km in both hot and cold subduction systems (Wada & Wang 2009; Syracuse *et al.* 2010). In colder subduction systems, olivine is stronger with respect to quartz (i.e. small *r*) and diapir growth is more sensitive to these changes in *L* (Fig. 6). At depths of ∼70 km in thermal models of Izu, *L* shortens rapidly as *r* increases with temperature (Fig. 9). This short-lived zone where both *L*/*h* and *r* are small enables rapid diapir growth while quartz and olivine viscosities are still relatively large, promoting diapirism on cold slabs.

Variability in sediment-layer geometry can produce large differences in diapir detachment depths (Fig. 10). As with other parameters influencing diapirism, instability growth is most sensitive to changes in these parameters in colder subduction systems where absolute viscosities are larger. For diapirs growing from a 30 per cent initial perturbation () in the Izu model of Syracuse *et al.* (2010), reducing the sediment layer thickness from 2000 to 1000 m deepens the depth for diapir detachment from ∼130 to ∼190 km. In the Cascadia models, the same reduction in layer thickness only produces a ∼20 km change in the depth of detachment. Because diapirs grow significantly in the exponential (background-strain-rate dominated) regime (Fig. 8), increasing the initial perturbation amplitude has a similar effect as increasing the layer thickness (eq. 19).

### Sensitivity of instability times to variability in viscosity

The growth of sediment diapirs is sensitive to uncertainties in rheologic parameters and background strain rates. Varying these parameters changes effective viscosity, and thus we assessed the sensitivity of instability times to uncertainty in viscosity by calculating timescales for different background strain rates in the mantle wedge (Fig. 11). As with other instability parameters, diapir growth is most sensitive to changes in viscosity in colder subduction systems, such as Izu, where changing from 10^{−13} to 10^{−12} s^{−1} reduces the minimum thickness required for diapir formation by up to ∼550 m. In the hot Cascadia models, changing by the same amount reduces this minimum thickness by only ∼100 m. This one order of magnitude change in background strain rate is equivalent to a one order of magnitude change in effective viscosity (Fig. 2).

We assessed the effect of the quartz-to-coesite transition on sediment rheology in our instability calculations. This phase change occurs at a pressure of ∼2.5 GPa (∼700 °C for Izu and ∼900 °C for Cascadia) (Fig. 2b) and would increase viscosity by up to two orders of magnitude at low temperatures (700 °C) (Fig 2a). To determine the maximum effect of this viscosity change on diapir growth, we calculated instability times for a dry-coesite layer underlying a wet-olivine half-space (Fig. 7). For exponential growth, instability times are similar to those for wet quartz, while for superexponential growth, instability times are up to an order of magnitude slower. Consequently, including the quartz-to-coesite transition would delay the onset of superexponential growth from ∼700 °C (for a quartz layer) to >950 °C (for a coesite layer). In both hot and cold subduction zones, however, a significant portion of diapir growth occurs in the exponential regime (i.e. at temperatures <700-950 °C) (Figs 8 and 9). Further, because vertical velocity increases either exponentially or superexponentially with time, this late-stage change in growth rate will have a minimal effect on the time-integrated instability times.

## Discussion

Growth rates calculated from numerical models predict that, for sufficiently thick sediment layers and/or large enough initial perturbations, sediment diapirs will form in a broad range of subduction zone settings (Fig. 10). For initial perturbations that are larger than ∼30 per cent of the layer thickness, diapirs form and detach from the slab within the mantle wedge (i.e. within ∼100 km of the slab depth below the arc) when sediment layers are as thin as 200 m in hot subduction systems (young slabs and/or slow subduction rates, for example, Cascadia) or ∼800 m in cold systems (old slabs and/or fast subduction rates, for example, Izu) (Fig. 11). Globally, estimates for the thickness of subducted sediment layers range from 150 to >1000 m (e.g. Clift & Vannucchi 2004), suggesting that subducted sediments commonly form diapirs. Our predicted minimum layer thickness for diapirism at the Izu arc is similar to the result of Currie *et al.* (2007), who found that, in trenches with old (>70 Ma), cold subducting slabs, diapirs forming in sediment layers >350 m will detach within the mantle wedge. However, Currie *et al.* (2007) required a much larger density contrast (<−400 kg m^{−3}) than the ∼−200 kg m^{−3} density contrast at which diapirs grow in our models.

Our calculations of instability times suggest that sediment diapir growth depends strongly on subduction zone thermal structure, which is largely influenced by the depth of the decoupling/coupling transition (e.g. Wada & Wang 2009; Syracuse *et al.* 2010). In general, the deeper the depth of this transition, the longer it takes subducting sediments to be heated by mantle flow from the back arc, resulting in absolute viscosities that remain larger to greater depths, increasing instability times (Fig. 9). At the decoupling/coupling transition, rapidly increasing temperatures and steepening mantle geotherms create a narrow zone in which the viscosity ratio between the sediments and mantle increases sharply and the length scale for viscous decay in the mantle shortens-both conditions that promote instability growth. These effects are strongest in cold subduction zones where the contrast in slab and mantle temperatures is more pronounced and diapir growth is most sensitive to changes in secondary controls (i.e. *r* and *L*).

A key result of our study is that that the ∼−200 kg m^{−3} density contrast between subducted sediments and the mantle wedge alone is sufficient to drive diapirism. This contrasts with previous studies in which sediments are entrained in cold diapirs that arise from hydration of the mantle above the subducting slab (Gerya & Yuen 2003; Gerya *et al.* 2006; Gorczyk *et al.* 2006, 2007). In such models, fluids supplied to the mantle wedge by slab dehydration at depth may serpentinize the mantle, reducing the mantle density by 100-300 kg m^{−3} in a thick (∼5-20 km) layer above the slab (e.g. Gerya & Yuen 2003). Hydration would also drive partial melting of the mantle, further reducing mantle density. These reductions in mantle density would inhibit the growth of sediment diapirs into the serpentine layer (by reducing the density contrast between the sediment layer and overlying mantle), but may result in the development of mixed plumes that rise off the top of the slab and entrain the underlying sediment layer. Numerical models have suggested that these hydrous, cold plumes initiate at depths of ∼50-200 km (Gerya & Yuen 2003), similar to the depths predicted for the detachment of sediment diapirs. However, because these plumes originate from a much thicker buoyant layer, they are typically larger than diapirs arising from the sediment layer alone, and thus have a much greater effect on the dynamics and thermal structure of the mantle wedge than the sediment diapirs that we consider here. In our models, sediment diapirs typically form at temperatures and pressures exceeding the stability field for serpentinite (e.g. Wada & Wang 2009), implying that the density reduction by serpentinization of the mantle would not strongly influence the growth of sediment diapirs.

In hot subduction zones, diapirs detach from the slab at temperatures between ∼850 and 1000 °C and at pressures between ∼2.5 and 5.5 GPa (Fig. 11). In cold subduction zones, diapirs detach from the slab at greater depths (pressures of ∼4-6 GPa), but at slab-top temperatures as low as ∼700 °C. Sediments at these temperatures and pressures are near the fluid-saturated solidus (Nichols *et al.* 1994; Schmidt *et al.* 2004), but still below 1050 °C-the temperature threshold associated with the depletion of key trace elements forming the sediment melt signature in high- and ultra-high pressure metasediments that have undergone subduction (Behn *et al.* 2011). Thus, our results support the hypothesis that the sediment-melt signature observed in many arc magmas is generated when subducting sediments detach from the slab as diapirs and rapidly rise into the mantle wedge, where they would undergo melting (Gerya & Yuen 2003). Both sediment-derived melts and sediment diapirs would rapidly (cm yr^{−1} to m yr^{−1}) traverse the mantle wedge (Kelemen *et al.* 1997 and references therein; Hall & Kincaid 2001; Gerya & Yuen 2003), satisfying U/Th isotope measurements from arc magmas that imply slab-to-surface transfer times on the order of a million years (e.g. Hawkesworth *et al.* 1997). Ascending sediments may also be re-laminated to the base of the upper plate (Hacker *et al.* 2011).

To assess how efficiently sediment diapirs transport sediments into the mantle wedge, we allowed several models to run past the instability time (Fig. 12). We found that diapirs form plume heads by ∼2.1*t _{b}* at lower temperatures (400 °C), where instability growth is exponential. At higher temperatures (1000 °C), where superexponential growth dominates, plume heads are nearly fully formed at

*t*. In both cases, plume heads have a width of . Since diapir accent rates eventually reach m yr

_{b}^{−1}, greater than the cm yr

^{−1}rate of mantle corner flow, these instabilities can efficiently transport sediments into the core of the mantle wedge near the location of detachment from the slab. We also note that, in both cases, the rigid bottom boundary condition causes a thin veneer of sediment to remain attached to the slab, which implies that diapir formation may not be 100 per cent efficient in recycling sediments into the mantle wedge.

The models considered here are 2-D with a horizontal bottom boundary, and thus our instability timescales correspond to diapirs forming in trench-parallel sheets. An important question for future work is how sediment diaper growth rates will change for the case of a 3-D, dipping slab. Theoretical and observational studies suggest that both finger-like and sheet-like instabilities form in three dimensions (Ribe 1998; Zhu *et al.* 2009) and have growth rates that are similar to those for 2-D diapirs (Kaus & Podladchikov 2001). Nonetheless, understanding the morphology of such instabilities in three dimensions will be important for determining how sediments released from the slab will eventually be incorporated into the melting region beneath the arc.

## Conclusions

Using instability growth rates calculated from numerical models of quartz diapirs forming in an olivine half-space, we predict that sediment diapirs detach from the downgoing slab and ascend into the mantle wedge in a broad range of subduction zone settings. The growth of sediment diapirs is largely controlled by layer buoyancy and the absolute viscosities of the sediments and mantle wedge. In hot subduction zones (young slabs and/or slow subduction rates), lower absolute viscosities enable rapid diapir growth and instabilities can form and detach from the slab in layers as thin as ∼200 m. Diapirs also form in arcs with cold subducting slabs, although larger sediment layer thicknesses (≥800 m) are required. These values are similar to estimates of the subducted sediment-layer thickness in a majority of subduction zones, suggesting that sediment diapirism is a common feature of arcs and may be responsible for the ‘sediment signature' in the chemistry of many arc magmas. The growth of diapirs is sensitive to the thermal structure of the mantle wedge, which, in subduction zone thermal models, is largely controlled by the depth of the transition from a mechanically decoupled slab and mantle to full slab/mantle coupling at depth. In models with a deeper decoupling/coupling transition, the advection of heat from the back arc to the slab by mantle flow is suppressed, lowering mantle wedge viscosities and slowing the growth of diapirs. This effect is most pronounced in cold subduction zones where the contrast between the cold slab and hot ambient mantle is more extreme and where cold quartz sediments are weaker with respect to mantle olivine. In both hot and cold subduction systems, we predict that sediment diapirs detach from the slab at temperatures below 1050 °C, the temperature associated with depletions in key trace elements in UHP rocks that endured subduction. This supports a hypothesis in which sediments are transported into the mantle wedge by viscous flow before they undergo significant degrees of melting.

## Acknowledgments

We thank Ikuko Wada and Ellen Syracuse for generously sharing their subduction-zone models, Daniel Lizarralde for his insights during development of the numerical model, and Taras Gerya for making his finite-difference codes freely available. We are also grateful to Peter Molnar, William Landuyt and two anonymous reviewers for their critiques of this manuscript, as well as to John Whitehead and Ikuko Wada for comments on an earlier version of this paper. This work was supported by NSF Grant EAR-0652707 and a WHOI Deep Ocean Exploration Institute Fellowship to MB.

## References

_{2}O subduction beyond arcs