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.
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.
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)Molnar et al. 1998):
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
For Newtonian materials, Chandrasekhar 1961) showed that maximum instability displacement Y increases exponentially from this initial deflection according to
To quantify the timescale for instability formation, we used the ‘instability time' (e.g. Houseman & Molnar 1997), which we refer to by the symbol tb. For the Newtonian (i.e. exponential growth) case, we define tb as the time at which diapir displacement reaches a height equal to the initial layer thickness 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 (BL and nL, respectively) as characteristic scales, we can define a characteristic timescale (Houseman & Molnar 1997)(6) and (7) is then
For cases with non-Newtonian rheologies in which BL and BH 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:Houseman & Molnar 1997), 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 byConrad & 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 timescaleeq. (12).
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% SiO2 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,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 (1020-1023 Pa s for strain rates of 10−12-10−16 s−1) and the viscosity ratio of quartz to olivine (i.e. r = Bqtz/Bol) at a given strain rate is <10−3. As temperature increases to ∼1000 °C, olivine viscosities decrease to ∼1018-1021 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% H2O in both compositions, which is representative of the amount of H2O 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 .
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 T0 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 T0 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, x2] and y=y2).
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 (T0= 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 (T0 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=BL/BH, and we thus included via the scaling for 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
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', ts) is tb(ts), then the time at which diapir displacement approaches the initial layer thickness is found by solving for when the sum of fractional displacement
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 (∼109 to 106 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 T0, 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≃ 1018-1024 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.
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.1tb 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 tb. In both cases, plume heads have a width of . Since diapir accent rates eventually reach m yr−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.
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.
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.