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.

Figure 1

Schematic diagrams showing the initial model setup and parametrization. (a) Model domain, boundary conditions, and initial particle layout. (b) Density structure and initial layer geometry. (c) Initial viscosity structure.

Figure 1

Schematic diagrams showing the initial model setup and parametrization. (a) Model domain, boundary conditions, and initial particle layout. (b) Density structure and initial layer geometry. (c) Initial viscosity structure.

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)  

1
formula
where A is an empirical constant, graphic 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 graphic is related to strain rate graphic by the constitutive relationship  
2
formula
we define an effective viscosity graphic as  
3
formula
where graphic is the second invariant of the strain-rate tensor graphic. The rheologic parameter B describes the sensitivity of viscosity to changes in strain rate and temperature and is defined (Molnar et al. 1998):  
4
formula
For n= 1 (Newtonian), graphic, 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  

5
formula
where Y0 is the perturbation amplitude and graphic 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  

6
formula
where t is time and q is a growth rate. Here, maximum vertical velocity w is given by  
7
formula
where w0=qY0 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 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  

8
formula

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 (graphic), 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)  

9
formula
and the following non-dimensional (primed) quantities  
10
formula
The exponential growth rate from eqs (6) and (7) is then  
11
formula
where r=BL/BH is the ratio of the layer to half-space (BH) rheologic parameters.

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:  

12
formula
and  
13
formula
respectively, where graphic is a growth rate analogous to graphic. In this case, the instability time tb is defined as the time at which maximum vertical velocity accelerates to infinity (Houseman & Molnar 1997), given by  
14
formula
As the velocity accelerates to infinity, 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  

15
formula
where T0 and graphic describe a linear, mantle-wedge geotherm, graphic. 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 graphic defined by the background strain rate. For this initial growth stage, Molnar et al. (1998) defined a second timescale  

16
formula
where graphic is the second invariant of the background strain rate tensor in the layer. The initial, exponential growth rate is then:  
17
formula
As long as the background strain rate is much faster than strain rates from local, diapiric flow, instability growth remains exponential. If local strain rates never exceed the background strain rate, the final instability time is simply  
18
formula
If, however, local strain rates become sufficiently large, growth will become super exponential and displacement will follow the form of eq. (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 (graphic3 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).

Table 2

Rheologic parameters.

Table 2

Rheologic parameters.

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.

Figure 2

(a) Effective viscosity as a function of temperature for hydrous quartz (Hirth et al. 2001) and olivine (Hirth & Kohlstedt 2003) at different strain rates. For comparison, shaded regions indicate the viscosity range of mica shist (strong T and graphic dependent rheology of Shea & Kronenberg 1992) and dry coesite (Renner et al. 2001). In subduction zone thermal models, both slab-surface and mantle wedge temperatures and strain rates vary from ∼200 to ∼1000 °C and from ∼10−16 to ∼10−12 s−1, respectively. (b) Density contrast between UHP rocks (sediments) and harzburgite (mantle wedge) with 2 wt% H2O as a function of temperature and pressure (Behn et al. 2011). Lines show the pressure-temperature trajectory of the slab surface in the Cascadia (red) and Izu (black) subduction zone models of Wada & Wang (2009).

Figure 2

(a) Effective viscosity as a function of temperature for hydrous quartz (Hirth et al. 2001) and olivine (Hirth & Kohlstedt 2003) at different strain rates. For comparison, shaded regions indicate the viscosity range of mica shist (strong T and graphic dependent rheology of Shea & Kronenberg 1992) and dry coesite (Renner et al. 2001). In subduction zone thermal models, both slab-surface and mantle wedge temperatures and strain rates vary from ∼200 to ∼1000 °C and from ∼10−16 to ∼10−12 s−1, respectively. (b) Density contrast between UHP rocks (sediments) and harzburgite (mantle wedge) with 2 wt% H2O as a function of temperature and pressure (Behn et al. 2011). Lines show the pressure-temperature trajectory of the slab surface in the Cascadia (red) and Izu (black) subduction zone models of Wada & Wang (2009).

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 (graphic) 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 graphic. That is,  

19
formula
where z is depth in the subduction zone model, Vslab is the slab velocity, and zDC 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. graphic 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 graphic 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 (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 graphic0.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.

Figure 3

Key parameters controlling instability growth from models of the Cascadia (red) and Izu (blue) subduction zones (Wada & Wang 2009) plotted as a function of slab depth. (a) Slab-surface temperature. Circles indicate increments of 1 Myr in subduction time from an arbitrary reference point. (b) Strain rate in the sediment layer and mantle wedge. (c) Density constrast of UHP rocks (sediments) and harzburgite (mantle wedge) (Behn et al. 2011). (d) Effective viscosity of quartz and olivine. (e) e-folding length scale for viscous decay in the mantle wedge (solid lines) and ratio of quartz-to-olivine rheologic parameters (dashed lines).

Figure 3

Key parameters controlling instability growth from models of the Cascadia (red) and Izu (blue) subduction zones (Wada & Wang 2009) plotted as a function of slab depth. (a) Slab-surface temperature. Circles indicate increments of 1 Myr in subduction time from an arbitrary reference point. (b) Strain rate in the sediment layer and mantle wedge. (c) Density constrast of UHP rocks (sediments) and harzburgite (mantle wedge) (Behn et al. 2011). (d) Effective viscosity of quartz and olivine. (e) e-folding length scale for viscous decay in the mantle wedge (solid lines) and ratio of quartz-to-olivine rheologic parameters (dashed lines).

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 graphic 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. graphic) 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 graphic) consider instability growth from initial perturbations of 10 and 30 per cent of the layer thickness (graphic).

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 graphic or graphic. For the sediment thicknesses that we consider (10-2000 m), wavelengths corresponding to graphic 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 graphic.

Numerical model

To calculate instability times graphic and graphic (eqs (18) and (14)), we first calculated the growth rates graphic and graphic. 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 graphic and graphic, 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 graphic 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 graphic 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 graphic, 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. graphic at x= [0, x2] and y=y2).

In all of our models, we used a 64 × 256-node grid (resolution of graphic 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 graphic as a function of graphic produces a linear relationship with a slope of graphic, verifying that the numerical instabilities grow according to eq. (7) (Fig. 4a). Similarly, plotting graphic versus graphic for models with non-Newtonian rheologies yields a linear relationship with slope graphic, indicating that the code can model the superexponential growth predicted by eq. (13) (Fig. 4b).

Figure 4

Plots of non-dimensionaled, maximum vertical velocity graphic as function of non-dimensional time graphic (top panels) in numerical models with (a) Newtonian and (b) non-Newtonian rheologies and different initial perturbation wavenumbers graphic. In the Newtonian case, velocity increases exponentially and plotting ln graphic versus graphic (lower panel) produces a linear function (thin line). Similarly for the non-Newtonian case, velocity evolves according to a power law and can be linearized by plotting graphic versus graphic. In both cases, the slope of linearized velocities is proportional to the growth rate factor (i.e. graphic or graphic).

Figure 4

Plots of non-dimensionaled, maximum vertical velocity graphic as function of non-dimensional time graphic (top panels) in numerical models with (a) Newtonian and (b) non-Newtonian rheologies and different initial perturbation wavenumbers graphic. In the Newtonian case, velocity increases exponentially and plotting ln graphic versus graphic (lower panel) produces a linear function (thin line). Similarly for the non-Newtonian case, velocity evolves according to a power law and can be linearized by plotting graphic versus graphic. In both cases, the slope of linearized velocities is proportional to the growth rate factor (i.e. graphic or graphic).

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 graphic 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.

Figure 5

Comparison of non-dimensional growth rates plotted as a function of the ratio of rheologic parameters r found numerically (dots) and analytically (lines) for Newtonian (n= 1) and non-Newtonian (n= 4) rheologies and different initial perturbation wavenumbers (graphic).

Figure 5

Comparison of non-dimensional growth rates plotted as a function of the ratio of rheologic parameters r found numerically (dots) and analytically (lines) for Newtonian (n= 1) and non-Newtonian (n= 4) rheologies and different initial perturbation wavenumbers (graphic).

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 (graphic) and superexponential (graphic) 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 (graphic 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.

Figure 6

Numerically calculated growth rate factors as a function of wavenumber graphic for quartz diapirs in an olivine half-space with a homogenous layer temperature of 600 and 1000 °C for different length scales for viscous decay in the half-space (L/h). In (a), viscosities are held constant, and the diapirs grow exponentially. In (b), viscosities decay with increasing strain rate and diapir growth is super-exponential. If L/h is small, the half-space viscosity decays rapidly with distance from the layer boundary, while if Lh, viscosity in the half-space is nearly homogenous. Changes in L/h have a much stronger effect at low temperatures where quartz is weak with respect to olivine (i.e. small r) than at higher temperatures where quartz and olivine viscosities are comparable.

Figure 6

Numerically calculated growth rate factors as a function of wavenumber graphic for quartz diapirs in an olivine half-space with a homogenous layer temperature of 600 and 1000 °C for different length scales for viscous decay in the half-space (L/h). In (a), viscosities are held constant, and the diapirs grow exponentially. In (b), viscosities decay with increasing strain rate and diapir growth is super-exponential. If L/h is small, the half-space viscosity decays rapidly with distance from the layer boundary, while if Lh, viscosity in the half-space is nearly homogenous. Changes in L/h have a much stronger effect at low temperatures where quartz is weak with respect to olivine (i.e. small r) than at higher temperatures where quartz and olivine viscosities are comparable.

Table 3

Non-dimensional, exponential and superexponential growth rates for quartz diapirs in an olivine half-space for various initial perturbationwavenumbers (k′). Models were run with fixed temperatures (T) and fixed relative length scales for viscous decay in the half-space (L/h).

Table 3

Non-dimensional, exponential and superexponential growth rates for quartz diapirs in an olivine half-space for various initial perturbationwavenumbers (k′). Models were run with fixed temperatures (T) and fixed relative length scales for viscous decay in the half-space (L/h).

We then used these growth rates to calculate dimensional instability times graphic and graphic using eqs (18) and (14), respectively, for a range of thermal conditions (T0 and graphic) and, for graphic, a range of background strain rates (graphic and graphic). In the calculations of exponential growth (eq. 18), graphic enters only through the timescale graphic (eq. 16). For a constant value of graphic, increasing graphic 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 graphic via the scaling for r. We found that graphic and graphic 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:  

20
formula

Figure 7

Example determination of diapir instability time. At low temperatures, diapirs grow fastest when Newtonian viscosities are calculated for a finite background strain rate (graphic s−1) (dashed line). At higher temperatures, diapirs grow faster when viscosities are non-Newtonian and decay with increasing local strain rates (dotted line). To determine an instability time as a function of temperature, we choose the minimum of these two curves (solid line). Instability times are shown for wet-quartz (black lines) and dry-coesite (grey lines) layers underlying a wet-olivine half-space.

Figure 7

Example determination of diapir instability time. At low temperatures, diapirs grow fastest when Newtonian viscosities are calculated for a finite background strain rate (graphic s−1) (dashed line). At higher temperatures, diapirs grow faster when viscosities are non-Newtonian and decay with increasing local strain rates (dotted line). To determine an instability time as a function of temperature, we choose the minimum of these two curves (solid line). Instability times are shown for wet-quartz (black lines) and dry-coesite (grey lines) layers underlying a wet-olivine half-space.

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, graphic and graphic are related to a difference in instability time by  

21
formula
where n= 1 (exponential growth) or n=nL (superexponential growth) and the subscripts 1 and 2 refer to the original and scaled values, respectively.

Figure 8

Scaled instability times as function of temperature for initial perturbation sizes of (a) 10 per cent and (b) 30 per cent for different layer thicknesses, background strain rates and half-space viscous decay length scales. Shaded regions show the range of instability times calculated as in Fig. 7 for background strain rates of 10−16and 10−16 s−1. In most arcs, sediments are subducted in a few million years, and thus diapirs are only likely to form and detach from the slab when tb < O(106) yr.

Figure 8

Scaled instability times as function of temperature for initial perturbation sizes of (a) 10 per cent and (b) 30 per cent for different layer thicknesses, background strain rates and half-space viscous decay length scales. Shaded regions show the range of instability times calculated as in Fig. 7 for background strain rates of 10−16and 10−16 s−1. In most arcs, sediments are subducted in a few million years, and thus diapirs are only likely to form and detach from the slab when tb < O(106) yr.

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  

22
formula
is graphic, where 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 ts.

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 (graphic800 °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 graphic 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 graphic. At these low temperatures, diapirs grow fastest at small wavenumbers (graphic) (Fig. 6). As temperature increases, olivine weakens with respect to quartz (r= 4.6 × 10−2 at 1000 °C), enabling faster diapir growth graphic. By 800 °C, the preferred wavenumber for diapir growth begins to converge on graphic (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 Lh to L=h increases graphic by a factor of ∼3 at 600 °C and graphic (Fig. 6); this same change in the decay length has only a minimal effect at hotter temperatures and/or short wavelengths (graphic). In terms of dimensional instability time, this same change in L/h results in a shortening of graphic 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 graphic 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 graphic 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.

Figure 9

(a-c) Subset of the model parameters shown in Fig. 3. (d) Instability times scaled for conditions at depth (lines) and depths for diapir detachment found by integrating these instability times (circles) for Cascadia (top) and Izu (bottom). W&W refers to models from Wada & Wang (2009) in which the decoupling/coupling transition is set to 80 km and the mantle-wedge temperature is allowed to evolve to steady state. Other models are from Syracuse et al. (2010) and have the following prescribed parameters: (D80) Depth of decoupling/coupling transition is set to 80 km. (T550) The decoupling/coupling transition occurs where the slab surface reaches the brittle/ductile transition at 550 °C. (W1300) The minimum temperature in the mantle wedge beneath the arc is held fixed at 1300 °C. (X25) The decoupling/coupling transition is placed 25 km from the arc. See references for additional model details.

Figure 9

(a-c) Subset of the model parameters shown in Fig. 3. (d) Instability times scaled for conditions at depth (lines) and depths for diapir detachment found by integrating these instability times (circles) for Cascadia (top) and Izu (bottom). W&W refers to models from Wada & Wang (2009) in which the decoupling/coupling transition is set to 80 km and the mantle-wedge temperature is allowed to evolve to steady state. Other models are from Syracuse et al. (2010) and have the following prescribed parameters: (D80) Depth of decoupling/coupling transition is set to 80 km. (T550) The decoupling/coupling transition occurs where the slab surface reaches the brittle/ductile transition at 550 °C. (W1300) The minimum temperature in the mantle wedge beneath the arc is held fixed at 1300 °C. (X25) The decoupling/coupling transition is placed 25 km from the arc. See references for additional model details.

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).

Figure 10

Locations of sediment diapir detachment in the D80 Syracuse et al. (2010) and the Wada & Wang (2009) models of the Cascadia and Izu subduction zones. Scale on slab surface indicates the location of diapir detachment for different layer thicknesses and a 30 per cent initial perturbation.

Figure 10

Locations of sediment diapir detachment in the D80 Syracuse et al. (2010) and the Wada & Wang (2009) models of the Cascadia and Izu subduction zones. Scale on slab surface indicates the location of diapir detachment for different layer thicknesses and a 30 per cent initial perturbation.

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 (graphic) 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 graphic 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 graphic 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).

Figure 11

Slab-surface pressures (top) and temperatures (bottom) at the location of diapir detachment plotted as a function of sediment layer thickness in thermal models for the (a) Cascadia and (b) Izu subduction zones. Values are for a 30 per cent initial perturbation with wavenumber graphic. Colour coding is the same as in Fig. 9.

Figure 11

Slab-surface pressures (top) and temperatures (bottom) at the location of diapir detachment plotted as a function of sediment layer thickness in thermal models for the (a) Cascadia and (b) Izu subduction zones. Values are for a 30 per cent initial perturbation with wavenumber graphic. Colour coding is the same as in Fig. 9.

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.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 graphic. 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.

Figure 12

Snapshots of quartz diapirs in an olivine half-space for an initial perturbation wavenumber k′ and amplitude graphic. In (a), viscosities are held constant but are calculated for a temperature of 400 °C and a background strain rate of 10−13 s−1, and growth is exponential. In (b), temperature is set to 1000 °C, viscosities are allowed to change with increasing local strain rates, and growth is superexponential.

Figure 12

Snapshots of quartz diapirs in an olivine half-space for an initial perturbation wavenumber k′ and amplitude graphic. In (a), viscosities are held constant but are calculated for a temperature of 400 °C and a background strain rate of 10−13 s−1, and growth is exponential. In (b), temperature is set to 1000 °C, viscosities are allowed to change with increasing local strain rates, and growth is superexponential.

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

Abers
G.A.
van Keken
P.E.
Kneller
E.A.
Ferris
A.
Stachnik
J.C.
,
2006
.
The thermal structure of subduction zones constrained by seismic imaging: implications for slab dehydration and wedge flow
,
Earth planet. Sci. Lett.
 ,
241
,
387
397
.
Arcay
D.
Tric
E.
Doin
M.P.
,
2007
.
Slab surface temperature in subduction zones: influence of the interplate decoupling depth and upper plate thinning processes
,
Earth planet. Sci. Lett.
 ,
255
,
324
338
.
Armstrong
R.L.
,
1968
.
A model for the evolution of strontium and lead isotopes in a dynamic earth
,
Rev. Geophys.
 ,
6
,
175
199
.
Armstrong
R.L.
,
1991
.
The persistent myth of crustal growth
,
Aust. J. Earth Sci.
 ,
38
,
613
630
.
Behn
M.D.
Kelemen
P.B.
Hirth
G.
Hacker
B.R.
Massone
H.
,
2011
.
Diapirs as the source of the sediment signature in arc lavas
,
Nature Geosci.
 ,
4
,
641
646
.
Behr
W.M.
Platt
J.P.
,
2011
.
A naturally constrained stress profile through the middle crust in an extensional terrane
,
Earth planet. Sci. Lett.
 ,
303
,
181
192
, doi:10.1016/j.epsl.2010.11.044.
Boettcher
M.S.
Hirth
G.
Evans
B.
,
2007
.
Olivine friction at the base of oceanic seismogenic zones
,
J. geophys. Res.
 ,
112
,
B01205
, doi:10.1029/2006JB004301.
Castro
A.
Gerya
T.V.
,
2008
.
Magmatic implications of mantle wedge plumes: experimental study
,
Lithos
 ,
103
,
138
148
.
Chandrasekhar
S.
,
1961
.
Hydrodynamic and hydromagnetic stability
 ,
illustrated edn
,
Courier Dover Publications
, New York, NY,
652
pp.
Chauhan
A.P.S.
Singh
S.C.
Hananto
N.D.
Carton
H.
Klingelhoefer
F.
Dessa
J.-X.
Permana
H.
White
N.J.
Graindorge
D.
,
Sumatra OBS Scientific Team
,
2009
.
Seismic imaging of forearc backthrusts at northern Sumatra subduction zone
,
Geophys. J. Int.
 ,
179
,
1772
1780
.
Clift
P.D.
Vannucchi
P.
,
2004
.
Controls on tectonic accretion versus erosion in subduction zones: implications for the origin and recycling of the continental crust
,
Rev. Geophys.
 ,
42
,
RG2001
, doi:10.1029/2003RG000127.
Conder
J.A.
,
2005
.
A case for hot slab surface temperatures in numerical viscous flow models of subduction zones with an improved fault zone parameterization
,
Phys. Earth planet. Inter.
 ,
149
,
155
164
.
Conrad
C.P.
Molnar
P.
,
1997
.
The growth of Rayleigh-Taylor-type instabilities in the lithosphere for various rheological and density structures
,
Geophys. J. Int.
 ,
129
,
95
112
.
Currie
C.A.
Beaumont
C.
Huismans
R.S.
,
2007
.
The fate of subducted sediments: a case for backarc intrusion and underplating
,
Geology
 ,
35
,
1111
1114
, doi:10.1130/G24098A.1.
de Bremond d’Ars
J.
Jaupart
C.
Sparks
R.S.J.
,
1995
.
Distribution of volcanoes in active margins
,
J. geophys. Res.
 ,
100
,
20 421–20 432
.
Dell'Angelo
L.N.
Tullis
J.
,
1996
.
Textural and mechanical evolution with progressive strain in experimentally deformed aplite
,
Tectonophysics
 ,
256
,
57–82
, doi:10.1016/0040-1951(95)00166-2.
Elliott
T.
Plank
T.
Zindler
A.
White
W.
Bourdon
B.
,
1997
.
Element transport from slab to volcanic front at the Mariana arc
,
J. geophys. Res.
 ,
102
,
14991
15019
, doi:10.1029/97JB00788.
England
P.C.
Katz
R.F.
,
2010
.
Melting above the anhydrous solidus controls the location of volcanic arcs
,
Nature
 ,
467
,
700
703
.
Fletcher
R.C.
Hallet
B.
,
1983
.
Unstable extension of the lithosphere: a mechanical model for Basin-and-Range structure
,
J. geophys. Res.
 ,
88
,
7457
7466
.
Gerya
T.
,
2010
.
Introduction to Numerical Geodynamic Modelling
 ,
illustrated edn
,
Cambridge University Press
, Cambridge,
358
pp.
Gerya
T.V.
Yuen
D.A.
,
2003
.
Rayleigh-Taylor instabilities from hydration and melting propel cold plumes’ at subduction zones
,
Earth planet. Sci. Lett.
 ,
212
,
47
62
.
Gerya
T.V.
Connolly
J.A.D.
Yuen
D.A.
Gorczyk
W.
Capel
A.M.
,
2006
.
Seismic implications of mantle wedge plumes
,
Phys. Earth planet. Inter.
 ,
156
,
59
74
.
Goff
J.A.
Jordan
T.H.
,
1988
.
Stochastic modeling of seafloor morphology: inversion of sea beam data for second-order statistics
,
J. geophys. Res.
 ,
93
,
13589
13513
.
Gorczyk
W.
Gerya
T.V.
Connolly
J.A.D.
Yuen
D.A.
,
2007
.
Growth and mixing dynamics of mantle wedge plumes
,
Geology
 ,
35
,
587
590
, doi:10.1130/G23485A.1.
Gorczyk
W.
Gerya
T.V.
Connolly
J.A.D.
Yuen
D.A.
Rudolph
M.
,
2006
.
Large-scale rigid-body rotation in the mantle wedge and its implications for seismic tomography
,
Geochem. Geophys. Geosys.
 ,
7
,
Q05018
, doi:10.1029/2005GC001075.
Hacker
B.R.
,
2008
.
H2O subduction beyond arcs
,
Geochem. Geophys. Geosys.
 ,
9
,
Q03001
, doi:10.1029/2007GC001707.
Hacker
B.R.
Kelemen
P.B.
Behn
M.D.
,
2011
.
Differentiation of the continental crust by relamination
,
Earth planet. Sci. Lett.
 ,
307
,
501
516
, doi:10.1016/j.epsl.2011.05.024.
Hall
P.S.
Kincaid
C.
,
2001
.
Diapiric flow at subduction zones: a recipe for rapid transport
,
Science
 ,
292
,
2472
2475
, doi:10.1126/science.1060488.
Hawkesworth
C.J.
Turner
S.P.
McDermott
F.
Peate
D.W.
Van Calsteren
P.
,
1997
.
U-Th isotopes in arc magmas: implications for element transfer from the subducted crust
,
Science
 ,
276
,
551
555
, doi:10.1126/science.276.5312.551.
Hirth
G.
Kohlstedt
D.
,
2003
.
Rheology of the upper mantle and the mantle wedge: a view from the experimentalists
, in
Inside the Subduction Factory
 , Geophys. Monogr. Vol. 138, pp.
83
105
,
American Geophysical Union
, Washington, DC.
Hirth
G.
Teyssier
C.
Dunlap
J.W.
,
2001
.
An evaluation of quartzite flow laws based on comparisons between experimentally and naturally deformed rocks
,
Int. J. Earth Sci.
 ,
90
,
77
87
.
Houseman
G.A.
Molnar
P.
,
1997
.
Gravitational (Rayleigh-Taylor) instability of a layer with non-linear viscosity and convective thinning of continental lithosphere
,
Geophys. J. Int.
 ,
128
,
125
150
.
Jull
M.
Kelemen
P.B.
,
2001
.
On the conditions for lower crustal convective instability
,
J. geophys. Res.
 ,
106
,
6423
6446
.
Karig
D.E.
Kay
R.W.
Brown
G.C.
Armstrong
R.L.
,
1981
.
Fate of sediments on the descending plate at convergent margins
,
Phil. Trans. R. Soc. Lond. A.
 ,
301
,
233
251
.
Kaus
B.J.P.
Podladchikov
Y.Y.
,
2001
.
Forward and reverse modeling of the three-dimensional viscous Rayleigh-Taylor instability
,
Geophys. Res. Lett
 ,
28
,
1095
1098
.
Kelemen
P.B.
Hirth
G.
Shimizu
N.
Spiegelman
M.
Dick
H.J.
,
1997
.
A review of melt migration processes in the adiabatically upwelling mantle beneath oceanic spreading ridges
,
Phil. Trans. R. Soc. Lond. A
 ,
355
,
283
318
, doi:10.1098/rsta.1997.0010.
Kelemen
P.B.
Rilling
J.L.
Parmentier
E.M.
Mehl
L.
Hacker
B.R.
,
2003
.
Thermal structure due to solid-state flow in the mantle wedge beneath arcs
, in
Inside the Subduction Factory
 , Geophys. Monogr. Vol. 138, pp.
293
311
,
American Geophysical Union
, Washington, DC.
Kirby
S.H.
,
1983
.
Rheology of the lithosphere
,
Rev. Geophys.
 ,
21
,
1458
1487
, doi:10.1029/RG021i006p01458.
Kneller
E.A.
Van Keken
P.E.
Karato
S.
Park
J.
,
2005
.
B-type olivine fabric in the mantle wedge: insights from high-resolution non-Newtonian subduction zone models
,
Earth planet. Sci. Lett.
 ,
237
,
781
797
.
Kneller
E.A.
Van Keken
P.E.
Katayama
I.
Karato
S.
,
2007
.
Stress, strain, and B-type olivine fabric in the fore-arc mantle: sensitivity tests using high-resolution steady-state subduction zone models
,
J. geophys. Res
 ,
112
,
B04406
, doi:10.1029/2006JB004544.
Marsh
B.D.
,
1979
.
Island arc development: some observations, experiments, and speculations
,
J. Geol.
 ,
87
,
687
713
.
McKenzie
D.P.
,
1984
.
The generation and compaction of partially molten rock
,
J. Petrol.
 ,
25
,
713
765
.
Molnar
P.
England
P.
,
1990
.
Temperatures, heat flux, and frictional stress near major thrust faults
,
J. geophys. Res.
 ,
95
,
4833
485
.
Molnar
P.
England
P.
,
1995
.
Temperatures in zones of steady-state underthrusting of young oceanic lithosphere
,
Earth planet. Sci. Lett
 ,
131
,
57
70
.
Molnar
P.
Houseman
G.A.
Conrad
C.P.
,
1998
.
Rayleigh–Taylor instability and convective thinning of mechanically thickened lithosphere: effects of non-linear viscosity decreasing exponentially with depth and of horizontal shortening of the layer
,
Geophys. J. Int.
 ,
133
,
568
584
, doi:10.1111/j.1365-246X.1998.00510.x.
Nichols
G.T.
Wyllie
P.J.
Stern
C.R.
,
1994
.
Subduction zone melting of pelagic sediments constrained by melting experiments
,
Nature
 ,
371
,
785
788
.
Plank
T.
Langmuir
C.H.
,
1993
.
Tracing trace elements from sediment input to volcanic output at subduction zones
,
Nature
 ,
362
,
739
743
.
Plank
T.
Langmuir
C.H.
,
1998
.
The chemical composition of subducting sediment and its consequences for the crust and mantle
,
Chem. Geol.
 ,
145
,
325
394
.
Plank
T.
Cooper
L.B.
Manning
C.E.
,
2009
.
Emerging geothermometers for estimating slab surface temperatures
,
Nature Geosci.
 ,
2
,
611
615
.
Rayleigh
L.
,
1883
.
Investigation of the character of the equilibrium of an incompressible heavy fluid of variable density
,
Proc. Lond. Math. Soc
 ,
14
,
170
177
.
Renner
J.
Stöckhert
B.
Zerbian
A.
Röller
K.
Rummel
F.
,
2001
.
An experimental study into the rheology of synthetic polycrystalline coesite aggregates
,
J. geophys. Res.
 ,
106
,
411
419
.
Ribe
N.M.
,
1998
.
Spouting and planform selection in the Rayleigh-Taylor instability of miscible viscous fluids
,
J. Fluid Mech.
 ,
377
,
27
45
.
Schmeling
H.
,
1987
.
On the relation between initial conditions and late stages of Rayleigh-Taylor instabilities
,
Tectonophysics
 ,
133
,
65
80
.
Schmidt
M.W.
Vielzeuf
D.
Auzanneau
E.
,
2004
.
Melting and dissolution of subducting crust at high pressures: the key role of white mica
,
Earth planet. Sci. Lett.
 ,
228
,
65
84
.
Shea
W.T.
Kronenberg
A.K.
,
1992
.
Rheology and deformation mechanisms of an isotropic mica schist
,
J. geophys. Res
 ,
97
,
15201
15237
.
Spiegelman
M.
McKenzie
D.
,
1987
.
Simple 2-D models for melt extraction at mid-ocean ridges and island arcs
,
Earth planet. Sci. Lett.
 ,
83
,
137
152
.
Stöckert
B.
Brix
M.R.
Kleinschrodt
R.
Hurford
A.J.
Wirth
R.
,
1999
.
Thermochronometry and microstructures of quartz-a comparison with experimental flow laws and predictions on the temperature of the brittle–plastic transition
,
J. Struct. Geol.
 ,
21
,
351
369
.
Sun
S.S.
,
1980
.
Lead isotopic study of young volcanic rocks from mid-ocean ridges, ocean islands and island arcs
,
Phil. Trans. R. Soc. Lond. A.
 ,
297
,
409
445
.
Syracuse
E.M.
van Keken
P.E.
Abers
G.A.
,
2010
.
The global range of subduction zone thermal models
,
Phys. Earth planet. Inter
 ,
183
,
73
90
, doi:10.1016/j.pepi.2010.02.004.
Taylor
G.
,
1950
.
The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. I
,
Phil. Trans. R. Soc. Lond. A.
 ,
201
,
192
196
.
Tera
F.
Brown
L.
Morris
J.
Sacks
I.S.
Klein
J.
Middleton
R.
,
1986
.
Sediment incorporation in island-arc magmas: inferences from 10Be
,
Geochim. Cosmochim. Acta.
 ,
50
,
535
550
.
van Keken
P.E.
Kiefer
B.
Peacock
S.M.
,
2002
.
High-resolution models of subduction zones: implications for mineral dehydration reactions and the transport of water into the deep mantle
,
Geochem. Geophys. Geosyst.
 ,
3
,
1056
, doi:10.1029/2001GC000256.
Vogt
P.R.
,
1974
.
Volcano spacing, fractures, and thickness of the lithosphere
,
Earth planet. Sci. Lett.
 ,
21
,
235
252
.
von Huene
R.
Scholl
D.W.
,
1991
.
Observations at convergent margins concerning sediment subduction, subduction erosion, and the growth of continental crust
,
Rev. Geophys.
 ,
29
,
279
316
.
Wada
I.
Wang
K.
,
2009
.
Common depth of slab-mantle decoupling: reconciling diversity and uniformity of subduction zones
,
Geochem. Geophys. Geosyst.
 ,
10
,
10009
, doi:10.1029/2009GC002570.
Wada
I.
Wang
K.
Hyndman
R.D.
,
2008
.
Weakening of the subduction interface and its effects on surface heat flow, slab dehydration, and mantle wedge serpentinization
,
J. geophys. Res.
 ,
113
,
B04402
, doi:10.1029/2007JB00519.
White
W.M.
Dupré
B.
Vidal
P.
,
1985
.
Isotope and trace element geochemistry of sediments from the Barbados Ridge-Demerara Plain region, Atlantic Ocean
,
Geochim. Cosmochim. Acta.
 ,
49
,
1875
1886
.
Whitehead
J.A.
Luther
D.S.
,
1975
.
Dynamics of laboratory diapir and plume models
,
J. geophys. Res
 ,
80
,
705
717
.
Zhu
G.
Gerya
T.V.
Yuen
D.A.
Honda
S.
Yoshida
T.
Connolly
J.A.D.
,
2009
.
Three-dimensional dynamics of hydrous thermal-chemical plumes in oceanic subduction zones
,
Geochem. Geophys. Geosyst.
 ,
10
,
Q11006
, doi:10.1029/2009GC002625.
Zimmerman
M.E.
Kohlstedt
D.L.
,
2004
.
Rheological properties of partially molten lherzolite
,
J. Petrol.
 ,
45
,
275
298
.