Advances in the thermal and petrologic modeling of subduction zones

In 1994, Gray Bebout, David Scholl, Steve Kirby, and John Platt organized a highly stimulating conference on subduction zones on Santa Catalina Island off the coast of Los Angeles, California (USA), resulting in the publication of the well-cited American Geophysical Union volume Subduction: Top to Bottom (Bebout et al., 1996). In this contribution, I review some of the major advances in understanding subduction zones made by the thermal-petrologic modeling community over the past 20–25 years. In particular, I focus on hybrid kinematic-dynamic thermal models where the flow of the subducting slab is prescribed kinematically and flow in the mantle wedge is driven by viscous coupling with the subducting slab (e.g., van Keken et al., 2002). In these models, the slab geometry, commonly defined by Wadati-Benioff seismicity, and convergence rate are treated as model inputs. Fully dynamical models, which calculate both the slab and mantle-wedge flow fields, include terms for thermal and compositional buoyancy (e.g., Gerya, 2011; Arcay, 2017).

In this paper, I adopt a petrologic point of view, choosing to represent and compare thermal models of the subducting plate through the use of metamorphic pressure-temperature (P-T) paths. These P-T paths, derived from analytical expressions and numerical models, represent the conditions “experienced” by subducting sediments, oceanic crust, and mantle. In steady state, the P-T conditions along the subduction interface coincide with the P-T path followed by a rock at the top of the subducting plate. Most kinematic-dynamic thermal models focus on steady-state solutions where the convergence rate and slab geometry remain constant, but these parameters, and therefore the thermal structure, vary with time (e.g., during subduction initiation, subduction zones are considerably warmer than at steady state). The petrologic evolution of a subducted rock depends primarily on the bulk composition and volatile (H₂O) content of the rock prior to subduction, and on the P-T path followed by the rock during subduction. For example, metabasalts of the upper oceanic crust transform to dense eclogite during subduction accompanied by significant H₂O release. The P-T paths followed by the subducting oceanic crust determine where the transformation occurs, and seismological observations can be used to test P-T paths generated by thermal models.

The physics and mathematics of heat transfer are well known (e.g., Carslaw and Jaeger, 1959) and were first applied to subduction zones 50 years ago (e.g., Oxburgh and Turcotte, 1970; Minear and Toksöz, 1970; Hasebe et al., 1970). The thermal structure of subduction zones can be approximated by using either analytical or numerical techniques to solve the equations governing heat transfer. In subduction zones, heat transfer involves conduction, advection (heat transported by the subducting plate and flowing mantle wedge), and a variety of heat sources (e.g., shear heating) and sinks (e.g., partial melting). Both analytical and numerical thermal models are simplified representations of the real world that strive to capture the essential physical processes occurring in subduction zones and can be tested against observations. Surface heat-flux measurements provide an important constraint on thermal models of subduction zones, in particular the low surface heat flux measured in the forearc and the high surface heat flux in the arc.

Analytical approximations (e.g., England, 2018) provide useful insight into the thermal structure of subduction zones and permit the systematic investigation of parameter space and uncertainties. The importance of different variables in determining subduction-zone temperatures can be readily visualized by inspecting the analytical expressions (equations). Analytical expressions are commonly used to estimate subduction-zone conditions at depths <50 km where temperatures are not affected by flow in the mantle wedge.

The three terms in the numerator of Equation 1 represent upward heat flow from the subducting lithosphere, shear (frictional or viscous) heating along the interface, and radiogenic heating in the upper plate, respectively. The denominator, S, reflects the cooling effect of subduction (advection) on interface temperatures, with faster convergence rates resulting in greater advective cooling and cooler interface temperatures.

Additional equations presented by Molnar and England (1990) and England (2018) allow the calculation of steady-state temperatures within the subducting plate as well as the initial temperatures experienced during the early stage of subduction. van Keken et al. (2019) presented slightly modified versions of Molnar and England’s (1990) equations that increase the accuracy of the analytical approximations. England and Wilkins (2004) and England and Katz (2010) presented analytical approximations for greater depths where the convecting mantle wedge strongly influences temperatures in the subducting slab.

Numerical approximations (e.g., Syracuse et al., 2010) allow calculation of the induced mantle-wedge flow field for different mantle rheologies and are commonly used to estimate temperatures in subduction zones at depths >50 km. Numerical models can be tailored to the geometry and structure of specific subduction zones, allowing the comparison of predicted temperatures with site-specific observations. In numerical models, appropriate forms of the equations describing the conservation of mass, momentum, and energy are solved numerically, most commonly using two-dimensional finite-element or finite-volume techniques. In hybrid (coupled) kinematic-dynamic numerical models, the overriding lithosphere is treated as fixed (rigid), the geometry of the subduction interface is defined by Wadati-Benioff seismicity, the motion of the subducting plate below the interface is prescribed kinematically, and viscous corner flow in the mantle wedge induced by the subducting plate is computed numerically. As discussed below, simulating the induced flow in the mantle wedge depends critically on the viscosity of the mantle wedge and the downdip transition along the subduction interface from uncoupled to fully coupled flow. van Keken et al. (2008) described a set of community benchmarks created to assess the different approaches used to construct numerical models of the dynamic and thermal structure of subduction zones. In order to accurately model the thermal structure of a subduction zone using a hybrid kinematic-dynamic approach, van Keken et al. (2008) found that the velocity discontinuity along the seismogenic zone must be carefully implemented, thermal boundary layers require high resolution, and the wedge-flow boundary condition must be modified to avoid a pressure singularity at the tip of the mantle wedge.

The accuracy of predicted temperatures at depth in subduction zones may be viewed as depending on two types of uncertainties: (A) uncertainties associated with the parameters in the equations above and (B) the extent to which the models capture the essential physical processes that control heat transfer in subduction zones.

From my perspective, the two most important type A uncertainties are both associated with the nature of the subduction plate interface—the magnitude of shear heating along the shallow subduction interface, and the depth at which the subducting slab becomes fully coupled to the overlying mantle wedge. Radiogenic heat production in the overriding continental crust also represents a significant uncertainty. In the following discussion section, I attempt to illustrate these uncertainties and how geophysical, geochemical, and petrological observations have been used to reduce these uncertainties.

Type B uncertainties are more difficult to assess. A number of researchers have constructed complex dynamical models of subduction zones that exhibit phenomena not captured by the relatively simple analytical expressions and hybrid kinematic-dynamic numerical models (e.g., Gerya et al., 2002, 2006; Honda and Saito, 2003; Gerya, 2011, 2019; Arcay, 2017). Dynamical modeling by Gerya et al. (2002) demonstrated that large-scale return flow of subducted materials may occur in a broad subduction channel composed of low-viscosity materials (subducted sediment and hydrated oceanic crust, and serpentinized mantle derived from the hanging wall). Well-documented geologic examples of subduction-zone return flow include the Central Belt mélange of the Franciscan Complex in northern California (e.g., Cloos, 1982) and the Pelona Schist of southern California (e.g., Platt et al., 2018). In my opinion, it remains to be determined whether large-scale return flow is the norm in subduction zones or whether such flow occurs only under unusual conditions such as when large amounts of sediment are subducted. The incorporation of compositional buoyancy into numerical models has revealed the potential importance of material and heat flow by diapirism including cold, hydrated-mantle plumes (Gerya and Yuen, 2003), sedimentary diapirs (Behn et al., 2011), and mélange diapirs (Marschall and Schumacher, 2012). Developing accurate thermal models for arc magmatism requires incorporating the transport of heat by magmatic flow, as well as the enthalpies of partial melting and crystallization. The thermal effect of single-pass aqueous fluid flow is likely of secondary importance, but hydrothermal circulation within subducted oceanic crust can alter the shallow thermal structure of subduction zones (e.g., Spinelli and Wang, 2008), and metamorphic dehydration and hydration reactions can change the compositional buoyancy and resulting solid flow field. A complete evaluation of fully dynamical models that include thermal and compositional buoyancy is beyond the scope of this paper, and I refer readers to Gerya (2011; 2019) and Arcay (2017).

The primary factors governing the thermal structure of subduction zones as determined from hybrid kinematic-dynamic models are depicted schematically in Figure 1.

Seafloor-depth and heat-flux measurements clearly demonstrate that oceanic lithosphere cools with age, and both data sets are fit well by cooling plate models (e.g., Parsons and Sclater, 1977; Stein and Stein, 1992; Hillier and Watts, 2005; Crosby et al., 2006). Subduction zones characterized by young incoming lithosphere are warmer than those with old incoming lithosphere, as reflected by the Q₀z_f term in Equation 1. Near spreading ridges, considerable scatter in heat-flux measurements documents the advection of heat by vigorous hydrothermal circulation within the oceanic crust driven by high thermal gradients (e.g., Davis and Lister, 1977; Stein and Stein, 1994). This hydrothermal circulation cools the oceanic crust and can reduce Q₀ out of very young oceanic lithosphere by up to 50% (Stein and Stein, 1994). England’s (2018) analytical expressions explicitly include the cooling effect of hydrothermal circulation for subduction zones where the incoming oceanic lithosphere is <67 Ma.

As oceanic lithosphere moves away from spreading ridges, hydrothermal circulation continues to redistribute heat within the upper oceanic crust (e.g., Davis et al., 1997; Fisher and Von Herzen, 2005). Additional hydrothermal circulation in oceanic lithosphere may occur along plate-bending normal faults in the trench–outer rise region (Grevemeyer et al., 2005) and along propagator wakes (Nedimović et al., 2009).

The upper-plate geotherm chosen for the arc-side boundary condition (Fig. 1) affects the thermal structure predicted for the convecting mantle wedge. Depending on whether one is modeling an oceanic or a continental subduction zone, the geotherm may include radiogenic heat production in the crust. For continental subduction zones, a number of hybrid kinematic-dynamic models (e.g., Peacock and Wang, 1999; Syracuse et al., 2010) define an upper-plate geotherm consistent with the average continental surface heat flux of 65 mW/m² (Pollack et al., 1993). Currie and Hyndman (2006) presented compelling arguments that many back-arc regions are significantly hotter than this, and proposed that vigorous thermal convection in the back-arc region, promoted by low viscosities, would result in significantly higher back-arc temperatures and a thinner lithosphere.

Global heat-flux measurements demonstrate that subduction-zone forearcs are cold relative to other tectonic environments, reflecting the subduction of cool lithosphere (e.g., Uyeda, 1977). In principle, surface heat-flux measurements can be used to constrain the temperature and rate of shear heating on the underlying subduction interface. In practice, surface heat-flux measurements vary considerably over relatively small distances, reflecting local sediment deposition and erosion and advection of heat by fluid flow (e.g., England, 2018), such that heat flux measurements do not provide firm constraints on subduction interface temperatures and the rate of shear heating. Nevertheless, the low heat flux observed in subduction-zone forearcs places an upper bound on the rate of shear heating along the interface, constraining μ′ in most cases to be <∼0.1.

Shear heating along the subduction interface refers to both frictional heating and viscous heating averaged over time scales spanning multiple earthquake cycles. Analytical approximations are well suited for calculating temperatures in the shallow forearc where the influence of mantle-wedge convection is minimal. The P-T paths depicted in Figure 2 for the northern Honshu subduction zone (Japan) illustrate the effect of shear heating on predicted slab-surface interface temperatures and the thermal structure of the shallow forearc. For the northern Honshu subduction zone, predicted interface temperatures at z_f = 50 km using Equation 1 range from ∼100 °C for μ′ = 0 to >1000 °C for μ′ = 0.1, where the calculated apparent shear stresses on the fault at z_f = 50 km are 0 and 150 MPa, respectively.

In a study of ten subduction zones, Gao and Wang (2014) combined surface heat-flux measurements with two-dimensional finite-element numerical models to estimate the strength of the fault interface in the seismogenic zone. They concluded that the subduction interface is relatively weak with μ′ ranging from 0.02 (northern Cascadia, western North America) to 0.13 (Hikurangi, New Zealand). For six of the studied subduction zones, Gao and Wang (2014) estimated μ′ to be ∼0.03, but the uncertainties are large, typically ±0.03, owing to the considerable scatter in heat-flux measurements. In a similar study, England (2018) used the analytical expressions described above to analyze nine subduction zones with surface heat-flux data, concluding that 0.05 < μ′ < 0.07 along the subduction interface; uncertainties in μ′ calculated for individual subduction zones exhibit a greater range.

For northern Japan (Honshu), Gao and Wang (2014) demonstrated that surface heat-flux measurements were best fit by μ′ = 0.025 (green curve in Fig. 2), but are consistent with 0 ≤ μ′ ≤ 0.06. England (2018) found that surface heat-flux measurements are best fit by μ′ = 0.05 (red curve in Fig. 2) with a range of 0.04–0.07 (uncertainties in heat production expand this range to 0.03–0.08). Estimates of μ′ for the northern Honshu subduction zone by Gao and Wang (2014) and England (2018) agree within stated uncertainties, but the difference in the best-fitting μ′ values correspond to significant differences in calculated subduction interface temperatures. For z_f = 50 km, Gao and Wang’s (2014) preferred value of μ′ = 0.025 yields T_f = 320 °C, whereas England’s (2018) preferred value of μ′ = 0.05 yields T_f = 555 °C. Note that the predicted interface P-T conditions calculated by Gao and Wang (2014) for μ′ = 0.025 agree very well with the P-T conditions calculated using England’s (2018) analytical approximations for μ′ = 0.02–0.03 (Fig. 2). With the exception of the Hikurangi subduction zone, Gao and Wang’s (2014) and England’s (2018) estimates for μ′ agree within uncertainties, with the differences in the “best-fitting” values largely reflecting the scatter in the heat-flux measurements, uncertainties regarding radiogenic heat production discussed below, and different assumptions underlying the analytical and numerical approximations.

Rather than being a simple planar feature, the subduction-zone plate boundary is a complex shear zone of finite thickness that may vary considerably in space and time. In the Barbados, Costa Rica, Nankai (southwest Japan), and Honshu (northeast Japan) subduction zones, Integrated Ocean Drilling Program drill cores reveal that the active subduction fault zone is several tens of meters thick near the trench (<1 km below the seafloor) (Rowe et al., 2013). Based on observations of exhumed faults, Rowe et al. (2013) proposed that the thickness of the active shear zone increases to ∼100–350 m at 1–2 km below the seafloor, a thickness that is maintained down to ∼15 km depth. Similarly, Agard et al. (2018) concluded that, at depths of 30–80 km, the regular plate interface is <300 m thick, but must increase dramatically during periods of interface down-stepping accompanying detachment and exhumation of kilometer-thick slices of the subducting plate. In the subduction-channel model (Shreve and Cloos, 1986; Cloos and Shreve, 1988), the thickness of the plate boundary interface may exceed 1 km where large amounts of sediment are being subducted. Mélange zones, such as those exposed in the Central Belt of the Franciscan Complex, can be several kilometers thick (e.g., Cloos, 1982; Bebout and Penniston-Dorland, 2016). Thus, rather than acting as a planar heat source, shear heating would be distributed across the active shear zone, which may be several hundred meters or more across.

Within the plate-boundary shear zone, deformation processes change with increasing depth from brittle (frictional) to ductile (viscous or plastic). Several different approaches have been used to estimate the depth of the brittle-ductile transition in subduction shear zones (e.g., Wang and Tréhu, 2016). For example, Tichelaar and Ruff (1993) and England (2018) equated the maximum depth of thrust faulting, defined by earthquake focal mechanisms, along the interface with the maximum depth of brittle behavior. Other researchers have defined the brittle-ductile transition based on temperature linked to an assumed rheology (e.g., 350 °C for quartz-bearing rocks; Hyndman and Wang, 1993). The concept of a relatively simple T- or P-dependent brittle-ductile transition along the subduction shear zone is complicated by (1) variable fluid pressure along the plate interface (e.g., Saffer and Tobin, 2011), (2) a transitional regime where both brittle (frictional) and ductile (viscous) processes operate (e.g., Holdsworth, 2004; Fagereng et al., 2014; Gao and Wang, 2017), and (3) the wide range of rock types that may be present in the shear zone (e.g., Fagereng and Sibson, 2010). Aqueous fluids can promote brittle behavior where fluid pressures approach lithostatic pressures, but such fluids also promote “ductile” flow by forming weak, platy alteration minerals (clays, micas) and enhancing dissolution-precipitation deformation mechanisms (Holdsworth, 2004).

In some thermal models (e.g., Peacock and Wang, 1999), shear heating along the megathrust is restricted to the brittle, frictional regime. In these models, shear stresses in the frictional regime are assumed to be constant or to increase linearly with depth (Equation 4) to a specified depth or temperature beyond which lies the ductile regime where shear stresses and shear heating drop to zero. Other thermal models incorporate shear heating along the deeper, ductile portion of the megathrust. For example, Gao and Wang (2014) calculated shear stresses and shear heating along the ductile portion of the megathrust by employing a wet granite rheology, a 500-m-thick shear zone, and a parameterization that smooths the frictional-viscous transition. Gao and Wang (2017) went a step further by incorporating a zone of high pore-fluid pressure (resulting in low shear stresses) downdip of the seismogenic zone in the region of episodic tremor and slow slip. The plate-boundary shear zone is composed of a variety of rock types, and may include metamorphosed sediments, basalt, and gabbro derived from the downgoing plate; granite, gabbro, and gneiss derived from the overriding crystalline crust; and metamorphosed ultramafic rocks derived from the overriding mantle wedge. Most of the shear heating occurs in the shallow frictional regime where shear stresses are highest, such that the exact rheology chosen for the ductile regime is less critical, as long as it is relatively weak.

The abundance and distribution of radioactive elements in the overriding plate are not well constrained, which adds considerable uncertainty to estimates of the rate of shear heating and μ′ based on surface heat-flux measurements. As can be seen from Equation 1, for a given surface heat flux, higher rates of radiogenic heat production result in lower calculated rates of shear heating and μ′, and vice versa. As part of their modeling efforts, Gao and Wang (2014) compiled existing heat-production data for eight subduction zones, finding values for A ranging from 1.8 µW/m³ for the northern Japan and Nankai subduction zones to 0.8–0.9 µW/m³ for Chile and 0.08 µW/m³ for Tonga. In Gao and Wang’s (2014) models, an uncertainty in A of ±0.2 µW/m³ (and D = 15 km) translates to an uncertainty of ±3 mW/m² in the calculated surface heat flux and ±0.01 in the calculated μ′. England (2018) considered estimates of AD to be uncertain by ±50% (and possibly more), which results in uncertainties in calculated surface heat flux of up to 20 mW/m². For continental subduction zones, England (2018) showed that the uncertainty in calculated μ′ introduced by uncertainties in the radiogenic heat production is comparable to, and in some cases greater than, the uncertainty owing to the scatter in heat-flux measurements.

During subduction, hydrothermal circulation may persist within the uppermost basaltic crust, transporting heat trenchward from deeper levels in the subduction zone (e.g., Spinelli and Wang, 2008). The thermal effects of hydrothermal circulation within oceanic crustal aquifers during subduction has been investigated by Spinelli and Wang (2008, 2009), Harris et al. (2010, 2017), Rotman and Spinelli (2013), and Spinelli et al. (2018), among others. In these models, the redistribution of heat by hydrothermal circulation is simulated by increasing the thermal conductivity of a thin (∼600 m) layer at the top of the oceanic crust representing the high-permeability (∼10⁻⁹ m²) oceanic crustal aquifer. Numerical models suggest that hydrothermal circulation within the crustal oceanic aquifer can reduce predicted interface temperatures at 50 km depth by >100 °C for subduction zones, like Nankai and Cascadia, characterized by slow subduction of young lithosphere (Rotman and Spinelli, 2013; Harris et al., 2017). In contrast, the predicted thermal effect is minor (<25 °C) for subduction zones characterized by rapid subduction of old lithosphere, like northeast Japan (Rotman and Spinelli, 2013). The clearest case for hydrothermal circulation in the oceanic crustal aquifer during shallow subduction comes from the Nankai subduction zone, where dense surface heat-flux measurements reveal very high fluxes in the trench of ∼200 mW/m², far greater than the 130 mW/m² expected for conductively cooled 15 Ma oceanic lithosphere (Spinelli and Wang, 2008). The permeability structure of the crustal aquifer and the depth to which it persists before porosity collapses represent significant uncertainties.

On geologic time scales, Earth’s mantle flows in a viscous fashion. Predicting the thermal structure of subduction zones at depths >∼50 km requires accurately simulating flow in the mantle wedge. Flow in the mantle wedge is driven primarily by mechanical coupling with the underlying subducting slab and is commonly referred to as corner flow. Batchelor (1967) derived a useful analytical expression for corner flow in an isoviscous fluid, which can be used to simulate mantle-wedge flow in subduction-zone thermal models (e.g., Peacock and Wang, 1999). These isoviscous corner-flow models demonstrated that the subducting slab warms considerably as it descends beneath the flowing mantle wedge, but the predicted cool slab-surface temperatures beneath the arc conflict with arc geochemical evidence for higher temperatures and partial melting of subducted sediments. This inconsistency was reconciled by employing more appropriate (non-Newtonian) rheologies for the mantle wedge, which result in more vigorous induced convection and higher predicted slab-surface temperatures. Many researchers have adopted Karato and Wu’s (1993) temperature- and stress-dependent rheology for wet olivine (e.g., van Keken et al., 2002; Currie et al., 2004; Wada et al., 2008). The effect of different mantle-wedge rheologies on the thermal structure of the mantle wedge can be seen in Figure 3, adapted from van Keken et al.’s (2008) subduction-zone benchmark paper. Modeling the mantle wedge as an isoviscous fluid results in flow lines and isotherms parallel to the boundaries of the mantle wedge (Fig. 3A). In contrast, using an olivine dislocation-creep rheology for the mantle wedge results in enhanced flow of material toward the wedge tip, higher wedge temperatures, and enhanced thermal gradients across the slab surface (Fig. 3B); the thermal boundary layer at the base of the mantle wedge is thinner and predicted temperatures along the slab surface are higher.

Between great megathrust earthquakes, the shallow (<30–50 km) plate interface is commonly considered to be “locked”. Viewed over time scales of thousands to millions of years involving many earthquake cycles, the subducting plate slides freely beneath the upper plate, and the shallow plate interface is effectively “decoupled”. Indeed, large-scale models of mantle convection commonly introduce an inclined plane of weak plate-boundary nodes extending from the surface down to depths of 100 km or more in order to simulate subducting plates (e.g., Zhong and Gurnis, 1998). Downdip of the seismogenic zone, viscous flow in the mantle wedge is driven, in part, by mechanical coupling with the subducting plate. Low surface heat fluxes and seismological observations require the forearc “nose” of the mantle wedge to be relatively cold and not participate in mantle flow (Fig. 1). Closer to the arc, viscous flow of the mantle wedge is required to reproduce the observed high surface heat fluxes and generate the high temperatures (>1200 °C) needed for partial melting beneath the volcanic arc. Thus, the plate-interface transition from decoupled to coupled must lie at depths greater than the forearc Moho (typically 30–40 km) and less than the sub-arc slab depth (∼100 km).

The depth and cause of the transition from a decoupled to a coupled plate interface has been the focus of considerable attention since the pioneering study of Furukawa (1993), who estimated that the transition occurred at ∼70 km depth based on surface heat-flow measurements above the northeast Japan subduction zone. Various approaches have been used to numerically simulate the decoupled-to-coupled transition along the subduction interface including: (1) imposing a no-flow condition for the near-trench portion of the mantle wedge (Peacock and Wang, 1999; van Keken et al., 2002; Currie et al., 2004), (2) prescribing a weak or T-dependent rheology for the subduction shear zone (e.g., Conder, 2005; Wada et al., 2008; van Dinther et al., 2013), and (3) defining a velocity discontinuity across the interface that gradually decreases with depth (Kneller et al., 2005, 2007; Abers et al., 2006). Wada et al. (2008) investigated the interface decoupling-to-coupling transition by constructing a two-dimensional finite-element model with a 100-m-thick uniform-viscosity layer along the plate interface. In numerical simulations where the viscous strength of the interface is less than that of the mantle wedge, the overlying forearc mantle wedge stagnates and effectively decouples from the subducting slab. The nonlinear rheology of the mantle results in a sharp downdip transition, over several kilometers, from decoupled to fully coupled along the subduction interface. For the relatively warm northern Cascadia subduction zone, the numerical simulations of Wada et al. (2008) suggest a maximum decoupling depth (MDD) of 70–80 km, consistent with surface heat-flow measurements and an extensively serpentinized forearc mantle wedge. Wada and Wang (2009) expanded the analysis of Wada et al. (2008) to 17 subduction zones, spanning a range in convergence rates and slab ages, and concluded that a common MDD of 70–80 km best predicted the temperatures and fluid fluxes required to generate arc magmas and the relatively uniform slab depth (∼100 km) beneath volcanic arcs.

Syracuse et al. (2010) presented two-dimensional finite-element models for 56 subduction-zone segments around the world that demonstrate the range in predicted thermal structures as a function of (1) subduction geometry, (2) convergence rate, (3) age of incoming lithosphere, (4) nature of the upper plate (e.g., continental versus oceanic), and (5) partial coupling between the subducting slab and overriding plate. The first four parameters were defined by observations specific to each subduction-zone segment. For the fifth parameter, four different partial-coupling models were considered, highlighting the importance of, and uncertainty regarding, convection in the mantle wedge. The range in predicted slab-surface P-T paths for the 56 subduction-zone segments and two of the four different models for slab-mantle coupling are depicted in Figure 4 (Syracuse et al., 2010). Predicted P-T paths are broadly similar for depths <60 km and depths >200 km; predicted slab-surface temperatures increase moderately from 0 to 60 km depth and increase slowly at depths >200 km. Between 60 and 200 km depth, predicted slab-surface temperatures increase rapidly, with the depth of rapid heating dependent on the choice of partial-coupling model. In the simplest case (Fig. 4A), with the MDD fixed at 80 km, subducting slab surfaces heat rapidly at 80 km depth and predicted slab-surface temperatures at 100 km depth (P = 3.3 GPa) range from 650 °C to 950 °C (Syracuse et al., 2010). Syracuse et al.’s (2010) models do not consider shear heating along the shallow plate interface and, as a result, predict relatively cool slab-surface temperatures beneath the forearc, with most models predicting T <300 °C at 50 km depth (Figs. 4 and 5). The highest temperatures are predicted for the Cascadia and Mexico subduction zones characterized by the subduction of young (7–11 Ma) oceanic lithosphere at moderate convergence rates (30–50 mm/yr). The coldest temperatures are predicted for subduction zones in the western Pacific, characterized by the subduction of old (>100 Ma) lithosphere and rapid convergence rates (>80 mm/yr). In the other three slab-mantle coupling models, the MDD varies from 50 to 175 km. Figure 4B shows the predicted slab-surface P-T paths for the same 56 subduction zone segments where the MDD is defined by the 550 °C isotherm. The MDD is found to vary from 54 to 89 km, and predicted slab-surface temperatures at 100 km depth (P = 3.3 GPa) range from 680 °C to 940 °C (Syracuse et al., 2010). Note that for the T-dependent coupling model, the rapid convergence and steep dip of the Tonga subduction zone results in anomalously cold predicted slab-surface temperatures (270 °C at 100 km depth) and a very deep MDD at 129 km (Syracuse et al., 2010).

Arcay (2017) constructed a two-dimensional thermal-mechanical model including a weak subduction channel to investigate how the maximum decoupling depth varied as a function of the viscosity contrast (parameterized as the activation energy contrast, ΔE_a) between the subducting crust and mantle. For subduction zones where the incoming lithosphere is >100 Ma, a low viscosity contrast leads to a shallow decoupling depth (e.g., MDD = 60 km depth for ΔE_a = 20 kJ/mol) and warmer forearc mantle. A high viscosity contrast leads to deeper decoupling (e.g., MDD = 135 km depth for ΔE_a = 120 kJ/mol) and a colder mantle wedge. Arcay (2017) concluded that constraints on slab-surface temperatures and temperatures of arc magmatism are best satisfied by an intermediate value for the activation energy contrast (80–120 kJ/mol) coupled with localized weakening of the mantle by H₂O.

During subduction, rocks in the downgoing plate undergo metamorphism as pressure and temperature increase. Constructing an accurate petrologic model of a subduction zone requires knowledge of the P-T structure and paths (derived from thermal models) and the bulk composition of the rocks in the subduction zone (sedimentary rocks, mafic oceanic crust, and ultramafic oceanic mantle). Metamorphosed mafic rocks (basalt, gabbro) contain minerals, such as amphibole and feldspar, that exhibit complex solid solutions, and metamorphic reactions in mafic bulk compositions are smeared out in P-T space. Thus, early petrologic models for metabasalt made use of the metamorphic facies concept (e.g., Eskola, 1920; Miyashiro, 1973), where characteristic mineral assemblages, rather than specific metamorphic reactions, defined regions in P-T space. In contrast to metabasalts, minerals in metasedimentary and meta-ultramafic bulk compositions tend to exhibit less solid solution, and constructing phase diagrams for these compositions is generally simpler.

In the past several decades, researchers have constructed phase diagrams for subduction-zone P-T conditions and bulk compositions based on experimental studies, mineral assemblages observed in the field, and thermodynamic calculations. Figure 6 shows phase diagrams and maximum H₂O contents constructed for metamorphosed mid-ocean ridge basalt using these three different approaches. Figure 6A is based on extensive experiments for T > 550 °C (Schmidt and Poli, 1998) and mineral assemblages observed in nature for T < 550 °C (Schmidt and Poli, 2014). Mineral metastability and sluggish reaction rates make it difficult to construct a phase diagram for metabasalts at T < ∼600 °C based solely on experiments (Hacker et al., 2003a). Phase diagrams constructed using mineral assemblages observed in the field need to be coupled to experiments (Fig. 6A) or thermodynamic calculations (Fig. 6B) that define key reactions bounding metamorphic facies (e.g., Hacker et al., 2003a). In experiments, the bulk composition of the system is specified by the experimentalist, whereas mineral assemblages observed in nature represent a range in bulk composition. Increasingly, thermodynamic calculations have been used to construct subduction-zone phase diagrams. For example, pseudosections (P-T phase diagrams for a specific bulk composition) can be constructed using Perple_X, a powerful thermodynamic software package based on minimization of Gibbs free energy (Connolly and Kerrick, 1987; Connolly, 2005). Figure 6C shows a phase diagram and maximum H₂O contents constructed by Hacker (2008) using the Perple_X package for a mid-ocean ridge basalt composition. One of the advantages of the Perple_X thermodynamic package is the prediction of mineral modes (percentages of minerals that make up the rock), which can be combined with experimental measurements of seismic properties (e.g., Christensen, 1996) to predict the seismological structure of a subduction zone that in turn can be tested against seismological observations. The Perple_X and other thermodynamic approaches benefit from an internally consistent thermodynamic database, but rely on thermodynamic properties of minerals, determinations of which vary in accuracy and precision.

The three phase diagrams for mid-ocean ridge basalt presented in Figure 6 are broadly similar, but differ in detail. The stability field for lawsonite (11 wt% bound H₂O) in metabasalts remains uncertain. Experiments demonstrate that lawsonite is stable to high pressures and temperatures (Fig. 6A), but lawsonite-bearing eclogites are rare in nature (Hacker et al., 2003a). Similarly, the experimentally based phase diagram (Fig. 6A) identifies chloritoid (7.5 wt% bound H₂O) as an important mineral in metamorphosed basalts, but chloritoid is rarely observed in naturally metamorphosed basalts (Hacker et al., 2003a) and is not included in the phase diagrams based on thermodynamic calculations (Figs. 6B and 6C).

Phase diagrams have also been constructed using these different approaches for sedimentary (e.g., Kerrick and Connolly, 2001; Poli and Schmidt, 2002; Hacker, 2008) and ultramafic (e.g., Ulmer and Trommsdorff, 1995; Schmidt and Poli, 1998; Hacker et al., 2003a; Rüpke et al., 2004; Fumagalli and Poli, 2005; Hacker, 2008) bulk compositions. For a detailed discussion of the differences between these phase diagrams, I refer the reader to Hacker et al. (2003a), Hacker (2008), and Schmidt and Poli (2014).

Hacker (2008) investigated where H₂O is released from the slab during subduction, and how much H₂O is subducted past the arc. Key findings include: (1) in warm subduction zones, metabasalts largely dehydrate at <100 km depth, whereas in cool subduction zones, metabasalts remain hydrated (∼3 wt% H₂O) to greater depths; (2) the amount of H₂O subducted by metabasalts in cold subduction zones depends on the initial H₂O content; (3) like metabasalts, meta-ultramafic rocks are predicted to undergo little dehydration in cold subduction zones at <100 km depth, but complete dehydration in warm slabs; (4) metamorphosed sediments can transport significant amounts of H₂O to depths >100 km, even in warm subduction zones; and (5) the amount of H₂O subducted in metasediments is proportional to the bulk-rock K₂O content, which stabilizes white mica (phengite).

The seismogenic zone of a subduction megathrust, which is the depth range inferred to be capable of producing large stick-slip earthquakes, may be defined in several ways (Wang and Tréhu, 2016). For subduction zones that have experienced recent great earthquakes, the extent of the seismogenic zone can be estimated by inverting geophysical observations. For example, Simons et al. (2011) inverted strong-motion, teleseismic, geodetic, and tsunami data sets for the great 2011 Tohoku earthquake (offshore Japan), concluding that most slip occurred at depths <40 km, but some slip occurred at 40–60 km depth. In subduction zones that have not experienced a recent great earthquake, the seismogenic zone is commonly defined by the extent of small to moderate thrust earthquakes. Small to moderate thrust earthquakes, including the aftershocks of great earthquakes, tend to define a slightly deeper downdip limit to the seismogenic zone.

Small to moderate thrust earthquakes occur at greater depths in cool subduction zones compared to warm subduction zones, suggesting that a critical temperature may control the maximum depth of thrust faulting on the subduction interface. This is also supported by recent seismo-thermo-mechanical models exploring T-dependent rheologies for the oceanic crust (van Dinther et al., 2013). If the downdip limit of the seismogenic zone could be defined by a critical temperature, then the downdip limit of thrust earthquakes could be used to constrain thermal models, and the critical temperature could help define the seismogenic zone of subduction zones, like Cascadia, that have not ruptured in historic times. Unfortunately, a number of studies suggest that temperatures at the downdip limit of the seismogenic zone vary considerably. For example, using analytical expressions to estimate temperatures along the subduction interface, Tichelaar and Ruff (1993) found that the maximum depth of interplate thrust earthquakes occurred at ∼400 °C where the overlying continental crust was thick, and ∼550 °C where the overlying continental crust was thin, assuming a constant coefficient of friction. If shear stresses are assumed to be constant with depth, then the maximum depth of thrust earthquakes could be explained by a single critical temperature of ∼250 °C (Tichelaar and Ruff, 1993). England (2018) found a wide variation in estimated temperatures at the downdip limit of small- to moderate-sized thrust earthquakes, ranging from 250 °C at 20 km depth for the Hikurangi subduction zone to ∼700 °C at 60 km depth for cold subduction zones like Honshu and Kermadec (southern Pacific). Rather than a constant interface temperature, England’s (2018) analysis of 91 subduction zones found the surface heat flux (corrected for radiogenic heating) above the maximum depth of thrust faulting to be relatively constant, at ∼40 mW/m², consistent with a constant μ′ = 0.07 ± 0.01 along the plate interface. Although the historical seismic record is relatively short, it seems unlikely that the wide variation in calculated temperatures is an artifact of the sparse record, suggesting that the downdip limit of the seismogenic zone is controlled by a number of factors in addition to temperature, including pore-fluid pressure and the type of rocks along the subduction interface.

There is very good agreement between thermal-petrologic models and a prominent low-velocity layer atop the subducting slab observed in seismological studies of subduction zones. Analysis of teleseismic scattered waves (receiver functions) in subduction zones reveals a landward-dipping, low seismic-velocity layer extending to 50–150 km depth (e.g., Yuan et al., 2000; Rondenay et al., 2001; Bostock et al., 2002; Ferris et al., 2003; Nicholson et al., 2005; Abers et al., 2006). The dipping low-velocity layer is also revealed by analysis of local P- to S-converted waves (e.g., Helffrich and Abers, 1997) and dispersion of seismic body waves (e.g., Abers and Sarker, 1996; Abers, 2005). This layer is on the order of 5 km thick and encompasses many (but not all) Wadati-Benioff zone earthquakes (Abers, 2005). It is commonly interpreted as metamorphosed oceanic crust and sediments at the top of the subducting plate (e.g., Hacker et al., 2003b) that have lower seismic velocities than mantle rocks (Christensen, 1996). Consistent with the metabasalt phase diagram and calculated slab P-T paths (Hacker et al., 2003b; Abers et al., 2013), the low-velocity layer ends at relatively shallow depth (∼40–50 km) in the warm Cascadia subduction zone (e.g., Nicholson et al., 2005) and at greater depth (∼120–130 km) in the cooler Alaska subduction zone (northernmost Pacific; Abers et al., 2006). The spatial resolution of seismological studies continues to increase rapidly, and we may look forward to resolving multiple layers in the subducting slab and adjacent mantle wedge.

The spatial distribution of earthquakes within the subducting slab correlates strongly with the location of continuous (multivariant) metamorphic dehydration reactions in the subducting crust and mantle, as predicted by thermal-petrologic models (e.g., Hacker et al., 2003b; Abers et al., 2013). In warm subduction zones, intraslab earthquakes cease by 100 km depth, whereas intraslab earthquakes extend to depths of 700 km in cold western Pacific subduction zones. Kirby et al. (1996) proposed that intermediate-depth earthquakes, which occur between 50 and 300 km depth, reflect dehydration embrittlement where H₂O-rich fluids released by metamorphic reactions promote brittle reactivation of preexisting faults. In addition to increasing pore pressure (Raleigh and Paterson, 1965), dehydration reactions generate fine-grained, weak reaction products that may help promote brittle behavior. Ferrand et al. (2017) proposed that dehydration-driven stress transfer, rather than fluid overpressure, causes embrittlement. Other hypotheses have been proposed for intermediate-depth earthquakes such as viscous instabilities (e.g., Kelemen and Hirth, 2007) and thermal runaway (e.g., John et al., 2009), but dehydration-related embrittlement remains the leading hypothesis for intermediate-depth earthquakes near the top of the subducting slab.

Thermal-petrologic modeling by Peacock and Wang (1999) and Hacker et al. (2003b) has predicted the location of hydrated metamorphic rocks in warm and cold subducting slabs, which correlates strongly with the location of intermediate-depth earthquakes. In order to evaluate the dehydration embrittlement hypothesis, it is important to determine the precise location of the intermediate-depth earthquakes with respect to the plate boundary. The thermal-kinematic models described in this paper commonly use Wadati-Benioff zone seismicity to define the top of the subducting slab, so the precise location of earthquakes is best determined independently using receiver functions (e.g., Abers et al., 2006). An important result of such analysis is that intermediate-depth earthquakes occur both in the subducting oceanic crust and in the uppermost subducting mantle (Abers et al., 2013). It is well accepted that spatially heterogeneous hydration of the oceanic crust occurs at spreading centers and transform faults. Earthquakes in subducted oceanic crust appear to be linked to continuous dehydration reactions that occur in metabasalt with increasing T and P (e.g., Kirby et al., 1996; Hacker et al., 2003b). The dehydration embrittlement model for slab mantle earthquakes requires infiltration and incorporation of H₂O into the upper mantle and subsequent dehydration during subduction. Recent studies have documented that substantial hydration of the oceanic crust and uppermost mantle can occur in the outer-rise region along normal faults associated with plate bending (e.g., Ranero et al., 2003, 2005), and intermediate-depth earthquakes within the uppermost slab mantle likely reflect dehydration of hydrated mantle rocks.

Earthquakes defining the lower plane of double Wadati-Benioff zones may also be caused by dehydration embrittlement linked to serpentine dehydration reactions (Peacock, 2001). The dehydration embrittlement hypothesis for these earthquakes remains controversial because it requires the infiltration of H₂O several tens of kilometers deep into the upper mantle along outer-rise faults. Numerical modeling by Faccenda et al. (2009) demonstrated that deep slab hydration may be induced by pressure gradients associated with plate bending. The recent 2017 Chiapas, Mexico, earthquake provides compelling evidence that normal faulting can extend several tens of kilometers into the slab mantle (Melgar et al., 2018; Zhang and Brudzinski, 2019). Low seismic velocities, interpreted to reflect partial serpentinization, extend to 24 ± 5 km depth below the incoming slab Moho in the Mariana subduction zone (western Pacific; Cai et al., 2018). On the basis of numerical modeling, Kelemen and Hirth (2007) proposed that lower-plane Wadati-Benioff zone earthquakes reflect shear instabilities (thermal runaway) associated with the onset of highly localized viscous creep in anhydrous mantle, rather than dehydration embrittlement. The shear-instability hypothesis is supported by seismic tomography studies consistent with anhydrous slab mantle (e.g., Reynard et al., 2010) and deformation experiments (Ohuchi et al., 2017). Perhaps both dehydration and thermal runaway mechanisms operate at these depths, as envisioned in the dehydration-driven stress model of Ferrand (2019). Adding to the list of possibilities, Faccenda et al. (2012) proposed lower-plane Wadati-Benioff zone earthquakes are triggered by hydrofracturing resulting from the upward fluid flow driven by slab unbending.

As discussed above, the low surface heat fluxes observed in the subduction-zone forearcs require the forearc mantle wedge to be cold (<800 °C), and petrologic models indicate that these ultramafic mantle rocks could be hydrated if H₂O is present. Depending on the P-T conditions and bulk composition, hydrous phases such as antigorite (serpentine), chlorite, brucite, and talc may be stable in the mantle wedge. Multiple geophysical and geological observations, including reduced P-wave velocities (<8 km/s), suggest that the forearc mantle wedges are partially serpentinized (Hyndman and Peacock, 2003). Recently, Abers et al. (2017) challenged this interpretation for cold subduction zones, suggesting that subducting plates in cold subduction zones release insufficient H₂O at shallow depths to hydrate the mantle wedge, even when integrated over the life span of the subduction zone. This suggestion depends, in part, on the assumption that subducted materials have uniform, homogenous H₂O contents. Hydration of oceanic crust and mantle is known to be spatially heterogeneous and concentrated along faults and fractures. If this heterogeneity is preserved during subduction, then H₂O-rich veins may dehydrate at shallower levels than predicted assuming homogenous bulk H₂O contents.

Numerous researchers have used the depth to the slab surface beneath arc volcanoes to constrain the thermal structure of subduction zones. Using the slab depth beneath the volcanic front to constrain subduction-zone temperatures requires assumptions about where partial melting occurs beneath the volcanic front—e.g., within the subducting slab, along the slab interface, and/or within the overlying mantle wedge. For example, based on the assumption that arc andesites are derived from partial melting of subducted basaltic crust, early thermal models of subduction zones relied upon high rates of shear heating to generate high slab-surface temperatures of ∼1000 °C beneath the volcanic arc (e.g., Turcotte and Schubert, 1973). In the 1980s and 1990s, most researchers tended to focus on the relatively constant depth to the slab beneath the arc volcanic front (∼110–125 km) (e.g., Gill, 1981; Tatsumi, 1986). The relatively constant slab depths suggested that arc magmatism was triggered by a critical P-dependent dehydration (or partial-melting) reaction in the subducting slab. More recent work has emphasized the observed global range in slab depths beneath volcanic arcs (∼65–170 km) (e.g., England et al., 2004; Syracuse and Abers, 2006), which argues against a simple P-dependent reaction controlling arc magmatism.

A growing consensus has developed that the location of the volcanic front is controlled by the thermal structure of the mantle wedge, specifically the closest approach to the trench of a particular isotherm or solidus. England et al. (2004) demonstrated that the global range in slab depths beneath the arc correlated inversely with slab descent velocity, suggesting that arc volcanism was controlled by the temperature at the top of the slab or in the mantle wedge. Later, England and Katz (2010) ruled out the slab-temperature hypothesis, noting that isotherms at the top of the slab were subparallel to the slab interface, which should act to smear out arc volcanism with depth if controlled by the slab-surface temperature. England and Katz (2010) concluded that the volcanic front was located above the mantle wedge where temperatures were 1250–1325 °C at pressures of 2–3.5 GPa, corresponding to the solidus of peridotite containing 200–500 ppm H₂O. Using a larger data set than England et al. (2004), which included subduction zones with steeper slabs, Syracuse and Abers (2006) found slab depths correlated with slab dip (and not descent velocity), concluding that the partial-melting region in the mantle wedge was displaced upward from the Wadati-Benioff zone by a constant-thickness boundary layer. On the basis of numerical modeling, Perrin et al. (2018) concluded that the trenchward location of the anhydrous peridotite solidus governed the position of the volcanic arc. Convection in the mantle wedge induced by the subducting slab exerts a primary control on the thermal structure of the wedge, but additional processes and parameters play an important role, including heat transported by rising magma (e.g., England and Katz, 2010), the thickness of the overriding lithosphere relative to the slab coupling depth (Perrin et al., 2018), and secondary convection (Currie and Hyndman, 2006).

The distinctive geochemistry of arc magmas appears to require a contribution to the arc magma source region of a partial melt derived from subducted sediments (e.g., Gill, 1981). In the past decade, fluid (melt) geothermometers for magmas have been developed based on the trace elements hosted in accessory minerals that remain stable during partial melting (Hermann and Spandler, 2008; Plank et al., 2009). For example, Cooper et al. (2012) measured H₂O/Ce ratios in melt inclusions in arc lavas from ten different subduction zones and calculated slab-surface temperatures of ∼730–900 °C. These high slab-surface temperatures are consistent with recent thermal models incorporating flow in a non-Newtonian mantle wedge driven by slab coupling (e.g., Syracuse et al., 2010). The inferred slab-surface temperatures assume that sediment melting occurs atop the subducting slab rather than in the mantle wedge. Based on analyses of sedimentary rocks metamorphosed at high pressures (2.7–5 GPa), Behn et al. (2011) concluded that diagnostic “slab signature” elements (e.g., Th, Sr, Pb, Nd) do not become mobile until temperatures exceed ∼1050 °C, significantly higher than the onset of partial melting at the wet solidus. Thermal modeling suggests that temperatures >1050 °C are difficult to achieve along the slab surface beneath volcanic arcs, but are readily achieved in the hot core of the convecting mantle wedge. The geochemical evidence for a high-T sediment signature supports dynamical subduction models that predict metasedimentary rocks detaching from the subducting slab form buoyant diapirs that intrude the overlying mantle wedge (e.g., Gerya and Yuen, 2003; Currie et al., 2007; Gerya, 2011; Behn et al., 2011; Marschall and Schumacher, 2012).

Metamorphic rocks form at depth in the Earth. The P-T conditions recorded by metamorphic rocks now exposed at the surface provide insight into the thermal structure of plate margins at depths greater than can be obtained by drilling. Blueschist- and eclogite-facies metamorphic rocks reflect relatively high-P, low-T conditions consistent with formation in subduction zones where observed forearc heat fluxes are low. In contrast, hornfels-, amphibolite-, and granulite-facies metamorphic rocks reflect relatively low-P, high-T conditions consistent with formation in magmatic arcs where surface heat fluxes are high. In the geologic record, the juxtaposition of these “paired” metamorphic belts is commonly used to infer the presence of an ancient subduction zone, with the direction of ancient subduction taken to be from the high-P, low-T belt (accretionary prism) toward the low-P, high-T belt (magmatic arc) (Ernst, 1971; Miyashiro, 1972).

In principle, we can invert the P-T conditions recorded by blueschists and eclogites to constrain important thermal parameters, such as the rate of shear heating along the subduction shear zone. In practice, however, there are significant challenges with using metamorphic P-T data to constrain thermal parameters including: (1) metamorphic rocks now at the surface formed millions of years ago when parameters like convergence rate, slab geometry, and slab age are not well constrained; and (2) there is considerable uncertainty as to when and where in the tectonic history of an orogenic belt the P-T conditions recorded by metamorphic rocks were achieved. A recent debate in the literature highlights the challenges of using metamorphic rocks to constrain subduction-zone temperatures.

Using a global data set, Penniston-Dorland et al. (2015) demonstrated that temperatures recorded by blueschists and eclogites at their maximum pressure were, on average, 100–300 °C warmer than predicted by thermal models presented by Gerya et al. (2002) and Syracuse et al. (2010). Penniston-Dorland et al. (2015) explored possible biases in both the petrologic data (e.g., peak temperatures may be achieved during exhumation) and in the thermal models (e.g., lack of shear heating). They noted that the prograde P-T paths recorded by blueschists and eclogites were also warmer than predicted by the thermal models. Kohn et al. (2018) extended the Penniston-Dorland et al. (2015) analysis to include surface heat-flow observations and used analytical expressions to argue that subduction-zone thermal models must incorporate shear heating (with μ′ = 0.05 ± 0.015) or else the models may underestimate slab-surface temperatures by 100–500 °C at 30–80 km depth.

Agard et al.’s (2018) global compilation of P-T conditions recorded by subducted oceanic rocks excluded continental fragments and yielded a geologic data set indicating significantly cooler conditions than the Penniston-Dorland et al. (2015) compilation. Agard et al.’s (2018),P-T data set is more consistent with the predictions of some thermal models (e.g., Gerya et al., 2002), but still significantly warmer than those of other thermal models (e.g., Syracuse et al., 2010).

van Keken et al. (2018) assessed two sets of thermal models with and without shear heating (Wada and Wang [2009] and van Keken et al. [2011], respectively) against a broadly similar global data set of P-T conditions recorded by exhumed subduction-zone rocks. van Keken et al. (2018) demonstrated that including reasonable amounts of shear heating in the thermal models increases slab-surface temperatures, but the average temperature of the 7-km-thick subducting oceanic crust increases by <50 °C. van Keken et al. (2018) concluded that blueschists and eclogites must be exhumed under relatively warm, anomalous conditions such as those experienced during subduction initiation or the subduction of very young lithosphere. van Keken et al. (2019) demonstrated that analytical approximations for frictional heating should not be extended past the brittle-ductile transition, reinforcing their earlier proposal (van Keken et al., 2018) that subduction-zone metamorphic rocks record anomalously warm conditions.

The P-T paths depicted in Figure 5 illustrate that models that neglect shear heating along the shallow plate interface (e.g., Syracuse et al., 2010) predict cooler slab-interface temperatures beneath the forearc than models that include modest shear heating (e.g., Wada and Wang, 2009). In contrast, both sets of models predict similar slab-interface temperatures at greater depth where slab heating is dominated by the overlying convecting mantle wedge. For P < 2 GPa, some of the discrepancy (perhaps 100–300 °C) between the temperatures recorded by blueschist and eclogites and those predicted by thermal models can be eliminated by including shear heating along the interface. For example, Wada and Wang’s (2009) model with μ′ = 0.03 yields slab surface temperatures that are 150 °C hotter at P = 1.5 GPa than the models of Syracuse et al. (2010). And, as discussed earlier, modest rates of shear heating appear to be required to fit the observed surface heat-flux measurements in the forearc.

Most hybrid kinematic-dynamic thermal models do not explicitly simulate the detachment and exhumation of rocks from the subducting slab, whereas select fully dynamical models, such as that of Gerya et al. (2002), simulate these processes. Penniston-Dorland et al. (2015) illustrated that “average” subduction-interface temperatures predicted by Gerya et al. (2002) are ∼100–150 °C warmer than the temperatures predicted by Syracuse et al. (2010) for P = 1–2 GPa, but are still lower than the temperatures recorded by most blueschists and eclogites. The higher temperatures predicted by the Gerya et al. (2002) models reflect, in part, the heat advected by rock underplating and circulation in the subduction channel; such processes would also affect the predicted forearc heat flow. More importantly, the high temperatures predicted by Gerya et al. (2002) for tectonic blocks in the circulating subduction channel appear to reflect peak metamorphic conditions attained during the early (warmer) stages of subduction. Rapid cooling of subduction zones after initiation is a fundamental prediction of virtually all subduction-zone thermal models and is consistent with numerous field-based studies of paleo–subduction zones that document that the highest-temperature metamorphic rocks form early in the history of the subduction zone (e.g., Platt, 1975; Cloos, 1985; Krebs et al., 2008; Agard et al., 2018).

Resolving the discrepancy between the P-T conditions recorded by blueschists and eclogites and the P-T conditions predicted by most thermal models is an important area for future research. In my opinion, there is compelling evidence that thermal models of subduction zones should incorporate modest amounts of shear heating along the brittle portion of the subduction interface, but shear heating alone is insufficient to explain the discrepancy between the rock record (especially high-P, low-T metagabbros) and the thermal models. Most blueschists and eclogites likely record anomalously warm subduction-zone conditions, such as those associated with subduction initiation. Following the important insights revealed by Gerya et al. (2002), work needs to be done to assess the thermal effects of underplating, exhumation, and subduction-channel circulation on the thermal structure of subduction-zone forearcs. I note that this has been an active area of research with respect to the exhumation of ultrahigh-pressure metamorphic rocks (e.g., Warren et al., 2008; Hacker et al., 2013).

Ideally, P-T conditions inferred from analysis of blueschist-facies rocks being brought to the surface today could be used to constrain the rate of shear heating because the convergence rate and slab geometry are known. In the Mariana subduction zone, active serpentine mud volcanoes occur on the forearc seafloor (e.g., Fryer et al., 1985, 1999). These mud volcanoes are fed by fluids derived from the subducting Pacific plate and entrain clasts of serpentinite and, less commonly, blueschist-facies rocks (e.g., Maekawa et al., 1993). Maekawa et al. (1993) estimated peak metamorphic conditions of T = 150–250 °C and P = 0.5–0.6 GPa for a blueschist clast derived from the Conical Seamount, which Peacock (1996) used to calculate that μ′ = 0.024–0.049 for the Mariana subduction zone. However, Tamblyn et al. (2019) recently presented U-Pb ages of 51.1 ± 1.2 Ma (zircon) and 47.5 ± 2.0 Ma (rutile) for a blueschist clast recovered from the South Chamarro Seamount, strongly suggesting that the Mariana blueschists formed during the Eocene initiation of subduction in the western Pacific and therefore cannot be used to infer P-T conditions at depth in the modern-day subduction zone.

Over the past two decades, advances in thermal and petrologic modeling, coupled with increasingly detailed seismological and geological observations, have improved our understanding of subduction zones. I expect thermal and petrologic modeling to continue to advance subduction-zone science as we further test model predictions against increasingly sophisticated geologic, seismological, and geophysical observations, and evaluate the importance of dynamical processes not captured in hybrid kinematic-dynamic models. Important questions that I hope will be answered in the near future include: What controls the depth at which the slab becomes fully coupled to the mantle wedge? What is the permeability structure of subduction zones and the spatial and temporal patterns of fluid flow? How does fluid infiltration and hydration affect the rheology and flow structure of the mantle wedge?

Subduction-zone science is inherently multidisciplinary, and advances in thermal and petrologic modeling will continue to be linked to advances in a wide range of observational and experimental geoscience disciplines. For example, experimental measurements of seismic velocities in natural rocks at elevated P and T (e.g., Christensen, 1996) have allowed researchers to invert seismological observations to image the basalt-to-eclogite transition in subducting oceanic crust. Further advances in this area will require the ability to differentiate among the competing effects of rock compositions and mineralogy, fluids, and rock and crack anisotropy on seismic velocities.

This paper has focused on hybrid kinematic-dynamic models of subduction zones. Looking ahead, we need to reconcile further the predictions of these models with the geologic record and with the predictions generated by fully dynamical models. Dynamical models of subduction zones have revealed a number of processes, like subduction-channel flow and diapirism, that are not represented in hybrid kinematic-dynamic models. If these processes commonly occur in subduction zones, then hybrid kinematic-dynamic models may have limited utility despite their success at explaining a number of important first-order observations like the spatial distribution of intermediate-depth earthquakes. In my opinion, it remains to be determined whether such dynamical processes are a common feature of most (all) subduction zones or whether they occur only under a restricted set of circumstances. For example, does the geologic record suggest that subduction-channel return flow is common or rare in paleo–subduction zones? Might we be able to observe subduction channels or diapirs in modern subduction zones using increasingly high-resolution seismological techniques?

Finally, to what extent are modern subduction zones in approximate thermal steady state, excluding those subduction zones that have formed in the last few million years? Dynamical models suggest that the geometry and rate of subduction can fluctuate rapidly over time scales of a few million years (e.g., Billen and Arredondo, 2018), and there are well-documented examples of changes in slab geometry in the geologic record, such as the shallowing of the slab during the Laramide orogeny (United States; e.g., Bird, 1998). Can we better use the geologic record, such as the location of the volcanic arc over time, to constrain the stability or instability of slab geometry over millions of years? And how might fluctuations in slab geometry affect our understanding of flow in the mantle wedge?

I thank Taras Gerya and Ikuko Wada for their constructive reviews of an earlier version of this manuscript. I wish to express my deep appreciation to Gray Bebout, Steve Kirby, John Platt, and David Scholl for organizing the 1994 SUBCON Interdisciplinary Conference on the Subduction Process, editing the original Subduction: Top to Bottom volume, and editing this special Geosphere volume. This research was supported by funding from the University of British Columbia.