Research ArticleEARTH SCIENCES

Monitoring transient changes within overpressured regions of subduction zones using ambient seismic noise

See allHide authors and affiliations

Science Advances  08 Jan 2016:
Vol. 2, no. 1, e1501289
DOI: 10.1126/sciadv.1501289


In subduction zones, elevated pore fluid pressure, generally linked to metamorphic dehydration reactions, has a profound influence on the mechanical behavior of the plate interface and forearc crust through its control on effective stress. We use seismic noise–based monitoring to characterize seismic velocity variations following the 2012 Nicoya Peninsula, Costa Rica earthquake [Mw (moment magnitude) 7.6] that we attribute to the presence of pressurized pore fluids. Our study reveals a strong velocity reduction (~0.6%) in a region where previous work identified high forearc pore fluid pressure. The depth of this velocity reduction is constrained to be below 5 km and therefore not the result of near-surface damage due to strong ground motions; rather, we posit that it is caused by fracturing of the fluid-pressurized weakened crust due to dynamic stresses. Although pressurized fluids have been implicated in causing coseismic velocity reductions beneath the Japanese volcanic arc, this is the first report of a similar phenomenon in a subduction zone setting. It demonstrates the potential to identify pressurized fluids in subduction zones using temporal variations of seismic velocity inferred from ambient seismic noise correlations.

  • Ambient seismic noise
  • Subduction zones
  • Pore fluid pressure


Processes accompanying earthquakes can produce changes in the mechanical properties and, therefore, seismic velocities of the surrounding crustal material. Dynamic and static stress changes, fluid migration, fault zone damage and healing, and nonlinear response of near-surface material to strong ground motion have all been implicated in causing changes to measured seismic velocities following earthquakes (13). These time-dependent velocity changes have been detected from measurement of time lags within repeating earthquake sequences (4) or from cross-correlation of continuous ambient seismic noise recordings (5). Although temporal variations in seismic velocity accompanying earthquakes have revealed information about the mechanical behavior of the shallow surface and near–fault zone environment, susceptibility of a region to coseismic velocity change has recently been used to characterize crustal rheology farther from the earthquake source. In particular, ambient noise cross-correlations showed that the strongest velocity reductions associated with the 2011 Tohoku earthquake occurred under the volcanic arc, rather than in the epicentral region (6). The velocity reduction was attributed to dynamic stress associated with seismic waves fracturing the fluid-pressurized weakened crust. We use a similar method to document seismic velocity reductions in the forearc crust of the Nicoya Peninsula following the 2012 Mw (moment magnitude) 7.6 underthrusting earthquake. Correspondingly, we ascribe this to an increase in fracture density in the high–pore fluid pressure, weakened crust due to dynamic stresses. Our observations extend the use of seismic velocity change susceptibility to identify pressurized fluids to subduction zones. The presence of such fluids is important in controlling plate boundary behavior and magma genesis; however, their distribution is poorly constrained by existing methods.

The Nicoya Peninsula in northern Costa Rica (Fig. 1) has been recognized as a natural subduction zone laboratory because of its unique occurrence directly above the plate interface, where the Cocos Plate subducts beneath the Caribbean Plate at a rate of 85 mm year−1 along the Middle American Trench (7). This location resulted in excellent seismic and continuous Global Positioning System coverage of the September 2012 (Mw 7.6) earthquake (Fig. 1), making it one of the best-recorded megathrust earthquakes ever (8). Several studies have shown that the seismogenic zone below the Nicoya Peninsula is characterized by strong lateral variations in thermal structure (9) and earthquake behavior due to differences in the origin of the subducting lithosphere (10). In addition to large earthquakes occurring every 50 to 60 years, aseismic, slow-slip events take place regularly both up- and down-dip of the seismogenic zone (11, 12). Water released from hydrous minerals within the subducted slab likely plays various roles in these geodynamic processes: it modifies the mechanical properties of the surrounding crustal material, decreases fault strength (reduces the effective normal stress), and controls the episodic seismic and aseismic activity of the fault (13, 14).

Fig. 1 Map of northern Costa Rica showing seismic velocity reductions following the 2012 Mw 7.6 Nicoya Peninsula, Costa Rica earthquake.

The white dashed line marks the position where the origin of the subducting oceanic lithosphere changes from warmer Cocos Nazca Spreading center (CNS) to colder East Pacific Rise (EPR) (28). (A) The focal mechanism (~21-km centroid depth) and area of concentrated slip for the mainshock [dashed black contour (29)] are indicated. Black inverted triangles correspond to the broadband seismic stations used in this study. Triangles show the location of strong ground motion sensors from the Universidad de Costa Rica with their color indicating peak ground velocity (PGV). Colored lines indicate coseismic velocity reductions obtained in the period band 1 to 3 s along the interstation path. (B) Colored lines indicate seismic velocity reductions obtained in the period band 3 to 10 s. Inverted triangles indicate the Vp/Vs values obtained for the forearc crust from receiver function analysis (30). Note that there were no relative coseismic velocity increases in the frequency ranges in our analysis.


Symmetric Green’s functions (GFs) with little amplitude variation as a function of time (indicating relatively stationary noise sources with little seasonal variation) were obtained from the cross-correlation of daily ambient seismic noise between most station pairs (Fig. 2). These were used to compute crustal velocity changes by measuring the time delay between the coda of the reference Green’s function (RGF, computed by stacking all GFs before the earthquake) and 30-day stacks of the daily GFs. The GFs and their coda are composed of surface waves with sensitivity to velocity that depends on period (short periods are sensitive to shallow crustal structures and longer periods are sensitive to deeper structures) (fig. S1). To differentiate velocity changes as a function of depth, we calculated time shifts relative to the corresponding RGF using different period bands: 1 to 10, 1 to 3, 3 to 10, and 5 to 10 s.

Fig. 2 Cross-correlation of ambient seismic noise.

Daily GFs filtered between 1 and 10 s for the station pair INDI-SAJU. The black trace represents the reference GF (RGF). The dotted line indicates the time of the Nicoya Peninsula earthquake. Note the high symmetry between the causal and acausal signals.

Our results reveal frequency-dependent coseismic velocity reductions that vary spatially. Analysis of the shortest period band (1 to 3 s) indicates that the largest velocity reduction of 0.4% occurs between station pairs that lie adjacent to the area of maximum slip associated with the 2012 Nicoya earthquake (Fig. 1A). In this period range, the velocity reductions are constrained to the upper 5 km of the crust (fig. S1). The precise timing of the velocity reduction onset is not well resolved because of the use of 30-day GF stacks to stabilize measurements of time shifts relative to the RGF for individual station pairs (Fig. 3A). However, the progressively earlier onset times obtained when using GF stacks of fewer days indicate that the velocity reduction occurs very shortly after the earthquake. On the basis of the timing, location, and shallow depth of the velocity reductions, we attribute them to near-surface damage produced by strong ground shaking accompanying the 2012 Nicoya earthquake. Static stress changes following the earthquake are also a possible cause of the coseismic velocity reduction. Dilatational volumetric strain changes are predicted at shallow depths in the region of the velocity reduction (fig. S2). Dilatational strain would favor the opening of fractures and the reduction of seismic velocity following the earthquake. Our results indicate that less than 50% of this coseismic velocity reduction was recovered a year after the mainshock (Fig. 3A). Several studies using ambient noise–based monitoring have reported reductions in seismic velocities following large earthquakes (15, 16); most of these reductions were relatively small (<0.4%), located at shallow depths (upper 2 km of the crust), and correlated with peak coseismic slip and/or peak ground velocity. A nonlinear response of near-surface material to strong ground motion was implicated in most of these studies. A similar magnitude and timing of velocity recovery to our results was obtained following many of these other earthquakes (17, 18).

Fig. 3 Seismic velocity reductions, GFs, and time delays associated with the 2012 Nicoya Peninsula earthquake for the station pair INDI-SAJU in the following period bands: 1 to 10, 1 to 3, and 3 to 10 s.

(A) The preseismic, coseismic, and postseismic velocity variations between April 2012 and October 2013. The dashed line indicates the time of the mainshock. Note the near-identical velocity reduction in the 1- to 10-s and 1- to 3-s period bands and the lack of a velocity variation in the 3- to 10-s band. (B) RGF, daily GFs, and time delays resulting from the cross-correlation of their coda filtered between 1 and 10 s: 5 days before and 15 days after the earthquake. The longer time period used after the earthquake was required to avoid mixing of postseismic signals or aftershocks.

We find even stronger coseismic velocity reductions (0.4 to 0.6%) in the 3- to 10-s period band that also show spatial variability (Fig. 1B). Rayleigh waves in this frequency range are sensitive to velocity structure in the 5- to 15-km depth range (fig. S1). Velocity reductions at longer periods are only found in the southeastern portion of the Nicoya Peninsula, between stations located farther from the coseismic slip patch than those having the shorter period velocity reduction (Fig. 1). Analysis of these stations in the 5- to 10-s period band shows the same magnitude velocity reduction, whereas the 1- to 3-s period band indicates no velocity change (Fig. 4A), confirming that the velocity reduction occurs at depths below 5 km. This deeper (5 to 15 km) velocity reduction is located in the region where a previous receiver function study found extremely high forearc Vp/Vs ratios that were interpreted as pore fluid overpressures. The Vp/Vs ratios determined at each station (Fig. 1B) reflect the average value for the entire forearc crust from the top of the subducting plate to the seismic station, a crustal section ranging from 16 km at the coast to about 30 km at the most inland station. Further support for the presence of fluids in the forearc in this region comes from results of an onshore-offshore magnetotelluric study conducted in 2007–2008 (19). This work imaged a low electrical resistivity anomaly (very conductive zone) at depths of 12 to 15 km in the southeast portion of the Nicoya forearc that was interpreted as the presence of fluids derived from clay mineral dehydration reactions in the subducted crust.

Fig. 4 Seismic velocity reductions, GFs, and time delays associated with the 2012 Nicoya Peninsula earthquake for the station pair INDI-LAFE in the following period bands: 1 to 10, 1 to 3, 3 to 10, and 5 to 10 s.

(A) The preseismic, coseismic, and postseismic velocity variations between February 2012 and June 2013. The dashed line indicates the time of the mainshock. Note the similarity in the velocity reductions in the 1- to 10-s, 3- to 10-s, and 5- to 10-s period bands, confirming that the physical processes responsible are occurring at depths of ~5 to 15 km. The delay in the velocity reduction apparent in the 5- to 10-s period band reflects the longer, 50-day stacks required to stabilize the analysis. (B) RGF, daily GFs, and time delays resulting from the cross-correlation of their coda filtered between 3 and 10 s: 5 days before and 15 days after the earthquake.

We propose that the deep coseismic velocity reduction results from overpressured fluids reducing the effective stress and weakening the forearc crust, making it more susceptible to damage (increasing fracture density) from dynamic stresses. Static stress changes or postseismic deformation following an earthquake may also be important in modifying seismic velocities at depth. Although it is unclear how crack opening and closure may operate under static stress (or strain) perturbations at depth, recent work reporting velocity changes associated with slow-slip events in Mexico (20, 21) indicates that slow stress changes are capable of affecting elastic properties of the deep crust. The negative volumetric strain (compression) in the southeastern portion of the Nicoya Peninsula at depths between 0 and 15 km, generated by the 2012 mainshock (fig. S2), is inconsistent with the velocity reductions we observed. To evaluate deformation related to afterslip as a possible source of the observed velocity reductions, we excluded from analysis the 70-day period following the earthquake that has previously been identified as the time constant associated with afterslip (22). Results in the longest period bands (3 to 10 and 5 to 10 s) still show a maximum reduction in velocity of ~0.5% with respect to the premainshock levels. This is comparable to results obtained using the entire postseismic time period and indicates that deformation associated with afterslip is not a main contributor to the deep velocity reduction. Therefore, we favor the interpretation that dynamic stresses from seismic waves increased crack damage and decreased seismic velocities in areas that were already highly weakened from overpressured fluids.

Unlike the shallower velocity reductions that attain their peak value shortly following the mainshock and recover only slightly in the following year (Fig. 3A), the onset of the deeper velocity reductions also closely follows the earthquake but such reductions take several months to attain their peak value. They also appear to recover more than 50% of their velocity reduction in the year following the earthquake (Fig. 4A). This corroborates our conclusion that different processes are responsible for the deep and shallow velocity reductions. It also indicates that some permanent deformation is likely occurring in the near-surface crust whereas a slow recovery process operates at depth. Our interpretation for the deeper velocity reductions under the Nicoya Peninsula is similar to that of a previous study implicating fluid-pressurized crustal weakening from seismic waves as the cause of velocity reduction beneath the Japanese volcanic arc. Ours is the first observation of this mechanism operating in a subduction zone environment and demonstrates the potential to identify pressurized fluids in subduction zones using temporal variations of ambient seismic noise.


Cross-correlation of ambient seismic noise

Crustal seismic velocity variations can be measured by the long-term cross-correlation of ambient seismic noise. A growing number of theoretical and experimental efforts have shown that the coherence of diffuse wave fields (for example, coda waves or ambient seismic noise) can be used to retrieve direct seismic waves (GFs) between two points on Earth’s surface. This reconstruction is related to the principles of reciprocity and time reversal of the random wave field and, therefore, directly linked to the mechanical properties of the medium.

Here, we used the cross-correlation of ambient seismic noise (23, 24) (diffuse wave field composed of surface waves) to characterize crustal seismic velocity variations following the September 5, 2012, Nicoya Peninsula earthquake. We analyzed 2 years of continuous seismic records, between April 2012 and April 2014, using all stations of a local broadband network that contained no significant data gaps during this period. We used MSNoise ( (25) to compute daily vertical cross-correlation functions or GFs (seismic response of Earth’s crust if an impulse is applied to each station) between all station pairs.

The daily seismic records were demeaned, tapered, and corrected for their instrument response. We applied a band-pass filter of 8 to 0.01 Hz and downsampled from 100 to 20 Hz as a preprocessing stage. To reduce the effects of earthquake signals, we applied 1-bit normalization and spectral whitening in the period ranges of our analysis: 1 to 10, 1 to 3, 3 to 10, and 5 to 10 s. The spectral normalization acts to enhance the band of ambient noise signal in the cross-correlations and reduce possible degradation caused by persistent monochromatic sources. We calculated daily cross-correlation functions for the 10 station pairs in the period ranges listed above to retrieve Rayleigh waves. For these period ranges and a local one-dimensional velocity model (26), Rayleigh waves are sensitive to velocity variations occurring from the uppermost crust to a depth of 15 km (fig. S1).

We defined an RGF by stacking all daily GFs before the earthquake to ensure that the RGF is representative of a background premainshock value while the daily GFs contain the information on the current state of Earth’s crust. To avoid mixing of preseismic and postseismic signals, we removed the days containing the mainshock and the largest aftershocks.

Crustal seismic velocity variations

To estimate the crustal velocity changes between all station pairs, we adopted the Moving-Window Cross-Spectral Analysis (MWCS) method (27). This method computes the time delay (∂t) between the RGF and 30-day GF stacks from the phase of their cross-spectrum. Our application of the MWCS uses a moving time window of 10 s with an overlap of 80%. The 30-day GF stack window length was used to reduce the uncertainty of our seismic velocity measurements for the individual station pairs. To ensure that our measurements are independent of variations in noise sources, we measured travel time delays within the coda of the GFs, which are primarily composed of surface waves. Assuming an isotropic medium, propagation time varies proportionally to distance between station pairs, and the phase shift between the records varies linearly with time. Accordingly, the relative crustal velocity change ∂v can be measured from the slope of the time delay ∂t as a function of lag time: ∂t/t = − ∂v/v.


Supplementary material for this article is available at

Fig. S1. Rayleigh wave velocity sensitivity kernels in our study region.

Fig. S2. Volumetric strain changes. (A to D) Volumetric strain changes generated by the Mw 7.6 2012 Nicoya Peninsula earthquake at depths of (A) 0 km, (B) 5 km, (C) 10 km, and (D) 15 km calculated using the finite fault model provided by the United States Geological Survey ( and Coulomb 3 software (31).

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.


Acknowledgments: We thank M. Protti and V. González from OVSICORI-UNA (Observatorio Vulcanológico y Sismológico de Costa Rica Universidad Nacional) for installing and operating the Nicoya seismic network with critical support from A. Newman and D. Sampson. A. Moya from LIS-UCR (Laboratorio de Ingeniería Sísmica, Universidad de Costa Rica) provided the peak ground velocity measurements. Discussions with P. Fulton, K. Wang, L. Screaton, J. F. Pacheco, and T. Lecocq as well as comments from two anonymous reviewers greatly improved the manuscript. Funding: This work was supported by the NSF under grant no. EAR-1321550 to S.Y.S. and by a UNAVCO COCONet (Continuously Operating Caribbean GPS Observational Network) Fellowship to E.J.C. Author contributions: E.J.C. performed all the computations; E.J.C. and S.Y.S. analyzed the data, interpreted the results, and wrote the paper. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All seismic data required to evaluate the conclusions in the paper are available at the IRIS DMC (Incorporated Research Institutions for Seismology Data Management Center; All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.
View Abstract

Stay Connected to Science Advances

Navigate This Article