## Abstract

Fault slip behavior during episodic tremor and slow slip (ETS) events, which occur at the deep extension of subduction zone megathrust faults, is believed to be related to cyclic fluid processes that necessitate fluctuations in pore-fluid pressures. In most subduction zones, a layer of anomalously low seismic wave velocities [low-velocity layer (LVL)] is observed in the vicinity of ETS and suggests high pore-fluid pressures that weaken the megathrust. Using repeated seismic scattering observations in the Cascadia subduction zone, we observe a change in the seismic velocity associated with the LVL after ETS events, which we interpret as a response to fluctuations in pore-fluid pressure. These results provide direct evidence of megathrust fault-valve processes during ETS.

## INTRODUCTION

Periodic phenomena of spatially and temporally correlated seismic tremor and aseismic slow slip are collectively referred to as episodic tremor and slip (ETS) and have been observed in several subduction zones (*1*–*3*). It is generally believed that ETS occurs along the subduction zone megathrust down dip of the seismogenic zone and that interaction between segments of the megathrust may be important for understanding increases in earthquake triggering (*4*). In most subduction zones, geophysical images reveal a landward-dipping seismically low-velocity layer (LVL) (*5*, *6*), with low electrical resistivity (*7*) in the vicinity of ETS that is interpreted as hydrated upper oceanic crust of the subducting plate and/or shear zone material. These studies, along with thermal and petrological models (*8*), provide evidence consistent with ETS occurrence in fluid-rich environments. Investigations on the mechanical properties of the LVL in subduction zones suggest near lithostatic pore-fluid (*P*_{f}) pressures facilitated by low permeability along the plate interface (*6*, *9*, *10*), which act to lower effective normal stress. This is supported by frictional modeling (*11*) and observations that suggest tremor sensitivity to small dynamic stress perturbations (*12*, *13*). At ETS depths (30 to 40 km), most fluids derive from prograde metamorphic dehydration reactions of subducting material (*8*). Studies of subduction zone forearcs suggest that fluids eventually migrate along the megathrust and into the overlying crust, near the mantle wedge corner (*14*–*16*). It has been suggested that these fluids may not only play a passive role in ETS by weakening faults via high *P*_{f} but that ETS may also be coupled with cyclic fluid processes (*14*–*18*). However, observations of these cyclic processes are limited and indirect (*17*, *18*). Characterizing temporal variations in physical properties in subduction zones in relation to ETS may help in understanding fluid processes and their effect on slip behavior along the plate interface.

Receiver functions (RFs), a type of teleseismic scattering data, have been used to investigate the seismic velocity structure of subduction zone forearcs in the ETS source region (*6*). In this study, we investigate temporal variations in RF data that may be associated with fluid-related *P*_{f} cycling during ETS. These data are sensitive to discontinuities in seismic shear-wave velocity (*v*_{s}) at vertical scale lengths of 1 to 10 km. RF data reveal the signature of the LVL characterized by double-polarity pulses of direct P-to-S–converted (Ps) and back-scattered (Pps and Pss) waves from the top and bottom boundaries of this layer (*9*). We use high-quality broadband seismic data from >1300 teleseismic earthquakes recorded between 1993 and 2019 at five long-running seismograph stations in the forearc of the Cascadia subduction zone (Fig. 1) and calculate vertically (*S*_{V}) and horizontally (*S*_{H}) polarized RFs. The raypaths of the scattered seismic waves that produce the RFs are nearly vertical. As a result, RFs are sensitive to the velocity structure directly beneath the recording station and sample the medium repeatedly. We use these repeated observations to constrain time-varying velocity structure in relation to ETS. The time resolution of RF data is low as they depend on the random and infrequent occurrence of large, distant earthquakes. However, combining RF attributes before and after 21 observed ETS events between 1994 and 2019 allows for observation of changes in material properties beneath the stations.

## RESULTS

The complexity of the structural environment can make it difficult to identify temporal changes in RF data. Dipping layers and/or anisotropy cause variations in the vertical incidence and back-azimuth angles of incoming waves, which results in periodicity of observed RF data with respect to the back azimuth of the recorded earthquakes. We account for this by performing a harmonic decomposition of all *S*_{V} and *S*_{H} data and produce a data-predicting model for a given input back-azimuth angle. We calculate the cross-correlation (*CC*) coefficients between observed and predicted RFs, remove outliers, and sort them chronologically. Using 15-day time bins before and after the midpoints of the ETS episodes determined using GPS data, we observe a statistically significant (95% confidence level) reduction of ~5 to 10% in the mean *CC* values of *S*_{V} data following ETS at three of the five stations (PGC, SNB, and VGZ; Fig. 2). This implies a systematic error in the harmonic decomposition model to accurately predict RFs following ETS and suggests that the velocity structure following ETS is different from the velocity structure averaged over the entire recording period. However, this analysis does not indicate where in the depth profile the change in velocity structure is occurring. No significant changes in *CC* values are observed for the *S*_{H} data, suggesting that changes occur predominantly in the background isotropic structure (see Materials and Methods).

To investigate the origin of the changes in *CC* values for *S*_{V} data, we convert the amplitudes of the direct Ps phases for each individual RF into estimates of shear-wave velocity contrasts (δ*v*_{s}) across the top and bottom of the LVL using a geometric model of the subducting slab (*19*). This accounts for effects of back-azimuth angle and slowness on slab-normal incidence (corrects for geometrical variations in the amplitude of RF data) and reduces the data information content to model parameters that reflect bulk structural properties of the LVL. After removing long-term amplitude variations in δ*v*_{s} estimates and outliers, we observe a statistically significant (95% confidence) increase in the mean δ*v*_{s} of ~0.05 to 0.1 km s^{−1} across the top boundary of the LVL following ETS for the same stations that revealed a change in *CC* values (Fig. 2). A reduction in the average δ*v*_{s} across the bottom boundary of the LVL is observed following ETS for these stations, although our confidence in rejecting the null hypothesis is lower in this case. We note that the two stations that do not reveal changes in velocity structure are those where the megathrust is at shallower depths. The lack of observation may be structural in nature, although it may also be a consequence of lower quality (and fewer) data.

Using longer duration time bins shows that the changes remain observable at statistically significant levels but at decaying magnitudes, suggesting that any change in the velocity structure is short lived (Fig. 2). However, we are unable to determine the absolute time scale for which the change persists with these data. Bootstrap analyses (by randomly resampling the dates of the events used to compute the RF data) of the *CC* and δ*v*_{s} estimates as a function of time bin width provide further evidence that the observed decaying change in these quantities is robust. The lack of scatter (relative to bootstrap samples) in the trend for estimated changes in δ*v*_{s} across the bottom boundary of the LVL suggests this observation may also be real. Considering changes in δ*v*_{s} across both the top and bottom boundaries would imply that the change is predominantly within the LVL. Last, we observe statistically significant changes in δ*v*_{s} estimates using higher frequency corners of RF data (up to 0.7 Hz), suggesting that they occur over a scale length of ~3.5 km.

## DISCUSSION

Recent studies that observe temporal changes in seismic velocities coincidentally with large earthquakes (*20*–*22*) interpret their results in terms of rock damage from strong ground motion, static stress redistribution, and/or fluid processes. Our observations of temporal changes in seismic velocities are associated with aseismic slow slip during ETS, which has a static stress drop several orders of magnitude smaller than regular earthquakes (*23*). ETS does not produce strong ground shaking, and seismically induced rock damage would only be expected in the uppermost part of the crust (*20*), far above the LVL. Furthermore, neither rock damage nor stress redistribution can explain the rapid recovery of our observed velocity changes. Alternatively, the “fault-valve” model (*24*, *25*) has been used to explain observations of changes in seismic velocities associated with large earthquakes along subduction megathrusts (*22*). This model describes breaching of low-permeability materials within the megathrust shear zone during slip, which may be trapping and pressurizing fluids. Breaching of these seals opens permeable pathways and allows for fluid migration within and across the megathrust. Low permeability is eventually re-established, and fluid pressurization re-accumulates. It has been suggested that aseismic slow slip (during ETS) may be linked to a similar model (*17*, *18*), which would necessitate periodic fluctuations in near lithostatic *P*_{f}.

We suggest that the observed *v*_{s} changes within the LVL are the result of temporary reductions in *P*_{f} following ETS (Fig. 3), consistent with the fault-valve model. This is supported by observations of temporal changes in bulk seismic attenuation coincident with inferred slow slip in Japan, which suggest the episodic fluid drainage of the plate interface (*18*). However, these observations are made in a colder subduction zone that does not display typical ETS behavior. Time-dependent measurements of the stress ratio within the subducting oceanic crust have recently been found to correlate with shallow slow slip events and suggest reductions in *P*_{f} on the order of several megapascals(MPa) (*17*). Assuming that all of our observed velocity change occurs within the LVL under high *P*_{f} conditions, we estimate that the increase in *v*_{s} of ~0.1 km/s corresponds to a *P*_{f} drop of ~1 to 10 MPa (*26*), consistent with estimates from stress ratio fluctuations (*17*). In addition, models of fluid compartmentalization within the fracture network of the subducting crust or shear zone may produce changes in *P*_{f} of several MPa when low-permeability barriers between compartments are ruptured (*27*, *28*). The fluid flux from the subducting slab at ETS depths is low given the short recurrence time of ETS (*8*). However, a fracture network model does not necessitate a high flux of fluids from the subducting plate, as high *P*_{f} can also be achieved through localized sealing of fluid compartments (*27*), which may be a result of fault healing or hydrothermal mineral precipitation (*14*–*16*).

Last, the observation that only stations located in the down dip extension of the slow slip source region experience a change in *v*_{s} correlates with greater fluid availability near the mantle wedge corner (*5*, *8*, *14*). Fluids may ultimately be mobilized into the forearc crust during this process (Fig. 3), consistent with Helium isotopic signatures from crustal fluids in Cascadia (*29*). Enhanced vertical fluid mobility following slow megathrust rupture would result in greater silica precipitation and explain the concentration of seismicity and reduced electrical resistivity in the forearc crust above the mantle wedge corner (*7*, *14*, *18*, *19*). These results lend support to models that favor hydrologic control of slow slip at the deep extension of megathrust faults.

## MATERIALS AND METHODS

### RF analysis

Data used in this study came from permanent broadband seismic stations of the Canadian National Seismograph Network located in the forearc of the Cascadia subduction zone (Fig. 1). We used 120-s-long, three-component recordings of the *P*-wave coda of teleseismic earthquakes with magnitudes above 5.8 that are located at epicentral distances between 30° and 90°. At each station, we compiled all available data between 1993 and 2019 with high signal-to-noise ratio (>5 dB) on the vertical component. Seismograms were decimated to 5 Hz, low-pass–filtered using a corner frequency of 1 Hz, rotated to the vertical-radial-transverse orientations, and then decomposed into up-going *P*, *S*_{V}, and *S*_{H} wave modes (*30*). The *P* mode was deconvolved from the *S* modes using Wiener spectral deconvolution (*31*) to obtain *S*_{V} and *S*_{H} component RFs. This deconvolution procedure has been shown to reduce the seasonal amplitude variations in the power spectral density function of RF data due to increased noise levels in the winter and is important in this context to avoid mapping seasonal variations into the ETS-related signal (*31*). RFs were then filtered at corner frequencies of 0.05 to 0.5 Hz. Figure S1 shows all *S*_{V} data for station PGC sorted by back azimuth, showing the signature of the LVL.

### Harmonic decomposition

In the presence of horizontal discontinuities in seismic velocity structure separating two isotropic layers, the timing and amplitude of the various phases depend only on the vertical incidence angle (characterized by the slowness vector) of incoming *P* waves. In this case, there is no energy on the *S*_{H} component of RFs; the signal is entirely contained on the *S*_{V} component and displays no variation with back azimuth. In the presence of either dipping discontinuities or anisotropy in seismic velocities, RF data will show variations in both timing and amplitude of the converted phases as a function of slowness and back azimuth (*32*–*34*). These patterns are generally characterized by periodic variations in amplitude with back azimuth observed in both *S*_{V} and *S*_{H} components, with a characteristic 90° phase shift between them. The harmonic decomposition method exploits this periodicity by decomposing RF data into a set of periodic functions with back azimuth that are described by sin(*k*ϕ) and cos(*k*ϕ), where *k* is the harmonic degree and ϕ is the back azimuth (*32*, *34*).

Before the harmonic decomposition, we converted all RF data from time to depth using an average seismic velocity profile under Vancouver Island (*35*). This step reduces the influence of small phase delays produced by variations in horizontal slowness (*34*). We specified a maximum depth of 60 km to avoid contamination by the free-surface reverberations and to include only the direct Ps conversions. Here, we limit the decomposition to *k* = 0, 1, 2, which captures most of the signal for simple models of dipping layers and low order anisotropy (*32*–*34*). The system of equations is written as

The term *A* captures the signal that is independent of back azimuth and corresponds to a simple mean of all *S*_{V} component data. The terms *B*_{1} and *B*_{2} describe RF signals that vary by one cycle over all back azimuths and mainly reflect either a dipping interface or hexagonal anisotropy with a plunging axis of symmetry. The terms *C*_{1} and *C*_{2} describe variations of two cycles over back azimuths and are normally interpreted as hexagonal anisotropy with a subhorizontal axis of symmetry.

We performed the harmonic decomposition using all available RF data by inverting Eq. 1 using singular value decomposition (fig. S2). The harmonic components combine to provide a data-predicting representation of the structure beneath each station, averaged over the station deployment time (>25 years). Once these components were obtained, we reversed the process and calculated the predicted *S*_{V} and *S*_{H} data for each teleseismic source in the RF dataset. This allowed us to quantitatively compare observed and predicted RFs by calculating *CC* coefficients as a function of teleseismic source chronologically ordered for both the *S*_{V} and *S*_{H} datasets separately (fig. S3). We then removed outliers within each dataset using a median absolute deviation threshold of 2.5.

### Interface velocity contrast

At an interface separating two materials with different seismic velocities, RF amplitudes vary according to the interface-normal incidence angle or interface-parallel slowness. For a dipping interface, corrections are therefore required to convert vertical incidence and back-azimuth angles into interface-normal incidence. We used a smooth slab interface model (*19*) and extracted the strike and dip angles of the slab surface to carry out the angle corrections. This procedure was performed by gridding the slab contours using a spline interpolation under tension with minimum curvature and measuring the directional derivative and slab dip at the locations of the seismic stations considered herein (fig. S4 and table S1). We then extracted the amplitudes of the main Ps phases on the *S*_{V} component from the top and bottom of the LVL to estimate the *S*-wave velocity contrast at those interfaces. We converted the Ps RF amplitudes to estimates of *S*-wave velocity contrasts for plane waves at near-normal incidence (*36*)*P*- and *S*-wave velocities, respectively, *p* is the slab-parallel slowness, *q*_{α} and *q*_{β} are the slab-normal slowness values for *P* and *S* waves, respectively, and ρ is density. Although the amplitudes are sensitive to both *S*-wave velocity and density contrasts, the equation is an order of magnitude more sensitive to the contrast in *S*-wave velocity (second term on right-hand side of Eq. 2). We therefore ignored variations in density contrast and simply inverted the amplitudes for Δβ/β, which we scaled to an absolute velocity contrast δ*v*_{s} = (Δβ/β)* using an average upper mantle *S*-wave velocity of β = 4.6 km s^{−1}. We estimated δ*v*_{s} for each RF and chronologically sorted these values. We then removed a third-order polynomial fit to the estimated δ*v*_{s} values to remove long-term variations and outliers using a median absolute deviation threshold of 2.5 (fig. S5).

### ETS times from GPS data

To estimate changes in the seismic velocity structure relative to ETS, we required precise estimates of ETS times from 1993 to 2019 obtained in a consistent way. In this case, however, we did not require the formal inversion of GPS data since we were only interested in the timing of slow slip in the vicinity of the seismic network. We estimated ETS start, middle, and end times using GPS data from station ALBH, located nearest station VGZ (Fig. 1). We used the NAM08 reference frame solution from UNAVCO, but similar results were derived using alternative reference models (e.g., IGS08). We then detrended and smoothed the raw east-component GPS times series using a rolling window average of 14 days, extracted shorter time series for dates surrounding known ETS events (*37*), and estimated ETS events from visual inspection (fig. S6). For each shorter ETS time series, we fit a fifth order polynomial curve and estimate its first and second derivatives. We accepted dates where the first derivative is approximately 0 to mark the beginning and the end of the GPS reversal (i.e., the start and end of the ETS event). Since these dates were obtained from smoothed and processed data, they are likely overestimates (i.e., errors of ~5 to 10 days), but with no consequence to our analysis. The date where the second derivative is approximately 0 was taken as the midpoint of the given event and represents the time of maximum slip rate during ETS that we used as reference in our seismic data analysis (fig. S6). These dates are summarized in table S2 and were also verified to correspond with times when tremors during slow slip are located within 50 km of station PGC (*38*, *39*). ETS start and end dates for the events in 1994 and 1995 could not be determined from the available GPS data, and the midpoint dates were determined independently (*37*).

### Temporal stacking

RF data sample a medium repeatedly beneath a seismic station; however, temporal sampling is highly irregular, and the results can be biased by the teleseismic source distribution (*31*, *40*). Furthermore, RFs may contain substantial noise that is difficult to quantify. Observing time-varying changes in seismic velocity structure using these data therefore requires stacking over multiple ETS events and evaluating the resulting stacks using statistical inference. Using the dates estimated from the GPS data (see above), we binned RF times into pre- and post-ETS events for all events in the catalog. We used the midpoint of ETS events as reference, and defined time bins of 15, 30, 45, 60, 90, 120, and 150 days. We then stacked the RF-derived quantities (i.e., either the *CC* or δ*v*_{s} values) to estimate the sample mean values and SEs. We performed a two-sample *t* test for the difference in sample mean values and calculated *P* values to evaluate the statistical significance (*P* values for station PGC are provided in table S3). Taking the midpoint of individual slow slip events also implies that any potential bias due to tremor activity on the individual waveforms will similarly affect the pre- and post-ETS time bins. We verified that randomly perturbing the ETS times using a normal distribution with SD of ±5 days does not change the statistical significance of our results.

### Boostrap analysis

In addition to the *t* test, we further tested the null hypothesis by randomly shuffling the dates of the RF data (without replacement) 20,000 times and reprocessing the randomized datasets for both *CC* and δ*v*_{s} estimates. This was performed for all time bins considered and shown in Fig. 2 for station PGC. We calculated the 95% confidence envelope from the randomized datasets and found that the difference in mean values is highly scattered when considering single randomized datasets as a function of time bin duration. This scattering is not observed with the real datasets, further confirming that our results are not due to chance. Figure S7 shows the bootstrap analysis for stations LZB, PFB, SNB, and VGZ for the change in δ*v*_{s} values. Among these, stations SNB and VGZ show trends that are consistent with those observed at station PGC (i.e., low scatter; comparable but lower change in δ*v*_{s}). These stations are located along similar (or slightly greater) plate interface contours. Results for station LZB and PFB are not statistically different than a purely random process.

## SUPPLEMENTARY MATERIALS

Supplementary material for this article is available at http://advances.sciencemag.org/cgi/content/full/6/4/eaay5174/DC1

Fig. S1. *S _{V}* RFs for station PGC.

Fig. S2. Harmonic components for station PGC.

Fig. S3. Cross-correlation coefficients between observed and predicted RFs.

Fig. S4. Example calculation of strike and dip angles at each station location.

Fig. S5. Estimated *S*-wave velocity contrast across the top and bottom of the LVL.

Fig. S6. GPS data used in the determination of ETS event times.

Fig. S7. Changes in δ*v*_{s} across the top and bottom boundaries of the LVL.

Table S1. Strike and dip angles of the plate interface at each station location.

Table S2. Estimated start, end, and midpoint dates of ETS events from 1994 to 2019.

Table S3. *P* values from the *t* tests for the change in *CC* coefficients and δ*v*_{s} values following ETS events measured at station PGC.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is **not** for commercial advantage and provided the original work is properly cited.

## REFERENCES AND NOTES

**Acknowledgments:**We thank H. Dragert and R. Bürgmann for early discussions on this work. We also thank M. Bostock for useful suggestions and for sharing tremor density and slab geometry data.

**Funding:**This work is supported by the Natural Science and Engineering Research Council (Canada) through a Vanier Canada Graduate Scholarship to J.M.G., a Discovery Grant to P.A., and a Canada Graduate Scholarship to S.G.M.

**Author contributions:**J.M.G. and P.A. conceived the hypothesis and scope of the project. P.A. performed RF analysis. C.E. calculated slab geometry from published sources. M.M. estimated ETS times from available GPS data and tremor catalogs. S.G.M. and A.J.S. considered tremor and seismicity catalogs. A.J.S. compiled seismic data. J.M.G. wrote the paper with substantial input from P.A. and additional input from all co-authors.

**Competing interests:**The authors declare that they have no competing interests.

**Data and materials availability:**Seismic waveform data for the Canadian National Seismograph Network are available at the Canadian National Data Centre at https://earthquakescanada.nrcan.gc.ca/stndon/CNDC/index-en.php. Continuous GPS data are available at www.unavco.org/data/data.html. Cascadia tremor locations are available at https://pnsn.org/tremor. Worldwide tremor locations are available at www-solid.eps.s.u-tokyo.ac.jp/idehara/wtd0/Welcome.html. All other data are available in the main text, Supplementary Materials, or may be requested from the authors.

- Copyright © 2020 The Authors, some rights reserved; exclusive licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. Distributed under a Creative Commons Attribution NonCommercial License 4.0 (CC BY-NC).