Fluid pressure and shear zone development over the locked to slow slip region in Cascadia.

At subduction zones, the deep seismogenic transition from a frictionally locked to steady sliding interface is thought to primarily reflect changes in rheology and fluid pressure and is generally located offshore. The development of fluid pressures within a seismic low-velocity layer (LVL) remains poorly constrained due to the scarcity of dense, continuous onshore-offshore broadband seismic arrays. We image the subducting Juan de Fuca oceanic plate in northern Cascadia using onshore-offshore teleseismic data and find that the signature of the LVL does not extend into the locked zone. Thickening of the LVL down dip where viscous creep dominates suggests that it represents the development of an increasingly thick and fluid-rich shear zone, enabled by fluid production in subducting oceanic crust. Further down dip, episodic tremor, and slip events occur in a region inferred to have locally increased fluid pressures, in agreement with electrical resistivity structure and numerical models of fault slip.

INTRODUCTION

Subduction megathrust faults produce damaging earthquakes in seismogenic zones through stick-slip processes (, ). The seismogenic zone itself is defined as the area of the fault that is partially to fully locked during the interseismic period and slips during earthquakes; several other definitions exist but most involve dominantly frictional over viscous processes (, ). Because earthquake magnitude is controlled partly by rupture area, knowledge of the down-dip limit of the seismogenic zone is key to modeling the upper limit of earthquake magnitude and the associated seismic hazard in coastal regions (). For the Cascadia subduction zone, thermal and plate coupling models based on heat flow and geodetic data predict that the seismogenic zone is dominantly located offshore (, ). Episodic tremor and slow slip events () occur further down dip and are separated from the seismogenic zone by a ~70-km-wide apparent aseismic gap (). A growing body of evidence indicates that tremors and low-frequency earthquakes (LFEs) occur near or within a dipping seismic low-velocity layer (LVL), characterized by anomalously high compressional-to-shear velocity ratio (vp/vs) values, which are commonly interpreted as near-lithostatic pore-fluid pressure (Pf) (). The inferred low effective fault-normal stress (σe) near the plate interface is consistent with numerical simulations of episodic slow slip and tremor-triggering properties that require a weak megathrust down dip of the seismogenic zone (). Numerical models of fault rupture require a decrease in σe from a strongly coupled interface (~50 MPa) in the seismogenic zone region to a weakly coupled interface (~3 MPa) down dip, presumably controlled by an increase in Pf (); however, observational evidence of this change remains sparse (–).

We use observations of converted teleseismic waves (that is, receiver functions) recorded over a combination of onshore and offshore broadband seismograph deployments in northern Cascadia (Fig. 1). The offshore component uses selected ocean-bottom seismograph (OBS) stations from the Cascadia Initiative (CI) deployed on the continental forearc shelf with an average station spacing of ~8 km. The land component uses broadband stations from the Cascadia Array for Earthscope (CAFE) experiment with station spacing of ~10 km. Together, these deployments form a dense linear array extending from the middle of the seismogenic zone, ~55-km down dip from the surface projection of the trench, to the end of the slow slip source region. On land, receiver functions are dominated by a set of oppositely polarized pulses [forward P-to-S (Ps), back-scattered P-to-S (Pps), and S-to-S (Pss) conversions] that characterizes the signature of the dipping LVL (Fig. 2), as imaged previously in this region (–). Offshore, the forward and back-scattered conversions appear to merge and interfere with each other, and the signature of a possible shallow LVL becomes difficult to identify in this data set (, ).

Fig. 1

Map of the Cascadia forearc region.

Onshore (CAFE) and offshore (CI) broadband seismic stations shown as inverted blue triangles and yellow squares, respectively. Plate interface contours are shown as dashed lines [from the study of Audet et al. ()] and dash-dotted lines [from the study of McCory et al. ()]; the inferred down-dip limit of the locked zone () is denoted by red (for 70% locked) and yellow (20% locked) lines. Earthquakes from the Northern California Earthquake Data Center catalog () are denoted by gray circles and LFEs () as colored diamonds.

Fig. 2

Offshore (CI) and onshore (CAFE) receiver functions and LVL signature in northern Cascadia.

Onshore data show the signature of an LVL as dipping negative-positive–converted (Ps) and reverberated (Pps and Pss) pulses from (A) the bottom (b) and (B) the top (t) of the LVL shown here for station S050 (C) and for each station across the linear profile (D). Vertical dashed line indicates the coastline. Offshore signals are more difficult to interpret. Gray-shaded area in (D) around station S050 highlights the signal shown in (C).

RESULTS

We apply common conversion point (CCP) stacking of the receiver functions projected along a line orthogonal (purple line; Fig. 1) to the trench to image the down-dip variability in LVL signature. CCP stacks are produced for each of the three main P- to S-converted phases (Ps, Pps, and Pss) (fig. S1). A weighted sum of these images reveals laterally continuous signals associated with the LVL (Fig. 3A). Greater weight is given to back-scattered (Pps and Pss) phases due to their higher coherency and larger time separation for the double pulse, as well as the difficulty in resolving the Ps phase at high frequency with ocean-bottom seismic data (, ). We further process receiver functions using Gaussian beam weighting () for each individual phase to account for P-wave sensitivity and recompute the CCP stacks (fig. S2). The individual Gaussian-weighted stacks are then combined using a phase-weighted sum (), which suppresses incoherent noise, as well as mismapped converted phases that inevitably contaminate the individual stacks. The resulting image (Fig. 3B) shows a semicontinuous LVL characterized by negative-positive impedance contrasts across the top and bottom layer boundaries, respectively, consistent with previous teleseismic images of the Cascadia forearc structure (, , ). The amplitude of the negative pulse associated with the top boundary increases down dip and appears to flatten out at a distance of ~200 km. At depths shallower than 20 km, the same negative signal decreases, whereas the amplitude of the positive pulse associated with the bottom LVL boundary reaches a maximum and the pulse becomes more horizontal (from 100-km trenchward).

Fig. 3

Receiver function images from CCP stacking of converted Ps and back-scattered Pps and Pss phases projected along the linear profile.

(A) Weighted sum and (B) Gaussian- and phase-weighted sum of the Ps, Pps, and Pss CCP phase stacks. Yellow and gray circles indicate low-frequency and regular seismicity, respectively. Receiver function (RF) and U.S. Geological Survey (USGS) slab models are taken from refs. () and (). (C) Same as (B) for synthetic data calculated for a model that incorporates an LVL (vp/vs, ~2) with thickness increasing from 2 km at the coastline to >5 km near the intersection with the mantle wedge corner (fig. S4 and section S1, model 8). Higher vp/vs (~3.5) near the source of LFEs [white-dashed area in (C)] suggests locally increased Pf. (C) also illustrates the approximate boundaries in 70 and 20% plate coupling (), the episodic tremor and slip (ETS) region, and static strength regimes (). These transitions correspond with the changes in fluid pressure (Pf) and shear-zone thickness across the coastline (red inverted triangle). Down dip of the viscous zone, locally increased Pf leads to embrittlement, where ETS occurs.

We model receiver functions using laterally varying seismic velocity models, simulating the geometry of the Cascadia subduction zone (fig. S4 and table S1) (, , ); synthetic receiver functions are further processed using the CCP algorithm. Quantitative comparison is performed through cross-correlation (CC) between the synthetic and observed images and from the mean data misfit calculated over the array (fig. S5). A suite of synthetic CCP images demonstrate that the disappearance of the LVL in the up-dip portion is an artifact of the processing; however, models that do not include a thick (1 to 3 km) offshore LVL with high vp/vs produce the lowest misfit, indicating that the offshore LVL is not required by these data. Our preferred model (Model 8, fig. S4, and section S1), for which the resulting synthetic CCP stack is illustrated in Fig. 3C, is characterized by a dipping LVL located above the plate interface [from the study of McCrory et al. ()] with vp/vs of ~2 and thickness increasing from 0 km (that is, no LVL) in the middle of the seismogenic zone offshore to >5 km at the intersection with the mantle wedge corner (), where vp/vs locally increases from ~2 to ≥3.5. The geometry of the LVL is consistent with the seismic reflection signature of the plate interface from active source seismic data (, ).

DISCUSSION

Similar to previous work, we interpret the low-velocity signature with high vp/vs as overpressured subducting material (), with fluids sourced from low-grade dehydration reactions in porous, altered oceanic crust and capped by a low-permeability seal (). The observed variations of this signature along dip coincide with important transitions in plate coupling and slip modes (Fig. 3C) (, ). Combined with evidence of water flux in the upper oceanic crust immediately up dip from refraction profiles (), the lack of a developed LVL with high vp/vs suggests that Pf are maintained at hydrostatic levels within the fully coupled region. Development of a thick (>2 km) LVL with high Pf coincides with the coastline, where the plate interface coupling is from transitional (15 to 70%) to steady sliding (<15%) (). Further down dip, LFEs within tremor () occur marginally trenchward, although deeper than the region where the LVL has the highest vp/vs values (Fig. 3C). Furthermore, this maximum vp/vs region also coincides with the intersection of the subducting plate with the mantle wedge corner (). Here, blueschist-to-eclogite reactions occurring within the oceanic crust locally enhance fluid production, resulting in upward migration and increased serpentinization of the mantle wedge corner (, ). The large permeability contrast between the juxtaposed gabbroic lower continental crust and underlying serpentinite may locally enhance the accumulation of fluids, further contributing to fluid overpressures and weakening the plate interface, promoting episodic slow slip (, ).

These results provide observational support to two classes of megathrust slip models: thermorheological models of static fault strength () and dynamic frictional models of fault rupture (). In the static strength model, the competing effects of σe (fluid-pressure–dependent) and temperature control the relative importance of frictional or viscous strength, respectively, and produce the transition in slip modes along dip (). Within the seismogenic zone, low plate interface temperatures and well-drained conditions allow plate locking and frictional slip. As temperature and fluid production increase down dip, viscous creep and stable sliding dominate. Episodic tremor and slow slip events occur where fluid pressure is locally increased, enabling the regime to regionally transition back to frictional slip (Fig. 3C) (). In the framework of dynamic friction, the competing rates of shear stress reduction (or weakening, proportional to σe) and elastic unloading (independent of σe) during slip determine fault stability (). At low Pf and high σe (~50 MPa) in the offshore region (), slip weakening occurs more rapidly, which leads to unstable slip and earthquakes. At high Pf and low σe, slip weakening is reduced and rupture is stopped, leading to stable creep. Near-lithostatic Pf (very high vp/vs) and very low σe (~3 MPa) promote conditionally stable and slow fault slip under strong dilatancy ().

Finally, the nature of the LVL remains ambiguous, at least in this part of Cascadia. Further north under southern Vancouver Island, LFEs within tremor occur within the LVL itself (). In that region, the layer is interpreted as overpressured upper oceanic crust whose low-velocity signature terminates at peak metamorphic dehydration conditions, consistent with observations worldwide (). Here, however, despite potential uncertainties in LFE depths, LFEs appear to cluster on the plate interface model of McCrory et al. () below the LVL (Fig. 3, A and B). The increasing thickness of the LVL could alternatively represent the progressive development of the megathrust into a viscous and fluid-rich ductile shear zone (), perhaps composed of underplated sediments () or a thick mélange of heterogeneously metamorphosed basaltic rocks derived from the oceanic crust (). In this framework, LFEs within tremor could represent slip along overpressured mylonitic fabrics () and/or heterogeneous brittle failure within a viscously deforming matrix (, , ). Our results suggest that the aseismic gap in Cascadia can be explained by the development of high Pf within a progressively thickened shear zone at levels that favor viscous creep and stable sliding. These findings are also consistent with lateral variations in coupling and structure along oceanic transform faults (), which indicates similar controls of fault slip modes from fluid pressure and shear zone development in different plate boundary fault types.

MATERIALS AND METHODS

Data selection and preprocessing

Broadband seismic data used in this study come from stations that are part of the CAFE () experiment that took place between 2005 and 2008 and the CI (), where OBS stations were deployed between 2011 and 2015. These stations followed a linear profile approximately perpendicular to the trench (Fig. 1). We only used data for the third year (2013 to 2014) of CI deployment due to significantly higher data quality (). At each selected station, we collected 360-s-long, three-component P-wave seismograms from all available teleseismic events (magnitude, M > 5.8) with epicentral distance between 30° and 90°. Waveforms from the two data sets were preprocessed separately due to nonoverlapping recording periods. We preprocessed seismograms by removing the mean value, detrending, downsampling to 10 Hz, and high-pass filtering using a zero-phase filter with corner frequency of 0.085 Hz. At those frequencies, OBS data were not affected by compliance noise even for stations deployed in shallow water (). For each event, we measured signal-to-noise ratio (SNR) on the vertical component for each station and calculated the median value over all stations within their respective network. We selected events for which all SNR values are higher than 0 dB and for which the network-wide median SNR is higher than 3.5 dB. We then discarded OBS stations for which fewer than five earthquakes were available, resulting in eight OBS stations with a sufficient number of recorded earthquakes for receiver function processing. This resulted in a range of 11 to 13 earthquakes per OBS station, and an average of 65 earthquakes for the stations deployed on land. Horizontal component OBS data were then rotated to the East-North coordinate system using the station orientations given in the study of Janiszewski and Abers ().

Receiver functions

Following data selection and preprocessing, data from both land and OBS stations were processed similarly. We first rotated horizontal components to the radial-transverse coordinate system. We did not further rotate the waveforms into P- and S-wave modes because this procedure is not advisable for OBS data (). The vertical components were then deconvolved from the radial components using Wiener spectral deconvolution to obtain individual receiver functions (). We then produced three sets of bandpass-filtered receiver functions differing only in their high-frequency corner: 0.085 to 0.75 Hz, 0.085 to 0.3 Hz, and 0.085 to 0.2 Hz, which highlight more clearly the forward-converted Ps and back-scattered Pps and Pss receiver function sets, respectively ().

CCP imaging

We used a CCP stacking technique [for example, the study of Tauzin et al. ()] to map receiver function amplitude to depth for each converted phase (Ps, Pps, and Pss) using an S-wave velocity model of the subsurface from a global tomography study (). The receiver function amplitudes were projected and averaged (or stacked) onto a trench-perpendicular profile, resulting in one CCP stack for each converted phase (fig. S1, A to C). We produced a first CCP image by taking a weighted sum of the three CCP stacks along the profile (fig. S1D). We selected a weighting scheme where weights for the back-scattered Pps and Pss phases are three times higher than those for the converted Ps phase because Pps and Pss phases are more coherent at those frequencies than Ps for OBS data. We produced a second set of CCP stacks by applying a Gaussian filter to the individual receiver functions in the horizontal direction using a width of 15 km to simulate a first-order P-wave sensitivity kernel (fig. S2, A to C) (). The three Gaussian-filtered CCP stacks were then combined into a second CCP image using a phase-weighted sum (). The phase weights here refer to the instantaneous phase of the analytic signal for each CCP stack to enhance coherent signals between the three CCP stacks and to suppress the inevitable effect of contamination from mismapped phases in the final image (fig. S2D).

Receiver function modeling

We modeled receiver functions using a reflectivity technique for stacks of uniform isotropic layers separated by sharp, flat interfaces and modified to take into account reverberations from the water column [see previous study () for details]. Although our approach does not consider dipping layers, these effects were expected to be minor because the LVL dip is low (~10°) (). For each station, we generated synthetic radial receiver functions using the same event parameters (back azimuth and slowness of incoming teleseismic plane P wave). Each set of synthetic receiver functions was further processed using the Gaussian- and phase-weighted CCP stacking algorithm. Models were evaluated on the basis of the CC between the synthetic and observed images, as well as the mean station misfit over the array, where the station misfit was calculated from the mean of squared residuals for all events at a given station. We also estimated a relative improvement factor from the difference between the percent increase in CC (positive for improvement) and percent increase in mean misfit (negative for improvement), calculated for two different seismic velocity models. The suite of models considered is described in the Supplementary Materials (section S1).

We noted that OBS receiver functions only display coherent arrivals between 0 and approximately 12 s after which time there is significant scatter in the data (fig. S3, A and B). Synthetic receiver functions calculated for realistic subduction zone models (see below) indicate that all coherent energy arrives within the first 10 to 15 s, depending on slab depth. The misfit for OBS stations will therefore only be sensitive to the early part of the receiver function waveforms.