The fate of the downgoing oceanic plate: Insight from the Northern Cascadia subduction zone

Earth and Planetary Science Letters 408 (2014) 237–251

Contents lists available at ScienceDirect

Earth and Planetary Science Letters www.elsevier.com/locate/epsl

The fate of the downgoing oceanic plate: Insight from the Northern Cascadia subduction zone Nicola Piana Agostinetti a,∗ , Meghan S. Miller b a b

Geophysics Section, School of Cosmic Physics, Dublin Istitute for Advanced Studies, 5 Merrion Square, Dublin 2, Ireland Department of Earth Sciences, University of Southern California, 3651 Trousdale Pkwy, MC 0740, Los Angeles, CA 90089-0740, USA

a r t i c l e

i n f o

Article history: Received 8 April 2014 Received in revised form 6 October 2014 Accepted 9 October 2014 Available online xxxx Editor: P. Shearer Keywords: subduction zones seismic anisotropy Cascadia receiver function

a b s t r a c t In this study, we use teleseismic receiver function analysis to image the seismic structure of the Juan de Fuca oceanic plate during its subduction beneath the North American plate. Seismic data have been recorded at 58 seismic stations deployed along the northern Cascadia subduction zone. Harmonic decomposition of the receiver function data-set along a trench-normal profile allows us to image both the isotropic and the anisotropic structure of the plate (slab). Our images highlight the presence of a highly anisotropic region at 40–70 km depths across the Cascadia subduction zone. The detected seismic anisotropy is interpreted to be related to both metamorphic facies (e.g. blueschists) and fluid released during the dehydration of the subducting mantle. The processes of dehydration and metamorphism produce the variations of the seismic properties within each lithologic unit that constitutes the subducted slab, i.e. basalts, gabbro layer and upper mantle, as the oceanic plate sinks in the upper mantle. Such variations make it almost impossible to recognize the “plate boundary” as a characteristic “velocity-jump” at depth (neither positive nor negative) along the Cascadia subduction zone. Based on the comparative interpretation of both the isotropic and the anisotropic structures retrieved, we propose a 4-stage model of the evolution of the Juan de Fuca oceanic plate during its subduction beneath the North American plate. © 2014 Elsevier B.V. All rights reserved.

1. Introduction The sinking of the oceanic plate at the trench, given by its gravitational instability, is considered the driving force of plate tectonics. Such process has been widely documented both from direct observations of subducted oceanic lithosphere (e.g. seismic tomography, Tsuji et al., 2008) and theoretical studies (e.g. numerical and analogues modeling, Royden and Husson, 2006; Funiciello et al., 2008). During subduction, the oceanic plate follows a pressure–temperature (PT) path that results in its dehydration and the outputs derived from such process (e.g. fluids) play a key-role in geophysical processes occurring along the plate interface (Wassmann and Stöckhert, 2013). For example, water released from the subducted plate can modify the thermal state of the lithosphere, extracting heat from the oceanic crust and affecting the downdip extent of the potential rupture area (i.e. the seismogenetic zone, Cozzens and Spinelli, 2012). A simplified scheme depicts the evolution of the subducted oceanic plate, and related fluid

*

Corresponding author. E-mail addresses: [email protected] (N. Piana Agostinetti), [email protected] (M.S. Miller). http://dx.doi.org/10.1016/j.epsl.2014.10.016 0012-821X/© 2014 Elsevier B.V. All rights reserved.

release, in three main processes: (1) shallow marine sediments compaction; (2) metamorphism of the crust; and (3) breakdown of antigorite minerals and other hydrous phases, such as chlorite and talc, in the oceanic upper mantle. The first process occurs at shallow depths, where the oceanic plate underthrusts beneath the overriding plate, and water stored in pore space is mostly expelled (cf. Faccenda, 2014). The other two processes occur at variable depths, depending on the thermal state of the subducted plate, which, in turn, depends on the plate age and motion. As the oceanic plate subducts, the PT conditions of the subducted materials change resulting in the eclogitization of the oceanic crust and the dehydration of the oceanic upper mantle (Hacker et al., 2003). While these two mechanisms of fluid release are widely accepted as key-processes for the diffuse distribution of the water along the overlying mantle wedge, the depth at which they occur is highly variable between different subduction zones (van Keken et al., 2011). Theoretical models of oceanic plate metamorphism and fluid flow have been presented, based on thermal modeling of the subduction zone, but a number of key-factors needs to be considered, such as weakening of the subduction interface (Wada et al., 2008), localized hydration in the incoming plate (Wada et al., 2012) and alteration to the thermal state due to hydrothermal

238

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

Fig. 1. Generalized map of study area. Plate convergence direction shown by a white vector, plate boundary indicated with a yellow line, and seismicity in the downgoing JdF plate from McCrory et al. (2006, 2012). Label “NAm” indicates the North American plate, “JdF” the Juan de Fuca plate. Red triangles indicate arc volcanoes. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

circulation in the subducting crustal aquifer (Rotman and Spinelli, 2013). Such models are usually tested against the distribution of intermediate-depth earthquakes, which are likely to be generated by both the dehydration of the oceanic crust (Tsuji et al., 2008) and the upper oceanic mantle (Reynard, 2013) and against the analysis of metamorphic rock in exhumed crustal sections (Breeding et al., 2004). In recent years, seismic velocity models of subduction zones have also been used to discriminate between different thermal models and to evaluate the influence of different factors to alter the thermal equilibrium of the oceanic plate (Cozzens and Spinelli, 2012). The high density of seismic stations and the development of new methods for seismic data analysis increased the resolution of the images of the subsurface, allowing to be directly compared to thermal and petrological models for subduction zones. The Cascadia subduction zone (Fig. 1) is an ideal location to explore the evolution and metamorphism of the subducted oceanic plate. Along the west coast of North America the young (<10 Ma) Juan de Fuca plate is subducting at approximately 4 cm/year (Gordon et al., 1990). Images of the Cascadia subduction zone structure at depth are primarily based on detailed seismic observations using various methodologies, even though there is a lack of prevalent Wadati–Benioff zone seismicity. Most recently results based on USArray Transportable and Flexible Array deployments have allowed for characterization of lithospheric structure of the convergent margin. Generalized models of the plate boundary (Fluck et al., 1997; McCrory et al., 2006, 2012) have utilized velocity models and hypocenter data to provide constraints on the slab morphology. Teleseismic tomography studies (e.g. Burdick et al., 2008, 2012; James et al., 2011; Obrebski et al., 2011; Roth et al., 2008; Schmandt and Humphreys, 2010, 2011; Sigloch et al., 2008; Sigloch and Mihalynuk, 2013) image the Juan de Fuca plate with fast P- and S-wave velocity perturbations down to at least the transition zone. The shallow slab and overriding plate structure have been imaged using various teleseismic converted wave imaging techniques (e.g. Receiver Function, RF)

with increasingly higher resolution (e.g. Rondenay et al., 2001; Bostock et al., 2002; Abers et al., 2009; Audet et al., 2009; Hansen et al., 2013). A distinctive dipping low S-wave velocity zone (LVZ), which is interpreted as the subducted oceanic crust was first noted by Langston (1977, 1981) and since then more detailed studies of this LVZ structure have been mapped by other groups. One of the initial detailed studies using RF migration clearly image the subducted basaltic crust of the Juan de Fuca plate beneath central Oregon down to ∼45 km depth, where its signal then disappeared due to the onset of eclogitization (Bostock et al., 2002; Hyndman and Peacock, 2003). Similar results have been recently obtained using seismic data collected in central Washington State (Abers et al., 2009). These images also show a dipping low velocity zone that extends to 45 km depth and is interpreted as the subducted crust of the Juan de Fuca plate. More studies at other subduction zones have focused on the dipping low-velocity zones and have begun to understand that these low-velocity zones are comprised of multiple layers within the subducted oceanic crust. For Cascadia, the low-velocity zone seems to be composed of three layers, such as sediments, hydrated pillow basalts and dykes, and gabbros and other ultra-mafic rocks (Hansen et al., 2013). Despite the large number of studies focused on the Cascadia subduction zone, there are a number of important open questions about the shallow structure in the 40–100 km depth range, as testified by the controversy over the depth of the plate boundary beneath northern Washington, a fundamental parameter that influences the assessment of seismic hazard potential in the area (Audet et al., 2010). We have used broadband seismic data from temporary and permanent network stations in Washington State to image the Swave velocity structure of the Cascadia subduction zone in order to better understand the morphology of the downgoing Juan de Fuca plate to depths of approximately 75 km. We have applied the harmonic decomposition of receiver functions to map seismic anisotropy and infer the structure and evolution of the Juan de Fuca subducted slab beneath Cascadia. Our observations suggest that the seismic properties of the Juan de Fuca plate change

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

239

Fig. 2. Map of seismic stations. Filled triangles indicate seismic stations used in this study. Filling colors show the different groups of stations: snow-white for station deployed along the coast; blue, stations where the RF data set is contaminated by near-surface reverberations from a thick sedimentary cover; dark-blue, stations with a RF dataset consistent with the hypothesis of a sealed plate-boundary (after Audet et al., 2009); light-blue, stations where the JdF upper interface is no more sealed; yellow, stations with clear signature of a continental crust; red, station with a RF dataset contaminated by reverberation and scattering from volcanic deposits; grey, low-quality dataset (for a detailed description of all groups of stations, see the Supporting On-line Materials, Section S2). Triangle size is proportional to RF dataset quality from smallest indicating the worst quality dataset, to largest indicating the best quality dataset. Triangles with a black border indicate stations used to draw Figs. 3 and 4. Crosses display piercing points at 40 km depth for teleseismic rays used in Fig. 4. Orange (light green) crosses indicate piercing points for station S040 (S050). A grey line shows the trace of the profile used in Figs. 3–8 (numbers display distance X along the profile for reference). Isobaths of oceanic Moho depth are also traced: dashed line for the model in Audet et al. (2010); unbroken line for the model in McCrory et al. (2006). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

significantly during the subduction process. The joint interpretation of the isotropic and anisotropic structures gives new insight on topical questions about the Cascadia subduction zone, like the nature of the dipping low-velocity zone and the geometry of the plate boundary, and allows us to propose a hypothesis for the fate of the Juan de Fuca slab that considers the dehydration of the subducted materials. 2. Data and methodology We analyzed more than 20 000 RF from approximately 1100 earthquakes recorded by a network of 58 broadband seismic stations (Fig. 2). Earthquakes between 2006 and 2009 with magnitudes greater than 5.5 and epicentral distances between 30 and 100◦ from the study area were selected. The 58 seismic stations used were a combination of Pacific Northwest Seismic Network (PNSN) permanent stations, USArray Transportable Array (TA) stations, and Cascadia Arrays for Earthscope (CAFE) Flexible Array stations (Fig. 2). The earthquakes recorded with these stations provided good azimuthal coverage, which is essential for the methods applied in this study. First, we discarded teleseismic records that displayed a low S/N ratio within the time window of the theoretical P-arrival. Seismograms were rotated into the ZRT reference system, where R is the Radial direction along the great circle path from the epicenter to the station, T is the Transverse direction perpendicular to R in the horizontal plane, and Z is the vertical component. Then, the RF were calculated using a frequency-domain deconvolution of

the vertical component from the horizontal (both R and T) components, following the approach in Di Bona (1998). Such method allows for estimation of the variance in a single RF, which is then used to weigh the RF in the analysis. We used a low-pass filter with a cutoff frequency of 1 Hz, limiting our vertical resolution of the RF to approximately 3–5 km. Finally, for each station, we visually inspected the RF dataset and discarded low-quality RF which displayed high variance and/or noisy signals (e.g. ringing). The final selection resulted in approximately 6700 high S/N ratio RF, with a mean of ∼100 RF for each station. To evaluate the quality of the single station RF dataset, we plotted both the R-RF and the T-RF dataset in 50% overlapping bins, 20◦ wide back-azimuth bins, as a function of the back-azimuth of the incoming wave. Each dataset was visually inspected and attributed to a different quality class. Visual inspection of the single station RF dataset allows for detection of coherent patterns in the R-RF and T-RF datasets between neighboring stations. Details on the classification of the RF dataset are given in the Supporting On-Line Materials. Both the isotropic and anisotropic structures are investigated in this study, where the isotropic structure represents the depthvariation of the S-wave velocity in the study area, and the anisotropic structures indicates the presence of anisotropic materials at depth. Isotropic structure affects the R-RF component and anisotropic and dipping structures affect both the R-RF and T-RF components. Thus, in this study, we apply two different analysis methods in order to exploit the information contained in both the R-RF and the T-RF data-set.

240

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

First, we compute the depth of the strongest positive S-velocity contrast beneath each single station, following the Zhu and Kanamori (2000) stacking method (ZK2000, hereinafter). Such approach gives first-order information on the strongest positive (i.e. velocity increase with depth) seismic velocity discontinuity at depth, and is routinely used to map the crust–mantle boundary beneath single seismic stations. In this study, we adopt the original stacking scheme using a mean crustal P-wave velocity value of 6.5 km/s, and weights for the single phases Ps, PpPs and PsPs + PpSs are 0.7, 0.2 and 0.1, respectively. Further details of this step of the analysis are given in the Supporting On-line Materials for two seismic stations, S040 and S050 (see Figs. S1(g)–(h) and S2(g)–(h)). This first method is based on a 1D assumption, and, thus, results can be biased by the presence of dipping structures beneath the seismic stations. However, in such cases, estimated errors on the depth of the velocity contrasts are quite large. Following the methodology in Piana Agostinetti et al. (2011), we mapped the anisotropy at depth in the Cascadia subduction zone, analyzing angular harmonics of the RF data-set. This method analyzes both R and T components simultaneously and extracts the back-azimuthal harmonics of the RF dataset as a function of the incoming P-wavefield direction, and allows for separation of the isotropic from the anisotropic component of the RF dataset (Bianchi et al., 2010; Piana Agostinetti et al., 2011). This analysis does not make use of any 1D assumption and it is suitable for tectonic settings where dipping structures are likely to be found (e.g. subduction zones). Here, we consider a profile that has the same strike as used in Abers et al. (2009) in order to directly compare our images to the results presented in their study (Fig. 2). Moreover, we select the same seismic stations used in Abers et al. (2009), which are the closest stations to the profile, to avoid the introduction of 3D artifacts given by the broad distribution of the CAFE stations along the subduction zone. The profile length is sampled with 50% overlapping, 80 × 20 km-wide bins, i.e. one location every 10 km. For each location (called “grid point”, hereinafter) the RF that have a piecing point at 40 km depth that falls within the correspondent bin area are analyzed together. First (k = 0) and second (k = 1) harmonics are extracted from the selected RF dataset. k = 0 harmonics contains information about the depth and the magnitude of the S-wave velocity contrasts beneath a seismic station, while the k = 1 harmonics can be used to constrain the presence of anisotropic materials with a dipping symmetry axis (Bianchi et al., 2010). k = 1 harmonics are projected along two normal axes, N–S and E–W and the sequence of traces along the profile are used to obtain a picture of the subsurface isotropic and anisotropic structure. To facilitate the interpretation, the traces are migrated to depth using the IASPEI91 model (Kennett and Engdahl, 1991). Uncertainties on the k = 0 and k = 1 harmonics are retrieved using a bootstrap approach, i.e. repeating one hundred times the computation of k = 0 and k = 1 harmonics using randomly resampled RF dataset. Finally, a Neighbourhood Algorithm is used to search the parameter space for acceptable seismic velocity models (both isotropic and anisotropic) along the given profile, following Piana Agostinetti et al. (2011). First, for each grid point we inverted the ensemble of selected RF to obtain a best-fit model. Then, the grid points are grouped in regions and the results from the first stage are used to put constrains on the isotropic velocity of the different regions, and the ensemble of RF is inverted for the anisotropic structure alone. For each inversion step and each grid points, a sample of 4 × 104 seismic models are tested against the observed harmonics. Errors on the inverted parameters are estimated from the best-fitting family (i.e. the ensemble of models which fit the observed wiggles with a misfit which is lower than 1.20 times the misfit of the best-fit model). Using this criterion, we obtain a conservative value for the depth-uncertainty of about ±2 km,

which is increased to ±4 km for the deeper interfaces. Errors on the thickness of the layers can be slightly overestimated due to frequency (1 Hz) used in this study (Gesret et al., 2010). Also, as we used only the Ps conversion, our layer thickness can be slightly overestimated with respect to other RF studies in the same area which make use of both Ps and multiples phases (e.g. Hansen et al., 2013). Uncertainties on S-velocity and anisotropic percentage parameters are estimated as large as ±0.2 km/s and 3%, respectively. 3. Results Using the methods depicted above, we can construct a model of the downgoing Juan de Fuca (JdF) slab below the North America margin along a trench-normal profile. This allows us to highlight the evolution (e.g. dehydration and metamorphisms) of the JdF crust and upper mantle, through interpreting the presence and character of seismic anisotropy at depth. To help the reader, we summarize all the observations presented below in Table 1. Hereinafter, the symbol X indicates the distance along the trenchnormal profile, increasing from west to east (cf. Figs. 3–8). 3.1. Analysis of the main S-velocity discontinuity at depth In Fig. 3, we present the results of the application of the ZK2000’s method to estimate the depth of the main S-velocity contrast, beneath each seismic station. Here, we show the results for stations along the transect (Fig. 2). Table S1 displays the results for all the stations in the study area. The application of the ZK2000 method displays coherent results for neighboring stations, indicating the high quality of the data used. Errors on the estimated depth of the S-velocity contrast are larger in the central portion of the profile, where more complex structures are expected (e.g. dipping interfaces). Overall, we observe a dipping structure seen as a S-velocity contrast beneath the western part of the profile and a more horizontal structure to the east. However, several other interesting features are observed in Fig. 3. First, the dip of the S-velocity contrast sharply changes from X = 40 to X = 50 km along the profile (from station S040 to S050, label A in Fig. 3), together with a downward shift of about 5 km. Between X = 50–80 km, the dip is steeper than the western portion of the profile. Yet, in the west, the S-velocity contrast seems to coincide to the top of the gabbroic part of the JdF oceanic crust, as reported in Abers et al. (2009). While east of X = 50 km, the S-velocity contrast follows the plate boundary defined by McCrory et al. (2006). The seismicity in the upper JdF mantle is located between 5 and 10 km below the observed S-velocity contrast, in the region west of S040 (label B in Fig. 3), while it is distributed above and below the velocity contrast east of S040. The “peak” of the occurrence of episodic tremor and slip (ETS) is confined to the region where the S-velocity contrast has shallow dip (label C in Fig. 3). Finally, we observe almost no seismicity above the S-velocity contrast between X = 0 and X = 40 km (label D in Fig. 3). 3.2. Harmonic decomposition Overall, harmonic decomposition of the RF dataset allows us to identify and map the prominent seismic velocity discontinuities at depth and to precisely locate sub-surface volumes where materials display a strong anisotropic behavior. In Fig. 4, we show the harmonic decomposition of the RF dataset along the profile. Panels on the left display the harmonic decomposition of the RFs, while panels on the right display our interpretation. The blue phases in the k = 0 harmonics are generated from positive velocity jump, i.e. where the S-velocity increases with

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

Table 1 Key-observations from the analysis of the results. The X coordinate refers to the distance along the trench-normal profile in Fig. 2, as specified in the text. Label

Brief description

Referred figure

A

Dip of the strongest, positive S-velocity discontinuity changes between station S040 and S050. Intraslab seismicity at 30–70 km depth (i.e. intermediate-depth seismicity, ID) is located 5–10 km beneath the estimates obtained using the ZK method westward from S040. Eastward of such station, it is located both below and above the estimates obtained using the ZK method. Peak density of ETS along the profile occurs between X = 0–40 km along the profile. Between X = 0–40 km, no crustal seismicity above the estimates obtained by the ZK method. Positive pulses dipping from 35 km depth, at X = 0 km, to 80 km depth at about X = 130 km (E-dipping). Between X = 0–90 km, negative pulses at about 30–40 km depth (E-dipping). Dip angle of the positive pulse in “E” changes three times toward East, at about X = 45 km, X = 75 km and X = 130 km. Amplitude of the pulse in “F” is constant (and very coherent between adjacent wiggles) between X = 0–40 km, and changes eastern of X = 45 km (broader toward E). Out-of-Plane energy (OoP) is confined between 25–50 km depth at X = 0–30 km, but it is deeper between X = 60–120 km (depth range 40–65 km). OoP disappears east of X = 120 km (depth range 25–85 km). OoP appears in the continental crust between X = 90–130 km. OoP has two main peaks between X = 0–50 km, but it has three peaks between X = 75–110 km. Two parallel, E-dipping pulses between X = 0–40 km, at about 30 km depth, separated by 5–6 km. Three pulses between X = 80–110 km. The deepest one is E-dipping, the other two are parallel and W-dipping. Depth range: 45–60 km. ETS and “ID seismicity below the pulse in E”, before X = 40 km. After, no ETS and “ID seismicity above the pulse in E”. OoP is confined between pulses in “E” and “F”, between X = 0–70 km. The vertical thickness of the anisotropic layer is uniform (6–8 km) between X = 0–30 km, and then increases to 16–18 km between X = 40–70 km. Between X = 80–110 km, pulses in “E” and the deepest pulses in “N” coincide. Between X = 90–110 km, middle pulse in “N” and red pulses in “k = 0” coincides. Between X = 90–130 km, the upper pulses in “N” and blue pulses in “k = 0” seem to be aligned, depicting an unique interface.

3

B

C D E F G

H

I

J K L M N

P

Q R

S T U

3

3 3 4ab 4ab 4ab

4ab

4cd

4cd 4cd 4cd 4ef 4ef

5

5 5

5 5 5

depth, while energy in the k = 1 harmonics is due to the conversion of a P-wave to an SV-wave and an SH-wave as it crosses into an anisotropic body. In panels (a) and (b) of Fig. 4, we observe the position of the main positive and the comparable negative velocity jumps. First, at the western termination of the profile ( X = 0 km) we identify a positive velocity jump at about 35 km depth (label E), just below a negative velocity jump at 30 km depth (label F). Both are dipping toward the east. The dip of the positive phases changes in three distinct locations along the profile, at X ∼ 45 km, X ∼ 75 km and X ∼ 130 km (label G). Finally, the amplitude of the negative pulse at 30 km depth is constant between X = 0–30 km along the profile, and it changes after X = 45 km becoming broader (label H). From panels (c) and (d), we infer the presence of highly anisotropic zones along the profile. Toward the western end of the profile, we observe some energy in the 25–50 km depth range, while it appears deeper (40–65 km depth)

241

between X = 60–120 km (label I). East of X = 120 km, the energy in the k = 1 harmonics completely disappears, while some energy is recorded in the shallow crust between X = 90 and X = 130 km (labels J and K). It is worth noting that the energy in the k = 1 harmonics seems to depict two parallel interfaces on the western end of the profile, while, we find three interfaces between X = 75–110 km (label L). The number of anisotropic interfaces can be best recognized in panels (e) and (f). There, we clearly identify two parallel, east-dipping anisotropic pulses between X = 0–40 km, at ∼30 km depth (labeled M). In the center of the profile, between X = 80–100 km, there are three main pulses, the deepest E-dipping and the other two, parallel and slightly W-dipping. The three pulses appear between 40 and 65 km depth (label N). To better compare the results with the observations in Fig. 4 (panels b–d–f), the seismicity and the location of the ETS, we plot the colored circles on a single profile (Fig. 5). The relative position of the main pulses in the harmonics and seismicity display some interesting features. First, the peak of ETS and a large number of intermediate-depth earthquakes (focal depths greater than 40 km) occur between X = 0–40 km along the profile. Here, the intraslab events occur below the blue circles (label P). The purple circles, which reflect the presence of anisotropic materials, are confined in between the red and blue circles, between X = 0 and X = 70 km (label Q). Such anisotropic materials appear to have a vertical thickness of about 6–8 km in the area X = 0–30 km, while its thickness increases to 16–18 km between X = 40–70 km (see the two labels R). In the center of our profile, intermediate-depth seismicity is no longer confined below the blue circles, but it appears in between the red and the blue pulses (label P). Finally, the blue circles and the lower branch of dark purple circles coincide between X = 90–110 km (label S), while the upper branch of dark purple circles and the light purple circles seem to follow the alignment of red and blue circles, respectively, west of between X = 90 km (labels T and U). 3.3. Inversion results From the analysis presented above, we draw a clear picture of the geometry of the structures and find evidence of the presence of anisotropic materials at depth. To quantify the seismic properties of such materials and give further details for the discussion, we perform a global inversion of the RF dataset, as described in Section 2. In Fig. 6, we plot the S-velocity profiles retrieved and highlight the presence of the anisotropic layers. In Fig. 7, we show the fit between the observed harmonic components of the RF dataset and the synthetic ones computed using the models in Fig. 6. Globally, the synthetic models are able to reproduce the main patterns in the harmonics which we highlight in Fig. 4. The model presented in Fig. 6 fails to reproduce the features in the 0–20 km depth range (e.g. a sharp blue pulse on the k = 1 at about 20 km depth between X = 0–40 km). However, such depth range is not considered in the analysis nor is essential to our interpretation, so we exclude it from the discussion. The S-wave velocity profiles depict two main features: (1) an east-dipping low-velocity layer with variable thickness, and (2) an anisotropic region which also dips toward the east, but less steeply with respect to the low-velocity layer. These two observations indicate that the bulk properties of the low-velocity layer changes with depth as the where materials undergo metamorphism, however, the inferred anisotropy cannot be solely related to metamorphism. If the anisotropy is only given by mineral orientation in metamorphic rocks, the isotropic and anisotropic structure should display the same geometry. In fact, in such a case, the P-to-S V and P-to-S H phases are generated from the same interface and the colored circles in Figs. 4b and 4f should overlap, which is not the

242

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

Fig. 3. Depth of the main S-velocity contrasts along the profile (profile trace is shown in Fig. 2). Stars display the depth of the S-velocity contrast obtained using the ZK2000’s method, with 1σ error bars. Star colors refer to different group of stations, as in Fig. 1, while star dimension refers to RF data-set quality (see SOM S3 for details). Lines represent plate boundary models: blue unbroken line, Audet et al. (2010) model; Dashed blue line for McCrory et al. (2006) model. Dashed orange lines show the gabbroic part of the oceanic crust in Abers et al. (2009), while the continuous orange line indicates the top of the metasediments/basalts in Abers et al. (2009) (see Fig. DR4 in Abers et al. (2009) for details). Crosses shows earthquakes occurred within the profile (from the catalogue in McCrory et al., 2006): light-(dark-)grey crosses are events beneath (above) the plate boundary, defined by McCrory et al. (2006). Stations along the profile are reported as grey triangles on top. Red triangles represent the stations described in the SOM, S040 and S050. A blue horizontal bar indicates the area along the profile where we have the maximum density (>100 per 5 km) of episodic tremor and slip (ETS) events, from Abers et al. (2009). A sketch of a volcano display the approximate position of the arc. Capital letters refer to observations in Table 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

case. In particular, we observe that, in the first three grid points ( X = 10–30 km), the low V S layer is composed of two parts where the upper layer has a lower V S with respect to the lower layer (V S = 3.2 km/s and V S = 3.8 km/s, respectively), and it is highly anisotropic (anisotropic percentage as high as 9%). The V S layer is coherently dipping to the east. Below this feature, we find a relatively high V S layer (V S = 4.2 km/s) with strong anisotropy (about 7%) that decreases in thickness toward the east. Shifting toward the center of our profile (between X = 40–60 km), the low V S layer is still composed of two parts, but now the total thickness increases and the inferred anisotropy is confined within the lower part of the layer. Below this feature, the high V S layer is still anisotropic, but its thickness is significantly reduced (from about 15 km at X = 0 km to about 8 km at X = 50 km), and the inferred anisotropy in the subducted upper mantle layer disappears east of X = 50 km. It is worth noting that where the high S-velocity layer below 60 km depth becomes isotropic, the low V S above reaches its maximum inferred anisotropy (about 14%). Between X = 70–90 km, the thickness of the low V S layer changes, and anisotropy is now apparent in both the upper and lower portions of this layer, although the upper layer appears to be more strongly anisotropic. In this part of the profile, the V S is about 3.5 km/s. In the central-eastern part of the profile (between X = 100–110 km), the thickness of the low V S layer is reduced (from more than 20 km at X = 60 km to less than 15 km at X = 100 km), and its V S increases to 3.7 km/s, while anisotropy is now apparent in the layer above it, in a zone of slightly higher V S . Finally, anisotropy disappears from the investigated depth range and a sub-horizontal, inverted V S jump is found at about 60 km. Below this feature, we observe a relatively low V S zone (V S is ∼3.7 km/s) of increas-

ing thickness toward east and a high V S lid that dips steeply to more than 90 km depth. Our results also give information about the orientation of the symmetry axis in the anisotropic materials, which, in turn, is related to the foliation in metamorphic facies or to the orientation of fluid-filled crack and fractures. However, the discussion of such parameters is beyond of the scope of the present study. 4. Discussion Our results describe the complexity of JdF plate subduction and suggest that the oceanic crust undergoes to a number of different processes during its subduction. This observation has two main consequences: (1) it is difficult to follow the evolution of a single structural unit of the plate during subduction, without considering metamorphism and interaction with the other components of the “subduction factory” (e.g. fluids migration between different parts), and (2) it makes mapping the plate boundary based on the recognition of a single feature along a subduction zone (e.g. the top of the low V S layer) difficult, due to the changes in the materials properties. We propose a qualitative hypothesis about the phases of oceanic crust subduction in the 30–80 km depth range, where the oceanic crust is expected to slowly change to eclogite. This hypothesis considers the complete process of metamorphism of the main structural units within the plate, including the layered oceanic crust and upper mantle, and the full path of the fluids from the subducted plate to the overriding plate. To strengthen our interpretations, we refer to the predictions of metamorphic facies for Cascadia computed from a high-resolution thermal model in van Keken et al. (2011). In such a warm subduction zone model, the oceanic crust passes the greenschist–

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

243

Fig. 4. (Left panels) Harmonic decomposition of the RF dataset along the profile (profile trace is shown in Fig. 2). (Right panels) Interpreted phases. (a–b) Isotropic (k = 0) component of the RF dataset. Positive (negative) pulses are show in blue (red). In (b), blue (red) dots indicate the strongest positive (negative) phases, used to constrain the structure at depth. (c–d) Sum of the squares of the cosine (N–S) and sine (E–W) contributions to the first harmonics (k = 1). Such quantity represents the “out-of-plane” energy. It should be zero for an isotropic medium where we do not expect any P-to-SH conversion, and can provide a precise location of the anisotropic material at depth. In (d) yellow dotted lines show the most energetic pulses, which testify the anisotropic behavior of the materials at depth. (e–f) Sine (E–W) contribution to the first harmonics (k = 1). In (f) light- and dark-purple dots indicated the strongest “out-of-plane” arrivals along the E–W direction. Capital letters refer to observations in Table 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

amphibolite–blueschists facies before entering the hydrous eclogite facies. The occurrence of blueschists is predicted in the ∼30 to 60 km depth range (cf. van Keken et al., 2011, Fig. 3b). 4.1. Principle-interpretations 4.1.1. Nature of the LVZ One of the most striking features we observe in our results is the presence of an east-dipping low-velocity zone (LVZ) along the profile. This LVZ is found at about 30 km depth at the westernmost end of the profile and disappears at about 70 km depth at X = 110 km. Its thickness is variable from 5-km-thick in the shallower part to 20–25 km in the center of the profile. Dipping LVZs have been recognized in many subduction zones (e.g. Ferris et al., 2003; Tsuji et al., 2008) and they are commonly interpreted as a keymarker of the sinking of the crust (or part of it), but a number

of different hypotheses can be considered (Rondenay et al., 2008). In Northern Cascadia, a LVZ in the upper mantle has been previously identified (Nicholson et al., 2005; Nikulin et al., 2009; Calvert et al., 2011), but its origin is still debated. While our results roughly agree with these previous studies about the location of such LVZ, we introduce a number of observations that can be used to discriminate between different hypotheses that describe its origin. Our results indicate that the LVZ displays variable thickness, from 6-km-thick in the shallower part, closer to the trench, to about 25-km-thick in the middle of the profile. Thickness variations of the LVZ along the profile mimic the shape-changes observed in Abers et al. (2009) along the same profile. Due to the large differences between the two extreme values of the LVZ thickness, we do not associate this LVZ to a single structural unit, such as the oceanic crust (Rondenay et al., 2008), under-plated

244

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

Fig. 5. Comparison of the most relevant pulses in Fig. 4. Red and blue circles report negative and positive amplitudes in the k = 0 plot, respectively. Light and dark purple circles show the position of the main pulses in the k = 1 (E–W) plot. Other symbols as in Fig. 3. Capital letters refer to observations in Table 1. Label “P” indicates the two regions where ETS and intermediate-depth seismicity display different behaviors. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 6. S-velocity models every 10 km along the profile. Colors indicate S-velocity. Over-plotted textures represent anisotropic zones within the profile. Highly (>10%) and low (<10%) anisotropic bodies are shown using dense and light textures, respectively. The upper continental crust (0–20 km depth) is not included. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

sediments (Calvert et al., 2011) or a slice of serpentinized upper mantle (Nikulin et al., 2009). We hypothesize that the LVZ imaged by our receiver functions could include a mixture of these structural units. In the shallower part (30–35 km depth) of the subduction zone, the LVZ has the lowest V S anywhere along the profile, which could be an indication of fluids trapped within the basaltic layer of the oceanic crust (Audet et al., 2009; Hansen et al., 2013; Bostock, 2013). However, its thickness is about 6 km, which is almost double the mean value for this layer found by Hansen et al. (2013). Although this difference can be attributed to the use of

solely Ps phases in our study, we suggest that fluid-filled, underplated sediments are required to match the observed thickness of the LVZ. In the center of the profile, where the thickness of the LVZ reaches its maximum, underplated sediments cannot explain the total thickness of the observed LVZ. Here, a portion of the LVZ can represent the lowermost section of the continental crust hydrated by the fluids coming from the “sealed” layer (see below and Audet et al., 2009). At the depth of 40–60 km, serpentinization of the continental upper mantle has been proposed to explain a 10-km-thick, low S-velocity layer on top of the oceanic crust

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

245

Fig. 7. Observed vs. synthetics wiggles. Synthetic wiggles are computed using the models in Fig. 6. On the left, panels (a)–(c)–(e), the observed k = 0, k = 1 (N–S) and k = 1 (E–W), respectively. On the right, panels (b)–(d)–(f), the synthetic harmonic coefficients, k = 0, k = 1 (N–S) and k = 1 (E–W), respectively.

(Park et al., 2004). The total thickness of the LVZ can represent both the oceanic crust and a slice of serpentinized upper mantle that undergoes retrograde motion along the slab surface due to its buoyancy (Nikulin et al., 2009). In cold subduction zones, such dipping LVZ have been found to extend as deep as 100–150 km, due to the combination of slab metamorphism and localized mantle serpentinization (e.g. Kawakatsu and Watada, 2007). However, in northern Cascadia, where the slab is young and warm, it cannot be identified at depths greater than 50–70 km (Audet et al., 2010). Our results confirm that the LVZ terminates at a relatively shallow depth and this fact implies that both JdF crust dehydration and metamorphism occur at shallower depths (above 70 km). Also, continental upper mantle does not undergo pervasive and localized serpentinization deeper than this depth (cf. the northeastern Japan subduction zone, Tsuji et al., 2008). The oceanic crust traverses the blueschist to hydrous eclogite facies at shallow depths (∼50 km) in the Cascadia subduction zone (cf. van Keken et al., 2011, Fig. 3b). Such depths for the eclogization requires relatively high perme-

ability, and, thus, fluid circulation, within the oceanic crust as predicted by thermal models that include the effects of fluid circulation in an oceanic crust aquifer (cf. Cozzens and Spinelli, 2012, Fig. 3). 4.1.2. The geometry of the main S-velocity jump at depth The Moho discontinuity is commonly associated with a sharp, positive seismic velocity jump at depth. A large number of studies identify dipping seismic velocity jump across subduction zone as the Moho of the subducted plate (e.g. Piana Agostinetti et al., 2011; Kim et al., 2010), while the Moho of the overriding plate can be quite elusive (Bostock, 2013) and, in some case, appear as reversed polarity reflection (Bostock et al., 2002). In our study, we find a strong positive S-velocity jump at depth, which marks the lower interface of the LVZ described above. The dip of this velocity jump has three distinct changes along the profile (at X = 30 km, X = 80 km and X = 110 km in Figs. 4 and 6). While the first change in dip angle at X = 30 km is less pronounced, so that

246

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

a single line can be used to fit, within the errors, the trend of the velocity jump between X = 0–80 km, the same line cannot be used between X = 60–110 km. This is partially caused by the previously mentioned variable thickness of the LVZ. In fact, the reported increase in the thickness of the LVZ occurs expanding the low S-velocity region at both its top and bottom, towards the east. The occurrence of these changes in the dip of the velocity jump suggests that we should not always consider the velocity jump as a marker of a single interface between two different structural units (e.g. the JdF crust and mantle, thus, the JdF Moho). Metamorphism of the oceanic crust and fluid flow from the JdF upper mantle can strongly affect the S-velocity of the subducted materials, so that different structural contacts can be the sharpest S-velocity contrast at different depths. In our case, we suggest that the velocity jump, between X = 0–30 km, represents the intra-crustal, basalt-to-gabbro, interface of the JdF oceanic crust (i.e. between Layer 2 and Layer 3 as defined by Le Pichon (1969)). Discriminating between both oceanic crust layers using RF analysis has been previously reported by Hansen et al. (2013) for the Cascadia subduction zone. Thus, due to the presence of dehydration fluids within the basalts layer, a strong, positive velocity jump is associated to the basalt-to-gabbro interface that results an even larger velocity jump than the one associated with the JdF Moho, especially if the JdF upper mantle is pervasively hydrated (see below). Nowack and Bostock (2013) also noted that a very large velocity jump (from ∼2.0 km/s to 3.4 km/s) between the basalt and the gabbro layers is needed to model the low-frequency earthquakes in the Cascadia area. Between X = 30–80 km, the dip of the velocity jump becomes steeper to the east and it seems to be steeper than the expected dip the JdF Moho (cf. the dip of the plate boundary in Fig. 3, which should be equal to the dip of the JdF Moho). Based on this dip angle for the JdF Moho, the S-velocity jump imaged here should cut the entire oceanic crust. Internal dipping, seismic interfaces can occur within the oceanic crust during its metamorphism (e.g. amphibole–blueschists interface in van Keken et al., 2011), however, the reduction of the velocity with depth is counter-intuitive. As previously mentioned, strong fluid circulation within the oceanic crust has been postulated for this subduction zone. Fluid circulation within the blueschists melange has been observed worldwide (e.g. Breeding et al., 2004, and references therein). We hypothesize that the velocity jump in this area represents the intra-crustal interface between the blueschists (above) and the not-metamorphosed JdF crust (below), where the low S-velocity within the blueschists arises from the presence of channelized fluid coming form the dehydration process (e.g. Ague, 2003). Moving westward, between X = 80–110 km, the almost horizontal trend of the velocity jump again does not match with expected dip of the JdF Moho. However, its dip strikingly resembles the low dip of the blueschist–eclogite interface in the model presented by van Keken et al. (2011) for the Cascadia subduction zone. These authors suggest that eclogization starts from the lower part of the JdF crust (gabbros) and ends with the upper layer that consists of volcanics or sediments. Also, here the velocity jump occurs at ∼70 km depth, which corresponds to the depth suggested by Cozzens and Spinelli (2012). So, we interpret the velocity jump to be the contact between blueschists and eclogite. At X = 110 km, after the complete dehydration and breakdown of the blueschists, the velocity jump would reflect the interface between an hydrated continental upper mantle and the eclogite and it can be defined as the plate boundary. Two main hypotheses can be considered for the hydration of the continental upper mantle. The mantle wedge corner is heavily hydrated and the serpentinized mantle can be dragged downward by the subduction-related flow, releasing water, with the breakdown of the serpentine, in a wide portion of the upper mantle (Kawakatsu and Watada, 2007). Another source of fluids in the upper mantle can be the final dehydration step of

the eclogite which releases fluids from 100 km depth (van Keken et al., 2011). Finally, it is important to notice that in our interpretation the velocity jump marks the JdF Moho along a small portion of our profile (at X ∼ 80 km), where the JdF crust in entirely composed of blueschists (van Keken et al., 2011). 4.1.3. Subduction zone anisotropy: contribution from metamorphism, hydrous phases and fluids flow in the subducted crust Seismic anisotropy is an indicator of many geodynamic processes, from mantle flow to metamorphisms (e.g. Ozacar and Zandt, 2004), as metamorphic materials often exhibit an anisotropic properties (e.g. Kim et al., 2013), even if anisotropy decreases with increasing metamorphic grades. Thus, observations of seismic anisotropy at depth can provide robust constrains of the many different components of the “subduction factory” (e.g. Hacker et al., 2003; Hacker and Abers, 2004). Our results highlight the presence of strong seismic anisotropy within the shallow part of the subduction zone (30–80 km depth), where asthenospheric mantle flow should be less prominent. At such depths anisotropy is also observed along the southern Cascadia subduction zone (Wagner et al., 2013). In our profile, seismic anisotropy can be linked to metamorphic facies, channelized fluid-flow within the subducted crust and/or hydrous phases in both the JdF and the supra-slab mantle. We present results that indicate that bulk seismic anisotropy develops with a dip angle that is less steep than the observed isotropic structures. Such a first-order observation suggests that fluid-flow plays a key-role in the development of anisotropy within the Cascadia subduction zone. Sealed plate boundary: from X = 0 to X = 40 km Along the western part of the profile, we identify two anisotropic regions at 30 km and 40–60 km depth. The upper one is a thin highly-anisotropic layer that coincides with JdF basaltic crust and, possibly, underplated sediments, as discussed above. Here, seismic anisotropy reflects the presence of a large amount of fluids sealed within this layer in cracks and/or fractures (Audet et al., 2009). Such fluids likely come from strong hydrothermal circulation at the oceanic ridge (Faccenda, 2014) and therefore are not of a mantle origin (i.e. fluids rising from the dehydrating JdF mantle) because we do not observe any effect of such fluids on the lower gabbro layer of the JdF crust. Evidence of strong anisotropy in the JdF oceanic crust at shallow depths in the Cascadia subduction zone has been previously observed by Hansen et al. (2013) by measuring the splitting parameters of teleseismic S-waves. The lower anisotropic zone, at 40–60 km depth, together with its relatively low S-velocity indicates that antigorite may be one of the main components of the JdF upper mantle (Mainprice and Ildefonse, 2009). This layer becomes thinner toward east, reflecting the continued breakdown of antigorite minerals and other hydrous phases. We suggest that the observation of no seismic anisotropy, and consequently the absence of mantle-derived fluids, within the gabbro layer may be due to the sealing of the JdF Moho (i.e. low permeability in the gabbro layer) in this portion of the profile, which facilitates the transportation of large amount of water, given by the dehydration of the JdF mantle, at 60 km depth (discussed below). The low-permeability of the gabbro layer has been previously postulated by Bostock (2013) to justify high pore-pressure in the basalts layer of the JdF crust. Permeable plate boundary: from X = 40 to X = 90 km In the central-western portion of the profile, seismic anisotropy is only observed beneath a single region within the 40–65 km depth range. This area partially overlaps with the lower portion of LVZ (where its maximum thickness is 25 km) and the JdF upper mantle, which is completely dehydrated at 60 km depth. In

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

addition, the intensity or amplitude of the anisotropy increases toward the west, reaching its maximum at X = 60 km, where the JdF upper mantle is completely dehydrated. Here, the LVZ has its maximum thickness and its lower part has a velocity decrease with respect to the values found to the east. We suggest that the anisotropy observed at this location is the result of the combination of the metamorphism of the gabbro layer, which gives rise to the bulk anisotropy due to the development of anisotropic minerals (i.e. schistosity), and a strong channelized fluid-flow from the JdF upper mantle to the JdF crust, which is needed to lower the S-velocity of the meta-gabbros. Specifically, the fluids coming from the breakdown of hydrous phases, such as antigorite, chlorite and talc, rising into the JdF crust are channelized in the blueschist melange through the lithologic contact (Breeding et al., 2004), possibly facilitated by (and facilitating) the occurrence of a large number of earthquakes. Channelized fluid-flow has been observed in blueschist melange exposed in subduction zones (Breeding et al., 2004), however, we cannot exclude that the fluids coming from the JdF mantle fill pre-existent fractures in the JdF crust. This process has its maximum intensity at about X = 60 km, where the percentage of the seismic anisotropy reaches its maximum and the largest number of intermediate-depth events occur within the JdF crust. Moving toward the east, seismic anisotropy is found in a region slightly deeper. The highly anisotropic part develops at the top of this region, following the pathways of the fluids from the JdF mantle. The lower part of the anisotropic region is bounded by the positive S-velocity jump described in the above subsection. This fact suggests that this region is composed of blueschists, where the fluids from the JdF mantle almost totally escaped upward toward the mantle wedge corner and the lower continental crust. The most popular tool to investigate subduction zone anisotropy, related to mantle flow, is the measure of the splitting parameters of SKS waves propagating through the upper mantle, which roughly indicate the fast-direction of lattice preferred orientation (LPO) in peridotite (e.g. Long, 2013, and references therein). However, SKS measurements can be affected by the anisotropic behavior of other components of the subduction zone (e.g. blueschists, Kim et al., 2013). In agreement with Wagner et al. (2013), we state that metamorphic rocks and hydrous phases need to be considered along the Northern Cascadia subduction zone in order to better evaluate the SKS measurements and, hence, the mantle flow driven by plate subduction. 4.1.4. Intermediate-depth seismicity Intermediate-depth (i.e. with a focal depth between 30 and 100 km) intraslab seismicity has been observed worldwide in subduction zones and its origin is still debated (Kirby et al., 2000; Brudzinski et al., 2007). Events occurring between 35 and 700 km depth have been attributed to metamorphic reactions occurring within the subducted plate through “dehydration embrittlement” (Barcheck et al., 2012). The presence of Double Benioff Zones (DBZ) has been generally ascribed to the occurrence of two metamorphic reactions: dehydration of metabasalts and meta-gabbros during eclogite formation, for the upper plane, and the breakdown of hydrous phases in the oceanic mantle for the lower plane (Chollet et al., 2011). While the presence of water within the subducted slab is necessary to generate such events, it is still unclear the relationship between the amount of water released by such reactions (i.e. the “dehydration flux”) and the seismicity rate at depth (Barcheck et al., 2012). In Northern Cascadia, Preston et al. (2003) relocated the intermediate-depth seismicity and suggested that intraslab events are relate to metamorphic reactions occurring within the slab. However, the paucity of intermediate-depth seismicity and the lack of defined DBZ prevented an accurate correlation between events and metamorphism. Here, comparing our

247

results to the catalogue from McCrory et al. (2006, 2012), we highlight a number of details about such reactions and help to depict a more complete picture of the subduction processes. Along the western end of our profile, intermediate-depth events are located at about 40–50 km depth, which is 5–10 km below the LVZ described above and within the lower of the two anisotropic regions identified there. Thus, in agreement with Preston et al. (2003), we associate the intermediate-depth events occurring in this portion of the profile to the breakdown of hydrous phases within the JdF upper mantle, and should form the lower plane of the DBZ. In this section of the profile, the paucity of events at relatively shallower depth (30–40 km depth), and the absence of anisotropy at the same depth level, suggests that fluids released from this reaction are confined within the top of the JdF uppermost mantle and do not enter the lower part of the JdF crust (i.e. there is no “hydration embrittlement”). Between X = 40–80 km, we observe that intermediate-depth events occur at slightly deeper depths (45–60 km). At about X = 40 km, the earthquakes are confined within the lower portion of the LVZ, where the anisotropy is observed. Progressively toward the east, the earthquakes extend throughout the whole LVZ, and, finally, at about X = 80 km, they shift to upper part of the LVZ and into the lower continental crust. It is worth noticing that no events occur in the high S-velocity below the LVZ in this part of the profile ( X = 80–100 km). From these observations, we suggest that in the portion of the profile between X = 40–80 km fluids coming form the JdF mantle pervasively hydrate the JdF crust and escape to the lower continental crust. A high-rate of fluid flux can be substained by a channelized fluid-flow through the blueschists melange (Breeding et al., 2004). Fluids trigger blueschist metamorphism, so that intermediate-depth events in this area should be considered part of the upper DBZ, while the disappearance of the lower DBZ plane indicates that JdF upper mantle is completely dehydrated. Complete dehydration of the JdF upper mantle at this shallow depth is possible in presence of a strong, heterogeneous hydration of the oceanic plate (Okazaki et al., 2013). 4.2. The fate of the JdF oceanic crust Our hypothesis is divided in four phases, where each phase reflects a different depth range reached by the JdF oceanic crust during subduction. Fig. 8, is a schematic diagram that represents this hypothesis, where each panel presents one of the four phases used to describe the fate of the JdF crust. The background is the seismic velocity model presented in Fig. 6, together with the anisotropic patterns. Full colors are used for the velocity model in the highlighted area. 4.2.1. Phase 1: JdF crust between 30 to 40 km depth The JdF crust appears to be composed of two layers. The uppermost layer of the JdF slab is composed of basalts in the JdF crust (also called “Layer 2”, e.g. Le Pichon, 1969). Due to the inferred thickness there may be a significant amount of sediment that are underplated beneath the NAM continent. The top boundary of the JdF slab has a very low permeability as postulated by Audet et al. (2010), thus this layer contains over-pressurized fluids, which come from the squeezing of the oceanic “Layer 1” and from lowgrade metamorphism of the basalts. Fluids strongly decrease the S-velocity in this layer, with respect to both continental lower crust and the “Layer 3” of the oceanic crust, and contribute to its anisotropy. The sharpest and strongest positive S-velocity contrast is associated to basalts–gabbro boundary (i.e. between Layer 2 and Layer 3). The JdF mantle is highly hydrated, has a similar velocity to the overlying gabbros, and is highly anisotropic. P-T conditions support the breakdown of hydrous phases and the intermediate-depth seismicity here represents the lower plane of

248

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

Fig. 8. Interpretative sketch. Each panel presents one of the four phases we used to describe the fate of the JdF crust. As background, we plot the seismic velocity model presented in Fig. 6, together with the anisotropic patterns. Full colors are used for the velocity model in the highlighted area. Dark (light) grey crosses represent seismicity within (outside) the highlighted area. Bars on top of each panel indicate different crustal seismicity behavior in the continental plate: orange – no crustal seismicity; yellow – crustal seismicity confined within the lower crust; green – crustal seismicity within the entire crustal volume; light green – crustal seismicity within the upper crust. On the right side of each panel, we reported the observations from Figs. 2–5, as coded in Table 1. A sketch of a volcano display the approximate position of the arc. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

the DBZ. However, at this depth, fluids coming from the JdF upper mantle do not propagate up into the oceanic crust, the JdF crust has not metamorphosed into blueschist (van Keken et al., 2011), and the upper plane of the DBZ is almost absent. This region coincides with the “sealed” region that was identified by Audet et al. (2010) and Bostock (2013) in which trapped fluids within the upper portion of the JdF crust have increased pore-fluid pressure in the subducting plate, facilitating the nucleation of ETS. 4.2.2. Phase 2: JdF crust between 40 and 60 km depth In this phase, the JdF plate uppermost layer increases its permeability and fluids from the JdF crust hydrate the overriding lower continental crust, which lowers its S-velocity and triggers earthquakes. Possibly, the permeability in the JdF plate is reduced due to the volume variation given by the metamorphic reaction occurring within the oceanic crust (Audet et al., 2009). We do not rule out the possibility that the low S-velocity is due to slices of serpentine moving updip along the plate contact from the mantle wedge (e.g. Nikulin et al., 2009). At this depth, basalts and gabbros start to turn into blueschists (van Keken et al., 2011), while the JdF mantle is almost completely dehydrated. Thus, intermediate depth events here represent the upper plane of the DBZ. Variability in the JdF crust volume increases the permeability of the JdF crust (Ahrens and Schubert, 1975). The occurrence of greenschists rather than blueschists cannot be excluded, however, blueschists turn out to be more stable at this depth level (40–60 km depth, Hacker et al., 2003) and display a significant anisotropic behavior (Kim et al., 2013). Dehydration of the upper oceanic mantle at shallow depth requires that the oceanic plate has been subject to heterogeneous hydration at the ridge, also for a warm subduction zone like Cascadia (Okazaki et al., 2013). Anisotropy in the blueschists is increased due to the presence of fluids from the JdF mantle dehydration, and such fluids also lower the S-velocity of the blueschists with respect to gabbro (i.e. the blueschists melange is wet, Breeding et al., 2004). 4.2.3. Phase 3: JdF crust between 60 to 70 km depth At this depth level (60–70 km), the JdF crust completes its final metamorphic phase to eclogite. The crust releases fluids via two mechanisms: (1) fluids coming from the JdF mantle, which have hydrated the JdF crust during Phase 2, which are almost immediately moved to the mantle wedge, creating a higher anisotropic layer above the JdF crust (i.e. the lithological contact act as fluid pathways, Breeding et al., 2004). The lower part of the blueschists melange starts to change to eclogite, hydrating the mantle wedge corner deeper than the fluids coming from the JdF upper mantle. Blueschist–eclogite metamorphic phase occurs at about 70 km depth defining a “dehydration level” which runs incident to the JdF plate interfaces (see van Keken et al., 2011). Fluids coming up to the corner of the mantle wedge can promote serpentinization of the continental upper mantle (Bostock et al., 2002). However, such serpentinization can be more localized to a thin layer standing above the JdF plate (Park et al., 2004). 4.2.4. Phase 4: JdF crust below 70 km depth At this depth level (>70 km), the JdF crust has been completely metamorphosed to hydrous eclogite (van Keken et al., 2011). Fluids rise up into the continental crust and intermediate depth seismicity disappears, together with any evidence of seismic anisotropy at this investigated depth range. Mantle wedge corner displays slightly lower seismic velocity with respect to lower continental crust (Bostock et al., 2002). High S-velocity in the uppermost portion of the continental mantle can represent a very dry portion of the mantle wedge where fluids are absent or not well connected (e.g. Matsuno et al., 2010).

249

5. Conclusion In this study, we present the results of the analysis of an extensive RF dataset recorded along the Northern Cascadia subduction zone. Our results allow us to reinterpret the fate of the metamorphism of the JdF crust and of the fluid migration from the oceanic upper mantle to the North America lower crust. Our findings comprise of 4 main components: 1. A highly anisotropic region at 40–70 km depth has been mapped in the Cascadia subduction zone using the harmonic decomposition of a P-RF dataset. This anisotropic region is shallower near the trench, and becomes deeper and thicker downdip, terminating almost below the volcanic arc. 2. Seismic anisotropy is found to be related to both metamorphic facies (e.g. blueschists) and fluid released during dehydration of the JdF mantle. Both processes, metamorphism and dehydration, are necessary to explain the variation of the S-wave velocity within the JdF crust and its anisotropic behaviour. 3. Seismic properties of the different structural units within the subducted slab, i.e. basalts, gabbro layer and JdF upper mantle, vary with depth, as expected. Blueschists-to-eclogite transition (i.e. eclogitization of the JdF crust) occurs at an almost constant depth of 70 km. 4. The plate boundary cannot be mapped through recognizing a single characteristic “velocity-jump” at depth (neither positive nor negative) beneath the Cascadia subduction zone. Acknowledgements The authors thank Stephen Rondenay and Goeffrey Abers for the seismic data of the CAFE experiment. The authors also thank the Editor Peter Shearer, Ikuko Wada and an anonymous reviewer for their useful comments. N.P.A. research is conducted with the financial support of Science Foundation Ireland & the Marie-Curie Action Cofund under Grant Number 11/SIRG/E2174. M.S.M. is supported in part by NSF-CAREER award EAR-1054638. GMT mapping tool has been used to plot the figures Wessel and Smith (1998). Appendix A. Supplementary material Supplementary material related to this article can be found online at http://dx.doi.org/10.1016/j.epsl.2014.10.016. References Abers, G., MacKenzie, L.S., Rondenay, S., Zhang, Z., Wech, A.G., Creager, K.C., 2009. Imaging the source region of Cascadia tremor and intermediate-depth earthquakes. Geology 37, 1119–1122. http://dx.doi.org/10.1130/G30143A.1. Ague, J., 2003. Treatise on Geochemistry, vol. 3. Elsevier, Amsterdam, pp. 195–228. Chapter ‘Fluid flow in the deep crust’. Ahrens, T.J., Schubert, G., 1975. Gabbro–eclogite reaction rate and its geophysical significance. Rev. Geophys. 13, 383–400. Audet, P., Bostock, M.G., Boyarko, D.C., Brudzinski, M.R., Allen, R.M., 2010. Slab morphology in the Cascadia fore arc and its relation to episodic tremor and slip. J. Geophys. Res. 115, B00A16. Audet, P., Bostock, M.G., Christensen, N.I., Peacock, S.M., 2009. Seismic evidence for overpressured subducted oceanic crust and megathrust fault sealing. Nature 457, 76–78. http://dx.doi.org/10.1038/nature07650. Barcheck, C.G., Wiens, D.A., van Keken, P.E., Hacker, B.R., 2012. The relationship of intermediate- and deep-focus seismicity to the hydration and dehydration of subducting slabs. Earth Planet. Sci. Lett. 349, 153–160. Bianchi, I., Park, J., Piana Agostinetti, N., Levin, V., 2010. Mapping seismic anisotropy using harmonic decomposition of receiver functions: an application to Northern Apennines, Italy. J. Geophys. Res. 115, B12317. http://dx.doi.org/10.1029/ 2009JB007061. Bostock, M., 2013. The Moho in subduction zones. Tectonophysics 609, 547–557.

250

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

Bostock, M.G., Hyndman, R.D., Rondenay, S., Peacock, S.M., 2002. An inverted continental Moho and serpentinization of the forearc mantle. Nature 417, 536–538. Breeding, C.M., Ague, J.J., Broecker, M., 2004. Fluid-metasedimentary rock interactions in subduction zone melange: implications for the chemical composition of arc magmas. Geology 32, 1041–1044. Brudzinski, M.R., Thurber, C.H., Hacker, B.R., Engdahl, R., 2007. Global prevalence of double Benioff zones. Science 316, 1472–1474. http://dx.doi.org/10.1126/ science.1139204. Burdick, S., Li, C., Martynov, V., Cox, T., Eakins, J., Astiz, L., Vernon, F., Pavlis, G.L., van der Hilst, R.D., 2008. Upper mantle heterogeneity beneath North America from travel time tomography with global and USArray transportable array data. Seismol. Res. Lett. 79, 384–392. Burdick, S., van der Hilst, R.D., Vernon, F.L., Martynov, V., Cox, T., Eakins, J., Karasu, G.H., Astiz, J.T.L., Pavlis G, G.L., 2012. Model update March 2011: upper mantle heterogeneity beneath North America from traveltime tomography with global and USArray transportable array data. Seismol. Res. Lett. 83, 23–28. Calvert, A.J., Preston, L.A., Farahbod, A.M., 2011. Sedimentary underplating at the Cascadia mantle-wedge corner revealed by seismic imaging. Nat. Geosci. 4, 545–548. http://dx.doi.org/10.1038/NGEO1195. Chollet, M., Daniel, I., Koga, K.T., Morard, G., van de Moortèle, B., 2011. Kinetics and mechanism of antigorite dehydration: implications for subduction zone seismicity. J. Geophys. Res. 116, B04203. http://dx.doi.org/10.1029/2010JB007739. Cozzens, B.D., Spinelli, G.A., 2012. A wider seismogenic zone at Cascadia due to fluid circulation in subducting oceanic crust. Geology 40, 899–902. Di Bona, M., 1998. Variance estimate in frequency-domain deconvolution for teleseismic receiver function computation. Geophys. J. Int. 134, 634–646. Faccenda, M., 2014. Water in the slab: a trilogy. Tectonophysics. http://dx.doi.org/ 10.1016/j.tecto.2013.12.020. Ferris, A., Abers, G.A., Christensen, D.H., Veenstra, E., 2003. High resolution image of the subducted Pacific (?) plate beneath central Alaska. Earth Planet. Sci. Lett. 214, 575–588. Fluck, P., Hyndman, R.D., Wang, K., 1997. Three-dimensional dislocation model for great earthquakes of the Cascadia subduction zone. J. Geophys. Res. 102 (B9), 20,539–20,550. http://dx.doi.org/10.1029/97JB01642. Funiciello, F., Faccenna, C., Heuret, A., Lallemand, S., Giuseppe, E.D., Becker, T.W., 2008. Trench migration, net rotation and slab-mantle coupling. Earth Planet. Sci. Lett. 271, 233–240. Gesret, A., Laigle, M., Diaz, J., Sachpazi, M., Hirn, A., 2010. The oceanic nature of the African slab subducted under Peloponnesus: thin-layer resolution from multiscale analysis of teleseismic P-to-S converted waves. Geophys. J. Int. 183 (2), 833–849. http://dx.doi.org/10.1111/j.1365-246X.2010.04738.x. Gordon, C.D.R.G., Argus, D.F., Stein, S., 1990. Current plate motions. Geophys. J. Int. 101, 425–478. Hacker, B., Abers, G., 2004. Subduction factory 3: an Excel worksheet and macro for calculating the densities, seismic wave speeds, and HO contents of minerals and rocks at pressure and temperature. Geochem. Geophys. Geosyst. 5. http://dx.doi.org/10.1029/2003GC000614. Hacker, B.R., Abers, G., Peacock, S., 2003. Subduction factory 1. Theoretical mineralogy, densities, seismic wave speeds, and H2 O contents. J. Geophys. Res. 108. http://dx.doi.org/10.1029/2001JB001127. Hansen, R.T.J., Bostock, M.G., Christensen, N.I., 2013. Nature of the low velocity zone in Cascadia from receiver function waveform inversion. Earth Planet. Sci. Lett. 337–338, 25–38. Hyndman, R.D., Peacock, S.M., 2003. Serpentinization of the forearc mantle. Earth Planet. Sci. Lett. 212, 417–432. James, D.E., Fouch, M.J., Carlson, R.W., Roth, J.B., 2011. Slab fragmentation, edge flow and the origin of the Yellowstone hotspot track. Earth Planet. Sci. Lett. 311, 124–135. Kawakatsu, H., Watada, S., 2007. Seismic evidence for deep-water transportation in the mantle. Science 306, 1468–1471. Kennett, B.L.N., Engdahl, E.R., 1991. Travel times for global earthquake location and phase identification. Geophys. J. Int. 105, 429–465. Kim, D., Katayama, I., Michibayashi, K., Tsujimori, T., 2013. Deformation fabrics of natural blueschists and implications for seismic anisotropy in subducting oceanic crust. Phys. Earth Planet. Inter. 222, 8–21. Kim, Y., Clayton, R.W., Jackson, J.M., 2010. Geometry and seismic properties of the subducting Cocos plate in central Mexico. J. Geophys. Res. 115, B06310. http://dx.doi.org/10.1029/2009JB006942. Kirby, S., England, R.E., Denlinger, R., 2000. Subduction: Top to Bottom. Geophysical Monograph, vol. 96. American Geophysical Union, pp. 195–214. Chapter ‘Intermediate-depth intraslab earthquakes and arc volcanism as physical expression of crustal and upper mantle metamorphism in subducting slab’. Langston, C.A., 1977. Corvallis, Oregon, crustal and upper mantle receiver structure from teleseismic P and S waves. Bull. Seismol. Soc. Am. 67, 713–724. Langston, C.A., 1981. Evidence for the subducting lithosphere under southern Vancouver Island and western Oregon from teleseismic P wave conversion. J. Geophys. Res. 86, 3857–3866. Le Pichon, X., 1969. Models and structure of the oceanic crust. Tectonophysics 7, 385–401.

Long, M., 2013. Constraints on subduction geodynamics from seismic anisotropy. Rev. Geophys. 51, 76–112. Mainprice, D., Ildefonse, B., 2009. Subduction Zone Dynamics. Springer-Verlag, Berlin, Heidelberg, pp. 61–85. Chapter ‘Seismic anisotropy of subduction zone: minerals – contribution of hydrous phases’. Matsuno, T., Seama, N., Evans, R.L., Chave, A.D., Baba, K., White, A., nori Goto, T., Heinson, G., Boren, G., Yoneda, A., Utada, H., 2010. Upper mantle electrical resistivity structure beneath the central Mariana subduction system. Geochem. Geophys. Geosyst. 11, Q09003. http://dx.doi.org/10.1029/2010GC003101. McCrory, P.A., Blair, J.L., Oppenheimer, D.H., Walter, S.R., 2006. Depth to the Juan de Fuca slab beneath the Cascadia subduction margin – a 3-D model sorting earthquakes. In: U.S. Geological Survey Data, vol. 91. V1.2. McCrory, P.A., Blair, J.L., Waldhauser, F., Oppenheimer, D.H., 2012. Juan de Fuca slab geometry and its relation to Wadati–Benioff zone seismicity. J. Geophys. Res. 117. http://dx.doi.org/10.1029/2012JB009407. Nicholson, T., Bostock, M., Cassidy, J.F., 2005. New constraints on subduction zone structure in northern Cascadia. Geophys. J. Int. 161, 849–859. Nikulin, A., Levin, V., Park, J., 2009. Receiver function study of the Cascadia megathrust: evidence for localized serpentinization. Geochem. Geophys. Geosyst. 10, Q07004. http://dx.doi.org/10.1029/2009GC002376. Nowack, R.L., Bostock, M.G., 2013. Scattered waves from low-frequency earthquakes and plate boundary structure in northern Cascadia. Geophys. Res. Lett. 40. http://dx.doi.org/10.1002/grl.50826. Obrebski, M., Allen, R.M., Pollitz, F., Hung, S.H., 2011. Lithosphere–asthenosphere interaction beneath the western United States from the joint inversion of body-wave traveltimes and surface-wave phase velocities. Geophys. J. Int. 185, 1003–1021. Okazaki, K., Katayama, I., Noda, H., 2013. Shear-induced permeability anisotropy of simulated serpentinite gouge produced by triaxial deformation experiments. Geophys. Res. Lett. 40, 1290–1294. http://dx.doi.org/10.1002/grl.50302. Ozacar, A.A., Zandt, G., 2004. Crustal seismic anisotropy in central Tibet: implication for deformational style and flow in the crust. Geophys. Res. Lett. 41, L23601. http://dx.doi.org/10.1029/2004GL021096. Park, J.J., Yuan, H., Levin, V., 2004. Subduction zone anisotropy beneath Corvallis, Oregon: a serpentinite skid mark of trench-parallel terrane migration? J. Geophys. Res. 109. http://dx.doi.org/10.1029/2003JB002718. Piana Agostinetti, N., Bianchi, I., Amato, A., Chiarabba, C., 2011. Fluid migration in continental subduction: the Northern Apennines case study. Earth Planet. Sci. Lett. 302. http://dx.doi.org/10.1016/j.epsl.2010.10.039. Preston, L.A., Creager, K.C., Crosson, R.S., Brocher, T.M., Trehu, A.M., 2003. Intraslab earthquakes: dehydration of the Cascadia slab. Science 302, 1197–1200. Reynard, B., 2013. Serpentine in active subduction zones. Lithos 178, 171–185. Rondenay, S., Abers, G.A., van Keken, P.E., 2008. Seismic imaging of subduction zone metamorphism. Geology 36, 275–278. http://dx.doi.org/10.1130/G24112A.1. Rondenay, S., Bostock, M., Shrag, J., 2001. Multiparameter two-dimensional inversion of scattered teleseismic body waves. J. Geophys. Res. 106, 30795–30807. Roth, J.B., Fouch, M.J., James, D.E., Carlson, R.W., 2008. Three-dimensional seismic velocity structure of the northwestern United States. Geophys. Res. Lett. 35. http://dx.doi.org/10.1029/2008GL034669. Rotman, H.M.M., Spinelli, G.A., 2013. Global analysis of the effect of fluid flow on subduction zone temperatures. Geochem. Geophys. Geosyst. 14, 3268–3281. http://dx.doi.org/10.1002/ggge.20205. Royden, L., Husson, L., 2006. Trench motion, slab geometry and viscous stresses in subduction systems. Geophys. J. Int. 167 (2), 881–905. http://dx.doi.org/ 10.1111/j.1365-246X.2006.03079.x. Schmandt, B., Humphreys, E.D., 2010. Complex subduction and small-scale convection revealed by body-wave tomography of the western United States upper mantle. Earth Planet. Sci. Lett. 297, 435–445. Schmandt, B., Humphreys, E.D., 2011. Seismically imaged relict slab from the 55 Ma Siletzia accretion to Northwest USA. Geology 39, 175–178. Sigloch, K., McQuarrie, N., Nolet, G., 2008. Two-stage subduction history under North America inferred from multiple-frequency tomography. Nat. Geosci.. http://dx.doi.org/10.1038/ngeo23. Sigloch, K., Mihalynuk, M.G., 2013. Intra-oceanic subduction shaped the assembly of Cordilleran North America. Nature 496, 50–56. http://dx.doi.org/10.1038/ nature12019. Tsuji, Y., Nakajima, J., Hasegawa, A., 2008. Tomographic evidence for hydrated oceanic crust of the Pacific slab beneath northeastern Japan: implications for water transportation in subduction zones. Geophys. Res. Lett. 35. http:// dx.doi.org/10.1029/2008GL034461. van Keken, P.E., Hacker, B.R., Syracuse, E.M., Abers, G.A., 2011. Subduction factory 4: depth-dependent flux of H2 O from subducting slabs worldwide. J. Geophys. Res. 116, B01401. http://dx.doi.org/10.1029/2010JB007922. Wada, I., Behn, M.D., Shaw, A.M., 2012. Effects of heterogeneous hydration in the incoming plate, slab rehydration, and mantle wedge hydration on slab-derived H2 O flux in subduction zones. Earth Planet. Sci. Lett. 353–354, 60–71.

N. Piana Agostinetti, M.S. Miller / Earth and Planetary Science Letters 408 (2014) 237–251

Wada, I., Wang, K., He, J., Hyndman, R.D., 2008. Weakening of the subduction interface and its effects on surface heat flow, slab dehydration, and mantle wedge serpentinization. J. Geophys. Res. 113, B04402. http://dx.doi.org/10.1029/ 2007JB005190. Wagner, L.S., Fouch, M.J., James, D.E., Long, M.D., 2013. The role of hydrous phases in the formation of trench parallel anisotropy: evidence from Rayleigh waves in Cascadia. Geophys. Res. Lett. 40, 2642–2646. http://dx.doi.org/10.1002/grl.50525.

251

Wassmann, S., Stöckhert, B., 2013. Rheology of the plate interface – dissolution precipitation creep in high pressure metamorphic rocks. Tectonophysics 608, 1–29. http://dx.doi.org/10.1016/j.tecto.2013.09.030. Wessel, P., Smith, W.H.F., 1998. New, improved version of the generic mapping tools released. EOS Trans. AGU 79, 579. Zhu, L., Kanamori, H., 2000. Moho depth variation in southern California from teleseismic receiver function. J. Geophys. Res. 105, 2969–2980.