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 (Pf) 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 Pf 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 Pf cycling during ETS. These data are sensitive to discontinuities in seismic shear-wave velocity (vs) 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 (SV) and horizontally (SH) 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.
Station locations considered for RF analysis are shown by triangles (solid triangles mark stations where changes in velocity structure are observed). The yellow square marks the location of the GPS station used to estimate ETS timing. Contours mark the depth of the plate interface inferred from the deformation front, up, and downdip limits of tremor and volcanic arc (19). The gray shaded region delineates the ETS source region interpreted from tremor density (19). The inset shows the tectonic setting of the Cascadia subduction zone.
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 SV and SH 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 SV 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 SH data, suggesting that changes occur predominantly in the background isotropic structure (see Materials and Methods).
Changes in mean CC values for 15-day time bins before and after ETS are shown for SV (A) and SH (B) data. Changes in mean δvs before and after ETS are shown for the top (C) and bottom (D) boundaries of the LVL. Mean values and estimated standard errors in (A) to (D) are shown by horizontal lines and shaded regions, respectively. Changes in RF data as a function of time bin duration are shown for CC analysis (E) and δvs analysis (F). Bootstrap resampling of RF times relative to ETS are shown by gray lines (500 of 20000 samples), and 95% confidence levels are shown by dashed lines in (E) to (F).
To investigate the origin of the changes in CC values for SV data, we convert the amplitudes of the direct Ps phases for each individual RF into estimates of shear-wave velocity contrasts (δvs) 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 δvs estimates and outliers, we observe a statistically significant (95% confidence) increase in the mean δvs 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 δvs 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 δvs 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 δvs across the bottom boundary of the LVL suggests this observation may also be real. Considering changes in δvs across both the top and bottom boundaries would imply that the change is predominantly within the LVL. Last, we observe statistically significant changes in δvs 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 Pf.
We suggest that the observed vs changes within the LVL are the result of temporary reductions in Pf 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 Pf on the order of several megapascals(MPa) (17). Assuming that all of our observed velocity change occurs within the LVL under high Pf conditions, we estimate that the increase in vs of ~0.1 km/s corresponds to a Pf 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 Pf 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 Pf 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).
During interslip periods (A), the deep extension of megathrust faults where ETS occurs is characterized by high Pf. During ETS (B), low-permeability seals within the plate interface shear zone are breached, which opens permeable pathways for fluid migration within and across the megathrust, causing a drop in Pf. Low-permeability barriers are re-established, and the cycle repeats. This figure is adapted from Peacock et al. (10).
Last, the observation that only stations located in the down dip extension of the slow slip source region experience a change in vs 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, SV, and SH wave modes (30). The P mode was deconvolved from the S modes using Wiener spectral deconvolution (31) to obtain SV and SH 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 SV 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 SH component of RFs; the signal is entirely contained on the SV 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 SV and SH 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 SV component data. The terms B1 and B2 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 C1 and C2 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 SV and SH 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 SV and SH 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 SV 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)
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 δvs 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 δvs 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 δvs 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 δvs). 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. SV 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 δvs 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 δvs 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
- 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).
