Diagenetic, metamorphic, and hydrogeologic consequences of hydrothermal circulation in subducting crust

Subsurface temperatures affect processes spanning all regions of subduction zones, including sediment alteration near the trench and slab alteration at subarc depths. For example, the thermal history of rocks and sediments passing through a subduction zone and/or residing near the plate-boundary fault controls the progress of the multiple reactions that affect frictional properties, fluid release, and potentially fault strength and deformation behavior (e.g., Hyndman and Wang, 1993; Moore and Saffer, 2001; Saffer et al., 2008). Deeper in the subduction zone, the magnitude and distribution of temperatures control metamorphic reaction progress, subduction dynamics, hydration of the mantle wedge, and melt generation below volcanic arcs (e.g., Hacker et al., 2003; van Keken and King, 2005; Peacock, 2009; Wada et al., 2012; Völker and Stipp, 2015).

Steady progress has been made in understanding the parameters that control the thermal state of subduction zones. The subducting plate age, convergence rate, and dip of the subducting slab have long been recognized as key parameters affecting subduction-zone temperatures (e.g., Cloos, 1985; Peacock, 1987; Dumitru, 1991; Kirby et al., 1991). The pattern of mantle-wedge flow, induced by interactions with the subducting slab, also exerts a strong control on the temperature distribution deeper in subduction zones (e.g., Wada and Wang, 2009).

More recently, scientists have recognized the importance of the thermal consequences of vigorous fluid circulation redistributing heat in subducting igneous crust (e.g., Kummer and Spinelli, 2008; Spinelli and Wang, 2008; Harris et al., 2017). This hydrothermal circulation in the igneous basement persists as the basement is subducted and can transport heat from the deeper part of the subducted basement to the trench; if sediment cover on the incoming plate is thick and laterally continuous, heat can be transported well seaward of the trench (Fig. 1). Evidence for persistent hydrothermal circulation and its main characteristics has been presented elsewhere (e.g., Spinelli and Wang, 2008; Wang et al., 2015; Hass and Harris, 2016) and summarized in this volume (Harris et al., 2017). Here, we review and further explore how this hydrothermal circulation can significantly impact the timing and spatial patterns of sediment diagenesis and the alteration of the upper parts of the subducting slab (i.e., sediments, basalt, gabbro, and upper few kilometers of oceanic mantle). We start with a brief description of petrologic and thermal models adopted in this study and then describe the impact of hydrothermal circulation on the temperature conditions of the incoming plate and the subducting slab beneath the forearc and the subarc regions. For the Nankai margin, we show that heat transported seaward of the trench accelerates alteration of sediments prior to subduction. For the Middle America Trench off Costa Rica and Mexico, we review how contrasting patterns of hydrothermal circulation impact sediment diagenesis differently. For the deeper part of the subduction zone, we first review observations that indicate delayed alteration of subducting basaltic crust in the Nankai margin and inhibited slab melting in southern Chile; then, we present new model results that tie many of the observations and inferences together. We conclude by discussing the hydrogeologic consequences of cooling subducting material via hydrothermal circulation, reviewing effects on fluid source distributions, fluid properties, and permeability of the basement aquifer in the oceanic crust.

Most fluid flow through oceanic lithosphere occurs in a high-permeability aquifer composed of pillow lavas and breccias that comprise the upper ∼600 m of the basement rock (Fisher, 1998; Becker and Davis, 2004). In generic lithologic models for oceanic lithosphere, the basaltic basement aquifer is underlain by a ∼1.4-km-thick layer of dikes and a ∼5-km-thick layer of gabbro (e.g., Hacker, 2008); these intrusive rocks have low permeability relative to the basaltic aquifer (e.g., Fisher, 1998). Away from mid-ocean ridges and other outcrops of basaltic basement rock (e.g., seamounts), sediments overlie the basement aquifer. These seafloor sediments have permeabilities ∼3–10 orders of magnitude lower than the basaltic basement (Spinelli et al., 2004), confining nearly all of the fluid flow to the basement aquifer. Within sediments, fault zones and sand bodies can have high permeability relative to the bulk of seafloor sediment; thus, they can focus fluid flow (e.g., Fisher and Hounslow, 1990; Saffer and Screaton, 2003; Saffer, 2010). Nonetheless, most thermally significant fluid flow in the oceanic crust occurs within the basement aquifer (e.g., Davis and Lister, 1977; Davis et al., 1997; Davis et al., 1999; Fisher, 2004).

Fluid circulation vigorous enough to redistribute large quantities of heat is pervasive in oceanic crust (Lister, 1972). In systems with open connections between the basement aquifer and the ocean, ventilated hydrothermal circulation occurs, extracting heat from the lithosphere (e.g., Davis et al., 1999; Fisher et al., 2003; Hutnak et al., 2008). Where sediment is thick and laterally continuous, insulated hydrothermal circulation redistributes heat within the basement aquifer. Sediment cover as little as ∼10 m thick is sufficient to inhibit thermally significant ventilation (Spinelli et al., 2004); thus, the lateral continuity of sediment cover is extremely important. On average, the global heat-flow data set indicates that ventilated hydrothermal circulation is present in crust as old as 60 Ma (Stein et al., 1995; Hasterok et al., 2011). Focused heat-flow surveys show that both ventilated and insulated hydrothermal circulation can extend to basement ages much greater than 60 Ma, indicating that in some situations, the permeability of upper basement remains large enough to host fluid circulation to great age and likely for the entire age span of oceanic basement (Fisher and Becker, 2000; Jarrard et al., 2003; Von Herzen, 2004; Fisher and Von Herzen, 2005; Harris and McNutt, 2007; Kawada et al., 2014; Yamano et al., 2014).

Vigorous fluid circulation tends to homogenize temperatures in the upper oceanic basement aquifer (e.g., Davis et al., 1997). The high permeability of the aquifer can persist as the crust is subducted, allowing fluid circulation to redistribute heat in the subducting crust (Spinelli and Wang, 2008). Thermal anomalies at some subduction zones (e.g., Nankai, Costa Rica, Chile, and Haida Gwaii) are most consistent with the subducting crust acting as a “heat exchanger,” in which fluid circulation mines heat from subducted crust and transports it to the crust seaward of the trench (Kummer and Spinelli, 2008; Harris et al., 2017). As a result, fluid circulation in subducting crust can simultaneously cool the subduction system landward of the trench and warm the system at and seaward of the trench (Fig. 1) (Spinelli and Wang, 2008).

For both the modeling case studies we review and new results we present, there are two main components: (1) modeling subduction-zone temperatures, which includes the thermal effects of hydrothermal circulation in the subducting basement; and (2) using modeled thermal histories for sediment or rock to estimate the progress of diagenetic or metamorphic reactions. There are key thermal connections between the incoming oceanic crust (i.e., seaward of the trench) and the subducted crust (e.g., Harris et al., 2017); thus, our subduction-zone thermal models include material seaward of the trench. Additionally, including the incoming plate in the thermal models allows us to examine alteration of incoming material resulting from heat transported by hydrothermal circulation in the subducting crust. We briefly review these models before reviewing studies from different subduction zones and presenting new results.

The base of the model is set to an isotherm, typically 1400 °C, the seaward side is set to an oceanic geotherm appropriate for the age of subducting lithosphere, and the landward side is set to a geotherm appropriate for a backarc. The velocity of the subducting plate is prescribed. The velocity of the mantle-wedge flow induced by the slab motion is calculated in PGCtherm2D but prescribed in the older code. In modeling the mantle-wedge flow, we use a temperature- and stress-dependent wet olivine rheology (Karato and Wu, 1993); the governing equation is not displayed here but can be found in van Keken et al. (2002) and Wada et al. (2008). Where the top of the subducting slab is <70 km deep, the overlying mantle wedge is stagnant; where the top of the slab is >70 km deep, the overlying viscous mantle wedge is permitted to flow (Wada and Wang, 2009). Heat production from radioactivity and frictional heating along a shallow portion of the plate-boundary fault can be incorporated. Frictional heating on the plate-boundary fault begins near the trench and extends downdip to the intersection of the continental Moho with the plate interface. The magnitude of frictional heating is calculated as the product of the slip rate and the shear stress τ =µ´σ_n, where µ´ is the effective friction coefficient and σ_n is the normal stress (Gao and Wang, 2014).

Because the calculation of Ra includes temperature-dependent fluid properties and the heat flux into the base of the aquifer, our simulations of the thermal effects of vigorous fluid circulation require an iterative procedure. We start by modeling subduction-zone temperatures with Nu = 1 in the aquifer (i.e., assuming no redistribution of heat by vigorous fluid circulation). We use the modeled temperatures to calculate fluid density, fluid viscosity, and the heat flux into the base of the aquifer for every aquifer element. We calculate Ra throughout the aquifer. We then use Ra to update Nu and increase the thermal conductivity of each aquifer element by a factor of Nu. With the new thermal conductivity values for the aquifer, we run the thermal model again, which can change the temperatures and fluid properties. We then use these temperatures to update Ra and Nu. This procedure is repeated until temperatures stabilize.

We use an example from the Nankai margin to show how heat transported from depth in subduction zones can warm the incoming sediment and accelerate diagenetic processes. Temperature-dependent diagenetic processes in these sediments include the opal-to-quartz transition and the smectite-to-illite transformation (e.g., Kastner et al., 1991). These diagenetic processes can strengthen the sediments and influence the position of the décollement and modulate fault-slip processes (Moore and Saffer, 2001; Trütner et al., 2015; Hüpers et al., 2017; Kameda et al., 2017; Lauer et al., 2017).

The Nankai margin of southern Japan is one of the most studied subduction zones, with data from >16 scientific drilling expeditions and associated activities generating a large volume of high-quality, heat-flux measurements. A first-order observation of the thermal state of the margin is that heat flux in the Nankai Trough is anomalously high. In some of the early detailed thermal studies of the margin (e.g., Yamano et al., 1984; Kinoshita and Yamano, 1986; Nagihara et al., 1989), the high heat flux from Nankai Trough was viewed as a curious observation, particularly after earlier observations of low heat flux from trenches around the Pacific (e.g., Watanabe et al., 1977). Surprise regarding this observation is encapsulated well in the title from Yamano et al. (1984), “Nankai Trough: A hot trench?”. A synthesis of thermal studies of the margin demonstrates that heat-flux measurements within 30 km of the deformation front average 20% above the expected values for conductively cooled lithosphere (Fig. 2) (Harris et al., 2013). The heat-flux anomaly in the trench extends for at least 300 km along strike (Fig. 2) (Yamano et al., 2003; Harris et al., 2013). The high heat flux in the trench is striking because heat-flux measurements >30 km seaward of the deformation front average 20% below predicted values for conductively cooled lithosphere (Fig. 2) (Harris et al., 2013). However, these values >30 km seaward of the deformation front are consistent with global averages (e.g., Stein and Stein, 1994; Hasterok et al., 2011), indicative of some hydrothermal removal of heat from the plate. Thus, the anomalously high heat flux in the trench requires an additional heat source to the oceanic lithosphere near the frontal subduction zone (Harris et al., 2013). Models of hydrothermal circulation in the subducting basement indicate that heat transported from depth within the subduction zone (i.e., updip) can explain the anomalously high heat flux at the trench (Spinelli and Wang, 2008; Harris et al., 2013; Harris et al., 2017).

Ocean Drilling Program (ODP) Site 1173 (Fig. 2) is ∼11 km seaward of the Nankai Trough and ∼35 km west of the inferred location of the fossil Shikoku Basin spreading center (Moore et al., 2001; Mahony et al., 2011). Magnetic anomalies in Shikoku Basin indicate the lithosphere at Site 1173 is ∼19 m.y. old (Okino et al., 1994; Sdrolias et al., 2004; Muller et al., 2008), and the oldest sediments are ∼15 m.y. old (Shipboard Scientific Party, 2001). Spreading ceased ∼15 m.y. ago (Okino et al., 1998). This site is within the region of anomalously high heat flux with a measured value of ∼180 mW m^–2, ∼1.5 times that of predicted values of 115–129 mW m^–2 for 19–15 Ma conductively cooled lithosphere (Yamano et al., 2003; Harris et al., 2013).

An early hypothesis for the cause of the anomalously high heat flux in the Nankai Trough was that the lithosphere is thermally young (i.e., it has only been cooling conductively for the past ∼5.5 m.y. [e.g., Yamano et al., 2003]). In this scenario, the lithosphere in the region is heated by volcanism that occurs after spreading along the Shikoku Basin spreading center ceased (∼15 m.y. ago). Based on the ages of basalt samples recovered from seamounts in Shikoku Basin, local postspreading volcanism continued until at least 7.7 m.y. ago (Fig. 2) (Sato et al., 2002; Ishizuka et al., 2009).

The current anomalously high heat flux at Site 1173 (or any particular site within ∼30 km of the deformation front in the Nankai Trough) is not necessarily indicative of anomalously high temperatures earlier in the history of the site. Spinelli and Underwood (2005) used depth profiles of reaction progress in mixed layer illite/smectite (I/S) to constrain the thermal history of the crust at Site 1173. Because the smectite-to-illite reaction progress is controlled by reaction kinetics and thermal history, depth profiles of changing I/S composition can be used with sediment accumulation rates to constrain thermal history (e.g., Huang et al., 1993; Elliott and Matisoff, 1996). Spinelli and Underwood (2005) found thermal histories most consistent with the clay-mineral observations closely follow the expectation for conductively cooling lithosphere for the first 10 m.y. of the site, with anomalously high temperatures only occurring in the most recent 5 m.y. Models that assume anomalously high temperatures early in the history of the site (e.g., Yamano et al., 2003) result in I/S reaction progress more advanced than what is observed (Spinelli and Underwood, 2005). Similarly, studies using vitrinite reflectance to constrain the thermal histories at ODP Sites 808 and 1174 (Fig. 2; ∼14–15 km landward of Site 1173, ∼2–3 km landward of the deformation front) indicate that anomalous heating of the sites is likely restricted to the most recent ∼2 m.y. (Horsfield et al., 2006; Ondrak et al., 2009). However, these solutions are not unique; many heat-flux scenarios yield results consistent with the observed depth profiles of I/S, provided that they do not result in excess heating early in the history of the accumulating sediment. For example, using a constant basal (bottom of the sediment column) heat flux at the present-day observed value (Steurer and Underwood, 2003) or a basal heat flux that is initially consistent with a plate-cooling model until it reaches the present-day observed value (Saffer and McKiernan, 2009) both result in modeled I/S profiles consistent with the observations. Thus, the cause of the recent heating of these sites remains enigmatic.

The short duration of anomalous heating for sediments in the Nankai Trough causes the sediments with anomalously advanced smectite-to-illite reaction progress to be spatially limited. At Site 1173, the amount of I/S reaction is small (<15% reduction in % smectite in the mixed layer phase) and only occurs in sediments buried to depths >500 mbsf. Other ODP and Integrated Ocean Drilling Program (IODP) sites in the Nankai Trough are on older oceanic lithosphere (i.e., farther from the fossil spreading center). At these colder sites, smectite-to-illite reaction progress is expected to be more limited (i.e., deeper and/or less advanced). For example, heat-flux observations around IODP Sites C0011 and C0012 (Fig. 2) require an addition of heat to the basement aquifer near the trench (Spinelli, 2014), likely mined from the subducted crust (Harris et al., 2013). However, on this >20 Ma lithosphere, the highest temperatures in the sediments (i.e., at the base of the sediment column) at Sites C0011 and C0012 are projected to be 79 °C and 64 °C, respectively (Expedition 333 Scientists, 2012a; Expedition 333 Scientists, 2012b). Consistent with those results, the clay-mineral assemblages at these sites show no significant smectite-to-illite reaction progress (Underwood and Guo, 2013; Underwood and Guo, 2017).

The clay-mineral assemblage at Site 1173 is not consistent with hypothesized substantial and widespread thermal rejuvenation of the lithosphere by postspreading volcanism. To generate the present-day heat flux in the trench near Site 1173 by conductive cooling, after a volcanically produced thermal rejuvenation, would require anomalously high heat flux into the sediment column for the past ∼5.5–7.5 m.y. (Fig. 3) (Yamano et al., 2003). Such a thermal history would advance I/S alteration in the sediment >500 mbsf more than observed (Fig. 4B). Thermal anomalies (and the associated acceleration of sediment alteration) resulting from postspreading volcanism may be localized around the young seamounts (Fig. 2); sites 808, 1173, and 1174 are likely unaffected by such thermal rejuvenation.

Finally, we note that alteration of smectite prior to subduction reduces the amount of smectite (and thus mineral-bound water) delivered into the subduction zone with the bulk sediment (Saffer et al., 2008). The limited amount of smectite-to-illite reaction progress deep in the sediment column at Site 1173 (resulting from the short duration of anomalous heating) allows more mineral-bound water to be subducted. If localized areas of more advanced I/S alteration exist around postspreading volcanic features, those could also be areas where less mineral-bound water is subducted.

As sediment on a subducting plate is underthrust beneath a margin wedge, its temperature history continues to impact diagenetic reaction progress. As a result, thermal anomalies within subduction zones affect the spatial distribution of reaction progress in subducting sediment. The distribution of common diagenetic reactions in underthrust sediment is important for the mechanical behavior of a subduction zone. Such reactions can: (1) release mineral-bound water to increase pore-fluid pressure (e.g., Moore and Vrolijk, 1992; Saffer and Bekins, 1998; Lauer and Saffer, 2015), weakening the megathrust and wall rock; and (2) change the sediment’s frictional properties (e.g., Saffer and Marone, 2003; Ikari et al., 2009). Additionally, the fluids from mineral dehydration reactions alter pore-fluid chemistry, providing useful tracers for fluid flow and affecting geochemical and biological processes (e.g., Kastner et al., 1991; Saffer and Bekins, 1998; Hensen et al., 2004; Solomon and Kastner, 2012).

The notion of linking the cooling of a subduction zone by hydrothermal circulation to spatial variability in fluid sources from mineral dehydration (and finally with spatial variability in excess fluid pressures and slip behavior on a subduction-zone megathrust) is most developed for the Costa Rican margin. Subduction-zone temperatures throughout the margin, from the Nicaragua and Costa Rica border in the northwest to the Osa Peninsula in the southeast are affected by hydrothermal circulation in the subducting crust (Harris et al., 2010). However, the nature of the fluid circulation varies between different segments of the margin. In the southeast, the hydrothermal circulation does not generate a regionally extensive heat-flow anomaly and is modeled as if it is insulated. Heat is redistributed laterally in the aquifer but not lost to the ocean. In the trench, this is manifested by anomalously high heat flux relative to the expectation with no hydrothermal circulation (similar to the Nankai margin) (Harris et al., 2010). In the northwest, the hydrothermal circulation is ventilated; in addition to redistributing heat, fluid circulation extracts heat from the oceanic crust by advecting it to the ocean (Harris et al., 2010). In this region, heat flux in the trench is extremely low (∼12 mW m⁻²) (Langseth and Silver, 1996; Hutnak et al., 2007). As a result, the subduction zone in the southeast is substantially hotter than that in the northwest. The effect on the thermal history of subducting sediment is dramatic, with subducting sediment in the southeast warming to 150 °C at ∼50–55 km landward of the trench but subducting sediment in the northwest reaching the same temperature at ∼70–75 km landward of the trench (Harris et al., 2010).

The differences in the thermal histories for subducting sediment along the margin result in substantial differences in the expected distribution of sediment alteration and dehydration. Spinelli and Saffer (2004) first examined the potential for differences in the distribution of fluid sources from mineral dehydration for cool versus warm segments of the margin, focusing primarily on how differences in heat content of the plate entering the subduction zone can affect the thermal history of subducting sediment. Subsequent thermal models for the margin (Harris et al., 2010) benefit from understanding the ability of the subducting aquifer to act as a heat exchanger and consider differences between ventilated and insulated hydrothermal circulation. Lauer et al. (2017) use this new generation of thermal models (Harris et al., 2010) and reaction kinetics for the smectite-to-illite transition (Pytte and Reynolds, 1989) to predict that the peak in fluid production from dehydration reactions occurs ∼30 km landward of the trench in the southeast and ∼50 km landward of the trench in the northwest.

The difference in fluid source locations can affect distributions of excess fluid pressure in the margin, as well as the likely pathways for drainage of fluid from the margin. The deeper fluid sources (in the cooler, northwestern portion of the margin) are more likely to flow up through the margin wedge along splay faults than the shallower fluid sources (in the warmer, southeastern portion of the margin), which are more likely to flow up the plate-boundary fault to the trench (Lauer and Saffer, 2015). Additionally, the locations of these fluid sources may help control the downdip location of nucleation for large earthquakes on the plate interface (Lauer et al., 2017).

The subduction zone in southwestern Mexico provides another example of the consequence of hydrothermally cooled underthrusting sediment delivering mineralogically bound water farther landward. These effects are most prominent in Guerrero, where the subducting slab is subhorizontal (i.e., flat-slab subduction occurs) from ∼150–290 km landward of the trench (Pardo and Suárez, 1995). In this portion of the margin, tectonic tremor and slow slip are observed on the plate interface (e.g., Kostoglodov et al., 2003; Husker et al., 2012); several models have proposed links between the generation of these events and excess fluid pressures (e.g., Shelly et al., 2006; Ide et al., 2007). Also in this region, Song et al. (2009) identify ultraslow seismic-velocity zones in the subducting crust, attributed to excess fluid pressures. Perry et al. (2016) generate thermal models of the margin and link them with petrologic models to predict the distribution of subducting slab alteration and dehydration. Models that include hydrothermal circulation predict regions of focused fluid release from subducting sediment in the areas of tremor, slow slip, and the ultraslow layer (Perry et al., 2016). In models without hydrothermal circulation, subducting sediment heats and dehydrates substantially before reaching the flat section of the slab, leaving very little fluid for release in the areas of tremor, slow slip, and the ultraslow layer (most dehydration of the subducting basement rock is predicted to occur deeper than the flat-slab section) (Perry et al., 2016).

Hydrothermally cooling a subduction zone can affect the pressure-temperature (P-T) history of a subducting slab throughout much of its path. The locations at which a subducting slab experiences substantial alteration is controlled by its P-T path (e.g., Hacker et al., 2003; Hacker, 2008; van Keken et al., 2011). Thus, temperature anomalies caused by hydrothermal circulation can affect the locations of slab alteration and dehydration.

Locations of slab alteration can provide additional constraints on the temperature distribution in the subduction zone. In subducting crust, the blueschist-to-eclogite transition results in large changes in P-wave velocity (Hacker et al., 2003); thus, this petrologic transition in subducting crust may be observed with seismic imaging (e.g., Hori, 1990; Abers and Sarker, 1996; Bostock et al., 2002; Nicholson et al., 2005; Rondenay et al., 2008; Abers et al., 2009; McGary et al., 2014). Constraining the location of a P-T–controlled transition deep in the system (i.e., far removed from the seafloor and land surface) can provide a useful complement to surface heat-flux observations (e.g., Spinelli and Wang, 2009; Cozzens and Spinelli, 2012).

For the Nankai margin, Spinelli and Wang (2009) use both surface heat-flux data and the seismically inferred location of the basalt-to-eclogite transition in the subducting slab to constrain thermal models for the system. In the Nankai margin, the seismically inferred basalt-to-eclogite transition (Hori, 1990) occurs ∼35 km farther landward (and >10 km deeper) than predicted in thermal models without hydrothermal circulation (Supplemental Fig. S1¹) (Spinelli and Wang, 2009). To sustain hydrothermal circulation sufficient to cool the subducting slab such that the modeled location of the basalt-to-eclogite transition is consistent with the interpretation of seismic observations, high permeability (≥10^–10 m²) in the basement aquifer must be maintained for ∼100 km into the subduction zone (Spinelli and Wang, 2009).

Fluid circulation in subducting crust generates larger thermal anomalies in hotter subduction zones (Rotman and Spinelli, 2013). In cooler subduction zones, the consequences of that circulation in affecting slab alteration are somewhat muted. Relative to the Nankai margin (above), the subduction zone along the Middle America Trench is cooler due to a much steeper slab dip and older lithosphere subducting (by ca. 5–10 Ma). Between the Costa Rica and Nicaragua segments of the subduction zone, there are differences in the degree of mantle hydration and volcanic arc geochemistry. Seismic tomography and attenuation studies indicate that the mantle beneath Nicaragua is more hydrated than that beneath Costa Rica (Rychert et al., 2008; Syracuse et al., 2008; Dinc et al., 2011). Arc geochemistry data (Ba/La and Ba/Th concentrations) are consistent with more slab dehydration under Nicaragua than Costa Rica (Patino et al., 2000; Carr et al., 2004). In order to investigate a possible explanation for these differences between the Nicaragua and Costa Rica segments, Rosas et al. (2016) developed thermal models for the margin and predicted the distribution of slab dehydration. In this system, nearly all of the slab dehydration is predicted to occur beneath the flowing portion of the mantle wedge (top of the slab at >70 km depth). In a transect through northwestern Costa Rica, hydrothermal circulation cools the subducting crust beneath the mantle wedge by ∼40 °C (Rosas et al., 2016). However, at the top of the subducting crust, the modeled location of the basalt-to-eclogite transition is essentially unaffected by hydrothermal circulation (Supplemental Fig. S1 [^{footnote 1}]). The basalt-to-eclogite transition at the top of the slab occurs shortly after it contacts the hot mantle wedge whether or not hydrothermal circulation cools the system (Rosas et al., 2016). However, the cooling of the system by hydrothermal circulation causes the basalt-to-eclogite transition at the base of the subducting crust to shift ∼15 km deeper than expected without the hydrothermal circulation. Nonetheless, the impact on the distribution of slab dehydration is expected to be small (Rosas et al., 2016). In this colder subduction zone (relative to Nankai and Cascadia), it is the mobile mantle wedge in contact with the subducting slab, not hydrothermal circulation in the slab, that controls the distribution of slab dehydration.

The cooling of the subducting slab caused by hydrothermal circulation has potential implications for slab melting, particularly in the hottest subduction zones (Rotman and Spinelli, 2013). At the expected P-T conditions in most subducting slabs at subarc depths (e.g., Syracuse et al., 2010), fluid-absent melting is unlikely (Hermann and Spandler, 2008; Bouilhol et al., 2015); and, in hot subduction zones, the subducting crust (including sediment, basalt, and gabbro) is likely dehydrated before it reaches subarc depths (van Keken et al., 2011). However, it may be possible to melt previously dehydrated subducting crust if the temperature is very high, especially if fluid sourced from dehydration deeper inside the slab (e.g., from dehydration of serpentinized oceanic mantle) migrates up into the overlying crust (e.g., Hermann et al., 2006; Cooper et al., 2012; Bouilhol et al., 2015; Walowski et al., 2015).

The geochemical signature of lavas from volcanoes in south-central Chile indicates their magmatic source includes a signature of subducting sediment, but not subducting basalt (Kilian and Behrmann, 2003; Shinjoe et al., 2013), providing a constraint on the P-T conditions of the subducting slab. Near the Chile Rise-Trench triple junction, the subducting Nazca plate is <2 m.y. old (Tebbens et al., 1997); as a result, the subduction zone there is expected to be very hot (e.g., Spinelli et al., 2015). Thermal models that do not include the effects of fluid circulation predict melting of both the subducting sediment and the subducting basalt (Spinelli et al., 2015), inconsistent with the geochemistry of the arc lavas (Kilian and Behrmann, 2003; Shinjoe et al., 2013). In contrast, thermal models that include the effects of fluid circulation predict a cooler subducting slab, such that subducted sediment is melted, but subducted basalt is not (Spinelli et al., 2015).

By cooling the hottest portions of the margin (i.e., those closest to the triple junction) the most, hydrothermal circulation reduces the along-arc variability in subduction-zone temperatures relative to what is expected in the absence of hydrothermal circulation. This can help to explain a lack of variation in arc lava composition (Shinjoe et al., 2013) over a region with a large along-margin variation in subducting plate age (<2–14 Ma), which, in the absence of hydrothermal circulation, is expected to result in a substantial variation in the subduction-zone thermal structure (Rotman and Spinelli, 2014).

In central Washington State, a transect (Fig. 5) approximately at the latitude of Mount Rainier has been a focus for both seismic and magnetotelluric studies that infer locations of alteration and dehydration within the subducting slab (Abers et al., 2009; McGary et al., 2014). The seismically inferred onset of eclogitization in the subducting crust occurs at ∼45 km depth (∼220 km landward of the trench) (McGary et al., 2014). However, a velocity contrast across the oceanic Moho persists to more than 80 km depth, indicating that some low-velocity (probably hydrated) crustal material persists in the subducting slab to this depth (McGary, 2013). In addition, low electrical resistivities in the mantle wedge (constrained by magnetotelluric data), consistent with ∼2% melt by volume (McGary et al., 2014), are likely tied to fluid sources from the subducting slab.

We use many of the ideas presented above to reassess the thermal state of the Cascadia margin along this transect. The models presented in this study are an advance over the previous models for the margin (Cozzens and Spinelli, 2012) in several ways. First, the geometry is derived from updated data for the margin. Second, the use of a temperature- and stress-dependent mantle-wedge rheology for calculating mantle-wedge flow (Cozzens and Spinelli [2012] used an isoviscous mantle wedge) allows for more accurate temperature estimates beneath the mantle wedge. This update is critical for allowing us to estimate the distribution of deep slab alteration and dehydration.

Our model transect extends from 100 km seaward of the trench to 400 km landward of the trench (Fig. 5). The depths to the top of the subducting slab are extracted from McCrory et al. (2012) and McGary et al. (2014). In the upper plate, the distribution of the Siletzia terrane (with thermal properties distinct from continental rocks) is from Wells et al. (2014). Continental Moho depths are from Bostock et al. (2002) and McGary et al. (2014). Thermal properties are summarized in Table 1. The convergence rate between the Juan de Fuca plate and the coastal blocks along the Mount Rainier transect is ∼36 mm yr^–1; the angle between the convergence vector and the transect is 20° (McCaffrey et al., 2007). Thus, the convergence rate in the direction of the transect is ∼32 mm yr^–1. At the seaward boundary (100 km seaward of the trench), we set temperatures for the incoming Juan de Fuca plate consistent with conductively cooled 5.6-m.y.-old lithosphere (Wang and Davis, 1992); at the given convergence rate, this yields 9-m.y.-old lithosphere at the trench (Wilson, 1993). We calculate the landward geotherm using the global average for backarc surface heat flux at the land surface, 80 mW m^–2 (Currie and Hyndman, 2006), and the assigned thermal conductivity and heat production values for the continental lithosphere. The top of the model is set to 0 °C; to define the temperature at the bottom of the model, we use a mantle potential temperature of 1420 °C that increases 0.5 °C per kilometer depth (Stein and Stein, 1994; Syracuse et al., 2010). The basement aquifer is assumed to be 600 m thick, and its permeability value seaward of the trench is set to k = 10^–9 m² (e.g., Becker and Davis, 2004) resulting in an average Nusselt number seaward of the trench of ∼6800; the permeability is assumed to decay logarithmically with increasing depth at a rate of log(k) = 0.07 per kilometer. We illustrate the sensitivity of the modeled surface heat-flux distribution and temperatures in the subducting crust to variations in aquifer permeability in Supplemental Figure S2 (^{footnote 1}).

We compare new thermal model results with surface heat-flow observations and locations of slab alteration and dehydration. We use the new thermal model with no hydrothermal circulation as a reference in order to examine the effects of hydrothermal circulation. The largest differences in modeled surface heat flux between the reference simulation and the simulation including the effects of hydrothermal circulation occur on the accretionary prism (Fig. 6). From ∼10–120 km landward of the trench, modeled heat flow in the reference simulation is ≥20 mW m^–2 higher than that for the simulation including the effects of hydrothermal circulation (Fig. 6). The simulation with hydrothermal circulation is more consistent with the observations. However, along much of the transect, the scatter in the surface heat-flux observations (e.g., Phrampus et al., 2017) spans the predicted differences in surface heat flux between the two models, thus making it difficult to definitively identify a preferred thermal model for this system based on the surface heat-flux distribution alone.

The presence or absence of hydrothermal circulation in the subducting crust has important implications for the P-T path of the subducting slab and the distribution of slab dehydration. In the simulation with hydrothermal circulation, the modeled P-T path for the subducting crust is ∼100 °C cooler than the reference simulation from the trench to ∼300 km landward of the trench (i.e., where the flowing mantle wedge impinges on the top of the subducting slab) (Fig. 7A). In the reference simulation, the peak in dehydration of the oceanic crust is predicted to occur at ∼250 km landward of the trench (top of the slab at 51 km depth) (Fig. 7B). In the simulation with hydrothermal circulation, peak oceanic crustal dehydration is predicted at ∼280 km landward of the trench (top of the slab at 63 km depth) (Fig. 7B). The seismic character of the subducting crust and oceanic Moho does not change appreciably at 250 km landward of the trench (Fig. 7B). Slightly landward of 280 km, the velocity contrast across the oceanic Moho (a likely indicator of low-velocity hydrated material in the oceanic crust) (McGary, 2013) begins to break down (Fig. 7B).

We use the modeled subduction-zone temperatures in combination with phase diagrams for subducting lithologies to estimate the distribution of fluid sources from dehydration reactions in the subducting slab following the approach of Wada et al. (2012). For these calculations, the subducting slab is composed of a 300-m-thick section of underthrusting sediment, 2 km of basalt, 5 km of gabbro, and 4 km of peridotite (we assume the oceanic mantle >4 km below the oceanic Moho is dry). Based on estimates derived from seismic velocities in the Juan de Fuca plate seaward of the trench (Horning et al., 2016), we set the initial water contents at 9.2 wt% for sediment, 4 wt% for basalt, 1.1 wt% for gabbro, and 0.7 wt% for peridotite. We use phase relations for each lithology calculated using Perple_X for H₂O saturated cases (Connolly, 2009). Bulk compositions for the lithologies are summarized in Supplemental Table S1 (^{footnote 1}). Our four calculated phase diagrams are similar to those presented in van Keken et al. (2011) and Wada et al. (2012) but with a lower maximum water content for peridotite based on the seismically derived estimate for the Juan de Fuca plate from Horning et al. (2016). For each 100 m × 100 m element in the subducting slab, we determine the pressure, temperature, and water content. The motion of the subducting plate controls the rate at which each lithology traverses its pressure-temperature path and is dehydrated. Upon a decrease in the mineralogically bound water content, we assume that fluid is completely released.

In hot subduction zones (e.g., Cascadia), melting in the mantle wedge and arc volcanism are often linked with dehydration of the subducting oceanic mantle (e.g., Cooper et al., 2012; Bouilhol et al., 2015). For example, in southern Cascadia, dehydration of subducting serpentinized mantle could provide a fluid source to drive melting in the upper slab (which had dehydrated at shallower depth) (Hermann et al., 2006; Walowski et al., 2015). Consistent with these ideas, our thermal models predict dehydration of the subducting basalt and gabbro primarily beneath the stagnant corner of the mantle wedge (i.e., seaward of the volcanic arc and the region of low electrical resistivity in the mantle wedge) (Fig. 8). In the reference simulation, the subducting mantle also dehydrates seaward of the region of inferred melting in the mantle wedge (Fig. 8). In the simulation with hydrothermal circulation, dehydration of the subducting mantle occurs under the region of low resistivity in the mantle wedge (∼20 km farther landward from the zone of mantle dehydration in the reference simulation). Thus, the cooling of the slab by hydrothermal circulation allows mineral-bound fluid in the subducting mantle to be transported deeper into the subduction zone, resulting in the location of subducting mantle dehydration (and, in turn, the flux of those fluids to the overlying subducting crust; e.g., Hermann et al., 2006; Walowski et al., 2015) that is more consistent with seismic and magnetotelluric observations and/or interpretations (McGary, 2013; McGary et al., 2014).

The relatively low mineral-bound water content in the warm subducting oceanic mantle along the Mount Rainier transect (Horning et al., 2016) limits the volume of water delivered below the mantle wedge. The initial mineral-bound water content for the upper 4 km of oceanic mantle used in our models, 0.7 wt% (from interpretation of seismic velocities in the Juan de Fuca plate near our study area; Horning et al., 2016), is quite low relative to estimates for some other subduction zones. Contreras-Reyes et al. (2007) estimate up to 3 wt% mineral-bound water in the upper mantle off Chile; Van Avendonk et al. (2011) estimate 3.5 wt% mineral-bound water in the upper mantle off Nicaragua. Thus, the potential volume flux of fluid released from the subducting oceanic mantle in the Mount Rainier transect is limited relative to other subduction zones because the incoming oceanic mantle is likely fairly dry (Horning et al., 2016); however, the predicted location of that fluid release in the Mount Rainer transect is below the region of low resistivity in the mantle wedge (Fig. 8). If we assume a higher initial water content in the subducting mantle (e.g., 2 wt%), the additional fluid release (i.e., dehydration of the oceanic mantle from 2 wt% to 0.7 wt% H₂O) would occur updip from the mantle dehydration shown in Figure 8 (cf. Walowski et al., 2015), while the location of the downdip end of subducting mantle dehydration would remain the same. One potential additional source of fluid to the flowing portion of the mantle wedge is serpentinized mantle from the wedge tip (e.g., Fig. 7B) being incorporated into a mélange and dragged to greater depth (e.g., Bebout, 1991; Bebout and Barton, 2002; Guillot et al., 2015). However, the expected P-T path for such material along the plate interface is much different from that for the subducting oceanic mantle (Supplemental Fig. S3 [^{footnote 1}]). Below the flowing mantle wedge, the top of the slab is heated much more rapidly (and to higher temperatures) than material below the oceanic Moho. For example, in our model with hydrothermal circulation, where the plate interface is at 70 km depth (i.e., the seaward limit for the flowing portion of the mantle wedge in this model), the predicted temperatures at the top of the slab and at the oceanic Moho are 695 °C and 605 °C, respectively (Supplemental Fig. S3 [^{footnote 1}]). In this case, any serpentinite along the plate interface would be expected to dehydrate near where the slab encounters the flowing mantle wedge; the cooler and more slowly heated subducting oceanic mantle could transport mineral-bound fluid farther landward.

A recent study highlighted an apparent disagreement in shallow (<2 GPa) subduction-zone P-T conditions between some numerical models and estimates from exhumed high- and ultrahigh-pressure metamorphic rocks. Penniston-Dorland et al. (2015) show that the temperature conditions along the plate interface estimated from exhumed metamorphic rocks are generally much hotter than those predicted by numerical models in two summary studies conducted by Gerya et al. (2002) and Syracuse et al. (2010). Penniston-Dorland et al. (2015) reason that because the models used in Gerya et al. (2002) and Syracuse et al. (2010) do not include frictional heating along the subduction thrust, they yield anomalously low temperatures (particularly for pressures <2 GPa). Further, Penniston-Dorland et al. (2015) contend that the P-T conditions estimated from exhumed metamorphic rocks describe the range of typically expected P-T conditions for subduction zones, and models that predict lower temperatures “misrepresent nature.” A counter argument made by Abers et al. (2017) is that the suite of exhumed metamorphic rocks may be biased toward hotter, more buoyant, forearcs; thus, the available rock record does not describe the full range of expected subduction-zone P-T conditions (i.e., not sampling cooler subduction zones as represented by this class of thermal models; Syracuse et al., 2010). Further, the numerical thermal models that incorporate arbitrarily high frictional heating to satisfy the P-T conditions from the exhumed rock record predict unrealistically high surface heat flow in the forearc and violate other geophysical observations, such as the relatively low seismic attenuation commonly observed in the mantle-wedge corner, indicating that such P-T conditions are not representative of present-day average subduction-zone temperature conditions (van Keken et al., 2018).

In the context of these ongoing discussions, it is interesting to consider how the results of thermal models that include the effects of both frictional heating and hydrothermal circulation compare with the P-T conditions estimated from exhumed metamorphic rocks. For systems that have been modeled with the effects of frictional heating and hydrothermal circulation, the temperature gradient in the subducting crust from 0 to 2 GPa is ∼10 °C km^–1 (i.e., through the heart of the data from exhumed rocks) for the Nankai margin (Spinelli and Wang, 2009), the Mount Rainier transect through the Cascadia margin (this study), and the flat-slab portion of the Mexican margin (Perry et al., 2016). Model results from southern Chile (Spinelli et al., 2015) and portions of the Mexican margin with a steeply dipping slab (Perry et al., 2016) predict temperature gradients in the subducting crust from 0 to 2 GPa of ∼5 °C km^–1 (i.e., cooler than the estimates from exhumed metamorphic rocks). In this volume, Harris et al. (2017) summarize the need to consider the competing effects of frictional heating and hydrothermal circulation in estimating subduction-zone temperatures; frictional heating increases the plate-boundary temperature (i.e., the “subduction geotherm”), while hydrothermal circulation decreases the subduction geotherm. Several studies have suggested that hydrothermal circulation transports a significant portion of the heat generated by frictional sliding toward the deformation front, thereby confounding attempts to use surface heat-flux observations to constrain the amount of frictional heating (Rotman and Spinelli, 2014; Wang et al., 2015; Hass and Harris, 2016). Harris et al. (2017, p. 1441) note: “in order to use heat flow observations above a subduction thrust to constrain the amount of frictional heating on the fault, independent knowledge of the vigor of hydrothermal circulation or the lack of hydrothermal circulation is required.” Thus, reconciling P-T estimates from the rock record and thermal models may not be as straightforward as simply increasing the effective coefficient of friction.

Hydrothermally cooling a subduction zone affects the hydrogeology of the system, impacting both the subducting aquifer itself and the surrounding material. Thermal anomalies in subduction zones (including those caused by hydrothermal circulation) can affect the hydrogeology of the system by controlling: (1) the distribution of fluid sources, (2) permeability, and (3) fluid properties. The examples from Nankai, Costa Rica, Mexico, and Cascadia summarized above highlight how hydrothermal circulation can influence the distribution of fluid sources from mineral dehydration in subduction zones. The spatial distribution of these fluid sources depends on the initial material composition and pressure-temperature (P-T) history (e.g., Bekins et al., 1995; Saffer and Bekins, 1998; Hyndman and Peacock, 2003). Below, we discuss how thermal anomalies can affect subduction-zone hydrogeology through their control on permeability and fluid properties.

The chemical and mechanical sealing of fractures and pore space with depth common in brittle rocks (e.g., Manning and Ingebritsen, 1999) is expected to reduce the permeability of the oceanic lithosphere during subduction. In volcanic rocks from oceanic crust and from ophiolites, prograde metamorphic reactions produce minerals (e.g., actinolite, prehnite, and pumpellyite) that fill voids and fractures (e.g., Alt et al., 1986; Harper et al., 1988). Laboratory experiments demonstrate that such authigenic mineral formation decreases permeability (e.g., Morrow et al., 2001; Kato et al., 2004; Tenthorey and Fitz Gerald, 2006; Yasuhara et al., 2006). However, for the basement aquifer of the oceanic crust, variations in permeability during subduction are only loosely estimated. Using thermal models constrained by surface heat-flux data and the seismically determined position of slab eclogitization for the Nankai margin, Spinelli and Wang (2009) estimate that permeability of the oceanic crustal aquifer is ∼10^–9 m² prior to subduction and is reduced by at least an order of magnitude by 20 km depth. This permeability decrease may be gradual or focused where the subducting crust passes through metamorphic facies boundaries (Spinelli and Wang, 2009). Hydrothermal circulation could affect the evolution of permeability in subducting crust by altering its P-T path during subduction. A positive feedback may occur, in which hydrothermal circulation reduces temperatures of the subducting crust, which delays prograde metamorphic reactions and allows the permeability of the subducting crustal aquifer (and thus the vigor of fluid circulation) to be maintained farther into the subduction zone.

Similar reductions in permeability with increasing temperature are likely important outside of the subducting aquifer; these too, will be impacted by hydrothermal circulation affecting subduction-zone temperatures. Kato et al. (2004) examine the effects of temperature on the permeability of subduction-zone fault rocks (i.e., hosted in sedimentary material overlying the basaltic basement). They find an approximately one order of magnitude decrease in permeability by increasing temperature from 30 °C to 250 °C (Kato et al., 2004). The cooling of a subducting aquifer cools surrounding material, including the subduction-zone megathrust. In hot subduction zones, the plate interface may be cooled by >100 °C (Rotman and Spinelli, 2013). This could result in higher fault-zone permeability than would be expected in the absence of hydrothermal circulation. Thus, hydrothermal circulation in the basement aquifer of the subducting crust can also affect the hydrogeology of broader regions of subduction zones (e.g., drainage of fluid along the plate interface and/or through the margin wedge; Lauer and Saffer, 2015).

Thermal anomalies generated by hydrothermal circulation in subducting crust affect fluid properties that control the ease with which fluid flows through a material. At the P-T conditions in the shallow portion of a subduction zone, temperature variations <100 °C have a large effect on fluid viscosity. For example, at 40 MPa, increasing temperature from 2 °C to 100 °C results in a >80% reduction in viscosity (Harvey et al., 2014). Relative to a system with no hydrothermal circulation, in a hydrothermally cooled subducting aquifer, fluid viscosity will be anomalously high and permeability will be anomalously high. The competing effects of high viscosity (decreasing hydraulic conductivity) and enhanced permeability (increasing hydraulic conductivity) will affect the vigor of fluid circulation in the system.

Spinelli and Saffer (2007) examined how thermally induced along-strike variations in fluid properties and permeability in a subduction zone could cause an along-strike component of fluid flow for water draining from the system updip along the plate-boundary fault. In the coldest systems examined, the anomalous viscosities dominate, driving fluid flow from the cold side of the system toward the hotter side (Spinelli and Saffer, 2007). Similarly, in subducting aquifers, anomalously high fluid viscosities will be manifested more often in cold systems due to the large variations in viscosity at low temperatures.

Thermal anomalies generated by hydrothermal circulation in the basement aquifer of subducting crust affect subduction-zone hydrogeology and the alteration of material approaching and within subduction zones. Vigorous fluid circulation redistributes heat within subducting crust, resulting in anomalously high temperatures in the crust near the trench and anomalously low temperatures farther into the subduction zone (relative to expected temperatures in the absence of hydrothermal circulation) (Fig. 1; Harris et al., 2017 and references within). The anomalously high temperatures seaward of the trench accelerate kinetically controlled diagenetic reactions in sediment prior to its entering the subduction zone. This affects both the amount of mineral-bound fluid delivered to the subduction zone (e.g., Saffer et al., 2008) and the frictional properties of the sediment (e.g., Ikari et al., 2011). Hydrothermal circulation in subducting crust causes subducting sediment to be cooler than expected in the absence of the circulation. As a result, the alteration of subducting sediments is delayed, occurring farther landward than it would in a hotter subduction zone (e.g., Spinelli and Saffer, 2004; Perry et al., 2016; Lauer et al., 2017). By shifting fluid sources from dehydration reactions farther landward, hydrothermal circulation in a subducting aquifer can affect the hydrogeology of the overlying sediments and margin wedge; fluid sources farther landward in the subduction zone are more likely to drain through the overlying wedge (Lauer et al., 2017). The cooling of the slab by lateral heat exchange in the basement aquifer in the shallow (less than ∼20-km-depth) subduction zone affects the P-T path of subducting material well landward of the cessation of the fluid circulation itself. This delays alteration of the subducting slab, shifting the locations of metamorphic dehydration reactions landward (e.g., Spinelli and Wang, 2009; Cozzens and Spinelli, 2012; Rosas et al., 2016) and affecting where (and if) subducting material may melt (e.g., Spinelli et al., 2015). The thermal legacy of hydrothermal circulation in the subducting slab is substantially diminished where the slab contacts the mobile mantle wedge; after abutting the mobile mantle wedge, the magnitude of legacy thermal anomalies (i.e., those from hydrothermal circulation farther seaward) in the slab are reduced by ∼50% (e.g., Spinelli et al., 2015; Rosas et al., 2016). In addition to affecting the distribution of fluid sources, hydrothermal circulation influences the hydrogeology of subduction zones by controlling fluid properties (particularly in the coldest portions of the system, less than ∼100 °C, where fluid viscosity varies greatly with temperature) and delaying metamorphic reactions that tend to reduce permeability. Because all of these consequences of hydrothermal circulation in subducting crust result from the thermal anomalies generated, they are most pronounced in the hottest subduction zones where the lateral heat exchange in the aquifer has its largest effects. However, the above inferences do not imply that vigorous hydrothermal circulation within the subducted crust must be occurring everywhere along warm-slab subduction zones. Permeability of the subducting crust is controlled by fractures, and the development of the fracture system depends on the deformation history of the oceanic plate and cannot be uniform along strike. It is expected that some segments of these subduction zones do not host vigorous circulation in the subducted crust. Conversely, it is conceivable that under very special conditions, extremely high permeability may be locally present in cold slabs and give rise to hydrothermal circulation in the subducted crust.

Spinelli was supported by National Science Foundation (NSF) grant OCE-1551587. Harris was supported by NSF grants OCE 1458211 and OCE 1355878. We thank Guest Associate Editor Gray Bebout, Patrick Fulton, and an anonymous reviewer for providing constructive comments that improved this study.