Constraints From Exhumed Rocks on the Seismic Signature of the Deep Subduction Interface.

Low Velocity Zones (LVZs) with anomalously high V p-V s ratios occur along the downdip extents of subduction megathrusts in most modern subduction zones and are collocated with complex seismic and transient deformation patterns. LVZs are attributed to high pore fluid pressures, but the spatial correlation between the LVZ and the subduction interface, as well as the rock types that define them, remain unclear. We characterize the seismic signature of a fossil subduction interface shear zone in northern California that is sourced from the same depth range as modern LVZs. Deformation was distributed across 3 km of dominantly metasedimentary rocks, with periodic strain localization to km-scale ultramafic lenses. We estimate seismic velocities accounting for mineral and fracture anisotropy, constrained by microstructural observations and field measurements, resulting in a V p/Vs of 2.0. Comparable thicknesses and velocities suggest that LVZs represent, at least in part, the subduction interface shear zone.

Introduction

Modern subduction zones exhibit a nearly ubiquitous Low Velocity Zone (LVZ) along the downdip extent of the megathrust that is 3–8 km thick and characterized by low velocities and high reflectivity, conductivity, Poisson's ratio (σ), and the corresponding V to V ratio (V /V ) (e.g., Audet & Bürgmann, 2014; Audet & Kim, 2016; Bostock, 2013; Y. Kim et al., 2014; Song et al., 2009; Toya et al., 2017) (Figures 1a and 1b), all consistent with near‐lithostatic pore fluid pressures (P ) (Audet et al., 2009; Bostock, 2013; Eberhart‐Phillips et al., 1989; Hansen et al., 2012; Peacock et al., 2011). Because near‐lithostatic P affects seismic velocities more than lithologic variations, the rock types that occupy the LVZ ‐ dominantly mafic (Audet & Schaeffer, 2018; Bostock, 2013; Hansen et al., 2012), dominantly sedimentary (Abers et al., 2009; Calvert et al., 2011; Delph et al., 2021), or a combination (Bostock, 2013; Delph et al., 2018) ‐ remain unclear from geophysical data. LVZs have been interpreted as overpressurized and relatively undeformed mafic crust sealed beneath a low‐permeability fault or narrow interface shear zone (Bostock, 2013; Calvert et al., 2020; Hansen et al., 2012; Kurashimo et al., 2013) (Figure 1c), or alternatively, as distributed viscous interface shear zones composed of mixed lithologies, including metasediments (Audet & Schaeffer, 2018; Calvert et al., 2020; Delph et al., 2018, 2021; Nedimović et al., 2003; Song & Kim, 2012) (Figure 1c).

Figure 1

The Low Velocity Zone (LVZ) (labelled margins in (a); example cross‐section at yellow circle of shear wave velocity structure from receiver functions in (b), after Delph et al., 2018) in modern subduction zones is collocated with transient seismic and aseismic slip (e.g., tremor frequency plot in (b), Delph et al., 2018). Base map produced in GPlates (Müller et al., 2018). (c) Schematics of competing LVZ interpretations.

Distinguishing between these endmember interpretations has important implications for rheological properties of deep subduction interfaces and associated seismic and transient deformation patterns. Transient seismic and aseismic slip ‐ for example, episodic tremor and slow slip, slow slip events, and low frequency earthquakes (LFEs) ‐ are very commonly collocated with LVZs (Audet & Kim, 2016; Calvert et al., 2020; Delph et al., 2018; Hirose et al., 2008; Song et al., 2009) (Figure 1b). Understanding these transient events, which factor into slip budgets and stress regimes related to megathrust earthquake probability (e.g., Rogers & Dragert, 2003; Wech et al., 2009), is crucial for hazard analysis, but both frictional slip along a discrete heterogeneous fault (e.g., Chestler & Creager, 2017; Ito et al., 2007; Lay et al., 2012; Luo & Ampuero, 2018; Shelly et al., 2006) or mixed brittle‐viscous deformation within a distributed shear zone (e.g., Beall et al., 2019; Behr et al., 2018; Hayman & Lavier, 2014; Kotowski & Behr, 2019; Tarling et al., 2019; Ujiie et al., 2018) (Figure 1c) are plausible sources. In addition, the composition and viscosity of the interface control coupling between the overriding and downgoing plates, contributing to, for example, slab velocities (Behr & Becker, 2018), upper plate topography (e.g., Delph et al., 2021), trench behavior (Čížková & Bina, 2013), underplating and recycling of material to the mantle (Bialas et al., 2011; Tewksbury‐Christle et al., 2021), and slab morphology (Čížková & Bina, 2013). The LVZ thus provides a possible window into the location and distribution of the subduction interface and the processes along it, which can be used to characterize modern subduction zones.

Investigations into subduction zone LVZs traditionally involve reflection seismology and/or receiver function waveform inversions. Here we take a complementary approach to investigating the LVZ by constraining the seismic signature of a shear zone that once occupied the subduction interface and is now exhumed. We focus on the Condrey Mountain Schist (CMS) in the Klamath Mountains of northern California/southern Oregon: a prograde, greenschist/epidote‐amphibolite to epidote‐blueschist facies, sediment‐dominated subduction complex exhumed from depths where the LVZ is recognized in modern subduction zones (Bostock, 2013; Helper, 1986; Tewksbury‐Christle et al., 2021). Previous work established the subduction context of the CMS and provided a structural framework for interpreting pulses of deformation and underplating through time (Helper, 1986; Tewksbury‐Christle et al., 2021). We use estimates of shear zone width, occupying rock types, and deformation styles to quantify the CMS' seismic properties during subduction. Our results suggest that the CMS interface shear zone was seismically anomalous due to mineral and fracture anisotropy, with elevated V consistent with modern LVZs.

An Exhumed Subduction Shear Zone in the Klamath Mountains

The CMS is a Late Jurassic to Early Cretaceous subduction complex on the Oregon‐California border that occupies a window through the older, overriding Klamath terranes (Helper, 1986; Snoke & Barnes, 2006) and sits inboard of the younger Franciscan Complex (Dumitru et al., 2010) (Figure 2a). The CMS comprises two main units with limited retrogression ‐ the upper CMS (greenschist to epidote‐amphibolite facies) and the lower CMS (epidote‐blueschist facies). The lower CMS is dominantly epidote‐blueschist facies schist intercalated with m‐ to km‐scale lenses of mafic epidote blueschist and serpentinized ultramafics; it was subducted to 400‐450°C and 0.7–1.1 GPa (∼12°C/km, 30–40 km) (Helper, 1986; Tewksbury‐Christle et al., 2021) (Figure 2b). This geothermal gradient is similar to estimated gradients for warm subduction zones, such as Cascadia and Mexico (Syracuse et al., 2010). The lower CMS schist protolith subducted along a sediment‐poor margin that was tectonically erosive updip of final CMS underplating depths (Tewksbury‐Christle et al., 2021), similar to the shallow erosion and deep underplating occurring along the modern Hikurangi margin (Bassett et al., 2010; Eberhart‐Phillips & Chadwick, 2002).

Figure 2

(a) Regional setting of the Klamath terranes (dark gray), Franciscan (white), and Condrey Mountain Schist (CMS) (blue) (after Snoke & Barnes, 2006). (b) Geologic map and cross section of the CMS (after Helper, 1986; Tewksbury‐Christle et al., 2021). Distributed deformation occurred between ductile thrust faults (red lines). (c) Schematic of the fossil subduction interface shear zone. (d and e) Quartz nodules in schist (d) and mafic blueschist (e) with tails (yellow arrows) elongated along foliation planes (teal lines) interpreted as prograde relict veins.

Neogene doming (Mortimer & Coleman, 1985) exposes 10+ km of lower CMS structural thickness, allowing for detailed characterization of interface shear zone deformation and occupying lithologies. Tewksbury‐Christle et al. (2021) identified three progressively underplated subduction interface shear zones (upper, middle, and lower sheets) in the lower CMS, of different thicknesses and formed at different times (Figure 2b). Here we focus on the middle sheet, for which both the upper and lower shear zone boundaries are preserved, allowing us to constrain shear zone thickness. Tewksbury‐Christle et al. (2021) documented two phases of strain localization within the middle sheet. An early stage of distributed deformation occurred over ∼3 km thickness of dominantly schist (94%) with minor mafic blueschist and serpentinite components (Figure 2c). Following this stage of distributed deformation, introduction of km‐scale serpentinite lenses to the subduction interface allowed for temporary strain localization in serpentinite to <10 m thickness proximal to the thrusts along which the lower CMS was assembled (Figure 2b) (Helper, 1986; Tewksbury‐Christle et al., 2021).

Distributed prograde ductile deformation in the CMS middle sheet resulted in a well‐developed foliation across the heterogeneous lithologies (Figures 2b, 2d and 2e). In the schist, a closely spaced cleavage‐microlithon fabric defined by alternating bands of quartz and graphite + aligned white mica is pervasively developed, consistent with pressure solution creep as the dominant deformation mechanism (e.g., Bell & Cuff, 1989; Durney, 1972; Passchier & Trouw, 2005). In mafic blueschists, sodic amphiboles are elongated within the foliation plane and define a stretching lineation. In addition to the ductile deformation, cm‐scale quartz nodules are common in both the schist and mafic blueschist and have elongated tails parallel to foliation (Figures 2d and 2e). We interpret these nodules as prograde dilational fractures/veins that were cyclically emplaced during progressive deformation, and variably transposed by subsequent ductile deformation, as part of the pressure solution process.

Methods

We estimated the CMS seismic properties for four different endmember assumptions, including: (a) isotropic (Abers & Hacker, 2016), (b) anisotropic (MATLAB Seismic Anisotropy Toolbox, MSAT, Walker & Wookey, 2012), (c) fractured isotropic (randomly oriented fractures, Peacock et al. (2011) and O'Connell and Budiansky (1974); oriented fractures, Hudson (1981) via MSAT; Text S1), and (d) fractured anisotropic (Hudson, 1981; Walker & Wookey, 2012) lithologies. These four scenarios bracket the predicted seismic signature of the CMS shear zone by characterizing the baseline velocities (a), as well as the independent (b and c) and cumulative effects (d) of mineral and fracture anisotropy.

We calculate seismic properties of individual lithologies and the bulk shear zone by averaging wave velocities and/or stiffness tensors weighted by mineral or rock volume fractions measured over thin section‐ to map‐scale (Tables S1 and S2). Results from different averaging schemes (Voigt, Reuss, and Voigt‐Reuss‐Hill) are presented in Section 4. Stiffness tensors are not corrected for pressure‐temperature (P‐T) conditions because of negligible effects (<0.05%, Table S3). For cases with mineral anisotropy, we assume interface‐parallel foliations with crystallographic preferred orientations (CPOs) for aligned minerals (c‐axis perpendicular to foliation: white mica; c‐axis parallel to lineation: glaucophane) based on observations from mineral fabrics in similar exhumed subduction complexes (Cao & Jung, 2016; Keppler et al., 2017; D. Kim et al., 2013; Kotowski & Behr, 2019) (Figure S1). We do not measure CPOs directly using electron backscatter diffraction (EBSD) due to poor indexing of the volumetrically predominant anisotropic mineral (white mica). For cases that include fracture anisotropy, we average porosity, calculated as measured vein area divided by total area, and aspect ratios over thin section, hand sample, and outcrop scales (Table S4). We assume that fracture distributions do not vary significantly in the third dimension and that primary fracture planes were perpendicular to foliation (i.e., σ was parallel to lineation), consistent with Mohr‐Coulomb theory for extensional fracturing (Sibson, 1998) and with similar observations in several other subduction complexes (e.g., Fagereng, 2011; Fisher & Brantley, 2014; Fisher & Byrne, 1987) (Figure S2). We assume lithostatic P (1.0 GPa) and calculate seismic properties using the thermodynamic properties of water (Burnham et al., 1969), although we also explore lower values of P representative of lower pore fluid factors. For cases with both mineral and fracture anisotropy, we merge MSAT's stiffness tensors derived for the mineral anisotropy and oriented fracture anisotropy cases and calculate velocities from the merged tensor (Text S2). For all MSAT velocities, we average the shear wave splitting velocities (V and V ) and calculate V /V and Poisson's ratio (σ) using V s (Figures 3 and 4). Table S5 presents V and V .

Figure 3

Pole figures showing the calculated Condrey Mountain Schist seismic signature accounting for mineral anisotropy (a), fracture anisotropy (b), and both (c). Color maps from Crameri (2018). Isotropic (white bars), maximum (squares), and minimum (triangles) values given for comparison. Assumed foliation (white dashed line) and lineation (white circles) given for orientation.

Figure 4

(a) Comparison of seismic signatures from modern subduction zones (Fukao et al., 1983 1; Matsubara et al., 2009 2; Kato et al., 2014 3; Kato et al., 2010 4; Tsuji et al., 2008 5; Toya et al., 2017 6; Kodaira et al., 2004 7; Audet & Bürgmann, 2014 8; Calkins et al., 2011 9; Delph et al., 2018 10; Peacock et al., 2011 11; Hansen et al., 2012 12; Audet & Schaeffer, 2018 13; Audet et al., 2009 14; Audet & Schwartz, 2013 15; Y. Kim et al., 2010 16; Hicks et al., 2014 17; Y. Kim et al., 2014 18) and calculated from the Condrey Mountain Schist (CMS). Red hexagons indicate CMS seismic properties during prograde deformation. (b) Schematic showing the CMS subduction interface shear zone as part of the Low Velocity Zone and the relationship to transient seismic and aseismic slip source regions.

Results

V /V , assuming isotropic lithologies for the 3‐km‐thick CMS interface shear zone, is low (Figure 4), consistent with experimental measurements of quartz at 1 GPa (Christensen, 1996). Introducing mineral or fracture anisotropy, or a combination, however, results in highly anisotropic V /V with maximums greater than isotropic values (Figure 3). Incidence angles that illuminate maximum and minimum V /V depend on the anisotropy assumptions. Maximum V /V for anisotropic lithologies is in the foliation plane at low angle to the lineation, and minimum V /V is at high angles to the lineation (Figure 3a). In contrast, maximum V /V for fractured isotropic lithologies with 10% porosity, as constrained from our vein measurements, is in a plane normal to the lineation, and minimum V /V is at low angles to the lineation (Figure 3b). Although assumed fracture orientation controls V /V anisotropy, ratios calculated for randomly oriented fractures at 10% porosity are also higher than for isotropic lithologies (Figure 4). The effect sums for fractured anisotropic lithologies, with maximum V /V occurring in the foliation plane but near‐perpendicular to the lineation, and minimum V /V occurring at high angles to the lineation (Figure 3c). Results are not significantly different at σ /P  = 0.8 (Table S6).

If we consider incidence angles perpendicular to the foliation, both anisotropic and fractured anisotropic lithologies have a small local V /V maximum perpendicular to the foliation that is higher than isotropic values (Figures 3a and 3c). Fractured isotropic lithologies are maximum for this incidence angle (Figure 3b).

Discussion

The preservation of strong mineral and fracture anisotropy in the CMS shear zone leads us to interpret our fractured anisotropic lithology results as the best‐constrained prediction of shear zone seismic velocities during prograde deformation. Estimated V /V for the fractured anisotropic case and foliation‐perpendicular arrivals is anomalously high (∼2.0, σ = 0.33, Figure 4a). To capture variations in incidence angle, slab dip, and foliation orientation (e.g., 20° to slab dip in modern subduction zones, Song & Kim, 2012), we average V /V over a ±20° cone centered on foliation‐perpendicular incidence angles (squares, Figure 4a). If foliations are also not perfectly parallel, velocities will approach fractured isotropic values. In the case of these compounding factors, the 20°‐averaged V /V for fractured isotropic lithologies is the best prediction, where V /V is reduced but still elevated (>1.8).

Modern subduction zones have LVZs with slow V and V (up to 70% slower for V , see Figure S3 for comparison of estimated CMS V and modern LVZs), and high V /V (1.8–3.3, 0.45– = 0.28σ) (Figures 1a, 1b and 4a) (e.g., Audet & Bürgmann, 2014; Audet & Kim, 2016; Bostock, 2013; Y. Kim et al., 2014; Song et al., 2009; Toya et al., 2017). Our results incorporating anisotropy and fracture porosity demonstrate that even metasedimentary rocks can reach the lower bounds of the high V /V values in modern subduction zones (Figure 4a). The very high V /V values (e.g., > 2.0, σ > 0.33) cannot be reproduced in our analysis, however. This may imply higher porosity in these regions or overestimated V /V and underestimated thicknesses due to the tradeoff between calculated thickness and V /V in receiver function studies (e.g., Bostock, 2013). Although the CMS is at the low end of modern LVZ thicknesses, thicker LVZs (e.g., Cascadia) may be explained by thicker incoming sediment packages or through contributions from previously underplated material/downgoing slab (Figure 4b) and may be thicker than previously estimated. In contrast, apparent absence of LVZs in margins like Izu‐Bonin could be due to a narrow interface or limited instrumentation.

Because the V /V range for modern LVZs is higher on average than values for isotropic rocks at LVZ depths (<2.0, Christensen, 1996), LVZs are typically attributed to high P (Audet et al., 2009; Bostock, 2013; Hansen et al., 2012; Peacock et al., 2011) based on experimental work correlating high V /V and high P (Christensen, 1996; Eberhart‐Phillips et al., 1989), where high P both maintains significant porosity at high confining pressures and increases V of water (e.g., Eberhart‐Phillips et al., 1989). This is consistent with our observations in the CMS shear zone of abundant quartz veins that were emplaced during brittle fracture associated with pressure solution creep and with our calculated high V /V assuming lithostatic P . Empirical relationships and magnetotelluric studies suggest 0.5%–4% porosity is needed to match LVZ velocities (Calvert et al., 2020; Peacock et al., 2011). Our estimates of up to 10% porosity are higher, but our calculated V /V is still compatible with V /V in modern environments because we also take into account mineral anisotropy. Porosities of up to 10% are compatible with vein exposure measurements in other subduction complexes exhumed from similar conditions, (e.g., 4%–11%, Muñoz‐Montecinos et al., 2020), so these slightly higher values may be more representative than existing experimental constraints.

Observations from the CMS fossil shear zone are also consistent with estimated LVZ thicknesses and some interpretations of the rock types that define the LVZ. Thickness estimates from modern LVZs range from ∼3–8 km (Abers et al., 2009; Audet & Kim, 2016; Audet & Schaeffer, 2018; Audet & Schwartz, 2013; Audet et al., 2009; Bostock, 2013; Delph et al., 2018; Hansen et al., 2012; Hirose et al., 2008; Y. Kim et al., 2014, 2010; Nedimović et al., 2003; Song et al., 2009; Toya et al., 2017) with along‐strike and down‐dip thickness variations (Audet & Schaeffer, 2018; Delph et al., 2018; D. Kim et al., 2019; Toya et al., 2017). The width of the CMS middle sheet was distributed over ∼3 km, comparable with the lower end of these LVZ thicknesses. Furthermore, the shear zone was dominated by metasedimentary protoliths, consistent with interpretations that the LVZ represents deforming and underplating sedimentary packages (e.g., Abers et al., 2009; Calvert et al., 2011; Delph et al., 2021), as opposed to relatively undeformed mafic crust. Sediment prevalence at depth in the CMS interface shear zone, despite subducting along a sediment‐poor, tectonically erosive margin, required stacking of thin incoming sediment packages through protracted underplating and entrainment (Tewksbury‐Christle et al., 2021). Down‐dip thickening observed in modern LVZs (Abers et al., 2009; Hansen et al., 2012; Toya et al., 2017) may be indicative of this progressive underplating process and may be independent of incoming sediment supply, contrary to previous assumptions (Hansen et al., 2012).

The fluid‐filled fracture anisotropy that we include in the CMS best estimate of seismic properties represents the subduction interface while it was actively deforming. However, once subducted material is detached from the downgoing slab and accreted to the upper plate via underplating, mineralization of fractures would change the fracture fill properties. To examine the potential seismic properties of this scenario, we averaged velocities over a ±20° cone with quartz‐filled fractures (Figure 4a). The resulting V (3.18 km/s) and 20°‐averaged V /V are anomalously low. Audet and Bürgmann (2014) previously interpreted low V /V at the base of the forearc crust in Japan and Mexico (Figure 4a) as silica enrichment. Delph et al. (2018, 2021) interpreted low V (<3.2 km/s) at the base of the forearc crust in Cascadia as hydrated underplated metasediments. Our results are consistent with these interpretations, for example, that these zones may represent previously underplated and abandoned metasedimentary material in the upper plate hanging wall or forearc region with mineralized quartz veins.

It is important to note that these predicted velocities are highly dependent on our assumptions. Varying lithologic proportions has limited impact compared to the primary effect of porosity. Decreasing porosity from our maximum value and/or varying mineral alignment will approach isotropic values. Furthermore, velocity behavior with incidence angle is strongly controlled by assumed fracture orientation (Figure 3). Foliation‐parallel veins observed in the Makimine mélange suggest that extreme fluid overpressure can transiently rotate σ 1 by 90° (Ujiie et al., 2018). In the case of our assumptions, this rotates the velocity anisotropy 90° such that V /V is anomalously low perpendicular to the foliation. Although vertical σ 1 is most common for underplated sediments in subduction zones based on rock record analyses (e.g., Fisher & Byrne, 1987), variations in the stress state could affect observed velocity patterns. Validity of these assumptions could therefore be tested by examining LVZ signatures with respect to incidence angle, which could help to deconvolve mineral and fracture anisotropy contributions, lithologic variations, and fracture orientations.

Our interpretation that LVZs in subduction zones may be consistent with a wide, sediment‐dominated shear zone deforming at high P also has implications for the source region and processes involved in slow slip and tremor. Transient seismic (e.g., LFEs) and aseismic slip (e.g., slow slip events) occur collocated with LVZs in modern subduction zones (e.g., Audet & Kim, 2016; Calvert et al., 2020; Delph et al., 2018; Hirose et al., 2008; Song et al., 2009). Temporal and spatial correlation of LFEs, tremor, and slow slip events suggest a genetic connection (Beroza & Ide, 2009; Obara & Hirose, 2006). Competing models for event sources invoke: (a) frictional slip on a heterogeneous fault (e.g., Chestler & Creager, 2017; Ito et al., 2007; Lay et al., 2012; Luo & Ampuero, 2018; Shelly et al., 2006) or (b) frictional failure of blocks or frictionally weak slip planes within a distributed ductile shear zone (e.g., Beall et al., 2019; Behr et al., 2018; Chestler & Creager, 2017; Hayman & Lavier, 2014; Kotowski & Behr, 2019; Tarling et al., 2019; Ujiie et al., 2018). Distinguishing between these two endmember models has important implications for estimating LFE and slow slip source properties, such as slip amount, stress drop, and recurrence (Behr & Bürgmann, 2021; Chestler & Creager, 2017; Frank et al., 2018). The observations that the CMS shear zone (a) accommodated subduction‐related deformation over a 3‐km‐thick zone, (b) records seismic properties that are compatible with modern LVZs, and (c) shows evidence for transient frictional vein emplacement during broader viscous deformation, all lend support to the latter model of distributed frictional‐viscous deformation dominating the deep subduction interface in the slow slip and tremor source region (Figure 4b).

Conclusions

We used estimates of deformation zone thickness, fabric anisotropy, and fracture porosity from a fossil subduction interface shear zone, now exposed at the surface, to calculate its seismic properties for comparison to LVZs in modern subduction zones. This fossilized subduction shear zone exhibits several features in common with modern LVZs, including (a) distributed deformation over a 3 km thick shear zone, compatible with observed LVZ thicknesses, (b) rock types that are consistent with low V and V velocities, and (c) mineral and fracture anisotropy that result in anomalously high V /V for near‐vertical incidence angles. These observations suggest that LVZs in modern subduction zones are compatible with a sediment‐dominated, distributed, subduction interface shear zone deforming under near‐lithostatic fluid pressures, rather than undeformed oceanic crust at near‐lithostatic fluid pressures located below the interface. This interpretation implies that zones of slow slip and tremor, commonly collocated with LVZs, record deformation within distributed frictional‐viscous shear zones rather than along discrete fault planes.

Supporting Information S1

Click here for additional data file.