Subduction interfaces are loci of interdependent seismic slip behavior, fluid flow, and mineral redistribution. Mineral redistribution leads to coupling between fluid flow and slip behavior through decreases in porosity/permeability and increases in cohesion during the interseismic period. We investigate this system from the perspective of ancient accretionary complexes with regional zones of mélange that record noncoaxial strain during underthrusting adjacent to the subduction interface. Deformation of weak mudstones is accompanied by low-grade metamorphic reactions, dissolution along scaly microfaults, and the removal of fluid-mobile chemical components, whereas stronger sandstone blocks preserve veins that contain chemical components depleted in mudstones. These observations support local diffusive mass transport from scaly fabrics to veins during interseismic viscous coupling. Underthrusting sediments record a crack porosity that fluctuates due to the interplay of cracking and precipitation. Permanent interseismic deformation involves pressure solution slip, strain hardening, and the development of new shears in undeformed material. In contrast, coseismic slip may be accommodated within observed narrow zones of cataclastic deformation at the top of many mélange terranes. A kinetic model implies interseismic changes in physical properties in less than hundreds of years, and a numerical model that couples an earthquake simulator with a fluid flow system depicts a subduction zone interface governed by feedbacks between fluid production, permeability, hydrofracturing, and aging via mineral precipitation. During an earthquake, interseismic permeability reduction is followed by coseismic rupture of low permeability seals and fluid pressure drop in the seismogenic zone. Updip of the seismogenic zone, there is a post-seismic wave of higher fluid pressure that propagates trenchward.

Large earthquakes in subduction zones are controlled by the rupture of asperities or patches of the plate interface with high shear strength (e.g., Thatcher, 1990) (Fig. 1). Some asperities may reflect roughness elements along the interface, such as seamounts on the subducting plate (Abercrombie et al., 2001; Husen et al., 2002; Bilek et al., 2003; Collot et al., 2017), but seamounts also likely play a role in limiting the area of some ruptures (Wang and Bilek, 2011). There is a positive correlation between great earthquakes and smooth plate interfaces (Scholl et al., 2015; Lallemand et al., 2018; van Rijsingen et al., 2019), so heterogeneity in strength along large, unsegmented stretches of the interface can develop during the interseismic period independent of roughness potentially due to geochemical processes in deforming sediments. Between large ruptures, there are interseismic changes in the physical properties of the plate boundary footwall due to deformation accommodated by mass redistribution and reduction of fracture porosity (Fisher et al., 2019a; Fisher et al., 2019b).

Mineral redistribution and fracturing along the interface impacts both fluid flow (by generating variations in crack porosity, permeability, and fluid pressure) and earthquake physics (by generating variations in strength). Therefore, there is potential for coupling between fluid and chemical processes with positive and negative feedbacks that modulate the mechanics of slip behavior through the seismic cycle (e.g., Bebout and Penniston-Dorland, 2016; Hooker and Fisher, 2021). In the seismogenic zone, these feedbacks may allow episodic movement of porosity waves composed of fluid generated at depth to pass through via fault-valve behavior (Sibson, 1992; Cox, 2005). This paper reviews the current state of knowledge of this complex system of fluid production, fluid flow, geochemical healing, mineral redistribution, and plate boundary slip from a field-based perspective.

The Subduction Top-to-Bottom American Geophysical Union Monograph in 1996 was published at a time when geophysical surveys defined the architecture of many convergent plate boundaries and their overriding accretionary prisms. In the context of Coulomb wedge mechanics, these observations argued for a weak basal décollement with a stronger wedge interior (Davis et al., 1983; Dahlen, 1984). By this time, many studies of on-land accretionary complexes had characterized the structures and microstructural fabrics of fossil subduction boundaries and provided evidence of critical roles for fluids, crack-seal deformation, and pressure solution (Fisher, 1996). It was apparent that an underthrusting sediment pile experiences a compactive strain path (Karig, 1990) and that fluids are released along this path through pore space reduction at shallow levels (e.g., at the deformation front) and thermally controlled dehydration reactions at greater depths (Moore and Vrolijk, 1992). Negative polarity reflections along the décollement in seismic profiles (Moore et al., 1990, 1995a), bulk density and fluid pressure anomalies for the material along the interface (Moore et al., 1995b), and the low taper of accretionary prisms (Davis et al., 1983) were consistent with locally high fluid pressure and low basal stress. Other observations such as low-chloride anomalies (Kastner et al., 1991; Bekins et al., 1995) and evidence for large fluid rock ratios (Bebout, 1991) argued for the décollement as a major conduit for fluids derived from metamorphic dehydration reactions. Thus, our understanding of the fluid-connectivity of the seismogenic system was rooted in the characteristics of time-averaged behavior, and it was not known how the spatial and temporal heterogeneity in the deformation-fluid flow system could play a role in slip instabilities associated with the seismic cycle.

Over the last two decades, there have been numerous advances that call for a reevaluation of our understanding of processes along the subduction interface. (1) From top to bottom, there is a range of plate boundary slip behaviors that shows some systematic variation with depth and requires a new understanding of earthquake physics that include tsunami earthquakes (Bilek and Lay, 2002), slow slip events (Ito et al., 2013; Saffer and Wallace, 2015; Araki et al., 2017), earthquake supercycles (Sieh et al., 2008; Goldfinger et al., 2013; Philibosian et al., 2017), creep, and episodic tremor and slip (Miller et al., 2002; Rogers and Dragert, 2003). Slip behavior varies downdip and along-strike at a given margin and across margins due to variations in upper plate rigidity (Sallarès and Ranero, 2019) and different thermal structure, incoming plate roughness, and subduction inputs (Beroza and Ide, 2011). (2) Great earthquakes have challenged our perspectives with regard to the seismogenic zone, including the occurrence of coseismic slip to the trench (Ide et al., 2011; Chester et al., 2013) due in part to the weakness of smectite-rich sediments at high velocities (Kameda et al., 2015). (3) We have developed a greater appreciation for the potential roles of incoming plate roughness on seismic behavior and plate boundary architecture from 3-D seismic observations and multibeam bathymetry (von Huene et al., 1995; Bangs et al., 2004; Bangs et al., 2016). (4) Deformation experiments have produced a range of observed slip behaviors and have documented the role of temperature on rate-state friction frameworks (Nakatani and Scholz, 2004; Nakatani and Scholz, 2006; den Hartog et al., 2012, den Hartog and Spiers, 2013). (5) There are now observations of the in situ stress state and the tectonic setting of deformation fabrics that form along the subduction interface from drilling along active plate margins (Kitajima et al., 2012; Huffman et al., 2016; Brodsky et al., 2017; Brodsky et al., 2020).

In this review, we focus on the Kodiak accretionary complex and the Shimanto belt of Japan because these terranes expose rocks that were accreted over tens of millions of years along convergent plate boundaries that are currently active. Ten regionally extensive mélanges from these areas were accreted at depths of seismogenesis during subduction of a sediment-rich plate (Vrolijk et al., 1988; Taira et al., 1988) similar to the Cascadia margin today. These observations are interpreted in the context of active convergent margins and the features of geophysical surveys, geodetic observations, and lab deformation experiments. The objective is to use geologic observations from ancient accretionary complexes to evaluate the different components of the fluid mechanical system and to consider where knowledge is needed to make quantitative predictions of slip behavior that can be compared with spatial and temporal distributions of earthquakes and results of monitoring along active convergent margins.

In the following sections, we first describe subduction mélange from ancient accretionary complexes of Japan and Alaska to demonstrate that these rocks provide a record of deformation kinematics in the seismogenic zone during underthrusting beneath the overriding plate. We then synthesize observations of microstructures in these rocks to develop an interpretation for the operative grain scale deformation mechanisms. These observations are used to define the rheology during underthrusting at temperatures typical of seismogenesis and to argue that most of this deformation occurred during the interseismic period as a consequence of diffusive mass transfer at rates that require a slip deficit. Based on this conclusion, a kinetic model of mineral redistribution is applied to estimate sealing rates for fracture networks in the footwall of the subduction interface. Given that the structural history of the underthrusting plate primarily reflects the interseismic period, we consider where coseismic slip is represented in the record of these accreted packages. Finally, we frame the discussion in terms of a conceptual framework that incorporates the processes that occur in subduction mélange, which include potential coupling of fluid flow and deformation. Here we highlight the need for a process-based understanding of rheology and porosity/permeability along the subduction interface that is rooted in observations of exhumed subduction interfaces; such an understanding can yield testable predictions for active margins.

Stratal disruption and development of a block-and-matrix fabric characteristic of mélange can occur through a variety of processes typical of convergent plate boundaries such as noncoaxial strain associated with faulting, diapirism, near-surface downslope creep, and olistostromal mass wasting (Cowan, 1985). However, there are some features of mélange belts that are diagnostic of underthrusting of sediments in the footwall of the subduction interface. We focus many of our observations on the Shimanto belt in Japan and the Kodiak accretionary complex in Alaska (Fig. 2)—rocks that include regional belts of mélange exhumed in the forearc of active subduction zones in the northern Pacific.

These subduction mélanges are regionally extensive, planar zones of disruption parallel to the boundaries of accreted packages. The Ghost Rocks mélange is continuous for tens of kilometers along the Kodiak archipelago (Fisher and Byrne, 1987), and the Uyak mélange extends parallel to the margin for hundreds of kilometers when linked with the McHugh Complex on Kenai Peninsula (Kusky and Bradley, 1999). Underthrust mélange units from the Shimanto belt also can be regionally correlated (Taira et al., 1988). In all of these cases, there is evidence in rock fabrics for a compactive strain path that includes an evolution from granular flow and stratal disruption of partially lithified sediments to localized deformation of lithified rocks that includes pressure solution in mudstones and growth of veins within mélange blocks (Fisher and Byrne, 1987).

The stratigraphy is pervasively disrupted on length scales of centimeters to meters, but the fabric of the mélange is structurally coherent with systematic variations in the shape and orientation of blocks within the mélange and a foliation that is defined by scaly fabric in mudstones and block shape-preferred orientations (Figs. 3A3C). The Ghost Rocks mélange has sandstone blocks with long axes that are parallel to the intersection lineation formed by an anastomosing web structure of cataclastic shears that break up the beds (Byrne, 1984; Fisher and Byrne, 1987). The Makimine (Kitamura and Kimura, 2012) and Uyak (Byrne and Fisher, 1990) mélanges have elongate blocks that are preferentially aligned with their long axes downdip and parallel to a stretching lineation in the adjacent mudstones.

Despite partial to complete disruption at the bed scale, there is typically a ghost stratigraphy that is diagnostic of oceanic crust from an underthrusting plate that is composed, from bottom to top, of altered, basalt-dominated mélange (oceanic crust); pelagic sediments and chert-dominated mélange (sedimentation in the deep ocean) and turbidites (trench-fill); or sand blocks in a mudstone matrix (Fisher and Byrne, 1987; Taira et al., 1988; Byrne and Fisher, 1990; Ujiie, 2002; Ikesawa et al., 2005; Kitamura et al., 2005; Kimura et al., 2012; Hashimoto et al., 2012). The ghost stratigraphy in these cases is imbricated during duplex accretion (e.g., Fig. 1, (Sample and Fisher, 1986; Byrne and Fisher, 1990; Kimura et al., 2012), which indicates that the mélange fabrics crosscut by these thrust faults were formed during underthrusting in the footwall of the low-angle subduction interface. There is macroscopic and microscopic evidence for noncoaxial strain in these zones with asymmetric tails on blocks and shear bands that crosscut the mélange fabric and a kinematic frame consistent with the subduction, which is typically layer-parallel shear and downdip or oblique underthrusting (Fisher and Byrne, 1987; Ujiie, 2002; Kimura et al., 2012; Kitamura and Kimura, 2012). Kimura et al. (2012) describe R1 and R2 Riedel shears in the Mugi mélange; these features are also observed in laboratory experiments on gouge (Haines et al., 2013).

Importantly, the subduction mélange belts line the footwall of major fault boundaries. In the Kodiak accretionary complex, the Uganik Thrust juxtaposes the Uyak Complex (Connelly, 1978) in the hanging wall with the Kodiak waterfall bay mélange in the footwall (Fisher and Byrne, 1987; Rowe et al., 2009). In the Shimanto belt, there are also regional faults with footwall mélange (Taira et al., 1988), which include: (1) the Goshikinohama Fault, which juxtaposes the coherent Susaki Formation in the hanging wall with the Yokonami mélange in the footwall (Hashimoto et al., 2012); (2) the Minami-Awa Fault (Kitamura et al., 2005; Kimura et al., 2012), which juxtaposes the coherent sedimentary rocks of the Hiwasa and Kainan Formations in the hanging wall and the Mugi mélange in the footwall; and (3) the “roof thrust” that juxtaposes the coherent Nonokawa Formation in the hanging wall with the Okitsu mélange in the footwall (Ikesawa et al., 2003).

Thermobarometric estimates derived from these mélange belts are variable, but they broadly agree with the temperatures inferred for the seismogenic zone of active convergent margins (150–350 °C) (Fig. 4). At this low metamorphic grade, temperature is constrained using stable isotope geothermometry (Brantley et al., 1997), vitrinite reflectance (Sample and Moore, 1987; Sakaguchi, 1996; Sakaguchi, 1999; Ikesawa et al., 2005; Kondo et al., 2005; Rowe et al., 2009), fluid inclusions (Vrolijk, 1987; Vrolijk et al., 1988; Brantley et al., 1997; Raimbourg et al., 2015; Hashimoto et al., 2003; Matsumura et al., 2003), and Raman spectroscopy of carbonaceous material (Rowe et al., 2009; Raimbourg et al., 2021). Fluid inclusions in veins are consistent with immiscibility and with H2O-rich and either CH4-rich or CO2-rich inclusions (Vrolijk et al., 1988; Brantley et al., 1997; Hashimoto et al., 2012). Salinity of H2O-rich inclusions varies between two endmembers: low-chloride inclusions and inclusions with salinity closer to seawater (Raimbourg et al., 2015). These end members occur in quartz with disparate cathodoluminescence, which suggests that fluid associated with precipitation within veins may be derived both locally and from external sources by advection (Raimbourg et al., 2015); this interpretation is consistent with some geochemical observations (Yamaguchi et al., 2012; Cerchiari et al., 2020). The veins in higher temperature mélange correspond with more homogeneous fluid compositions (Raimbourg et al., 2021). Thus, fluid composition and source can vary through time and with depth along the interface, but there is likely fluid derived both externally and internally due to temporal variations in the fluid flow system through the seismic cycle (Cerchiari et al., 2020) due to local redistribution during the interseismic period, and advective flow of fluid over larger distances in the footwall following coseismic ruptures that break down seals and link fluid reservoirs (e.g., Sibson, 1992; Cox, 2005).

It is difficult to estimate the active thickness of the subduction mélange zone (i.e., the incremental strain distribution) because only the finite strain is observed, but Rowe et al. (2013) argue, based on the variation in mélange thickness with increasing accretion temperature, that thickness of the active zone of shear increases to ∼100–350 m at ∼1–2 km below seafloor, and this thickness is maintained down to a depth of ∼15 km. These thickness estimates are in accordance with those derived from seismic surveys of active margins (Abers, 2005).

Microstructures: What Are the Dominant Deformation Mechanisms?

Deformation fabrics within the mélange reflect rheological heterogeneity between mudstones and isolated blocks of altered basalt, chert, and sandstone (Fagereng and Sibson, 2010), and this heterogeneity could play an important role in the deformation behavior of rocks along the interface (Fagereng and Beall, 2021). In the mudstones, there is a pervasive, scaly fabric of anastomosing slip surfaces that, combined with the preferred orientation of blocks, defines the mélange foliation (Figs. 5A5B). Individual scaly folia are polished and commonly striated. The noncoaxial strain that characterizes the kinematics of the mélange is largely accomplished through slip that is distributed within this web-like network of microfaults (Fig. 5A). Scaly fabrics are ubiquitous features of subduction mélange worldwide (Vannucchi, 2019) and are found in faults as shallow as the toe of the wedge (Maltmann et al., 1993). At shallow depths, the development of each scaly micro fault involves collapse of the phyllosilicate fabric; this process leads to strain hardening of scaly folia relative to the undeformed mudstone and densification and widening of this pervasive network (Moore and Byrne, 1987).

In the subduction mélanges studied that were underthrust to temperatures characteristic of the seismogenic zone, the scaly fabric is characterized by narrow bands of silica depletion relative to wall-rock lithologies (Kawabata et al., 2007; Fisher et al., 2019a; Ramirez et al., 2021) (Fig. 6A). Detrital feldspars throughout the matrix underwent dissolution within scaly folia, and there is evidence for the growth of illite, and in some cases, additional Fe-Ti phases (Ramirez et al., 2021). Shear failure of scaly folia at temperatures typical of the seismogenic zone is followed by low-grade metamorphic reactions, dissolution along scaly microfaults, and depletion of fluid-mobile chemical components. Elemental analysis of Mugi and Makimine scaly fabrics elucidates a record of leaching of fluid-mobile elements (Si, Na, Ca, large ion lithophile elements [LILEs]) and relative enrichment in fluid-immobile elements such as the rare earth elements (REEs), high field strength elements (HFSE), and Ti (Fig. 6B), which is consistent with slip dissolution, compactive strain, and volume loss in mudstones (Ramirez et al., 2021). Increasing temperature within the underthrusting sediment pile, from ∼150–300 °C, is depicted in X-ray diffraction patterns of mudstones by a prograde decrease in the proportion of smectite, changes in the proportion of illite polytypes, and the increasing abundance of chlorite (Ramirez et al., 2021). These observations show that noncoaxial strain in subduction mélange is accomplished through a combination of quartz and feldspar dissolution and metamorphic equilibration of phyllosilicates.

Undeformed mudstone phacoids in the mélange contain quartz and albite, whereas scaly folia are microgranular aggregates of phyllosilicates that include smectite, illite, and chlorite. These phase relations indicate the formation of sheet silicates during the dissolution of albite with development of the scaly fabric through reactions such as:

as proposed by Beach (1979) and Fry (1982) and documented in Ramirez et al. (2021).

The coherent blocks within the subduction mélange, which are stronger than the mudstones, are deformed by opening mode cracking accompanied by mineral precipitation, which forms veins with fractures oriented at high angles to the mélange fabric (Fig. 5B). The orientation of veins is consistent with low shear stress parallel to the mélange fabric (Fagereng et al., 2010; Yamaguchi et al., 2011). Vein infill forms up to 30–75% of block volumes, and the concentration of veins is significantly greater in mélange blocks than in coherent rocks outside of the mélange (e.g., Rowe et al., 2009). Veins are predominantly composed of quartz and subordinate calcite and albite; these minerals contain fluid-mobile chemical components such as silica, sodium, and calcium that are depleted from scaly folia during noncoaxial strain. Therefore, a large proportion of the precipitated minerals within blocks of the mélange are derived from the dissolution and local transport of chemical components from the surrounding mudstone matrix. Characteristic distances between scaly fabric (sites of mineral dissolution) and veins (sites of precipitation) are on the order of 10−3 m to 10−1 m.

Vein microstructures record a variety of crack behaviors: (1) composite quartz-albite-calcite veins devoid of planar inclusion bands (Fig. 7A), (2) quartz with planar syntaxial fluid inclusion bands oriented parallel to the walls of the vein and perpendicular to the long axis of the blocks (Fisher et al., 2019a; Fisher et al., 2019b) (Fig. 7B), and (3) quartz with antitaxial solid inclusion bands—either of wall rock or chlorite seeded from the wall rock margin—that can be continuous but also show spotty bands that are consistent with partial healing (Fisher and Brantley, 2014) (Fig. 7C). Vein microstructures 1 and 2 are common in the subduction mélanges that formed at lower temperature within the seismogenic zone (150–280 °C), whereas examples of the antitaxial microstructures are mostly observed in subduction mélange accreted near the downdip end of the seismogenic zone (280–350 °C) (Fisher and Brantley, 2014; Ujiie at al., 2018; Raimbourg et al., 2021). Microstructure 1 is interpreted as evidence for persistence of a substantial fracture porosity over the lifetime of veins, whereas microstructure 2 indicates growth of quartz from the walls of fractures inward, with repeated cracking and sealing, particularly in the necks of blocks that are separating. The necks of these attenuated blocks can record hundreds of crack-seal events (Fisher et al., 2019a). The repeated cracking and sealing within the seismogenic zone may record the passage of fracture porosity waves that heal in their wake (Cox, 2005). Such porosity waves may be expected during the rupturing of seals and the rehealing that accompanies slip events within the seismogenic zone.

The microstructural observations are significant because they indicate that the footwall of the subduction interface at depths of the seismogenic zone has a crack porosity that fluctuates due to the interplay of cracking and precipitation. Moreover, the pervasive dissolution along scaly microfaults in the mudstones, contemporaneous with cracking and precipitation within the stronger blocks, implies a system of local redistribution of mobile chemical components as the mechanism for accommodating deformation. This mechanism of deformation depicts a rheology defined by plastic shear failure followed by pressure solution slip, strain hardening, and the development of new shears in adjacent undeformed material.

The plate boundary deformation in subduction mélange does not account for coseismic slip and an abrupt reduction in slip deficit. Instead, slow deformation in the subduction mélange is related to unrecoverable deformation in the rocks adjacent to the plate interface in the interseismic period. Given a typical slab-top geothermal gradient for a warm subduction zone like Cascadia (van Keken et al., 2018) and a displacement rate of 30 km/m.y., it would take ∼1 m.y. to underthrust from 250 °C to 350 °C along the interface. Such a time interval could include thousands of large earthquakes with several meters of slip. If the active shear zone in the footwall were 100 m wide, a steady ductile strain rate of ∼10−11/s would be required to accommodate relative plate motions. A strain rate in this shear zone of 10−13/s would produce a γ of ∼3 over 1 m.y. In other words, a large visible strain can occur at rates that lead to a slip deficit but are difficult to parse from the plate coupling percentages inferred from GPS campaigns. This back-of-the-envelope calculation highlights the need for quantitative assessments of strain related to diffusive mass transfer from fossil subduction zones in assessing the processes that lead to slip instabilities along active margins.

Because the predominant deformation mechanism of diffusive mass transfer may be inadequate for facilitating plate boundary shear at rates that are compatible with relative plate motions, there is a slip deficit that increases until failure of the interface. Between earthquakes, the viscous coupling and dissolution during non-coaxial strain, combined with volume loss in the mudstones and precipitation in the sandstone blocks, are also likely to impact the permeability. Compactive strain in the mudstones and precipitation in cracks within the blocks are typical of Ruina healing in rate-state frictional mechanics (Marone, 1998); strength increases (e.g., cohesion) with log time due to increasing contact area that is facilitated by very slow deformation between slip instabilities.

Kinetics of Silica Redistribution: How Fast Do Cracks Heal in the Footwall of the Interface?

There are a variety of different models for the kinetics of crack seal deformation, but they generally differ in the driving force and transport mechanism for silica from sites of dissolution—in this case scaly fabric or from dissolved silica from deeper in the slab—to sites of precipitation (i.e., in cracked blocks). The kinetics are described as being driven by: (1) a drop in crack fluid pressure with diffusion from a local source (Fisher and Brantley, 1992; Renard et al., 2000; Fisher and Brantley, 2014; Fisher et al., 2019a) or advection down a fluid pressure gradient from an external source (e.g., Ujiie et al., 2018) or (2) a difference in mean stress between the weak, scaly mudstone and the stronger blocks (Robin, 1979; Fisher et al., 2019a).

In the case of a drop in fluid pressure, estimates for the magnitude of this drop range from 10 MPa (Fisher and Brantley, 1992) to the difference between lithostatic and hydrostatic pressure (Renard et al., 2000) to 100 MPa (Fisher and Brantley, 2014) to a few hundred MPa (Ujiie et al., 2018). The larger values are in some cases based on the variability in estimates of fluid pressure from fluid inclusions in veins (Vrolijk, 1987; Ujiie et al., 2018). Dilation and increased effective stress associated with hydrofracturing could stabilize a slip instability (e.g., Rubin, 2008), so cracking and sealing could plausibly reflect the passage of porosity waves at the front of a slow slip instability with healing in the wake of the stress pulse. This has led some workers to argue for healing at a time scale of days to weeks (Fisher and Brantley, 2014; Ujiie et al., 2018), which is compatible with slow slip events, albeit with large pressure drops that should cause collapse of the permeability. Given the potential for a large increase in effective stress, it is unclear how cracks with apertures of tens of microns (i.e., the spacing of crack-seal bands; e.g., Fisher et al., 1995; Fagereng et al., 2011; Remitti et al., 2012) can remain open if the fluid pressure drops by hundreds of MPa as this amount is larger than the tensile strength of the blocks (< 25 MPa; Zhang, 2002). Moreover, geophysical observations are consistent with fluid pressure drops that are one to two orders of magnitude smaller during slow slip events (Gosselin et al., 2020; Warren-Smith et al., 2019). Finally, the models based on large pressure drops do not account for the transience of the pressure changes and the restoration of fluid pressure as quartz precipitates in the crack, which may lead to a diminishing driving force and sealing rate through time (Fisher et al., 2019a). The size of pressure drops associated with slip instabilities remains an outstanding question, but it is difficult to explain drops in excess of a few tens of megapascals.

An alternative driving force for the kinetics of silica redistribution is related to the heterogeneity in strength between the weak mudstone and the stronger blocks, which can result in differences in mean stress of tens of MPa (Fisher et al., 2019a). This driving force is constant and does not diminish as the fracture aperture decreases; the model predicts healing of a 10-micron-wide fracture over hundreds of years. This number is relevant as a healing rate for the subduction interface particularly because partial healing of fracture porosity could have a large impact on the strength of the fault zone. In such a model, the rate of crack healing is determined by the competing rates of mineral interface kinetics (both dissolution and precipitation) and grain boundary diffusion. Whereas experimental constraints exist for the former process (e.g., Dove, 1994; Rimstidt and Barnes, 1980), there are no experimental or empirical measurements of silica diffusivity along scaly microfaults. Experimental constraints on the diffusivity of silica through an H2O solution and a quartz mylonite differ by four orders of magnitude at temperatures relevant to the seismogenic zone (Watson and Wark, 1997; Farver and Yund, 1999). Future studies should focus on the energetics of silica diffusion through fault zone rocks.

For the model of Fisher et al. (2019a), the boundary between interface-controlled (i.e., mediated by rates of dissolution/precipitation) and diffusion-controlled (i.e., mediated by the grain-boundary diffusion rate) behaviors overlaps with the pressure-temperature conditions of most rocks exhumed from subduction zones (Penniston-Dorland et al., 2015; Raimbourg et al., 2018) and coincides with the thermal predictions of models for warm subduction zones such as Cascadia (van Keken et al., 2018).

Given the number of poorly constrained variables in the calculations, the rates of healing, which will vary systematically with depth and the thermal structure of the subduction zone, could be faster. The silica kinetics model of Fisher et al. (2019a) also does not account for the hastening effects of salinity (Dove, 1994), organic material (Bennett, 1991), phyllosilicates (Renard et al., 1997), or silica gels (e.g., Faber et al., 2014; Kirkpatrick et al., 2013) on crack-seal deformation. These factors could be important if the rates are controlled by interface kinetics, which is likely for many of the world's subduction zones with a warmer thermal structure (Fisher et al., 2019a). Independent of these influences, rates of thermally activated diffusive mass transfer will be enhanced with increasing depth and temperature along all subduction interfaces and along convergent margins that involve the slow subduction of young oceanic lithosphere. This dependence of subduction interface healing on geochemical processes is implied by conceptual views of a seismogenic zone that has a thermally activated updip limit across a narrow temperature range of 140–160 °C (Moore and Saffer, 2001).

Coseismic Slip along the Interface: Where Is the Record?

If interseismic deformation controlled by diffusive mass transport across an ∼100-m-wide zone occurs at rates inadequate for the accommodation of plate motions, where is coseismic slip accommodated in the rock record (e.g., see Rowe et al., 2011)? A logical argument in the Shimanto locales is that coseismic slip occurs within the narrow zone of cataclasite at the top of the mélange (e.g., Minami-Awa fault at the top of the Mugi mélange (Kitamura et al., 2005; Ujiie et al., 2007; Hashimoto et al., 2009), the Goshikinohama fault at the top of the Yokonami mélange (Fig. 8) (Hashimoto et al., 2012, 2014), and the Okitsu mélange with its associated roof thrust (Ikesawa et al., 2003). In the Kodiak case, fault rocks are observed along a crosscutting fault within the mélange (e.g., the Ghost rocks mélange, Rowe et al., 2005; Meneghini et al., 2010). These fault zones from Shimanto and Kodiak are narrow (<15 m thick) and display cataclastic rock fabrics and pseudotachylites that are indicative of frictional melting (Ikesawa et al., 2003; Hashimoto et al., 2012; Ujiie et al., 2007; Rowe et al., 2005; Meneghini et al., 2010). JFAST drilling through the plate boundary interval at the updip end of the 2011 Tohoku rupture shows a 5–15-m-thick zone of scaly fabric in a smectite-rich layer, but coseismic slip is likely localized on bounding faults (Kirkpatrick et al., 2015).

If the subduction mélange and the brittle fault zone that caps the mélange represent plate boundary deformation during different phases of the seismic cycle, then the different structures should be broadly contemporaneous at the resolution of modern geochronological techniques. In the case of ancient exhumed examples of Kodiak and Shimanto, the simple interpretation that these faults record the coseismic slip associated with underthrusting is not supported by K-Ar illite dating analysis of clay-sized particles in the fault zones; there is evidence for authigenic illite growth tens of millions of years after the timing of subduction (Tonai et al., 2016; Fisher et al., 2019c). One possible explanation for the young ages of illite growth is that the subduction interface, which is a discrete boundary between underthrusting sediments and an older forearc, is a weakness reactivated post-accretion as an out-of-sequence fault in response to later collisional encounters of the upper plate with incoming bathymetric features (Raimbourg et al., 2018). Such collisions have occurred with triple junction migration in the Miocene related to the Izu-Bonin arc in Japan (Kimura et al., 2014) and the Kula-Resurrection ridge in the Paleocene in the Kodiak accretionary complex (Moore et al., 1983). Under these circumstances, it could be difficult to see the evidence of coseismic underthrusting through the overprinting microstructure. Pseudotachylite from the earlier record of subduction may also be rare because of the factors that hamper preservation (Kirkpatrick and Rowe, 2013) or reduce frictional heating such as low shear stress and thermal pressurization (Schmitt et al., 2011).

Fluid Flow and Kinetic Healing: A Coupled System

The simple kinetic model for geochemical healing (Fisher et al., 2019a) implies that footwall deformation along the subduction interface can produce changes in rheology on interseismic timescales of hundreds of years and potentially much faster as a consequence of kinetically controlled mass redistribution. The changes of state include: (1) increases in cohesion through tectonically driven compaction and crack healing and (2) decreases in permeability through a reduction in fracture porosity and interconnectivity. Thus, the slip behavior is envisaged as the response of a coupled system with feedbacks between different components such as fluid production, fluid flow, hydrofracturing, and silica kinetics. When fluid is produced from metamorphic dehydration reactions faster than it can escape, the fluid pressure rises to a level that is ultimately limited by the tensile strength of the rock and the occurrence of hydrofracturing with rupture of seals, updip migration of a fracture porosity wave (Sibson, 1992; Cox, 2005), and precipitation of externally derived silica (e.g., Ujiie et al., 2018). Interseismic, slow noncoaxial strain in mudstones and local silica redistribution leads to the reestablishment of fluid overpressures and re-strengthening of the interface. This conceptual model leads to two different sources of silica and provides a framework for numerical models of the fluid flow system.

Recent numerical models have shown that the coupled hydromechanical models can be an important factor in explaining slip behavior. Sun et al. (2020) show that such a model can explain patterns of megathrust slip in relation to subducting seamounts. The observed variety of slip behavior along active subduction interfaces also can emerge in numerical models that consider interseismic compaction of fluid-filled pore space within a visco-elastic-plastic system (Petrini et al., 2020). We emphasize here that the observations of fossil subduction interfaces can be used to constrain potential feedbacks between mechanics and hydrology. These feedbacks have implications for subduction slip behavior.

Hooker and Fisher (2021) developed a hydromechanical model that is based on processes observed in the footwall of fossil interfaces such as hydrofracturing/crack sealing and the influence of these processes on cohesion and slip behavior within the seismogenic zone. This model (Hooker and Fisher, 2021) explores the impact of increases in cohesion on slip behavior by coupling a simple earthquake simulator with a numerical model for fluid flow to evaluate the impact of different fluid flow parameters on the temporal and spatial distribution of earthquakes. The numerical model for slip behavior is a cellular automaton that is composed of cells held in place against an overlying plate interface by friction, with static frictional strength exceeding the dynamic strength (e.g., Brown et al., 1991; Huang et al., 1992). Cells are connected updip and downdip with coil springs and laterally with leaf springs (Fig. 9A). In cross sectional view, leaf springs connect all of the cells to a rigid plate that moves steadily downdip at the rate of relative plate motion (Fig. 9B). Cells remain stationary as the spring-loaded system accumulates elastic strain until a critical failure stress, or the static strength, is reached. The failure of a cell imparts stress on its neighbors, which has the potential to cause failures of surrounding cells. This cascade of cell failures continues until there is no further failure and a rupture is complete. The passage of time is fixed in the model by the motion of the rigid subducting plate, so the ruptures are considered instantaneous in that reference frame.

Fisher et al. (2019b) allowed for healing of the interface in this 2-D, block-slider model. Healing was treated as a process of stochastic nucleation of asperity cells that increase in cohesion linearly with log time similar to the increase in strength of gouge in slide-hold-slide experiments (Nakatani and Scholz, 2006). The probability that a geochemical asperity will nucleate in the model is related to the activation energy for silica redistribution as is the rate of increase in cohesion of the asperity after it nucleates. The healing of the interface is thus described by a population balance equation with nucleation, strengthening, and failure of asperity cells in the model. The length scale in the model is imparted by thermal structure, and the distance downdip corresponds to the slab-top geotherm in thermal models of active subduction zones (van Keken et al., 2018).

The model leads to behaviors that are not observed in the absence of geochemically controlled healing such as an updip limit to the seismogenic zone, clustering of earthquakes during supercycles of buildup and release of elastic strain, and a Gutenberg-Richter size distribution of earthquakes. Moreover, the slope of cumulative frequency plots under these circumstances depends on the temperature structure of the interface in the models, and larger earthquakes are observed in warmer margins, which is also seen in comparisons between active margins with different thermal structure (Nishikawa and Ide, 2014; Fisher et al., 2019b). In the model, power law size distributions only occur within a range of activation energies for nucleation and strengthening (Fisher et al., 2019b).

Silica kinetics can influence earthquake patterns through thermally activated increases in strength, but what is the impact on the fluid flow system and what are the feedbacks between mineral redistribution, permeability, effective stress, and slip behavior? Hooker and Fisher (2021) address this question by combining the model for slip behavior along the subduction interface with a fluid flow system governed by the distribution of fluid production in combination with temporal variations in either the porosity or permeability within the cellular network. The model also responds elastically to temporal and spatial variations in effective stress on the plate interface. The block slider and the fluid flow models are coupled through the effect of cementation; increases in strength are accompanied by decreases in porosity and permeability. Rupturing of the interface resets the strength and the permeability to background levels.

The numerical simulations depict a subduction zone interface where the slip behavior is governed by the negative feedback between silica redistribution and permeability and the positive feedback between silica redistribution and cohesion. The reduction of permeability relative to fluid production (e.g., Saffer and Tobin, 2011) leads to a rise in fluid pressure that is limited by the criterion for the opening of hydrofractures, which increases permeability and lowers fluid pressure. The model results illustrate how these feedbacks may be manifest in naturally observed phenomena. Specifically, modeling changes to porosity and permeability that are imparted by cementation, along with varying fluid production at depth, result in changes to earthquake magnitude distributions, recurrence intervals, and fluid flow behavior. The coupled slip-fluid flow model (MEFISTO, Hooker and Fisher, 2021) suggests that features of the plate boundary plumbing system can control the rupture characteristics, with high hydraulic conductivity and low fluid production leading to a lower b value (i.e., a greater proportion of larger earthquakes) and a longer average time between earthquakes. This result is capable of explaining some of the statistical differences in earthquake records of segments of the Aleutian plate boundary with varied locking characteristics (Hooker and Fisher, 2021). Analysis of earthquake patterns along active convergent boundaries can thus provide insight into the interplay between silica kinetics, fluid production, and fluid flow.

In addition to the potential controls on slip behavior, coupled hydromechanical models (e.g., Sun et al., 2020; Petrini et al., 2020; Hooker and Fisher, 2021) also have implications for observables associated with the plate boundary fluid flow system. For example, Figure 10 shows 17 points at different positions downdip along a plan view of a MEFISTO model simulation for the subduction interface immediately following a rupture. Fluid pressure within these cells is monitored through time from before to after the development of the rupture (Fig. 10). In these simulations, we include the reduction in permeability due to the sealing of cracks, but we do not include the increase in fluid pressure that may occur when a cell experiences a reduction in porosity. Histories are red or blue depending on whether the cell has nucleated as an “asperity cell” or has the default values of static strength and permeability, respectively. The plotted cells form a transect along dip, and the brightness of each cell's time series is proportional to that cell's depth along the interface (Fig. 11). At depths within the seismogenic zone, the seals that hold back fluid from deeper levels are opened up, and there is a drop in fluid pressure (Figs. 11 and 12A). Cells at shallower levels are less likely to heal and experience lower pressure through the interseismic period, but these cells are subjected to a post-seismic rise in fluid pressure to a maximum pressure that decays through time (Figs. 11 and 12B). This post-seismic shallow pulse in fluid pressure reflects a wave of high fluid pressure that propagates updip after the seals of the seismogenic zone are ruptured and fluids generated at depth migrate into the shallow part of the interface.

The magnitudes of fluid pressure are highly dependent on factors that are poorly constrained such as the default permeability and the permeability change that accompanies healing, but this general pattern of behavior (i.e., greater pressure at depth, drop in fluid pressure at depth, and upward migration of a pulse of fluid pressure at shallow levels) is not sensitive to these values. Therefore, the record of fluid pressure differs from the seismogenic zone to the shallow parts of the interface (Fig. 12), and the monitoring of fluid pressure changes at shallow levels related to a rupture at depth holds the potential to calibrate parts of this system.

Ancient accretionary complexes expose regionally extensive zones of subduction mélange with systematic fabrics that record noncoaxial strain during underthrusting adjacent to the subduction interface at temperatures associated with the seismogenic zone. The subduction mélange is characterized by a block-in-matrix fabric with an anastomosing scaly fabric of polished or striated slip surfaces in mudstones and pervasive veining of sandstone, chert, or greenstone blocks.

The scaly fabric shows evidence for dissolution compaction during slip, with removal of mobile elements such as Si, Na, Ca, and LILEs and concentration of immobile elements such as REEs, HFSEs, and Ti. Veins in blocks are composed of quartz, calcite, albite, and chlorite—minerals that contain the mobile elements removed from scaly fabrics. Together, these observations indicate local dissolution-precipitation kinetics and linear viscous flow associated with unrecoverable deformation during the interseismic period. Models for silica kinetics indicate that this process of dissolution compaction and associated crack-seal deformation has the ability to modulate the fracture porosity and cohesion at time scales relevant to seismogenesis and potentially slow slip processes. Currently, there are few quantitative constraints on important parameters that control this system such as the diffusion rate and the length scale of local transport. Important questions about the coupled system of geochemical redistribution and fluid flow can only be answered through geologic investigations of fossil subduction interfaces. These questions include: what is the driving force for local redistribution of mobile elements from scaly fabric to cracks? What is the magnitude of volume loss in the scaly fabric? What is the mass balance between the volume loss in the mudstone and the volume gain through precipitation in stronger blocks in the mélange? Finally, what is the proportion of vein material derived through local dissolution-diffusion-precipitation relative to vein minerals derived through episodic fluid advection down geochemical potential gradients as a consequence of ruptures of seals?

Coseismic slip is likely restricted to <15-m-thick cataclastic fault zones at the top of the mélange that show rare examples of pseudotachylite. Rupture of these areas during earthquakes may allow the upward migration of a porosity wave that impacts the fracture porosity that persists tens of meters deep into ductilely deforming footwall. The cataclastic fault zones at the inferred interface typically produce K-Ar dates that are younger than the timing of subduction from other constraints, which we interpret as being due to reactivation within the overriding wedge long after accretion during collision with incoming bathymetric features.

A model for fluid flow and thermally activated strengthening along the subduction interface shows that coupling between processes of dissolution-compaction/crack healing and fluid production/fluid flow has an impact on observable features along active margins such as earthquake size distributions, interseismic time intervals, and fluid pressure variations related to seismic events. Thermally activated strengthening can produce supercycles of buildup and the release of elastic strain energy as well as temperature-dependent, Gutenberg-Richter size distributions. High permeability and low fluid production favor a greater proportion of larger events and longer time intervals between slip events. Reduction in permeability within the seismogenic zone due to interseismic deformation leads to a post-seismic rise in fluid pressure updip of the seismogenic zone that propagates trenchward due to the rupture of low permeability seals.

We thank Gaku Kimura for introducing D.M. Fisher to the Shimanto belt and Rina Fukuchi, Yoshi Hashimoto, and Asuka Yamaguchi for contributions to the field work at the Shimanto belt. D.M. Fisher is supported by National Science Foundation grant EAR-1524530 from the Tectonics Program. J.N. Hooker is supported as a fellow of the GDL Foundation. A.J. Smye acknowledges the support of the Sauermann family through receipt of the Slingerland Early Career endowment; A.J. Smye also acknowledges support from National Science Foundation grant OISE1545903. Whitney Behr and an anonymous reviewer contributed constructive comments. The paper benefited from group discussions with H. Savage, G. Hirth, D. Saffer, R. Skarbek, and C. Marone.

Science Editor: Shanaka de Silva
Guest Associate Editor: Gray E. Bebout
Gold Open Access: This paper is published under the terms of the CC-BY-NC license.