Fluid-related deformation processes at the up- and downdip limits of the subduction thrust seismogenic zone: What do the rocks tell us?

The subduction zone setting involves oceanic lithosphere, including crystalline rocks in the upper mantle and oceanic crust, covered by marine and trench-fill sediments, thrust below a hanging wall composed of a forearc of sedimentary rock assemblages seaward of an oceanic or continental arc. This setting creates a plate-boundary interface that is generally not a simple discrete plane with a homogeneous composition; rather, the subduction megathrust is a mélange of different rock compositions with a range of rheological properties that evolve with progressive deformation and metamorphism (Shreve and Cloos, 1986; Fisher and Byrne, 1987; von Huene and Scholl, 1991; Bebout and Barton, 2002; Bachmann et al., 2009; Fagereng and Sibson, 2010; Angiboust et al., 2011; Rowe et al. 2013). During this time-progressive evolution, fluid-saturated oceanic crust and its sedimentary cover dehydrate in response to increasing pressure (P) and temperature (T). The release, migration, and absorption of these fluids have a significant effect on fault deformation.

Hubbert and Rubey (1959) established the importance of fluid pressure in lowering effective normal stress, allowing slip along very low-angle thrusts, to explain how a gently dipping subduction thrust remains active. However, low effective normal stress is not the only mechanical effect of fluids along faults (see, for example, reviews by Etheridge et al., 1983; Hickman et al., 1995; Faulkner et al., 2010). Fluids also lead to reaction weakening through facilitating growth of new minerals that may be weaker or finer than reactants (Wintsch et al., 1995), and they enable fluid-assisted diffusion and accommodation of strain by pressure solution creep at low driving stresses (Brodie and Rutter, 1987; Gratier et al., 2011; Fagereng and den Hartog, 2017), particularly in the presence of fine-grained reaction products (White and Knipe, 1978).

Deformation of the accretionary prism demonstrates the interplay of strain accumulation and release along the mega-thrust and the overlying accreted sedimentary rocks (Wang and Hu, 2006; Kimura et al., 2007; Hu and Wang, 2008; Cubas et al., 2013). On larger time scales and assuming a steady influx of materials, the accretionary prism deforms in an attempt to maintain a critical taper angle, a mechanical equilibrium where the sum of the accretionary prism surface slope and the dip of the subduction thrust interface is a function of the relative strength of the fault and the prism materials (Davis et al., 1983; Dahlen, 1984). At smaller scales, subducting and accreted sediments record the style and conditions of deformation, leading to modifications in their strength and geometry that affect further deformation (Karig and Sharman, 1975; Cowan, 1985; Sample and Moore, 1987; Fagereng and Sibson, 2010; Kimura et al., 2012a; Rowe et al., 2013). Here, we summarize some recent work investigating the progressive deformation of subducting sediments and oceanic crust, and the interplay between deformation and fluid production and migration. We focus on the fault-rock scale, but we emphasize how mineral-scale deformation and dehydration processes are reflected in the larger-scale geometry of subduction margins.

Subduction megathrusts are extremely weak relative to most rocks, with shear strengths of no more than tens of megapascals, as inferred from the small taper angles of overlying accretionary wedges (Davis et al., 1983; Dahlen, 1990; Fagereng, 2011a; Cubas et al., 2013), drilling in the shallow parts of décollements (Moran et al., 1993; Housen et al., 1996; Chester et al., 2013), low surface heat flow (Wang et al., 1995a), exhumed ancient examples (Byrne and Fisher, 1990; Fagereng et al., 2010), the distribution and nature of aftershock focal mechanisms (Magee and Zoback, 1993; Hasegawa et al., 2011), and geophysical measurements such as high Vp/Vs ratios indicating high fluid content (Eberhart-Phillips and Reyners, 1999; Kodaira et al., 2002; Audet et al., 2009; Bassett et al., 2014). Although there are arguments for low-friction clay minerals coating the fault plane and leading to low shear strength (Brown et al., 2003; Saffer and Marone, 2003; Faulkner et al., 2011; Ujiie et al., 2013), it is difficult to reconcile laboratory deformation data and field observations of active faults without invoking fluid pressures that are significantly above hydrostatic, and, in many cases, fluid pressure must be near, or temporarily exceed, lithostatic pressure. This statement arises from the Coulomb failure criterion:

Here, shear strength (τ) is equal to cohesive strength (C₀) plus the product of the frictional coefficient (µ) and effective normal stress (σ_n′ = σ_n[1 – λ]) of the megathrust fault, where λ is the ratio of pore fluid pressure (P_f) to the vertical stress (σ_v) corrected for weight of overlying seawater. For a depth of 10 km with an average rock density of 2500 kg/m³ (Dahlen., 1984), and neglecting overlying seawater, σ_v (assumed equal to the weight of the overburden) is 245 MPa. The low clay friction values measured in the laboratory are in the range 0.1 < µ < 0.4. Assuming a sub-horizontal megathrust so that σ_n = σ_v, and negligible cohesion, only the lower part of this range in clay friction allows for the weakness (τ equal to tens of megapascals) inferred from observations, and at greater depths, even the weakest clays are not sufficiently weak. Also note that the lowest friction values involve smectite (Brown et al., 2003; Saffer and Marone, 2003; Ujiie et al., 2013), and smectite clays are only stable at T < 150 °C. Therefore, near-lithostatic or higher fluid pressures are required for frictional failure at “seismogenic” depths (T > 150 °C), where weak clay minerals are no longer stable (cf. Sibson, 2014). Talc-bearing serpentinite may be sufficiently weak even at some seismogenic conditions (Moore and Rymer, 2007); however, it is uncertain whether talc is sufficiently abundant to control strength at depth in subduction zones.

If subduction megathrusts are fluid overpressured, two conclusions follow: (1) A fluid source provides sufficient fluid production to maintain a fluid overpressure; and (2) permeability is sufficiently low to contain overpressured fluids within the fault zone (Neuzil, 1995; Saffer and Tobin, 2011; Sibson, 2014; Saffer, 2015). As P and T increase in subducting crust and sediments, this P-T path will involve a prograde, progressive increase in both P and T, leading to mechanical breakdown of porosity with increasing pressure, and a series of mineral reactions to maintain thermodynamic equilibrium. Thus, the potential fluid sources are the free and mineral-bound waters present within the subducting rocks, but the exact mechanism of fluid release and the specific pore space or mineral source will vary with depth and with the composition of the subducting sediments. Consequently, a question arises as to whether fluid release is continuous or episodic during progressive subduction. Similarly, the permeability and rheology of the subduction thrust interface evolve with depth as porosity, mineral assemblage, and fault zone structure evolve with progressive deformation and fluid loss, and they can therefore vary in both time and space (cf. Sibson, 1996).

Many contributions have reviewed the interplay between fluids and deformation in subduction zones (e.g., Moore and Vrolijk, 1992; Saffer and Tobin, 2011; Saffer and Wallace, 2015). Our emphasis here is the progressive but potentially episodic nature of fluid production and rock deformation that is recorded in exhumed subduction thrusts (Fisher and Byrne, 1987; Fisher and Brantley, 1992; Ujiie et al., 2003; Fagereng et al., 2011). We emphasize two transitions visible in the rock record: (1) a shallow transition from ductile to brittle, and (2) a deep transition from brittle to ductile (cf. Saffer and Tobin, 2011). Under the hypothesis that these structural transitions are analogous to the limits of the seismogenic zone, which require at least local frictional and velocity-weakening behavior, we can discuss the controls on brittle-ductile deformation and associated controls on the depth limits of the seismogenic zone. First, however, we review the definitions of “brittle” and “ductile.”

From a geological perspective, “brittle” refers to permanent discontinuous deformation, where failure occurs after an initial but small elastic deformation. “Ductile” is a broader term, encompassing all deformations that are macroscopically continuous, without implying a particular microscale mechanism. Thus, the term “ductile” describes structures that are continuous at the scale of observation, but may have been accommodated, at the microscale, by one or more of the following: (1) independent particulate flow, (2) pressure-dependent cataclastic flow, (3) thermally activated dislocation creep, or (4) grain-size–sensitive diffusive mass transfer (Sibson, 1977; Hadizadeh and Rutter, 1983; Knipe, 1989; Dragoni, 1993). The deep transition from brittle to ductile may in some cases be better referred to as a change from discrete frictional sliding along one or more faults to distributed, thermally activated viscous flow in shear zones, observed in rocks as a transition from cataclasites to mylonites (e.g., Handy et al., 2007). However, for consistency throughout this manuscript, and to emphasize our observational approach, we will use the terms “brittle” and “ductile” to describe discontinuous and continuous deformation, respectively, except when explicitly discussing microscale mechanisms of deformation.

The temporal, as well as spatial, scale of observation is critical when discussing transitions between brittle and ductile behavior. Ductile flow with permanent finite strain is generally a long-term process, and at short time scales, rocks tend to behave elastically, as illustrated by the transmission of seismic waves through most rocks. Thus, processes typical of brittle rocks may occur in the ductile regime on short time scales, as happens in the up- and downdip propagation of earthquake ruptures that initiated within the brittle regime (Iio et al., 2002; Kodaira et al., 2012), and aftershock distributions including areas typically deforming by long-term creep (Rolandone et al., 2004). In particular, mutually crosscutting relationships between mylonite and pseudotachylyte (Sibson, 1980; Hobbs et al., 1986; Moecher and Steltenpohl, 2011; Price et al., 2012; White, 2012), formed by aseismic creep and seismic slip, respectively, illustrate strain rate fluctuations and local, transient downward movement of the deep brittle-to-ductile (or frictional to viscous) transition.

The shallow part of the subduction megathrust hosts several changes leading to progressive lithification of seafloor sediments at increasing P-T conditions, resulting in a shift from ductile to brittle deformation (e.g., Moore and Saffer, 2001). Since “ductile” is a descriptive term meaning continuous at the scale of observation, and “brittle” is a term for deformation involving macroscopic discontinuities, ductile processes may occur even at extremely shallow depths in subduction zones. This is reflected in several reported observations of continuous structures crosscut by discontinuities in rocks that have only reached subgreenschist metamorphic conditions or less (Fig. 1A; Fisher and Byrne, 1987; Labaume et al., 1991; Maltman, 1994; Bettelli and Vannucchi, 2003; Remitti et al., 2007, 2011; Fagereng, 2011b; Hashimoto and Yamano, 2014).

At low P and T, typically prior to reaching temperatures where quartz or carbonate mobility allow chemical cementation (T < 100 °C) and at pressures too low to mechanically destroy pore space, ductile deformation commonly occurs by independent particulate flow (Maltman, 1994), generally defined as frictional grain boundary sliding with little modification of single particles (Borradaile, 1981). Because this flow mechanism requires rigid grains to rotate and flow past each other, such that the aggregate must dilate, it is promoted by elevated fluid pressures and suppressed by high confining stress. As sediments are compacted during burial, particulate flow typically accommodates volume loss by ductile deformation, creating a clay fabric parallel to bedding, as observed in the Nankai accretionary prism (Morgan and Karig, 1993).

Mélanges commonly preserve features of early shear deformation without breaking of grains, implying independent particulate flow accommodating shear as well as compaction, illustrated here by an example from the Chrystalls Beach accretionary mélange, New Zealand (Fig. 1B). A complication arises, however, in that different rock types evolve along different rheology-time paths as consolidation occurs (e.g., Maltman, 1994; Bettelli and Vannucchi, 2003; Remitti et al., 2007), depending, for instance, on the starting composition, grain size, and permeability of the sediment, and the flux and chemistry of circulating fluids. As an example, consider a simple mudstone-sandstone sequence, and assume failure is governed by the Mohr-Coulomb criterion (Fig. 2). In this example, the clay-rich mudstone develops a foliation and cohesive strength at an early stage of burial, whereas sand only gains cohesion once a cementation process occurs (e.g., Dott, 1966; Hillier and Cosgrove, 2002). There is, therefore, an early stage of burial where mudstones have cohesive strength, and sand does not (Fig. 2A; e.g., Hillier and Cosgrove, 2002), so that mudstone may be stronger than sand at low effective normal stress. Considering that clay-rich mudstone likely has a low friction coefficient, say ≤0.4 (Underwood, 2002; Saffer and Marone, 2003), while cohesionless sand may have a typical Byerlee coefficient of internal friction ≥0.6 (e.g., Jaeger and Cook, 1979), the sand becomes stronger than mudstone as frictional resistance suppresses particle flow at elevated normal stress (Fig. 2A). The result of this strength reversal is that at shallow depth in wet rocks, sand injections into mudstones occur (Fig. 1C). Because optimally oriented faults in mudstones are rarely observed, whereas flow structures in sandstones are, prelithification deformation is likely governed by independent particulate flow in sand, and common observations of boudins, folds, and asymmetric phacoids of sandstone imply that effective stress is likely low at shallow depths; a consequence of these observations is that mélange fabric characterized by independent particulate flow in sandstones likely developed at shallow (less than a few kilometers) depths.

As cementation takes place, sandstone gains cohesive strength and tensile strength and becomes stronger than mudstone at all levels of effective normal stress (Fig. 2B). In this case, hydrofractures can occur in relatively high-tensile-strength sandstone units (Figs. 1C and 1D), while mudstone may fail in shear even if differential stress is low, particularly if foliation defined by aligned clays is well oriented for shear failure. Note that the development of tensile strength in sandstones, and the associated emergence of tensile cracks, implies that a greater fluid overpressure can be contained, sustaining lower effective normal stresses (Fig. 2).

The observational evidence for a transition in properties from weak, ductile low-cohesion sand to brittle sandstone with a tensile strength includes sandstone boudins showing both a macroscopically ductile pinch-and-swell structure and superimposed tensile fractures with the same extension direction (Figs. 1C and 1D; Byrne, 1984). This progression from early, high-porosity, bulk deformation of sandstone blocks to formation of quartz-filled tension gashes in boudins necks was discussed by Hashimoto et al. (2006) in relation to the Shimanto complex of Japan, and it was ascribed to the development of a quartz-calcite-chlorite cement in the sandstone pores at a depth of ~3 km. A similar progression, mainly related to carbonate cementation, has been described in the clay-rich basal sequence of the Ligurian units of the Northern Apennines, interpreted as the frontal part of an accretionary prism (Vannucchi and Bettelli, 2002; Bettelli and Vannucchi, 2003). Here, these sequences are complexly folded and sheared and have developed a pervasive block-in-matrix fabric (Fig. 1A). The deformation began in poorly lithified sediments, as recorded by mud injection, hydroplastic fracturing, and pinch-and-swell boudinage, which strengthened progressively, as testified by pervasive sets of tensile veins (Figs. 1A and 1D) still remaining in the smectite stability field. Even if the “ductile” deformation pervasively affected the whole volume of the sedimentary sequence, sharp shear surfaces without veins, intersecting but coeval to the folded units, are present (Fig. 1A; Bettelli and Vannucchi, 2003).

Hashimoto et al. (2006) recorded thrust-sense, discrete shearing in the Shimanto complex that was inferred to begin at T > 150 °C and overprint more distributed structures accommodating layer-parallel extension, such as shear bands and boudinage. More recently, Hashimoto and Yamano (2014) determined the P-T conditions of tensile and shear veins in the Shimanto complex and found a transition in deformation style at 175 °C < T < 210 °C. They suggested this as the range where lithification of sandstone leads to a change from dominantly ductile to dominantly brittle behavior. Similarly, in a shear zone interpreted as an analogue to a shallow subduction plate interface cropping out in the Northern Apennines (the Sestola Vidiciatico tectonic unit; Vannucchi et al. 2008), a transition from diffuse deformation involving partly lithified sediments to discrete faulting characterized by calcite shear veining was reported (Labaume et al., 1991; Remitti et al., 2007), with veins both oblique (Fig. 1E) and subparallel (Fig. 1F) to clay foliation. Similar calcite shear veins were described by Dielforder et al. (2015) in a Paleogene accretionary complex of the central European Alps. These veins were inferred to have formed along bedding planes during reverse faulting in the outer wedge, documenting horizontal contraction within the wedge and precipitation at 40–70 °C (Dielforder et al., 2015). This low temperature of vein formation suggests that in a carbonate-rich environment, cementation can lead to brittle behavior at much cooler conditions than the 150 °C boundary inferred in silica-rich systems, adding complexity to the interpretation and prediction of the shallow ductile-to-brittle transition in subducting sediments.

In the Chrystalls Beach complex, Otago Schist, New Zealand, incrementally developed quartz-calcite veins coat shear surfaces within mélange (Fig. 1G; Fagereng et al., 2010, 2011). Shear discontinuities have typically developed along phyllosilicate foliations and, in places, at a high angle to associated tensile veins (Fagereng et al., 2010). Tensile veins, both in dilational step-overs and in rigid sandstone lenses, imply that, at least locally and episodically, P_f was in excess of the least compressive stress (Secor, 1965), and differential stress was less than four times the tensile strength, on the order of tens of megapascals at the most (Etheridge, 1983). Similar, incrementally developed slickenfibers have been reported in calcite shear veins also in Sestola Vidiciatico tectonic units (T < 120 °C; Labaume et al., 1991; Vannucchi et al., 2010). Ductile deformation, illustrated by diffuse pressure solution and affecting the shear veins, is recorded where deformation has occurred at T > 150 °C in the Chrystalls Beach complex (Fagereng, 2011b) and in the central European Alps (Dielforder et al., 2015).

The transition from shallow, diffuse deformation of weakly consolidated sediments to sharp discrete faults coated by shear veins cutting lithified and cemented material may coincide with the boundary from aseismic to seismic deformation, inferred to be thermally controlled at 100–150 °C (e.g., Vrolijk, 1990; Oleskevich et al., 1999). For instance, Moore and Saffer (2001) concluded that decreased fluid production beyond the smectite-illite transition, coupled to sediment cementation by clays, carbonates, zeolites, and quartz, which becomes efficient at ~150 °C, can explain the updip limit of seismicity in southwest Japan, because the resultant increases in effective normal stress and rigidity change rock friction to velocity-weakening regimes. Bangs et al. (2004) reported a seismic reflection that decreases in amplitude downdip of the inferred updip limit of the Nankai Trough seismogenic zone, which is also indicative of a decrease in fluid pressure at this depth. Moore et al. (2007) provided a further argument that because quartz is velocity weakening at T > 150 °C (Blanpied et al., 1995), and 150 °C is also the temperature where quartz mobility by pressure solution becomes efficient (Rimstidt and Barnes, 1980), quartz cements allow velocity weakening and therefore seismic behavior above 150 °C. This final point implies that the transition from unlithified to lithified sediments may coincide with the updip limit of seismicity, at least in quartz-rich subducting sediment sections.

A thermally controlled, updip limit of seismicity at 100–150 °C is, however, quite simplistic. Several lines of evidence from active subduction zones imply a more complex and diverse behavior of the shallow part of the subduction megathrust. For instance, the updip limit of seismicity in Costa Rica, as defined by aftershock events of the Mw 6.9 Quepos earthquake, occurs at temperatures cooler than the 100–150 °C temperature typically predicted (Hass and Harris, 2016). Slow slip events, which have been suggested to be indicators of unstable frictional properties (Liu and Rice, 2005) and thought to be associated with transitions from seismic to aseismic behavior (Saffer and Wallace, 2015), propagate to near the seafloor in both Costa Rica (Davis et al., 2011) and northern Hikurangi, New Zealand (Wallace et al., 2016). These near-trench regions therefore allow slow slip, and possibly propagation of seismic slip, although no earthquake nucleation is observed here. Along the Nankai Trough, very low-frequency earthquakes have been observed within a seismic low-velocity zone below the shallow accretionary prism, indicating low-stress-drop frictional failure shallower than 10 km deep (Ito and Obara, 2006; Sugioka et al., 2012). Recently, repeating slow slip events have also been reported near the Nankai Trough (Araki et al., 2017). This varied, and conditionally unstable, near-trench behavior is not consistent with a simple, thermally controlled, updip limit of unstable slip along the interface.

The 2011 Tohoku-Oki earthquake further highlighted uncertainties around the mechanical properties of the shallow megathrust interface. This earthquake involved coseismic slip of ~60 m near the trench (Fujiwara et al., 2011; Kodaira et al., 2012), demanding a change in the paradigm that seismic slip cannot occur seaward from the updip limit of seismicity (near the 150 °C geotherm), and, if earthquakes occur on discrete planes, implying the existence of discrete slip surfaces as shallow as at the seafloor. J-FAST drill core yielded a thin (<5 m) zone of scaly clays, interpreted as the plate boundary, within which seismic slip is inferred to have propagated into unlithified smectite-rich rocks (Chester et al., 2013). This observation is analogous to discrete faults cutting smectite-rich sediments in the Apennines exhumed accretionary complex (Fig. 1A). It therefore seems that the evolution from (1) unlithified to partially lithified sediments deformed dominantly by distributed soft-sediment deformation to (2) the development of discrete faults as rocks are buried and lithified is a simplified view. What the simple view does not consider is the effect of propagating shear displacement that starts in velocity-weakening, cemented rocks, and the effect of a weak, conditionally stable, clay foliation. What the Tohoku earthquake has shown is that conditionally stable, fluid-rich, and frictionally weak unlithified sediments, although incapable of nucleating seismic slip, may nevertheless slip seismically (Chester et al., 2013; Ujiie et al., 2013).

Further complexity is given by the presence of fluids rich in carbonate, allowing cementation and development of brittle, incrementally developed faults at temperatures well below 150 °C (Dielforder et al., 2015). We also note that in fine-grained, phyllosilicate-rich rocks, distributed pressure solution creep may be favored over localized brittle shear, causing ductile deformation at temperatures beyond the inferred updip limit of seismicity at 150 °C, as reported in several examples (e.g., Schwarz and Stöckhert, 1996; Fagereng, 2011b; Dielforder et al., 2015; Fagereng and den Hartog, 2017). A likely scenario is that in a geologically complex plate-boundary interface, with disrupted sandy layers in a mudstone matrix, there is a time-dependent cyclicity between interseismic, slow pressure solution creep and episodic, fast earthquake slip triggered at or near stress risers such as rigid lenses or layers of stronger rock (Collettini et al., 2011; Fagereng, 2013), or subducting seamounts (Kodaira et al., 2000; Yang et al., 2012).

The previous section discussed the onset of nucleation of brittle deformation when sediments become lithified and obtain a tensile strength, which also leads to a transition from distributed to localized deformation in exhumed subduction thrust rock assemblages. At depths beyond this transition, the megathrust is capable of nucleating seismic slip; however, rock deformation experiments indicate a likely time-dependent competition between localized frictional sliding and distributed pressure solution creep (e.g., Niemeijer and Spiers, 2006; den Hartog et al., 2012; den Hartog and Spiers, 2013; Gratier et al., 2014), and seismic style may be determined by the volumetric proportions of competent and incompetent materials (Fagereng and Sibson, 2010; Collettini et al., 2011). With analogy to the frictional-viscous transition in continental crust, this zone of possible nucleation of seismic events and dominant brittle deformation in competent rocks is commonly inferred to extend to the depth of the 350 °C isotherm (Hyndman et al., 1997), where quartz becomes velocity strengthening (Blanpied et al., 1995) and deforms efficiently by crystal plastic flow (Hirth et al., 2001).

In inferred megathrust fault rocks exhumed from temperatures exceeding 300 °C, brittle features, such as quartz veins, are pervasively overprinted by ductile deformation (Fig. 3A; Angiboust et al., 2015; Bachmann et al., 2009; Fagereng et al., 2014; Meneghini et al., 2014; Stöckhert, 2002; Wassmann and Stöckhert, 2012). The presence of veins in boudin necks (Fig. 3B), and viscously deformed veins crosscutting crystal plastic foliation (Fig. 3C), however, implies continuation of brittle deformation, at least locally and temporarily, within the dominantly ductile regime (Bachmann et al., 2009; Fagereng et al., 2014). We take the Damara belt accretionary prism as an example. These rocks contain a foliation developed at 540–560 °C and ~1 GPa (Cross et al., 2015), and they are inferred to have deformed in an accretionary prism (Kukla and Stanistreet, 1991; Fagereng et al., 2014; Meneghini et al., 2014). The foliation is schistose and defined by muscovite, biotite, chlorite, and flattened quartz aggregates, but like the clay foliation seen at much lower metamorphic grade, the foliation is subparallel to lithological layering (Figs. 3D and 3E; Fagereng et al., 2014; Hartnady, 2014). In places, dismembered bedding is still recognizable, and a greenschist-grade, chlorite-muscovite-biotite foliation has developed within both psammitic and pelitic rocks (Fig. 3D).

The point to make here is that the foliation has a range of geometrical relationships with quartz veins. On the one hand, veins are commonly subparallel to the foliation, and boudinaged (Figs. 3A, 3D, and 3E). It is therefore possible that some veins were inherited from early in the deformation sequence and formed parallel to an early clay foliation, as suggested in our discussion on the shallow ductile-to-brittle transition. Some veins may also have formed along higher-metamorphic-grade foliations, which would have represented low-tensile-strength planes (e.g., Stöckhert et al., 1999; Fagereng et al., 2010). Some veins are isoclinally folded with hinge lines subparallel to the regional stretching lineation (Fig. 3A), and these represent cross-foliation fractures that were passively folded and rotated in a dominantly ductile regime. The fact that some veins crosscut the ductile foliation, but are themselves ductilely deformed (Fig. 3C), supports the interpretation that at least some brittle fractures occurred during dominantly ductile deformation (Bachmann et al., 2009; Fagereng et al., 2014). Veins are also common in boudin necks, where competent lenses deformed by crystal plastic deformation are separated by local, vein-filled fractures (Fig. 3B). Overall, we suggest that veins still develop when rocks are subducted to, or beyond, the depth of the 350 °C isotherm, where quartz may deform plastically. This implies that the hydrofracture criteria, i.e., fluid pressure in excess of the least principal stress (Secor, 1965) and differential stress less than four times the tensile strength (Etheridge, 1983), can be locally achieved at T > 350 °C. The veins are filled with blocky quartz, implying relatively fast growth in open space, indicating that the veins are indeed records of open fluid-filled fractures (Fisher and Brantley, 1992; Bons et al., 2012), and not records of vein growth by force of crystallization (Wiltschko and Morse, 2001) or other diffusive, local, growth mechanisms not requiring open fracture. An interpretation is that veins form by episodic brittle failure within the dominantly ductile regime below the base of the brittle crust, and thus deformation cycles between fast and slow strain rates (Bachmann et al., 2009; Wassmann and Stöckhert, 2012; Fagereng et al., 2014). The veins also provide evidence for a sink for precipitation of dissolved silica, consistent with dissolution-precipitation creep as an important mineral deformation mechanism along the deep subduction thrust interface (Stöckhert, 2002; Wassmann and Stöckhert, 2013).

A critical point arising from boudinaged quartz veins is that the quartz is clearly a relatively competent lithology; otherwise, the veins would not form boudins. On the other hand, the phyllosilicates that form the foliation are clearly highly strained, as evidenced by their tightly spaced foliation and the observation that the phyllosilicate foliation wraps around more rigid clasts of quartz or metapsammite (Figs. 3A, 3D, and 3E). A final point here is that, unless the phyllosilicates have broken down through dehydration reactions, they are both a potential source of fluids and likely the minerals that control bulk rheology.

The transitions from dominantly brittle to dominantly ductile behavior depend primarily on four parameters: composition, temperature, strain rate, and fluid pressure (Sibson, 1984; Tse and Rice, 1986). Of these, strain rate and fluid pressure may be highly variable in time and space, as mentioned earlier, and are the topic of this paper. Composition is also a space- and time-dependent variable, as various rock types reach the trench in spatially variable proportions (Underwood, 2007), and metamorphic reactions lead to depth-dependent changes in mineralogy. At shallow levels in particular, chemical reactions and precipitation of cements change both composition and microstructure of subducting sediments (e.g., Moore and Saffer, 2001; Moore et al., 2007). Thermal structure, on the other hand, varies significantly between subduction zones (e.g., Peacock, 1996; Syracuse et al., 2010), but if we assume shear heating is negligible, thermal structure is primarily determined by the age of subducting lithosphere and the vertical component of the trench-normal convergence (e.g., Molnar and England, 1990; Peacock, 1996), and therefore it is reasonably constant in time and space for a given margin.

Here, we are interested in the rheology of the subduction thrust interface at depths shallower than the mantle wedge. Thus, although numerous numerical models have been presented for subduction zone thermal structure (Peacock, 1996; van Keken et al., 2002; Wada and Wang, 2009; Syracuse et al., 2010; Rotman and Spinelli, 2013), for internal consistency, we use the analytical formulation of Molnar and England (1990) to compare some end-member, well-studied, subduction zones (Fig. 4A). Global parameters were chosen as in Table 1, including an effective coefficient of friction, defined as µ′ = µ(1 – λ). Following Wada et al. (2008), µ′ was set to 0.03, consistent with a weak subduction thrust interface and minor shear heating (e.g., Wang et al., 1995a; Gao and Wang, 2014). To show the global range in subduction zone thermal structure, we chose six subduction margins of variable age and convergence velocity (Table 2) and used margin-specific age, convergence rate, and slab dip values as presented in Syracuse and Abers (2006) and Syracuse et al. (2010). We did not include the Nankai Trough, despite it being exceptionally well studied, because ridge subduction at 15 Ma and the subsequent increase in subducting plate age from zero to 15 Ma make the analytical thermal model inappropriate (Wang et al., 1995b). For the Hikurangi margin, the dip is much gentler at depths less than 50 km than it is at deeper levels, as a result of the subduction of the Hikurangi Plateau, and so we used an average of 15° dip as calculated from the slab model of Williams et al. (2013).

The Cascadia margin, characterized by relatively slow subduction of young (10 Ma) oceanic lithosphere, stands out as a warm subduction margin (Fig. 4A). North Honshu (site of the 2011 Tohoku-oki earthquake) is particularly cold, characterized by relatively fast subduction of an old (older than 100 Ma) slab. Costa Rica is also relatively cold, because of fast convergence, despite the subducting slab being relatively young; Alaska, Hikurangi, and Chile represent intermediate but still quite cool thermal structures (Fig. 4A). Strikingly, and similar to numerical models by Syracuse et al. (2010), most considered margins remain at temperatures below 350 °C to depths of 50 km. This is important because the paradigm has long been that rocks are brittle and seismogenic to the depth of the 350 °C isotherm or the hanging-wall Moho, whichever is shallower (Hyndman et al., 1997). In this paradigm, of the subduction margins considered in Figure 4, only Cascadia would have a thermally controlled seismogenic zone base at a depth shallower than 40 km. The Nankai Trough, as modeled by Hyndman et al. (1995), is also relatively warm and fits a paradigm of a thermally controlled seismogenic zone.

We do not wish to embark on a detailed review of subduction thrust thermal structure. However, we make the point that it appears common for subduction thrusts to be actively deforming at temperatures of less than 350 °C to depths well below the hanging-wall Moho, when assuming that these faults are weak, and therefore shear heating is negligible. Wassmann and Stöckhert (2012, 2013) have made the point that dissolution-precipitation creep allows for very low effective viscosity along the subduction thrust interface in the blueschist stability field. We emphasize this point, and the fact that fluid flow along weak foliation enables element mobility and continued fault weakening, as illustrated by the development of an intense foliation in phyllosilicate-rich horizons and abundance of ductilely deformed quartz veins (Figs. 3A, 3D, and 3E).

In the exhumed examples we presented earlier in this paper, strain is generally localized within phyllosilicates, although the type of phyllosilicate varies with depth. At depths around the shallow ductile-to-brittle transition, the dominant structures are soft sediment folds, foliation defined by preferentially oriented clays, and veins parallel and oblique to this foliation. At depths around the downdip brittle-to-ductile transition, foliations defined by metamorphic phyllosilicates and flattened quartz and amphiboles wrap around boudinaged and folded quartz veins that both are parallel to and crosscut foliation. Quartz cements, and calcite at shallow levels, also play a rheological role in coating existing slip surfaces, and potentially forming load-bearing networks created by progressive and incremental vein formation. To gain a picture of shear strength as a function of depth, and potential depths of a downdip frictional to viscous transition, we calculated shear strength of phyllosilicates for the same six subduction zones as in Figure 4A. For the frictional regime, we continued the assumption of an effective coefficient of friction of 0.03 and that normal stress is approximated by the vertical stress; for the viscous regime, we used the muscovite flow law of Mares and Kronenberg (1993), adapted to the power law by Bukovská et al. (2016), and the biotite flow law of Kronenberg et al. (1990), both of the form:

For muscovite and biotite, respectively, the power-law exponent (n) is 11 and 18, the activation energy (Q) is 31 and 51 kJ mol^–1, and the material parameter A is 4 × 10^–20 and 1.2 × 10^–30 (Kronenberg et al., 1990; Mares and Kronenberg, 1993; Bukovská et al., 2016); R is the universal gas constant, and

The outcome of these strength calculations is that the warmest margin considered, Cascadia, has a relatively shallow predicted base of the seismogenic zone at 10–15 km and temperature of less than 300 °C (Figs. 4B and 4C). This is equivalent to the current geodetically determined locked zone at 10–15 km, but it is shallower than the depth of episodic tremor and slow slip, which is deeper than 30 km (McCrory et al., 2014). For the other subduction zones considered, the predicted deep frictional to viscous transition lies between 10 and 25 km depth (Figs. 4B and 4C), but at a temperature much less than 300 °C (Fig. 4A), assuming that phyllosilicates determine the behavior of the subduction thrust interface and that laboratory flow laws can be extrapolated to lower temperatures and slower strain rates than those used in the laboratory. If we make our calculations with a quartz flow law, we arrive at frictional-viscous transitions deeper than 40 km for all margins.

The calculations presented in Figures 4B and 4C assume a very low effective coefficient of friction; if fluid pressure is not near lithostatic, frictional strength will increase to more than the few tens of megapascals predicted here. Although we have not included flow strength of dissolution-precipitation creep in quartz, or quartz-phyllosilicate mixtures, in our calculations, the low flow strengths predicted are consistent with low viscosities inferred by Wassmann and Stöckhert (2013) and Fagereng and den Hartog (2017) assuming a dissolution-precipitation creep rheology. Similarly, the abundant presence of tensile veins implies differential stress of no more than a few tens of megapascals (Etheridge, 1983), consistent with the calculations. Boudinage of these veins implies that quartz has greater viscosity than the surrounding phyllosilicate-rich metasediments, supporting the assumption of a phyllosilicate-controlled bulk viscous rheology.

The observations from exhumed subduction thrust rock assemblages, including the inference of independent particulate flow at low effective stress at shallow levels, and pervasive vein growth at both the updip and downdip limits of the seismogenic zone, imply that the hydrofracture criterion is commonly met, particularly locally and transiently in zones of mixed brittle-ductile deformation. This observation implies that fluid overpressured conditions are common in these zones, but it also demands that a fluid source is available to supply the fluids and that these fluids do not escape through existing permeability pathways. We next model (1) the fluid release from shallow porosity loss and progressive very low-grade reactions as sediments enter the trench and are transported into the seismogenic zone, and (2) the higher-grade metamorphic reactions expected as oceanic crust and overlying sediments are subducted to depths around the deep frictional-viscous transition.

At shallow depths, where a clay fabric and early stage mélange structures develop, fluids are released as sediments lithify into rock and fluid-filled pores are squeezed, reducing primary porosity. We know from Ocean Drilling Program (ODP) and Integrated Ocean Drilling Program (IODP) drilling that sediments are typically underconsolidated below the décollement, and thus they carry fluids deeper than most, if not all, other settings (e.g., Morgan and Karig, 1993; Housen et al., 1996; Screaton et al., 2002; Tobin and Saffer, 2009). The feedbacks among porosity, permeability, effective rock stress, and fluid production can lead to a complex evolution in fluid pressure and rock properties that is tightly coupled to the mechanical behavior of the wedge (Saffer and Bekins, 1998; Ellis et al., 2015). This can be illustrated using a simple example generated from a fluid-thermomechanical numerical model (Fig. 5). Figure 5 examines the evolution of the outermost part of an accretionary wedge with a basal décollement and a layer of sediment that subducts beneath the wedge. We prescribe a uniformly tapered wedge with a horizontal décollement (Fig. 5A), although in real accretionary wedges, the base will dip gently landward. The model is not a realistic representation of a “real” subduction margin, but it serves to illustrate some of the important feedbacks among porosity, fluid pressure, and mechanical behavior in an accretionary wedge.

Figure 5B shows the updated wedge geometry and velocity field for flow of material into the wedge after 500 k.y. of deformation. The weak frictional décollement has allowed most of the wedge to decouple from the subducting oceanic crust and underthrust sediments (as in Fisher and Byrne, 1987), with active accretion occurring at the wedge toe. The porosity field at this step is also shown in Figure 5B. Porosity is computed at each step as a function of effective pressure, that is, the mean rock stress from gravitational loading and tectonics minus the fluid pressure ( Appendix 1, Eq. A2). Once porosity has been reduced due to loading in this manner, it cannot recover, although secondary (fracture) porosity may develop. Figure 5B shows that porosities near the base of the wedge are reduced to <10% by lithostatic loading. However, higher porosity is maintained below the décollement in the subducting sediments. This is because of the decoupling in stress across the décollement, with greater (near-horizontal) tectonic loading within the wedge that decreases in the underthrust sediments, so that maximum compressive stress in the lower plate is near-vertical. As mentioned already, this decoupling has been observed in drilling studies, for example, in Nankai, Barbados, and the Japan Trench (e.g., Moran et al., 1993; Housen et al., 1996; Lin et al., 2013). It is a key point that allows transportation of significant fluid volumes to depth below the accretionary wedge.

Close inspection of Figure 5B also reveals subtle variations in porosity across the wedge, which result from changes in dynamic stresses as wedge material deforms. Figure 6A shows the fluid release from porosity loss calculated for the model after 500 k.y. of deformation, assuming the material flow and porosity fields from Figure 5B (see  Appendix 1 for details of this calculation). This reveals a general trend of fluid release at the toe of the wedge (the zone over which horizontal tectonic stresses are high), focused at the deformation front and landward of this in shear bands where tectonic compaction is occurring. In contrast, the fluid release caused by porosity loss in the subducting sediments extends further landward owing to the higher porosities maintained there. Peak fluid release rates of ~2 × 10^–15 m³ m^–2 s^–1 are comparable to those estimated in fluid budgets for real subduction margins such as Hikurangi and Nankai (e.g., Pecher et al., 2010; Moore et al., 2011; Ellis et al., 2015).

In addition to porosity loss, low-grade mineral reactions occur within sediments as they are accreted or subducted beyond the trench. For example, smectite clays progressively transform to illite in the temperature range from ~65 °C to 130 °C (Pytte and Reynolds, 1989). The transition from opal to quartz is another diagenetic dehydration reaction occurring around T = 100 °C (Kameda et al., 2012), and oceanic crust also contains hydrous phases such as saponite and zeolites formed by seafloor hydration and alteration. Zeolites, however, remain stable through the shallow transition from ductile to brittle deformation and only dehydrate to form greenschist-facies minerals at 300 °C or more (e.g., Kerrick and Connolly, 2001; Hacker et al., 2003; Fagereng and Diener, 2011). Saponite also carries a significant crystal-bound water content into the subduction zone (Kameda et al., 2011), but it predominantly dehydrates at temperatures in excess of 200 °C (Hillier, 1993), and then it creates chlorite, which itself is more hydrous than illite-muscovite. Dehydration of the oceanic crust through the breakdown of hydrous minerals, therefore, predominantly contributes to the fluid release at depths greater than 15–20 km, but maybe also at shallower levels if the released fluids flow along the décollement. We discuss these deeper processes later.

Figure 6B shows color contours of associated dehydration fluid release in units of cubic meters per square meter per second (m³ m^–2 s^–1) calculated for the smectite to illite reaction taking place in sediments and the slab (details given in  Appendix 1), assuming an approximate thermal steady state that includes rock advection. We have assigned 20% initial smectite fraction in incoming sediments at low temperatures, and 10% smectite content in the slab and hemipelagic sediments; in both, we assume a volumetric fraction of water bound in the hydrated smectite of 0.35 that is available to be released upon transition to illite. Results indicate that the transition from smectite to illite is gradual. For the thermal field in the model, the transition mostly occurs within subducting sediment and begins at the toe of the wedge, extending ~40 km landward, but it is an order of magnitude smaller than the fluid release from porosity loss (Fig. 6A). Wedges that have a low geothermal gradient experience dehydration fluid release at greater depths, and if the fraction of clay in incoming sediment is significant, it may still produce a significant effect on fluid pressures, allowing porosity to be maintained in subducting sediment and thus interacting in a positive way with fluid release from porosity loss in the subducting sediments.

The combined fluid sources as shown in Figures 6A and 6B generate excess fluid pressures below the outer accretionary prism. This is comparable to calculations by Kitajima and Saffer (2012), who showed elevated pore pressures, coinciding with very low-frequency earthquakes, in the outer 50 km of the Nankai forearc. The extent to which this outer prism fluid pressure will affect wedge strength and porosity depends on the assumed permeability function in wedge and décollement sediments (e.g., Ellis et al., 2015). For the model example shown in Figure 6C, we have assumed the porosity-permeability relationship from Saffer and Bekins (1998) for all sediment and for the subducting slab. Since the décollement experiences high brittle strains, and we have prescribed permeability to change with strain, permeability there increases by a factor of 10 after a short time. For this simplified assumption, and without significant fault permeability that could drain the wedge, moderate fluid overpressures develop (λ ~ 0.6; Fig. 6C) after 500 k.y. Although some fluid overpressure is present seaward of the toe of the wedge within accreting sediment, most occurs within the subducting sediments and underlying slab, where localized pockets of higher fluid pressures develop where porosity loss drives fluid release into low-permeability sediment. For a lower décollement permeability and/or permeability that seals at increasing temperature, higher overpressures may develop, significantly affecting wedge mechanics and weakening the wedge. In contrast, Figure 6D shows fluid flow and overpressure predictions for a similar model but where slab permeability is a constant 10^–18 m², where considerable fluid flow within the slab limits the overpressure that can develop in the subducted sediments. The enhanced fluid flow along the high-permeability décollement and diffuse fluid flow in the moderately permeable slab bleed off sufficient fluid so that extreme overpressures do not develop in this example.

Note that in some subduction zones, markedly the Japan Trench and Costa Rica, incoming pelagic sediments are dominated by siliceous ooze, making the opal to quartz transition more important than what is modeled here (e.g., Tsuru et al., 2002; Spinelli and Underwood, 2004; Kameda et al., 2012). The Sumatra margin may also deviate significantly from our calculations, given that the thick incoming sediment pile may be largely dehydrated, and may have passed through the smectite-illite transition before entering the subduction zone (Hüpers et al., 2017). Thus, noting the composite fluid sources within real subducting sediments (e.g., Moore and Saffer, 2001), our models are a simplification to understand the dominant processes, but local variations are clearly present.

In the previous section, it was demonstrated that the main fluid source at shallow depths is porosity loss (Fig. 6A), with fluids released from the breakdown of smectite, and opal in some margins, providing a minor contribution (Fig. 6B). However, once the rocks have been completely cemented and lithified, the crystal-bound water in hydrous minerals becomes the main fluid reservoir. This fluid will only be released at the conditions where the relevant hydrous mineral breaks down, and consequently fluid production in and around the megathrust becomes a function of the mineral stability and mineralogical evolution of the different lithologies involved.

Mineral equilibria modeling calculations can be used to determine the pressure-temperature stability of mineral assemblages, including the conditions and amounts of fluid that will be produced when these minerals are destabilized (e.g., Peacock, 1990; Kerrick and Connolly, 2001; Fagereng and Diener, 2011). These calculations rely on the availability of appropriate activity-composition relations for the relevant minerals, and unfortunately models for many of the low-grade phases such as smectite, opal, saponite, illite, and the various zeolites have not been formulated. Consequently, the calculated equilibria at temperatures below ~300 °C are likely to be metastable to these phases. Because of this, we cannot model and describe the low- to high-temperature fluid release processes with perfect continuity.

Pseudosections showing the stability of mineral assemblages as a function of P and T were calculated for the composition of the average global subducting sediment (GLOSS of Plank and Langmuir, 1998) and a typical mid-oceanic-ridge basalt (the MORB composition of Sun and McDonough, 1989). These diagrams are presented in Figures 7A and 7B, and the details of the calculation procedure are presented in  Appendix 2. The calculated amount of fluid that is cumulatively held in the hydrous minerals and, conversely, the amount of fluid released when these minerals break down are shown on Figures 7C and 7D. What is immediately apparent from these diagrams is that fluid production is episodic, rather than continuous, and that there are large parts of a rock’s P-T history where little or no fluid is produced. In fact, many dehydration reactions occur at relatively confined P-T conditions, such that the majority of a rock’s fluid release occurs in a small number of discrete events, each releasing a relatively large volume of crystal-bound fluids over a small P-T interval.

The MORB composition is capable of holding much more water because of its greater abundance of hydrous minerals. The maximum H₂O content calculated for the MORB is ~16.5 mol % H₂O, whereas for GLOSS, it is only 10.4 mol % (Figs. 7C and 7D). The main minerals that are reservoirs of fluid, and the breakdown of which therefore controls fluid production in both GLOSS and MORB, are lawsonite and chlorite. Other hydrous minerals such as white mica (muscovite or phengite), amphibole, and epidote are stable over the entire P-T range, and their abundance does not change greatly, such that they do not contribute significantly to the fluid budget in subduction zones over this P-T range.

The breakdown of lawsonite, in particular, causes the release of large volumes of fluid over a narrow P-T band at conditions between 0.25 GPa at 200 °C and 1.6 GPa at 450 °C (Figs. 7C and 7D). In GLOSS, which is calculated to contain ~10–12 vol% lawsonite, this results in the release of 4–5 vol% H₂O. By contrast, MORB, which is calculated to contain more than 20 vol% lawsonite, produces more than 7 vol% H₂O when the lawsonite breaks down.

GLOSS contains only 5% chlorite that breaks down at temperatures ~50 °C higher than lawsonite breakdown, but only at P below 0.8 GPa (Figs. 7A and 7C). Chlorite breakdown produces an additional 1.5 vol% fluid, and the GLOSS composition does not release significant amounts of additional fluids at higher T or P. On the other hand, MORB contains 20–25 vol% chlorite, which releases an additional 4–5 vol% H₂O as it breaks down. There are a number of equilibria reactions that consume chlorite; the first is associated with the destabilization of omphacite in favor of actinolite over a narrow window from 0.8 GPa at 400 °C to 450 °C at 1.6 GPa. This reaction produces 1.5 vol% fluid. Then, at 450–500 °C and P < 1.1 GPa, the breakdown of chlorite is focused and corresponds to the transition from the greenschist or blueschist facies to the amphibolite facies. This transition can release between 3 and 5 vol% fluid (Figs. 7B and 7D).

The thermal gradients of many subduction zones are within the stability field of lawsonite, and they are subparallel to the lawsonite breakdown reaction, but at temperatures that are at least 50 °C colder (Figs. 7C and 7D). These subduction zones are therefore modeled to not produce much (if any) internal fluids at depths within and downdip of the seismogenic zone (Fig. 7E). Rather, the fluid in these subduction zones is retained in hydrous minerals (notably lawsonite) and carried to depths beyond 50 km, where it is likely released to the mantle wedge. However, a notable exception is Cascadia, which, because of its warmer geotherm, intersects the lawsonite breakdown reaction at 20–25 km depth (Figs. 7C, 7D, and 7E). At these conditions, we calculate that almost 7 vol% fluid is released from the MORB, whereas 5.5 vol% is produced by GLOSS. This one event accounts for the release of about half of the entire fluid budget held in both rock types (Fig. 7E). Following this dramatic fluid release event, GLOSS is modeled to not produce any additional fluid with further burial, but MORB is modeled to release an additional 1.5 vol% H₂O through chlorite breakdown at a depth of around 38–40 km (Fig. 7E).

The degree to which fluid production can lead to fluid overpressures depends on permeability, and subduction zone permeability is time- and space-dependent. High permeability can be expected along the décollement zone if the fault is not well cemented (Moore, 1989), and along subducting sandy layers that have remained continuous rather than dismembered layers (Carson and Screaton, 1998). At shallow depths, a variably fractured wedge with similarly variable porosity may allow fluids to escape from the fault zone and into the wedge (Fig. 6; Screaton et al., 2002). Thus, at depths where sediments are relatively permeable (i.e., cementation is not complete), we envisage that fluid pressure along the interface is limited to moderate values (Fig. 6), as a consequence of leakage through the overlying wedge and along the megathrust. On a small scale, drainage can be envisaged to occur through steep fractures in sandstone lenses, and along clay foliation (Fig. 1C), emphasizing the anisotropy of the megathrust permeability structure. However, as highlighted by Saffer and Tobin (2011), permeability anisotropy may only account for a factor of 10 variability in permeability (Kwon et al., 2004), compared to at least four orders of magnitude in fault-parallel permeability reported in a range of subduction zones, as reviewed by Saffer (2015). This variation is a testament to the range of materials present along the subduction thrust interface, as well as the dynamic nature of its deformation. Saffer (2015) concluded that transiently high permeability, orders of magnitude greater than that for the sedimentary matrix, is required to explain observations such as thermal anomalies across faults in Barbados (Cutillo et al., 2003) and pore-water freshening in Costa Rica borehole profiles (Spinelli et al., 2006). Such high permeability may be related to interconnected faults and fractures making a zone of elevated permeability (Sibson, 1996); however, this permeability may be created episodically through incremental brittle deformation, and it may be destroyed by subsequent vein precipitation (Fisher et al., 1995; Sibson, 1996). A fault-fracture mesh, which crosscuts all the components of the mélange after complete lithification (Figs. 1F and 1G), is present in several fossil examples where deformation occurred at T ≥ 100–150 °C (Fisher et al., 1995; Kondo et al., 2005; Rowe et al., 2009; Vannucchi et al., 2010; Fagereng et al., 2011; Fisher and Brantley, 2014); however, these permeability pathways are identified by veins and cements in exhumed accretionary complexes, and thus they are unlikely to have been long-lived before chemical equilibria required their sealing, or movement of a finite volume of pressurized fluid changed the fluid pressure and permeability distribution in the fault zone (Bekins and Screaton, 2007).

Starting at depths of around 10 km, where porosity becomes negligible in cemented sedimentary rocks, permeability surrounding the megathrust is likely negligible, and a greater permeability gradient may exist between the subduction thrust interface and surrounding rocks (Moore, 1989; Kato et al., 2004). Fisher and Byrne (1987) showed that veining and stratal disruption in the Kodiak tectonic mélange, Alaska, is restricted to the footwall, bounded by a late discrete fault at its upper boundary. In a study on the Uganik thrust in the Kodiak complex, Rowe et al. (2009) reported that a heterogeneous vein distribution occupies a greater thickness in the subduction thrust footwall than in the hanging wall. The fact that stratal disruption, fracture initiation, and vein generation dominantly occur in the footwall rocks of exhumed subduction thrusts has now also been demonstrated in several other locations (e.g., Meneghini and Moore, 2007; Fagereng, 2011b). This localization of veins to the footwall implies that the fault is a permeability barrier, as may also be inferred from the abundance of foliation-parallel veins (Fig. 3), indicating fluid flow parallel to, but rarely across, fault-related fabrics. The fluid pressure redistribution that may occur along a relatively permeable (at least transiently) fault at the upper ductile-to-brittle transition is therefore less likely at the downdip end of the seismogenic zone, where permeability is more uniform in the range of low-porosity rock types. The presence of extensive foliation-parallel veins implies elevated fluid pressures at this depth; on the other hand, relatively short and rarely interconnected veins (Fig. 3) may imply that fluids at the downdip end of the brittle regime are less likely to be interconnected. This interpretation does not, however, preclude transient events of significant fluid redistribution. For example, Husen and Kissling (2001) reported postseismic fluid flow after the 1995 M_w 8.0 Antofagasta earthquake, based on time series of Vp/Vs ratios. There is, therefore, a possibility that fluids are stored in low-permeability rocks at the downdip end of the seismogenic zone, and they are only able to escape through the impermeable subduction thrust fault zone and into the over-lying plate when fault zone damage and postseismic stress changes cause a transient permeability increase. A consideration here is that in an Andersonian stress regime favoring reverse faulting, the least principal stress is vertical, causing tensile fractures to be horizontal, or to form along near-horizontal weak foliation (Behrmann, 1991). Thus, even if the hydrofracture criterion is achieved, fluids trapped below a low-permeability fault at the downdip end of the seismogenic zone will be redistributed within fractures subparallel to the fault, rather than allowed passage for significant distances into the upper plate.

Saffer and Tobin (2011) postulated that there are two zones of extreme overpressure, roughly coincident with the transition zones in frictional behavior at the up- and downdip limits of the seismogenic zone. In between, they suggested that the seismogenic zone is moderately overpressured. This model is supported at the first order by the observations that hydrothermal vein systems appear to be generated, preferentially although not exclusively, at the ductile-brittle and brittle-ductile transitions marking the approximate up- and downdip limits of the frictional regime. We do, however, suggest some modifications to the model based on our observations and calculations (Figs. 6 and 7).

The message from calculations at shallow depths is that substantial fluid volumes are carried within subducting sediments, and they are released as pore spaces break down (Saffer, 2003; Moore et al., 2011; Ellis et al., 2015). Although the bulk of the pore fluids may be expelled during early accretion (Moore et al., 2011), it is significant that low effective stresses beneath the décollement, and a strain reversal from horizontal shortening in the prism to vertical shortening within underthrust, mostly low permeability, sediments, allow for subduction of underconsolidated, fluid-rich sediments beneath the accretionary wedge. As a result, pore fluids can be carried to depths of at least a few kilometers within the underthrust sediment pile (Fig. 6A). As a consequence, a gradual increase in fluid pressure is likely to occur at shallow depth where pore fluids are released, with a peak at very shallow levels, but continued fluid production for tens of kilometers landward of the trench. Furthermore, the smectite-illite reaction releases additional fluids, although, depending on the thickness of subducting sediments and their clay content, this is generally a smaller volume than what is released through porosity loss (Fig. 6B). The smectite breakdown reaction occurs at similar depth to the final porosity loss, although because it depends on temperature, smectite breaks down deeper in colder subduction zones (Fig. 8), where opal dehydration may also be important (Kameda et al., 2012). Overall, fluid release at shallow levels is continuous over a depth range reaching down to final smectite breakdown at ~150 °C.

Fluid pressure under the sedimentary, outer accretionary wedge is limited by the permeability of the décollement and the overlying prism. As attested by common fluid seeps at the seafloor above shallow décollements (e.g., Saffer and Bekins, 1998; Lewis and Marshall, 1996; Suess et al., 1998; Bohrmann et al., 2002; Sahling et al., 2008; Barnes et al., 2010; Pecher et al., 2010; Saffer, 2015), secondary porosity through hanging-wall faults is likely to keep fluid pressures at moderate levels, as is fluid escape along the décollement (Fig. 6C). However, soft-sediment deformation structures implying flow of fluidized sediments at low effective stress (Figs. 1A and 1B) are widespread, and low-frequency earthquakes also indicate low effective stress under the outer prism in Nankai (Ito and Obara, 2006; Kitajima and Saffer, 2012). Quartz and calcite vein generation (Figs. 1C–1G) attest to local fluid overpressures, both in subvertical cracks in high-tensile-strength lenses present from very shallow levels (Fig. 1C), and along faults (Figs. 1E and 1F), implying that fluid pressure overcomes a subhorizontal least compressive stress or is locally elevated along a weak fault surface. We note that the onset of vein precipitation, and associated cementation of sediments, also decreases permeability, and there is likely an increase in the maximum contained overpressure around the depth of the shallow ductile-to-brittle transition. This transition in permeability occurs because brittle deformation implies increased cohesion, but it also enhances fluid flow and cementation through cracks, transiently and locally increasing and subsequently decreasing permeability through fracture and vein generation.

Further downdip, not much fluid is released until specific dehydration reactions, particularly breakdown of lawsonite, lead to a very significant increase in fluid production (Figs. 7 and 8). In warm subduction zones, such as Cascadia, this fluid release occurs at depths of 20–30 km, i.e., above the hanging-wall Moho, but below the observed downdip limit of seismicity (McCrory et al., 2014). The fluid release does, however, coincide with the depth of tremor and slow slip (Fagereng and Diener, 2011). In the other subduction zones we considered here, the modeled geotherms are too cold for lawsonite breakdown at depths above the Moho (Fig. 7). Lawsonite and chlorite therefore remain stable through the seismogenic zone, and the slabs carry their crystal-bound fluids into the mantle wedge. Numerical thermal models predict a significant increase in slab and megathrust temperatures below the Moho, when the subducting slab comes into contact with the upper-plate lithospheric mantle (e.g., Syracuse et al., 2010). This temperature rise is likely to trigger a release of the water held in the lawsonite and chlorite crystal structures in short succession (Fig. 8B). This fluid release may explain the presence of tremors near the mantle wedge corner in many subduction zones (e.g., Brown et al., 2009; Katayama et al., 2012; Hyndman et al., 2015).

Permeability is highly anisotropic at depths beyond the final significant porosity loss, i.e., where the fault zone–parallel permeability can become much greater than the permeability of surrounding, low-porosity rocks (Kato et al., 2004), unless the fault is efficiently sealed (Audet et al., 2009). Fluids generated below the seismogenic zone may therefore migrate upward, but only if the megathrust fault zone is sufficiently permeable. Swarms of foliation-parallel, ductilely deformed quartz veins are common in subduction thrust interfaces exhumed from depths around or below the deep brittle-ductile transition (Fig. 3; Fagereng et al., 2014; Bachmann et al., 2009). These veins are testament to localized fluid overpressure in excess of a steeply plunging least compressive stress, and sealed fractures aligned parallel to, and therefore not allowing fluid flow through, the low-permeability subducting metasedimentary rocks. We therefore suggest fluid pressure will locally increase dramatically as a result of metamorphic fluid release within low-permeability rocks, as fluids can only travel up the subduction interface or through the interface into upper-plate rocks if there is a transient increase in permeability, for example, after a major earthquake (Husen and Kissling, 2001). The locally increased fluid pressure will generally, at least in cold subduction zones, occur below the brittle-ductile transition, and it may be observed geophysically as tremor and slow slip, commonly observed at, or more usually below, the interseismically locked zone (Dragert et al., 2001; Brown et al., 2009; Beroza and Ide, 2011).

We predict that fluid pressure generation is episodic in space (Figs. 7 and 8); consequently, fluid pressure is likely also compartmentalized, but the peaks predicted from fluid release may be altered by the subduction thrust permeability structure, where the generated fluids can flow from sites of production to sites of high permeability, if these locations do not coincide. At shallow depths, fluid production is relatively continuous, and permeability is heterogeneous, but probably, on average, relatively high, at least transiently through hanging-wall faults or along the décollement (Fig. 5B; Saffer, 2015). As a result, maximum sustainable overpressure is moderate and likely occurs in patches where permeability is locally (and maybe transiently) low (Fig. 6C). In contrast, the deep subduction thrust interface is surrounded by low-porosity rocks, and it is likely cemented and fine grained itself, leading to a low time-averaged permeability, particularly oblique to the intensely foliated fault zone. Fluids are therefore generally trapped below or within the megathrust fault zone (Byrne and Fisher, 1990; Meneghini and Moore, 2007; Rowe et al., 2009), leading to near lithostatic fluid overpressures sufficient for vein generation along gently dipping foliation (Fig. 3).

At shallow levels, porosity loss is an approximately continuous source of fluids, maintaining moderate fluid overpressures that allow soft-sediment deformation at low effective stresses, such as independent particulate flow, in the accretionary prism (Fig. 9). On the other hand, mineral dehydration occurs in focused, transient episodes as a subducting volume of rock passes through a phase boundary (Fig. 7), but major dehydration reactions are not intersected by typical subduction geotherms until the slab reaches depths deeper than the interseismically locked zone, below the inferred brittle-to-ductile transition (Fig. 9).

At shallow depths, we have also made the point that although there is a ready fluid source, permeability along faults and through poorly compacted and cemented sediments is such that fluid pressures are elevated but limited (Fig. 6C). Similarly, observations in rocks deformed around the shallow ductile-to-brittle transition highlight distributed deformation reflected by continuous structures such as folds and boudins, crosscut by discrete faults (Figs. 1A, 1B, and 9). The soft-sediment deformation structures are easier to form if effective stress is low, but they do not show evidence for spatially fluctuating stress or fluid pressure. Thus, in the ductile regime, permeability may be sufficient for a broadly constant fluid pressure in time, although local changes may occur as a function of time-variable permeability as sediments turn to rock through diagenetic processes.

Downdip of the shallow ductile-to-brittle transition, particularly below the subduction thrust interface, fluids are likely trapped by fine-grained, impermeable, preferentially oriented phyllosilicates making up a fault-parallel barrier to fluid flow. Incrementally developed, fault-parallel veins are common within foliated subduction thrust fault rocks (Figs. 1F and 1G; Fagereng et al., 2010), and they provide evidence for intermittent failure followed by mineral precipitation. These veins also require fluid pressure fluctuation, which implies a fluctuation in the failure stress of the fault surface, linked to fluctuations in fault effective normal stress (Eq. 1; Fig. 9E). These veins thereby highlight the intermittent, time-dependent nature of local fluid pressure and effective stress along the megathrust (Fig. 9E).

Local fluid overpressure has the effect of lowering the shear stress required for frictional failure, as effective normal stress is reduced (Hubbert and Rubey, 1959). An effect of this shear strength decrease is the embrittlement of faults that are surrounded by ductilely deforming host rocks. This embrittlement has been invoked to explain the presence of pseudotachylyte within mylonites (e.g., Sibson, 1980; Hobbs et al., 1986; White, 2012); however, because very low shear stresses are required for frictional failure to be preferred over ductile flow at elevated temperatures (e.g., Figs. 4B and 4C; Brace and Kohlstedt, 1980), it is likely that frictional failure below, or even within, the seismogenic zone occurs at shear stress values too low to generate friction melt. Exceptions are, however, illustrated by the presence of pseudotachylyte in the Shimanto (Ikesawa et al., 2003) and Kodiak (Rowe et al., 2005) exhumed accretionary prisms. Although slip rates along shear veins and discrete slip surfaces observed within rocks exhumed from below the brittle-to-ductile transition are not directly measurable, we speculate that transient brittle deformation and/or frictional failure, recorded as earthquakes, low-frequency earthquakes, and tectonic tremors below the interseismically locked zone at temperatures in excess of 300 °C, may be preserved in the rock record as discontinuous structures within dominantly ductilely deformed rocks. Similarly, brittle structures within ductilely deformed rocks at very shallow levels, where smectite is still preserved (Fig. 1A; Bettelli and Vannucchi, 2003), may also record local and transient embrittlement. For example, Kimura et al. (2012b) suggested that the shallow propagation of the Tohoku-Oki 2011 earthquake was aided by propagation into fault segments where fluid overpressures were contained by an impermeable clay fabric, and thus the rocks behaved brittlely when a rupture front arrived. It is an important caveat that these interpretations of local embrittlement and brittle failure involve secondary porosity generation, such that, by definition, areas where increased fluid pressure allows local embrittlement will be transient, and no longer fluid overpressured after failure, at least until a new overpressure can be generated by fluid influx from a fluid source.

We can, overall, raise a few questions and hypotheses from the interpretation that fluid pressure peaks are local and transient. Because fluid overpressure lowers effective normal stress, it increases the likelihood of slip nucleation; therefore, if “asperity” is defined as a point where the interface is interseismically locked and inferred to slip episodically in earthquakes, do zones of fluid overpressure represent asperities on the subduction thrust interface? The upper limit of the seismogenic zone likely coincides with the downdip limit of dominantly ductile deformation of low-cohesion sediments, and thus also the transition from moderate to possibly higher fluid overpressures. As a consequence, is there a fluid pressure peak at the updip limit of the seismogenic zone (e.g., Saffer and Tobin, 2011)? Similarly, the downdip limit of the seismogenic zone represents a transition from episodic slip to bulk shearing flow that is continuous in time, where the latter may maintain some dilatation and associated fault-parallel permeability. Is there, therefore, a fluid pressure peak maintained by a low-permeability seal at the downdip end of the seismogenic zone (Saffer and Tobin, 2011)? If these last two hypotheses are correct, the seismogenic zone is, in a time-averaged sense, surrounded by elevated fluid pressures. Finally, if we are correct in positing that the bulk of fluid release from metamorphic reactions occurs some distance downdip of the predicted lower limit of the brittle subduction thrust interface, this fluid release may coincide with the source area of tremors and slow slip. If so, it is likely that the episodic failure and healing as fluid pressure reaches a failure stress, driven by relatively constant fluid release at a point in space, are related to the near-periodicity of tremors on geophysical observational time scales.

Fluids are not only involved in weakening through lowering effective normal stress. Ongoing reactions may lead to additional weakening, such as through the creation of interconnected phyllosilicate networks. Analogous to continental crust, it has been inferred that the frictional to viscous transition is controlled by the onset of crystal plasticity in quartz at 350 °C (e.g., Hyndman et al., 1997). However, observations along subduction thrust faults imply that phyllosilicates are more likely to control rheology (Figs. 1 and 3), leading to much weaker rocks and a shallower frictional viscous transition coinciding with the depth of the 150–300 °C isotherm. The creation of the interconnected phyllosilicate network required for these minerals to control rheology is in part a function of the smearing capabilities of simple shear; however, ongoing reactions forming new phyllosilicates are efficient at creating preferentially oriented and connected phyllosilicate foliation (e.g., Wintsch et al., 1995).

Fluid flow, such as along anisotropic fault-parallel fabrics, also allows dissolution-precipitation creep, which is particularly efficient in fine-grained fault rocks (e.g., Niemeijer and Spiers, 2006; Gratier et al., 2011; den Hartog et al., 2012). This deformation mechanism allows low effective viscosities and deformation at very low shear stresses, as calculated from quartz microstructures in blueschist-facies rocks (Stöckhert, 2002; Wassmann and Stöckhert, 2013). Thus, in addition to the brittle deformation emphasized here, fluid movement down pressure gradients may allow deformation by dissolution-precipitation creep at low driving stresses and slow strain rates, for example, in the interseismic period, or in a creeping subduction thrust where no major seismic events are required to accommodate plate convergence (Fagereng and den Hartog, 2017).

An interesting caveat is whether fluids can also provide a mechanism for fault zone strengthening. For example, near the updip ductile-to-brittle transition, dissolution-precipitation creep allows cements to form, creating load-bearing strong frameworks (Moore et al., 2007). This enhanced diagenesis by fluid-driven cementation leads to a strength increase of the fault and a transition from independent particulate flow in cohesionless sediments to cataclasis of cohesive rocks. Similarly, at both the up- and downdip limits of the brittle subduction thrust interface, there is a process of episodic fluid pressure weakening resulting from repeated events of intermittent fluid pressure accumulation and release—this may be geophysically recorded as tremors or microseismicity. Geologically, this process can incrementally form quartz veins (Fagereng et al., 2011). These quartz veins can form strong layers, as shown by formation of boudins (Fig. 3); networks of veins that form in zones of elevated fluid pressure can, therefore, over time, strengthen these regions.

At the updip limit of seismicity, the process of fluid-related strengthening by cementing the fault rock has been discussed previously by Moore and Saffer (2001) and Moore et al. (2007). We add that the pressure-temperature range where this occurs may be broad, because fluid release at the updip edge of the seismogenic zone affects a large region. This large footprint of fluid release occurs because: (1) calcite dissolution and reprecipitation are active at very low temperatures, below the ~100 °C isotherm where quartz solubility becomes sufficient for effective pressure solution and quartz precipitation, whereas (2) underconsolidated, high-porosity sediments can carry fluids deeper than the 100 °C isotherm. Thus, cementation starts at T << 100 °C and carries on to T >> 100 °C. At the downdip limit, we suggest that a fluid-flow barrier exists, and therefore a similar process occurs, where quartz veins preferentially develop as fluids flow upward from the depth of fluid generation along a fault zone that may have some fault-parallel but little fault-oblique permeability. These veins preferentially form along foliation and may lead to long-term strengthening by breaking down the interconnected, weak foliations required for fault weakening, and replacing this foliation with a stronger, load-bearing framework (e.g., Handy et al., 2007).

The subduction thrust interface evolves from a mixture of fluid-saturated, unconsolidated, high-porosity trench-fill sediments deposited on top of oceanic crust covered in clay-rich pelagic sediments to a mélange of low-porosity, foliated metamorphic rocks (Fig. 9A). Through this evolution, the initial fluid content of the subducting crust and sediments is expelled, gradually and intermittently (Fig. 9B). Similarly, composition, structure, and strength also evolve within the subducting rocks, such that within the subduction thrust interface, a number of deformation mechanisms coexist, but their relative importance varies in time and space (Fig. 9). We have, in agreement with Saffer and Tobin (2011), identified two main fluid source regions:

• (1) Shallow porosity loss and clay dehydration: As sediments are accreted or subducted, pore spaces close, and fluids are expelled. This occurs over a large zone, largely because stress decoupling across the décollement allows underthrust sediments to retain their pore fluids in a regime of light effective tectonic loading and a vertical greatest compressive stress. Similarly, clay dehydration is a gradual process as low-temperature reaction kinetics are slow, and it occurs over a similar P-T range to pore fluid expulsion. The result in terms of megathrust deformation is a shallow zone, at T ≤ 150 °C, where soft-sediment, ductile flow dominates, but discrete faults develop as sediments become cemented and pore fluids are expelled (Fig. 9).

• (2) Deep, T-dependent release of crystal-bound waters: As low-porosity oceanic crust and underthrust sedimentary rocks approach and enter the ductile regime, where rocks may deform plastically, metamorphic dehydration reactions release fluids from hydrous mineral phases. This dehydration occurs at specific P-T conditions where mineral reactions are reached (Fig. 7). Because permeability is low within these rocks, which are buried to tens of kilometers depth at temperatures where crystal plastic flow can close pore spaces, fluid pressures are elevated as fluids are released. Thus, ductilely deformed vein swarms are common features in rocks from the deep brittle-to-ductile transition, recording fluid pressure fluctuations occurring in relation to metamorphic fluid release, fluid pressure buildup, fracture, and vein growth (Fig. 9E).

This review has summarized the observational evidence from a number of areas (Figs. 1 and 3) and compared these observations to some relatively simple models of fluid generation (Figs. 5–7). Clearly, much scope exists to refine, confirm, or disprove the hypotheses we have put forward. For example, before the recent discovery of tectonic tremor and slow slip, we would not have thought of interpreting ductilely deformed, foliation-parallel veins as records of episodic swarms of very small, fluid-enabled, seismically detectable events. With increasing resolution of geophysical data sets, new discoveries will likely open other avenues for geological interpretations, and improve the geometrical constraints to the models and ideas arising from studying the complexity revealed by exhumed subduction thrusts in outcrop. In closing, we therefore emphasize that the rocks point to a complex time-evolution of subduction thrust rheology arising from the interplay between solid-state reactions and the release and flow of fluids, in an environment where permeability creation and destruction compete (Fig. 9).

This work was supported by a Marsden grant GNS1204 (to Ellis), European Research Council Starting Grant 715836 “MICA” (to Fagereng), and separate University of Capetown Research Development Grants to Fagereng and Diener. We acknowledge comments from two anonymous reviewers and the editorial assistance of Asuka Yamaguchi and Tim Byrne, which improved the final version of the paper.

This work was supported by a Marsden grant GNS1204 (to Ellis), European Research Council Starting Grant 715836 “MICA” (to Fagereng), and separate University of Capetown Research Development Grants to Fagereng and Diener. We acknowledge comments from two anonymous reviewers and the editorial assistance of Asuka Yamaguchi and Tim Byrne, which improved the final version of the paper.

This work was supported by a Marsden grant GNS1204 (to Ellis), European Research Council Starting Grant 715836 “MICA” (to Fagereng), and separate University of Capetown Research Development Grants to Fagereng and Diener. We acknowledge comments from two anonymous reviewers and the editorial assistance of Asuka Yamaguchi and Tim Byrne, which improved the final version of the paper.

The two-dimensional thermo-fluid-mechanical model has a total length of 220 km (though only the left-hand 125 km are shown in Figs. 5 and 6). The finite-element code SULEC was used to model frictional deformation for a submarine wedge (Ellis et al., 2011, 2015; Buiter and Ellis, 2012). Pressure-sensitive frictional plastic flow was based on a Coulomb yield criterion and velocities solved on an Eulerian calculation mesh, while materials, strain, and other Lagrangian quantities were tracked using tracer particles that moved through the mesh. Boundary conditions representing the subduction velocity component normal to the margin (1 cm/yr) were applied at the right-hand boundary, along the base, and along the lowermost 12.6 km of the left-hand boundary. Above this, the left-hand boundary horizontal velocity was prescribed to be zero, while vertical velocity was free, as indicated by circles on Figure 5. These boundary conditions mean that outward material flux was 80% of incoming flux, so that the top 3.2 km of material accreted to the wedge.

The weight of overlying seawater was imposed on the model as a force boundary condition on the top surface. Fluid pressure was initially assumed to be hydrostatic throughout the wedge, but over time, high fluid pressures could develop, which we show in Figure 6 by plotting the fluid pressure ratio λ,

which is the modified Hubbert-Rubey fluid pressure ratio corrected for water depth, D. Here, p_f is fluid pressure, σ_z is the vertical normal stress, g is gravity, and ρ_w is fluid density (Hubbert and Rubey, 1959; Davis et al., 1983; Dahlen, 1984).

Primary porosity (i.e., from initial pore space in the sediment, which reduces due to loading and compaction) decreases exponentially as a function of the effective mean stress (pressure P_eff = rock pressure – fluid pressure) computed at every time and location in the model, where primary porosity n is given by:

and n₀ is the initial porosity for zero effective stress (given here as 0.4), with P_ref being a reference pressure (here prescribed as 100 MPa). A lower limit of 0.01 was imposed on porosities calculated this way. Note that effective pressure includes the effects of tectonic stress, so high differential stress can lead to an increase in effective pressure and a decrease in porosity, while high fluid pressures can maintain high porosity. We distinguish between primary and secondary porosity, where secondary porosity results from frictional deformation and/or hydrofracturing of the rock. Primary porosity as defined here is not allowed to recover; that is, if high effective pressures reduce primary porosity, and then later these high pressures are reduced, primary porosity retains memory of its lowest value and does not increase again.

In order to compute release of fluid from clay dehydration in the subducting and accreting sediment, the numerical code also solves the steady-state heat-flow equation, taking material advection into account. No radiogenic heating or shear heating is applied, and heat advection by fluids is neglected. The assumption of negligible shear heating is consistent with constraints requiring low strength along the detachment (as noted earlier). We assume a temperature of 160 °C at 20 km depth. Thermal conductivity was computed from the porosity using the relationship:

where n is porosity, and k is the geometrically averaged thermal conductivity derived from thermal conductivity of seawater (k_s = 0.67 W m^–1 K^–1) and of rock grains (k_r = 2.8 W m^–1 K^–1; Pecher et al., 2010). The seafloor is prescribed as a constant of 4 °C. Because we modeled the shallow wedge and not the entire subduction system, these calculations only encompass the simplest approximation of a subduction thermal structure. The thermal gradient is ~8 °C/km, consistent with Figure 4A, giving an interface temperature that increases from ~30 °C outboard to the wedge to 76 °C at the right-hand edge of the model.

Material frictional strength is a function of effective pressure, where the oceanic slab is represented by a strong frictional layer with Byerlee-type friction coefficient µ = 0.84 and a cohesion of 1 MPa. Subducting sediment and wedge sediment are prescribed as µ = 0.53 and a cohesion of 0.1 MPa, while the frictionally weak décollement is given µ = 0.09 and a cohesion of 0.01 MPa. The low value of décollement friction produces thickening at the toe of the accretionary wedge. The resulting material flow and porosity field are shown in Figure 5B. This model predicts a mostly stable wedge with thickening and accretion concentrated at the toe.

Accreting and subducting sediment loses water by compaction, because net reduction in porosity leads to expulsion of fluids. We combined the spatial gradient in the porosity field with the material flow field from the mechanical model (Fig. 5B), assuming that the local time derivative of porosity is negligible (e.g., Bethke, 1986; Bekins and Dreiss, 1992; Wang, 1994). We calculated fluid release from porosity reduction as:

(Bekins and Dreiss, 1992), where n is porosity at each point, u and v are horizontal and vertical velocities derived from the mechanical model, and

Following Saffer et al. (2008), we used the kinetic reaction model of Pytte and Reynolds (1989) to compute the transformation of smectite in the wedge and subducting sediments, using the steady-state temperature field, and the material flow field of the mechanical model (Fig. 5B). The material derivative (the change in smectite content of a tracked material point with time, assuming constant potassium content) is given by:

where S is the mole fraction of smectite in the mixed-layer illite-smectite (I/S) clay component of the sediment, T is temperature in Kelvin, and C₁ and C₂ are constants (3.9 × 10⁹ s^–1 and –19,089 K, respectively; derived from Pytte and Reynolds, 1989, and incorporating the activity of potassium into constant C₁ assuming equilibrium between albite and K-feldspar). If the wedge geometry changes very slowly, we can assume an approximate steady state, so that the local (Eulerian) time derivative is zero, so that the material derivative is balanced by the advective rate of smectite change:

where u and v are the two-dimensional velocities. We solved this equation iteratively for molar fraction S and the associated steady-state fluid release. We assumed an initial smectite weight fraction F at shallow depths of 20% in incoming sediments (Pecher et al., 2010) and 10% in the slab. We took an initial mole fraction of S in mixed I/S of 1.0, which decreased with depth (temperature). That is, we assumed that the smectite-clay molar fraction in incoming sediment and slab had equilibrated with the local temperature field for 10 m.y. prior to entering the margin, so that the smectite-illite reaction was already under way (as a function of temperature) when material entered the trench, yielding a reduction in the smectite mole fraction with depth in the incoming section.

Pseudosections and fluid contents were calculated using the THERMOCALC software (version 3.45; Powell and Holland, 1988) and an updated version of the Holland and Powell (2011) thermodynamic data set (file tc-ds62.txt, created 2 June 2012). The activity-composition models used were those presented by Holland and Powell (2011) for epidote; White et al. (2014) for biotite, chlorite, and muscovite; and Green et al. (2016) for actinolite, glaucophane, hornblende, and omphacite. Albite, lawsonite, aqueous fluid (H₂O), quartz, and sphene were taken as pure phases. Calculations were performed in the idealized Na₂O-CaO-K₂O-FeO-MgO-Al₂O₃-SiO₂-TiO₂-O (NCKFMASHTO) chemical system, with the compositions of GLOSS (Plank and Langmuir, 1998) and MORB (Sun and McDonough, 1989) converted to (in mol%): SiO₂: 52.47, TiO₂: 1.05, Al₂O₃: 9.10, Fe₂O₃: 1.47, FeO: 5.21, MgO: 12.71, CaO: 12.21, Na₂O: 2.61, K₂O: 0.23, and SiO₂: 68.57, TiO₂: 0.55, Al₂O₃: 8.22, Fe₂O₃: 1.48, FeO: 2.14, MgO: 4.33, CaO: 7.47, Na₂O: 2.76, and K₂O: 1.52, respectively. The calculated fluid content in each mineral assemblage was output by THERMOCALC as mole fractions, with each phase normalized to one oxide sum total to approximate vol% fluid released.