Abstract

The plate interface undergoes two transitions between seismogenic depths and subarc depths. A brittle-ductile transition at 20–50 km depth is followed by a transition to full viscous coupling to the overlying mantle wedge at ∼80 km depth. We review evidence for both transitions, focusing on heat-flow and seismic-attenuation constraints on the deeper transition. The intervening ductile shear zone likely weakens considerably as temperature increases, such that its rheology exerts a stronger control on subduction-zone thermal structure than does frictional shear heating. We evaluate its role through analytic approximations and two-dimensional finite-element models for both idealized subduction geometries and those resembling real subduction zones. We show that a temperature-buffering process exists in the shear zone that results in temperatures being tightly controlled by the rheological strength of that shear zone’s material for a wide range of shear-heating behaviors of the shallower brittle region. Higher temperatures result in weaker shear zones and hence less heat generation, so temperatures stop increasing and shear zones stop weakening. The net result for many rheologies are temperatures limited to ≤350–420 °C along the plate interface below the cold forearc of most subduction zones until the hot coupled mantle is approached. Very young incoming plates are the exception. This rheological buffering desensitizes subduction-zone thermal structure to many parameters and may help explain the global constancy of the 80 km coupling limit. We recalculate water fluxes to the forearc wedge and deep mantle and find that shear heating has little effect on global water circulation.

INTRODUCTION

In subduction zones, shallow deformation is largely localized along the plate boundary. Earthquakes rupture at seismogenic depths to account for plate motion in some but not all subduction zones (Scholz and Campos, 2012; Wang and Bilek, 2014), and slow-slip phenomena demonstrate that this localized deformation persists to depths of 40–80 km (e.g., Lay et al., 2012; Bürgmann, 2018). Geological observations of high-pressure metamorphic rocks document narrow shear zones along the plate interface to 2.8 GPa or ∼80 km depth, exhibiting a variety of localized mechanisms of brittle and viscous deformation (Agard et al., 2018; Tulley et al., 2020). In this region, the overlying upper-plate mantle does not participate in large-scale mantle flow, forming a cold and highly viscous “nose” of stagnant mantle (e.g., Kincaid and Sacks, 1997; Stachnik et al., 2004; Abers et al., 2006). It is cold enough to allow the large-scale stability of hydrous phases such as serpentine and chlorite (Hyndman and Peacock, 2003). Seismic imaging suggests that a pervasive low-wave-speed channel exists at the plate boundary at these depths, which is a likely candidate for the narrow shear zone (e.g., Abers, 2005; Bostock, 2013).

By contrast, upper-mantle temperatures beneath arcs must exceed 1200 °C for generation of typical arc basalts (e.g., Lee et al., 2009). To sustain such temperatures, the downgoing plate must be viscously coupled to the overriding mantle at subarc depths, driving the corner flow that advects hot asthenosphere to the region. Geophysical evidence discussed below suggests that the transition from a hot subarc mantle to cold forearc is abrupt and occurs where the slab surface reaches depths near 80 km. This depth appears to be roughly constant (within 10–20 km) over a wide range of subduction zone thermal conditions (e.g., Wada and Wang, 2009). Models that produce the abrupt change in temperature for realistic mantle rheologies only do so if the transition from localized shear to distributed deformation also occurs near 80 km depth (Wada and Wang, 2009; Syracuse et al., 2010).

Thus, three deformational regimes exist along the plate interface, denoted I, II, and III on Figure 1. These are separated by two transitions: (1) the brittle-ductile transition (BDT); and (2) the transition from mechanically decoupled plates to full viscous decoupling at the coupling depth (Dc). (Fig. 1). To better understand both transitions, we investigate the temperature and stress in the shear zone for a range of shear heating in region I (brittle fault) and plausible rheologies for region II (ductile shear zone) and focus on their interplay. For temperature-dependent creep laws, we find that temperatures deeper than the BDT are buffered by the creep law, such that shear-zone stress and temperature are relatively constant over a wide range of conditions updip. The net effect of shear heating is production of a narrow range of temperatures and viscosities in the ductile shear zone, with temperatures largely controlled by the rheology of this zone. For that reason, a temperature-dependent change in the shear zone itself seems unlikely to control Dc. Instead, the increasing thickness of overlying mantle and other pressure-dependent processes facilitate the transition to coupled flow deeper than Dc. Consequences are that shear heating plays a very limited role on the temperatures of subduction zones at forearc mantle depths and that temperatures are largely controlled by shear-zone rheology.

OBSERVATIONAL CONSTRAINTS ON PLATE-INTERFACE TEMPERATURE

Heat Flow

At frictional failure, the shear stress along the fault surface (τf) is described by: 
graphic
(e.g., Scholz, 2019), where μ′ is the effective friction coefficient including pore-pressure effects (assumed to be constant here), σn is the normal stress on the fault, and zf is the depth to a point along the fault. For weak faults and shallow dips, σn ≈ ρgzf, where ρ is the mean density of overlying rock and g is gravitational acceleration. Motion along a fault generates frictional shear heating at a rate: 
graphic

where Vf is the differential rate of motion across the fault and equal to the convergence rate V to depth Dc for simple systems.

However, heat flow above thrust faults is reduced by advection of cold material downward from the Earth’s surface, counteracting the shear heating (e.g., Molnar and England, 1990). Heat-flow measurements in forearcs show some of the lowest regional values on the planet, roughly 20–30 mW/m2 in northern Honshu (Japan; Furukawa, 1993) and 25–50 mW/m2 in Cascadia (offshore northwestern North America; Lewis et al., 1988; Blackwell et al., 1990a, 1990b). Such measurements are abundant in only a few subduction zones, primarily the Japanese segments, Cascadia, Hikurangi (New Zealand), and Kermadec (southwestern Pacific Ocean) (as reviewed by Wada and Wang, 2009). Less-dense data show that consistently cold forearcs (20–40 mW/m2) are observed in many subduction zones, including Bolivia and Peru (Henry and Pollack, 1988), central Mexico (as summarized in Manea et al., 2004), and Kyushu (Japan) (as reviewed by Suenaga et al., 2018). These observations place strong limits on frictional heating along the thrust zone, limiting maximum effective friction coefficients (μ′ in Equation 1) mostly to between 0.03 and 0.04 (Gao and Wang, 2014, 2017; van Keken et al., 2018) or possibly up to 0.06 ± 0.01 (England, 2018).

These low values near the trench abruptly increase within a few tens of kilometers of the volcanic arc, generally corresponding to where the subducting plate interface reaches depths of 70–80 km (Wada and Wang, 2009; Syracuse et al., 2010). Heat flow in arcs shows considerable scatter due to complex hydrologic circulation in permeable volcanic aquifers (Saar and Manga, 2004), but generally averages closer to 100 mW/m2, or two to four times higher than in the forearc (Furukawa, 1993; van Keken et al., 2002; Wada and Wang, 2009). The sharp boundary between low and high heat flow is clear in arc-crossing studies and generally lies 5–30 km trenchward of the arc. To explain the location of this sharp transition, thermal models with temperature-dependent viscosities invoke an abrupt transition from a decoupled to a coupled plate boundary (Dc) near 70–80 km depth (e.g., Furukawa, 1993; Peacock and Wang, 1999; van Keken et al., 2002; Conder, 2005; Wada and Wang, 2009; Syracuse et al., 2010; Maunder et al., 2019). This depth is significantly greater than the maximum depth of thrust-zone earthquakes (see the Thrust-Zone Earthquake Depths section) so requires a deep region of aseismic localized deformation along the plate boundary.

Seismic Attenuation

It has been known for a considerable time that seismic signals recorded over the cold forearc are greatly enriched in high-frequency energy (described by a high quality factor Q), while signals recorded nearer to the volcanic arc show notable absence of high-frequency energy indicating a significant change in seismic attenuation (1/Q) between earthquake source and receiver (Utsu, 1966; Barazangi and Isacks, 1971; Sacks and Okada, 1974; Umino and Hasegawa, 1984). This effect may explain why the earliest recordings of earthquakes in the Sea of Japan were recorded strongly on the Pacific coast of Honshu (Hasegawa, 1918). Measured attenuation commonly varies by an order of magnitude between the cold forearc and the hot subarc mantle and tends to have a sharp near-vertical boundary (Abers et al., 2006, 2014). A number of local-earthquake tomographic studies show that the 75–80 km slab-surface depth commonly represents this boundary (Fig. 2; Table S4 in the Supplemental Material1). The transition depth is remarkably constant between arcs. The sharp boundary requires an abrupt transition in rheological properties of the plate interface (Wada et al., 2011).

Attenuation can be measured from P waves (1/QP) or S waves (1/QS). High attenuation makes high-frequency S waves difficult to observe, so QP observations are more abundant, but QS has more direct connections to temperature through shear-modulus relaxation (e.g., Faul and Jackson, 2005; Takei, 2017). In the laboratory, shear-wave attenuation becomes observable at temperatures >900 °C in olivine-dominated rocks and produces order-of-magnitude variations in QS for 200 °C changes at 1100–1300 °C (Jackson and Faul, 2010; McCarthy et al., 2011). This order-of-magnitude temperature sensitivity means that even noisy QS measurements provide a tight constraint on thermal structure. Hence, observations of QS provide a complementary constraint to heat flow that has depth resolution.

Although studies vary in the details of their methodology, the body-wave studies discussed here fit amplitude spectra to a function with the form forumla (Knopoff, 1964), where f is frequency, T is the travel time, forumla is a path-averaged attenuation parameter, and S(f) describes the earthquake source. There should be an exponential decay of amplitude with f with the decay rate scaled to 1/forumla once source and site effects are removed. Some studies include a frequency dependence to Q within this framework, approximating forumla as Q0fα, where Q0 is the Q at 1 Hz, with α of ∼0.1–0.3 supported by many experimental studies (e.g., Karato and Spetzler, 1990). Those studies then report Q0 as the effective Q at 1 Hz. Local earthquake studies typically rely on signals in the 1–10 Hz band and path lengths of 50–500 km (T ∼10–100 s for S waves). Several important examples are discussed below.

In central Alaska (USA; Fig. 2F), earthquakes within the subducting plate show variations in QS of 10×–20× between cold nose and hot subarc (Stachnik et al., 2004). The boundary between the high- and low-Q regions is tightly constrained by seismometers 10–20 km apart, which show drastic amplitude variation. The same methodology has been used to show a >20× contrast across a similar subvertical boundary in Nicaragua and Costa Rica (Figs. 2C, 2D; Rychert et al., 2008). Using different methods, the central Andes at 21°–24°S shows a very similar pattern in QP, with a 5×–10× increase beneath and behind the backarc compared to the forearc with a boundary that is intersecting slab seismicity at 85–110 km depth (Schurr et al., 2003, 2006). The Hikurangi subduction zone exhibits a sharp vertical boundary between high and low attenuation 10–20 km trenchward of the Taupo arc front intersecting the slab at 80–90 km depth in both QP (Eberhart-Phillips et al., 2008) and QS (Eberhart-Phillips et al., 2020). As in many places, QP values in the Hikurangi cold nose are intermediate between those in the hot wedge and the cold subducting plate. Similarly, the boundary between high and low QP beneath the Tonga-Tofua arc intersects the slab at ∼80 km depth (Wei and Wiens, 2018, 2020). The cold nose is well imaged in northern Honshu, although the geometry is somewhat different from other examples; the boundary dips steeply arcward a couple tens of kilometers trenchward of the arc at the upper-plate Moho and intersects the slab surface near 100 km depth (Fig. 2E; Nakajima et al., 2013; Liu et al., 2014). Beneath Hokkaido, complex three-dimensional variations show that the cold forearc corner is restricted to mantle-wedge depths <80 km (Kita et al., 2014). Beneath central Java, the high-QP forearc is well imaged above where slab seismicity is <80–100 km deep, and hints of low-QP mantle are seen directly beneath the volcanic arc, although resolution beneath the arc is poor (Bohm et al., 2013). A boundary between high-Q forearc and low-Q subarc mantle is evident in images of the Aegean subduction zone, 30–40 km arcward of Santorini, intersecting the main belt of slab seismicity at 75 ± 10 km depth (Ventouzi et al., 2018).

Some exceptions exist to this common pattern. In the northern Marianas arc (Fig. 2B; Pozgay et al., 2009), the forearc shows moderate rather than low attenuation, perhaps because of unusual hydration of the forearc (Cai et al., 2018). Unusual hydration is also evidenced by venting in shallow seamounts (Fryer et al., 1999). The remoteness of the region and limitations of temporary seafloor seismic deployments may make it difficult to resolve the small mantle wedge there. The low-attenuation (high-Q) corner appears to be absent beneath Kyushu (Saita et al., 2015; Liu and Zhao, 2015), which is attributed to be the result of a high degree of serpentinization. Forearc heat flow there is similar to that in other forearcs, indicating that the relatively low Q is indeed not a temperature effect. A complex tectonic history following Neogene subduction of the Shikoku Basin may also complicate the setting here (Nakajima, 2019). Southern Peru shows little P- or S-wave attenuation in a flat-slab segment, and there is no evidence for a cold forearc in the normal-dipping segment (Jang et al., 2019). Chen and Clayton (2009) imaged 1/QP across the trans-Mexican volcanic belt, but the unusual flat-slab geometry there makes comparison with other regions difficult.

Many of these regions do not show a clear change in seismic velocity corresponding to the sharp 1/Q boundary. Seismic velocity is strongly influenced by compositional heterogeneity, variable hydration, and the geometry of the Moho, and therefore provides a less-useful constraint on temperature in this regime.

Thrust-Zone Earthquake Depths

Several previous reviews of megathrust earthquake depths (e.g., Tichelaar and Ruff, 1993; Ekström and Engdahl, 1989) have shown that earthquakes with thrust-fault mechanisms nucleate along the plate boundary at depths ranging from 10–15 km to ∼50–55 km, depending upon the subduction zone. Such behavior is illustrated by foreshocks and aftershocks with thrust-fault mechanisms associated with recent great earthquakes. For example, those of the 2011 Mw 9.0 Tohoku earthquake (Japan) extended to 50–55 km depth (Asano et al., 2011; Nettles et al., 2011) and those of the 2010 Mw 8.8 Maule earthquake (Chile) to 40–50 km depth (Rietbrock et al., 2012; Lange et al., 2012). These depths are resolvably less than the 80 km decoupling depth, consistent with the conceptual framework of Figure 1. The physical controls on the maximum depth are probably related to the rheological transition to velocity-strengthening frictional behavior (Scholz, 1998), although the wide variety of rupture and locking behaviors suggest more-complex relationships (e.g., Wang and Tréhu, 2016; Bürgmann, 2018). These depths commonly extend beyond the upper-plate Moho depth (as compiled by Abers et al., 2017), indicating that upper-plate composition does not have a primary control on frictional properties. One recent study found maximum interplate thrust earthquake depths of 60–80 km in several subduction zone segments, based on a comparison of global catalogs and a global slab-surface compilation (England, 2018). It is difficult to reconcile these great depths with the many studies using detailed regional data or broadband body-wave modeling, all of which show that maximum depths of interplate thrust earthquakes rarely extend past 50–55 km depth. Without understanding the cause of this discrepancy, we do not consider these results further.

ESTIMATING BRITTLE-DUCTILE TRANSITION TEMPERATURES

Shear Stresses on the Fault

A feedback exists between temperature along a shear zone (Tf) and viscous dissipation in region II (Fig. 1). Because many flow laws have a strong temperature dependence, at a constant speed and shear-zone thickness colder shear zones would move at higher stress (τ), thus generating more heat. That heat would in turn weaken the shear zone, reducing the shear stress and reducing the heat produced. This negative feedback buffers the heating that occurs along the subduction plate boundary, reducing the effects of shear heating. We demonstrate this behavior in the following sections.

In region I, shear stresses follow Equation 1. Deeper, along a ductile shear zone of width w, we assume flow follows a power-law relationship for differential stress (e.g., Bürgmann and Dresen, 2008): 
graphic

where T is temperature, P is pressure, fH2O is water fugacity, R is the ideal gas constant, Vf / w is the strain rate averaged across the shear zone, and n, r, E, V*, and B are laboratory-derived flow parameters (Table S3). This model has several simplifications, such as an assumed constant shear-zone width, but provides a useful framework for thinking about the influence of rheology on temperature.

Across the brittle-ductile transition, the weaker of these two mechanisms (friction or ductile flow) limits the shear stress. Gao and Wang (2017) suggested that shear stress across the BDT can be approximated as: 
graphic

where τ is the total shear stress. We approximate τv, the boundary-parallel shear stress in the ductile regime, by Δσv/2 in Equation 3; although the principal stress orientations are uncertain, the resolved shear stress reaches this value on planes 45° from them. Shear heating is then generated following Equation 2, substituting τ from Equation 4 for τf.

Approximate Equations for Slab-Surface Temperature and Flow-Law Feedback

The temperature along the shear zone increases due to shear heating, heat conducting from below, and radiogenic heat in the overlying plate. Precise estimates of temperature require numerical methods as discussed below, but some insight can be gained from analytical approximations. Molnar and England (1990) derived approximate expressions for steady-state temperature along a planar, shallow-dip thrust zone. We followed their approach with modifications to more accurately account for non-negligible dip, as described in van Keken et al. (2019). In the case of constant dip Δ, zero crustal heat production, a constant basal heat flux Q0, and shear heating via Equations 1 and 2, the interface temperature is approximately: 
graphic
where the last term is the addition of the adiabat assuming an adiabatic temperature gradient (dT/dz)a = 0.3 °C/km. The expression: 
graphic
varies with heating mode, with forumla for basal heating and forumla for shear heating (van Keken et al., 2019); κ is thermal diffusivity. For young plates, the assumption of constant Q0 along the seismogenic zone is invalid (Molnar and England, 1995; van Keken et al., 2019), and a better approximation replaces Q0 in Equation 5 with: 
graphic

where k is thermal conductivity, A is the age of the incoming plate at the trench (Molnar and England, 1995) and Tm is the mantle potential temperature, assumed to be 1350 °C here.

The derivation of these equations requires that temperature increases along the constant-dip plate interface as a function of zfν/2, where ν is a positive integer (Carslaw and Jaeger, 1959). For the case of pure basal heating, ν = 1, and for the case of shear-heating increasing linearly with depth, ν = 3. As shown in van Keken et al. (2019), this requirement is not satisfied if radiogenic heating in the overlying crust were to be included in Equation 5 although the error that is introduced is relatively small at depths shallower than the BDT. A significant error is introduced at depths deeper than the BDT, where shear heating has a potentially arbitrary functional form and so cannot increase monotonically with zf, which is a requirement for this functional form of Tf to be valid. Any precise estimate of the slab-surface temperature deeper than the BDT requires modeling that solves the full two-dimensional equations (as we will do below), but we can predict some of the influence of the viscous rheology deeper than the BDT on slab-surface temperature by assuming that the frictional shear stress τf(zf) in Equation 5 can be replaced by the aggregate shear stress τf(zf, Tf) defined by Equation 4.

In the ductile regime, the shear stress across the shear zone τ and therefore the shear heating τVf decreases as Tf increases because of its Arrhenius dependence on temperature. As a result, the shear heating contribution to Tf in Equation 5 decreases, so Tf increases more slowly. The net effect is a reduced change in temperature once the BDT is crossed. The critical temperature of this transition depends on the flow parameters in Equation 3 and is only weakly dependent on other parameters.

There are several effects that are not included in our approach. For example, hydrothermal circulation within the shallow subducted oceanic crust near the trench may reduce slab-surface temperatures by several tens of degrees Celsius (Spinelli and Wang, 2009). We suspect that shear heating and the temperature buffering described below may minimize this effect within and deeper than the ductile shear zone. Time-dependent variations in the incoming plate can create thermal transients that could be significant at the level of tens of degrees, as discussed below. We assume a simple model to take into account sediment erosion (as described in Syracuse et al., 2010). Additional perturbations to temperature could include subduction of heterogeneities, mass transfer between upper and lower plate, and upper-plate deformation. All of these effects lead to uncertainties of a few tens of degrees, particularly near the plate boundary, but none should alter the first-order patterns described here, in large part due to the self-regulating mechanism of shear zone rheology in region II.

Numerical Modeling

While the approximate equations tend to be quite accurate along the frictional portion of a constant-dip slab, particularly when the corrections to Molnar and England (1990) that were implemented by Molnar and England (1995) and van Keken et al. (2019) are considered, the failure to correctly approximate the temperature deeper than the BDT requires the use of numerical methods even before asthenospheric corner flow is modeled. Numerical modeling is also necessary to accurately account for radiogenic heating and curved slab geometries.

We will compare the predicted slab-surface temperature from the approximate equation described in the previous section with those obtained by numerical methods similar to those presented in Syracuse et al. (2010), Abers et al. (2017), and van Keken et al. (2018, 2019). The equations and model assumptions are summarized in detail in van Keken et al. (2018) and in the Supplemental Material (footnote 1). We will describe critical or differing assumptions here.

The model geometry is that of either a straight or a curved slab descending below a rigid overriding crust that extends to the Moho depth. Between the slab and overriding crust is a dynamic wedge. We prescribe the slab motion as an internal boundary condition on the slab surface and solve the Stokes equations with dry olivine rheology in both the slab and dynamic mantle wedge. On the wedge side we assume the velocity boundary condition has velocity of 0 down to a depth of 80 km and then ramps up linearly to velocity V over an interval of 2.5 km.

For shear heating, we assume Equation 4 and apply Equation 2 along the slab surface (with Vf = V) to a depth of 80 km. Deeper than this, the differential motion across the shear zone decreases linearly (Vf → 0) until the slab is fully mechanically coupled to the wedge at 82.5 km depth, below which no more shear heating is applied. The shear-heating term renders the heat equation nonlinear. The heat equation is also nonlinearly coupled with the Stokes equation because the viscosity of the asthenospheric wedge is temperature dependent. The Stokes equation itself is nonlinear due to the stress dependence of viscosity. The governing equations are spatially discretized with a high-resolution finite-element method with a minimum temperature nodal spacing of 1 km along the slab surface from the trench to the coupling depth. We use an outer Picard iteration to resolve the nonlinearities between the Stokes equation and temperature, and an inner Newton iteration to solve the heat equation. We apply boundary conditions as in Syracuse et al. (2010) except that we use a half-space cooling model on the trench side of the domain with a potential mantle temperature of 1350 °C. For the arc-side temperature boundary condition, we used a continental geotherm that is based on a surface heat flow of 65 mW/m2 and the analytical solution of the one-dimensional heat equation with our assumed heat production until it reaches the mantle potential temperature, below which it is assumed constant. We add an adiabatic gradient of 0.3 °C/km a posteriori to the temperature used in the viscosity and differential stress in Equation 3 and in all figures shown here. Other model details and parameters are described in the Supplemental Material (footnote 1). While we include shear heating along the plate interface, we ignore viscous dissipation in the mantle wedge. Kneller et al. (2007) evaluated the effect of viscous dissipation and found that it had negligible effect in the warm mantle wedge or at the plate interface, although it could increase temperatures in the core of the cold corner by ∼100 °C.

In all cases, we assume steady-state behavior. A comparison with time-dependent models (not shown) indicates that the asthenospheric corner flow establishes quickly (within a few million years) and warms the slab surface below 80 km depth. The cold corner establishes itself somewhat slowly but tends to develop to a steady-state structure over a period of ∼60 m.y., which was assumed in van Keken et al. (2018). The longer integration time causes slight changes from those published in Syracuse et al. (2010) in which the models were time-integrated to 20 or 40 m.y. The results presented below are, assuming identical model parameters (and μ′ = 0), therefore a little colder near the slab surface in the 40–80 km depth range than the models presented in Syracuse et al. (2010), van Keken et al. (2011), and Abers et al. (2017).

A full description of the governing equations, model formulation, and solution methods is provided in the Supplemental Material (footnote 1). We used the TerraFERMA finite-element software (Wilson et al., 2017) for all results presented here. We have benchmarked this software against the subduction-zone benchmark of van Keken et al. (2008) and by detailed reproduction of the results obtained with the licensed SEPRAN software (van den Berg et al., 2015) that we have used extensively in previous work. This direct comparison as well as internal convergence tests show that our predictions of slab-surface temperature are precise to within 1 °C. The Supplemental Material (footnote 1) also provides access to the TerraFERMA modeling software to allow interested readers to reproduce the results presented below.

RESULTS

Effect of Varying Frictional Shear Heating on the Ductile Shear Zone

The finite-element solutions provide numerous insights into the consequences of feedbacks between shear-zone rheology and temperature. We compare the finite element predictions for temperature to those from the approximate equations and show that they agree well for models with straight slabs at depths shallower than the BDT. Figures 3 and 4 show the effect of varying μ′ and hence the magnitude of shear heating along the plate interface. As a reference condition (Figs. 3, 4A, 4C), we assume a wet-quartz rheology for the ductile shear zone (Hirth et al., 2001) with water fugacity estimated following the approximation of Shinevar et al. (2015) for an activity of water of 1. The temperature in the brittle regime varies with μ′ such that at 20 km depth, increasing μ′ from 0.03 to 0.15 corresponds to a 139 °C increase in temperature (136 °C for the analytical approximation) for the parameters described in the caption for Figure 3. The highest increase is at the BDT, and in the ductile regime the range of temperatures is much smaller. For a depth of 65 km, that same range of μ′ gives a temperature range of 32 °C (1.3 °C for the analytical approximation). For μ′ >0.05, the slab surface remains nearly isothermal from the BDT to the decoupling depth Dc. In this example, the limiting temperature is 390 °C and is independent of μ′. An important observation is therefore that arbitrary increases in shear heating by increasing μ′ do not result in comparably hotter plate surfaces. The ductile shear rheology ultimately limits temperature. The temperatures and stresses inferred here are also consistent with those inferred for some exhumed faults from near the BDT, for instance the exposed megathrust shear zones of Kyushu, Japan (Tulley et al., 2020).

The right-side panels of Figure 3 show the difference between models with and without shear heating. This clearly shows how shear heating primarily heats the overriding mantle wedge because diffusion into a near-stationary medium is more efficient than that into a moving medium such as the slab below it. The advection of heat causes the slab to be somewhat warmer even below the coupling point, but the effect of shear heating on the warm part of the mantle wedge is minimal. As addressed in van Keken et al. (2018), the high-μ′ models create unrealistically high heat flow, although here this effect is partly mitigated by the depth limit of the BDT due to the thermal feedback on rheology.

The effect is similar for a much weaker published serpentine flow law (Hilairet et al., 2007) (Figs. 4B, 4D). This flow law features a very low activation energy of 8.9 kJ/mol, so τ does not change with depth once it reaches a critical value, which is ∼8 MPa for this example. As a result, all calculations with any shear heating show τ = 7.8 ± 0.2 MPa at 20–70 km depth and result in essentially identical Tf(zf). Temperatures do not exceed 300 °C until >70 km depth, where the hot wedge is approached, and the predicted BDT lies at <20 km depth. This rheology may be inappropriate because complex frictional properties likely dominate serpentine over much of this depth range (Proctor and Hirth, 2016). Still, the calculation shows that the very weak rheologies sometimes postulated to be abundant along the plate boundary (e.g., Peacock and Hyndman, 1999) would result in low temperatures regardless of μ′.

Comparison to numerical solutions confirms that the approximate equations (Eqs. 5 and 6) do well in the brittle field or when τ (zf) has a simple form, as in Figure 4D. The approximations do not adequately replicate Tf when shear heating is not simply linearly dependent on depth (for reasons discussed by van Keken et al., 2019), which is the case where the yield strength is strongly temperature dependent (Figs. 4A, 4C). The analytical approximations consistently and significantly underestimate Tf and overestimate τ along the ductile shear zone.

One consequence of the narrow range in shear heating is that the temperature structure near the base of the megathrust is relatively insensitive to brittle frictional strength. For example, based on surface heat flow, England (2018) estimated that μ′ = 0.06 ± 0.01, while Gao and Wang (2014, 2017) favored lower estimates of between 0.03 and 0.04. While surface heat flow may have sensitivity to the value, Tf at depths below the BDT should be fairly insensitive to such small differences in μ′. We have not been able to precisely determine why the μ′ values of England (2018) are systematically higher than those of Gao and Wang (2014, 2017), but we note that (1) we can precisely reproduce the numerical models of Gao and Wang (2014, 2017); and (2) we have shown in van Keken et al. (2019) that the approximate equation of England (2018) used to predict the thermal structure of the forearc cannot be accurately used for curved fault geometries and that radiogenic heating cannot be accurately approximated by the semi-analytical approach. The BDT effects discussed here complicate the interpretation of heat flow over the downdip portions of the thrust zone.

Sensitivity of Tf across the BDT to Plate Input Parameters

We test a range of planar fault models with varying incoming plate age (A), convergence rate (V), plate dip (Δ), and decoupling depth (Dc) for the wet-quartz shear zone (Fig. 5). With modest shear heating, Tf remains strongly sensitive to age of the incoming plate, but variations in V, Δ, and Dc cause the predicted Tf in the ductile regime to show less variation, and in most cases by significantly less than 100 °C, even though parameters are varied over large geologically plausible ranges. In most cases, when μ′ >0, a maximum temperature of 350–400 °C is reached. In the case of zero shear stress (μ′ = 0), the fault is colder than 350 °C and an increase in V results in lower Tf due to the more efficient advection of heat by the slab. In contrast, with even modest shear stress (μ′ = 0.03), the temperatures are 350–400 °C and do not vary significantly with V, because the increase in frictional heating due to the increase in V compensates for the cooling effect of the increased downward advection. A similar variation in Tf for zero shear stress has been observed in other similar calculations (for example, Maunder et al., 2019, their figure 3). Variations in Dc have little effect on the deeper Tf because the temperature immediately shallower than Dc remains roughly constant (Fig. 5I).

The age of the incoming plate can significantly influence temperature in the ductile regime. The only models that show T > 400 °C at 30–65 km depths are those with very young (10 Ma) incoming plates. These plates rapidly heat up because they are so thin. Still, even in these cases, a small amount of shear heating greatly reduces the range of temperatures seen. The temperature range between 10 and 100 Ma plates for μ′ = 0.03 is only about one-third the range for μ′ = 0.

Sensitivity of Tf to Shear-Zone Rheology

To understand the effects of variable composition of the shear zone, we test several alternatives (Fig. 6). The wet-quartz rheology (Hirth et al., 2001) provides a reasonable approximation for many high-silica rocks because for a significant majority of subduction zones, clastic terrigenous sediments or silicic pelagic muds likely dominate the plate-boundary input sedimentary section (Hacker, 2008) and may modulate the behavior of gabbroic compositions (Wang and Tréhu, 2016). Some studies (e.g., Gao and Wang, 2014, 2017) parameterized the ductile shear zone with a wet-granite rheology (Hansen and Carter, 1983), which behaves similarly to wet quartz (Fig. 6, compare WESTERLY to WETQZ). At a wide range of depths, white micas can be stable in metasediments and could control strength (Hacker, 2008), for which limited muscovite experiments (MUSC in Fig. 6; Mariani et al., 2006) provide some insight. At oceanic subduction zones, serpentine might be abundant near the plate interface (SERP in Fig. 6; Hilairet et al., 2007). Finally, nominally anhydrous mineralogies may be relevant if the plate interface is stripped of all metasediments and hydrous material; we show calculations for a wet-olivine dislocation-creep flow law (WETOLV in Fig. 6; Hirth and Kohlstedt, 2003). Details of the assumed rheological laws are in the Supplemental Material (footnote 1).

Several of these flow laws (granite, muscovite) have similar behavior to that of wet quartz. The serpentine and biotite laws have much lower activation energy than the quartz law and much higher stress exponent, resulting in a weak, near-constant stress along the ductile shear zone (Fig. 6). The net effect is that Tf is roughly 100–150 °C lower than for wet quartz, irrespective of μ′. The wet-olivine rheology and a similar wet-anorthite rheology (not shown; Rybacki et al., 2006) are significantly stronger. With those rheologies, the plate interface remains brittle everywhere shallower than Dc. The resulting stresses are quite high, and given such high strengths, it is not clear why strain would be localized in these materials at depth. For the other flow laws we tested, at depths >60 km, stresses are generally <30 MPa and temperatures are 300–420 °C, regardless of μ′.

Realistic Geometry—Alaska Peninsula Example

Real subduction-zone faults are not planar, and realistic increases in dip with zf produce significant effects on thermal structure (van Keken et al., 2019). To explore these effects, we repeat sensitivity tests to μ′ for a geometry and other parameters resembling the Alaska Peninsula subduction zone. Figure 7 shows the resulting thermal models, and Figure 8 shows the temperature in selected layers, along with τ and surface heat flux for a range of μ′ between 0.00 and 0.15. These tests resemble those for a planar case, although temperatures are 20–40 °C higher because a small amount of upper-plate radiogenic heat production is included (see the Supplemental Material [footnote 1]). As with simpler geometries, shear heating beyond μ′ = 0.03–0.05 makes little additional change to the temperature structure deeper than 40 km. Temperatures at the Moho of the downgoing plate (Fig. 8B) are less affected by shear heating and never exceed 400 °C at depths <75 km. Surface heat flow (Fig. 8D) shows somewhat more sensitivity to shear heating, particularly above the shallow brittle portion where each increase of μ′ by 0.05 results in a 10 mW/m2 increase in heat flow. Surface heat flow should be much more sensitive to shear heating than Tf is, particularly at depths where the slab surface is shallower than the BDT.

DISCUSSION

The Decoupling Point

One consequence of the τ-Tf buffering is that the temperature of the plate interface is relatively constant as the decoupling depth Dc is approached. Within ∼10–15 km of that depth, Tf increases more rapidly than at shallower depths due to lateral conduction from the hot wedge where isotherms are nearly perpendicular to the slab surface (e.g., Figs. 3, 7). Overall, the temperature of the interface at Dc is relatively unaffected by the earlier thermal history, but is affected by interface rheology and, to a lesser extent, incoming plate age. This may contribute to the relative invariance of Dc between arcs (see Observational Constraints on Plate-Interface Temperature section), in that the temperature at Dc is relatively insensitive to that depth, so temperature-dependent processes are probably unimportant. In other words, downdip variations in Tf are small in the ductile regime, so it is unlikely that a single critical isotherm is crossed in a narrow depth range for most subduction zones.

A pressure-dependent process seems more likely, and several have been proposed, although none have found general acceptance. For example, Wada and Wang (2009) suggested that the thickness of the overlying wedge governs upper-plate strength, because the vertically averaged shear strains (and hence basal stresses) should decrease as the thickness of the wedge increases. A nearly constant-strength shear zone, as demonstrated here and generated by near-constant Tf, would then become stronger than the overlying mantle wedge once a critical depth is passed, promoting flow there.

Other explanations rely upon the temperature and/or pressure dependence of mineral stability. The breakdown of hydrous phases such as antigorite or chlorite would strengthen the plate interface and increase coupling to the overriding mantle wedge. The dehydration reactions are almost isothermal in the 2–3 GPa range of relevance, with breakdown of antigorite at ∼650 °C (e.g., Reynard, 2013) and that of chlorite at ∼850 °C (Till et al., 2012). It is unlikely, however, that this mechanism can explain the near-uniform depth of Dc because progressive cooling of the cold corner would over time strengthen the plate interface to greater depths, and with that, increase Dc over time (Conder, 2005; Syracuse et al., 2010). Amphiboles break down near typical Dc depths over a wide range of temperatures (Schmidt and Poli, 1998), but it is not clear how that would significantly alter the shear-zone rheology. The quartz-coesite transition is also expected at these depths (Bohlen and Boettcher, 1982). This reaction can perhaps strengthen the plate interface, but it seems unlikely that a metasedimentary package in contact with presumably serpentinized overlying mantle would control rheology of the primary shear zone. Overall, while several candidate mineral transitions exist, all show inconsistences. While the causes of globally near-constant Dc remain unclear, the constancy of temperature between the BDT and Dc provides important context for this process.

Relevance to Earthquake Depths

Depths of thrust-zone earthquakes show that brittle behavior is common at 20–50 km depth along the megathrust (Tichelaar and Ruff, 1993). The observed maximum earthquake depth range corresponds closely with the maximum stress and maximum extent of brittle behavior for several flow laws with moderate strength and high activation energy (quartz, muscovite, and granite in Fig. 6). By contrast, the weak serpentine or biotite laws cross the brittle-ductile transition at depths <20 km. If these materials control the plate rheology over its entire contact, it is unclear how megathrust earthquakes would nucleate at >20 km depth as is commonly observed. On the other hand, the very strong olivine (and anorthite) flow laws show no ductile flow at any depths above Dc. If such rheologies controlled ductile flow, there should be no gap between the downdip end of seismicity and Dc. The distribution of earthquakes in many subduction zones is best represented by the moderate-strength flow laws tested here.

The downdip limit of earthquake rupture is commonly considered to be a transition from velocity-weakening to velocity-strengthening frictional behavior (Scholz, 1998), although physical processes affecting frictional stability, geodetic locking, and rupture of either great, moderate, or slow earthquakes seem to result in different inferred depths of the seismogenic zone (Wang and Tréhu, 2016). Regardless, with increased depth below that of earthquake rupture, some weakening along the plate interface seems likely. Over the 10 m.y. order-of-magnitude time scale over which thermal structure forms, the simple smoothed BDT considered here should provide a reasonable approximation to stresses and behavior, somewhere within the range of flow laws tested. There is likely to be some variation between regions as well; several oceanic arcs in the western Pacific (Izu-Bonin, Marianas, south Tonga, Kermadec) show little interplate thrust seismicity (Scholz and Campos, 2012) and scant siliciclastic sediment input. These thrust zones may be controlled by weak or frictionally stable serpentine rheology, and hence never show velocity-weakening behavior.

Off the western Alaska Peninsula, from where the geometry in Figures 6 and 7 is based, the deepest thrust-zone earthquakes reach depths of 45 km (Abers et al., 1995), consistent with the deeper BDT depth estimates. If ductile behavior resembling the quartz flow law is appropriate, then μ′ ≤0.03 results in brittle behavior to depths of 40 km or more. Greater shear heating (higher μ′) results in a shallower BDT so is difficult to reconcile with earthquake nucleation at >40 km depth. The maximum depth for thrust-zone earthquakes somewhat exceeds the upper-plate Moho depth, indicating that earthquakes can nucleate where the upper plate is peridotite that is potentially hydrated. Upper-plate crustal thickness is 35–40 km near the arc (Fliedner and Klemperer, 2000; Janiszewski et al., 2013), but the upper-plate crust appears to thin southward, reaching <25 km where the upper-plate Moho contacts the plate interface (Lizarralde et al., 2002). Serpentinization of the upper plate, if it occurs, does not seem to have an obvious effect on seismogenesis, supporting the inference that other slightly stronger rheologies are relevant.

Effect of Shear Heating on H2O Fluxes

Shear heating might affect dehydration of the subducting plate by changing the depths over which hydrous phases are stable. To test this possibility, we regenerate the calculations of van Keken et al. (2011) for H2O fluxes to great depth and those of Abers et al. (2017) for H2O fluxes into the wedge, for a range of shear heating and a subset of subduction zones representing the spectrum of behavior. We use the same incoming plate stratigraphy, compositions, and water contents as van Keken et al. (2011), but with updated H2O estimates generated by updated phase-equilibria calculations and a newer thermodynamic database (see the Supplemental Material [footnote 1]), leading to water contents shown in Figure 9 and dehydration paths shown in Figure 10. These calculations show slightly smoother dehydration paths than previously determined for some compositions (compare thick solid with dotted line in Fig. 10), and some phase boundaries move slightly. Differences between the μ′ = 0 case and previous estimates largely reflect changes in phase-equilibria estimates, while μ′ = 0.03 and μ′ = 0.05 cases illustrate the effect of moderate and high shear heating; we assume the wet-quartz rheology below the BDT. Results are summarized in Table 1 for four subduction zones: the Alaska Peninsula represents average ages and geometries; Nicaragua represents steep, fast subduction of a moderate-age plate; northern Honshu typifies cold subduction; and the Washington Cascades (northwestern USA) typify hot subduction zones. The first three have geometry and other parameters as described by Syracuse et al. (2010) and van Keken et al. (2011), while the Washington Cascades model is tied to more recent geophysical constraints (van Keken et al., 2018). Water flux to the cold mantle wedge is calculated as the amount of water released from the slab where its surface lies between the upper-plate Moho and a point just shallower than Dc = 80 km depth, chosen to be three model grid nodes shallower than Dc (∼3 km updip) to avoid the strong influence of sharp temperature gradients near Dc.

The effect of shear heating on dehydration is modest, although for some pressure-temperature trajectories that are very close to dehydration reactions, small temperature changes can have noticeable consequences (compare paths for μ′ = 0 with μ′ = 0.05 in Figs. 9, 10). Quantitative H2O loss at most subduction zones is slight to nonexistent at depths less than Dc (80 km here), with most dehydration immediately deeper where the slab encounters the hot mantle. Only the Cascades (and other hot subduction zones) show significant dehydration at shallower depths (Abers et al., 2017). The addition of shear heating (μ′ = 0.05) increases H2O loss beneath the forearc mantle wedge by 10% for Cascadia and <2.5% for other arcs, relative to the total H2O entering the subduction zone. Even with this additional input, it would take several hundred million years to fully hydrate most mantle wedges, much longer than the typical time over which subduction system geometries are stable so that wedge hydration is an inefficient process (Abers et al., 2017).

In all cases, the total H2O transported by the slab at 150 km depth is within 2% of that calculated previously for μ′ = 0 (van Keken et al., 2011). At 230 km depth, the new calculations for Alaska and Nicaragua predict several tens of percent change in the net H2O retained in the uppermost subducted mantle lithosphere (Table 1) caused by variations in H2O loss at >150 km depth (Fig. 10). The variations result from the Tf(zf) trajectory lying very close to a subparallel phase-A dehydration boundary, so that increase in Tf of a few tens of degrees Celsius with shear heating causes the boundary to be traversed; also, small changes in the calibration for phase-A dehydration produce changes from the results presented in van Keken et al. (2011). Because the depth and degree of hydration of incoming mantle are poorly known in most subduction zones, this uncertainty has only marginal significance. Ultimately, addition of even large shear heating (μ′ = 0.05) generally increases water loss from the slab by a small fraction of total input by 150 km depth and does not alter global flux calculations nearly as much as other uncertainties.

CONCLUSIONS

It has become increasingly evident that the plate interface in subduction zones passes through at least two distinct rheological transitions between the seismogenic zone and subarc depths. At depths of 20–50 km, a BDT separates a shallow region of earthquake nucleation from a deeper ductile shear zone. Near 80 km, the localized ductile shear zone of plate decoupling transitions to a viscously coupled mantle flowing wedge. The depths of earthquakes provide clear indication of the shallower transition, while heat flow and seismic attenuation show that updip of the deeper decoupling transition, the mantle wedge is relatively cold.

The temperature dependence of interface rheology reduces the range of Tf on the downdip end of the interface, and this places strong limits on the effect that shear heating has on Tf. This effect depends upon the rheology of the shear zone, but for a wide range of rheologies, temperatures are limited to <420 °C in the 40–70 km depth range. In the presence of even modest (μ′ = 0.03) shear heating, many common subduction parameters such as V, Δ, and Dc have little effect on interface temperature, either in region II (the ductile shear zone) or region III (the fully coupled slab). For larger values of μ′, the effect of varying subduction zone parameters can be larger. Only a few rheologies predict earthquake depth distributions resembling those actually observed.

We conclude that shear heating is important in determining surface heat flow as well as the location of the BDT but has only modest effect on temperature structure below the BDT. The exact value of μ′ is less critical than the ductile-shear-zone rheology in determining deeper temperatures. Notably, fluxes of H2O delivered by the slab to either the cold mantle wedge or the deep mantle are negligibly changed by shear heating; factors such as incoming plate age are far more important.

Several issues remain outstanding. First, we note that most careful global surveys of earthquake depths along the thrust zone are several decades old, despite a rapid growth in the quantity and quality of the relevant seismic data. Second, the ultimate controls on Dc remain poorly understood. The models presented here (and in virtually all other studies) do not treat the shear zone and flowing wedge in a consistent, fully dynamic manner that would allow coupling to be understood. This study does show that a simple increase in temperature along the shear zone may be an unlikely cause of coupling, because the slab-top temperature changes very little between the several models tested here. We suspect that the mechanics of driving wedge flow are a more critical factor in determining decoupling depth.

ACKNOWLEDGMENTS

We thank Kelin Wang and Saskia Goes for detailed and constructive reviews. We thank Yoko Teruahara (Earthquake Research Institute, University of Tokyo, Japan) and Miki Nakajima (University of Rochester, New York, USA) for help in translating and providing context to Hasegawa (1918). This project was supported by a Carnegie Tuve Fellowship to GAA and by the U.S. National Science Foundation through GeoPRISMS grant EAR-1850634 to PEvK. We thank Nathan Sime for help in prototyping our models in the early phases of this project.

1Supplemental Material. Details of the numerical methods, equations, code and model availability, full list of parameters used in the models, and a summary of previous attenuation studies. Please visit https://doi.org/10.1130/GEOS.S.12931637 to access the supplemental material, and contact editing@geosociety.org with any questions.
Science Editor: Shanaka de Silva
Guest Associate Editor: Gray E. Bebout
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.