Investigations of Ambient Noise Velocity Variations in a Region of Induced Seismicity near Greeley, Colorado

The Denver basin lies east of the southern Rocky Mountains in Colorado, and extends from the Rocky Mountain foothills to eastern Kansas and southern Wyoming (Higley and Cox, 2007). Oil‐and‐gas‐derived wastewater has been injected into the combined Denver disposal zone, a 500 m thick interval mainly composed of Pennsylvanian to Permian sandstones, for many years. The lowest basin unit—the Fountain Formation lies in contact with the Proterozoic Boulder Creek granodiorite of the Routt Plutonic suite—a metamorphic batholith. It has been hypothesized that the bottom of the Fountain Formation may be highly fractured, allowing hydraulic communication between the injection interval and the underlying basement faults (Yeck et al., 2016). Induced earthquakes in the Denver basin include the Rocky Mountain Arsenal earthquakes from the 1960s and 1970s (Healy et al., 1968), and a sequence of induced earthquakes east of Greeley, in Weld County, Colorado, starting in 2014 (Yeck et al., 2016). Pore‐pressure diffusion modeling indicates that 56% of the local pore pressure increase (totaling 0.15 MPa as of 2016) can be attributed to seven local injection wells (Brown et al., 2017). During 2014–2019, over 1000 low‐magnitude (M ≤3.2, with only 12 M >2.5) earthquakes were recorded in the region (Sheehan, 2020).

In this study, we use ambient noise interferometry to search for subsurface seismic velocity changes associated with wastewater injection in Northeastern Colorado. The motivation is to use the velocity changes as a proxy for poroelastic or other material property variations in the subsurface that can be linked to the onset of induced seismicity. We additionally apply the same analysis to a set of synthetic waveforms for several modeled velocity perturbation scenarios.

Poupinet et al. (1984) first used surface and coda waves of repeating earthquakes to demonstrate that seismic velocities can change slightly over time, evident in the minor offset in arrival of phases from identical source–receiver locations. Shapiro and Campillo (2004) demonstrated that cross correlations of ambient noise at two stations could approximate the Green’s function, which describes the seismic wave propagation between two seismic stations. Coda waves can also be reconstructed with cross correlation and have the benefit of sampling more of the subsurface than direct arrivals due to their scattered paths, and thus are more sensitive to material changes (Aki, 1957).

Ambient noise interferometry has been used in a variety of contexts and scales to detect geologic changes that would otherwise require costly active sources to measure. Sens‐Schönfelder and Wegler (2006) found an inverse correlation between precipitation and groundwater levels and seismic velocities. Pre‐ and postseismic changes have been observed, and have been attributed to healing and ground compaction following large earthquakes (e.g., Vidale and Li, 2003; Heckels et al., 2018). Interferometry has also detected changes associated with volcanic fluctuations (e.g., Duputel et al., 2009; Mordret et al., 2010; Yates et al., 2019).

In a system in which the medium between two stations has not experienced any change in velocity, the retrieved virtual seismic sources should show no difference. Velocity variations in most studies are typically on the order of 0.01%–1% dv/v, with higher values (∼10 s% dv/v) only apparent in shallow surveys in environments experiencing rapid change, such as permafrost thawing (James et al., 2017) and in which time resolution can be increased to detect short‐lived seismic velocity using converted phases (e.g., Lu and Ben‐Zion, 2022).

We analyze five years of continuous vertical‐component seismic data from June 2014 to July 2019 from the 16 broadband, intermediate, and short‐period seismic stations of the Weld County seismic network (Sheehan, 2016; Fig. 1). Waveform preprocessing included removal of instrument response, demeaning, detrending, filtering between 0.01 and 8 Hz, and spectral whitening.

We expect changes associated with injection to occur at depths of roughly 2–7 km, in which Brown et al. (2017) modeled the highest pore pressure increases. Based on the depth sensitivity of seismic surface waves, which are thought to form the dominant part of the ambient seismic noise field, we focus our analysis on the frequency band of 0.17–0.25 Hz of the cross‐correlation function (S5). We also analyze frequency bands of 1–4 Hz and 0.5–1 Hz, which have greater sensitivity to shallower changes (S6).

Once the seismograms have been preprocessed, stations are cross correlated and autocorrelated in 1 hr segments, stacked into daily average functions, and saved as 4 min long waveforms. We use the MSNoise software (Lecocq et al., 2014) for computing the noise cross‐correlation functions and measuring the seismic velocity changes. Following Clarke et al. (2011), we excluded cross correlations with an average signal‐to‐noise ratio below two. Single‐day cross correlations are generally unstable, so it is common to stack several tens of days; in this study we stack functions for 90 days. Large stacks limit the temporal resolution of seismic velocity changes, but we only seek a coarse relationship between seismic speeds and monthly injection volumes over a period of several years. We do not expect to see substantial short‐term changes in velocity associated with earthquake rupture, because the induced seismicity in the region is of small magnitudes and correspondingly low peak dynamic strains (e.g., Boschelli et al., 2021).

We measure travel‐time offsets in coda waves using the moving‐window cross‐spectral (MWCS) method (Clarke et al., 2011) as implemented in MSNoise. Short, overlapping time segments of the daily cross‐correlation function (short‐term average) are compared to the reference function (long‐term average). We use the 90‐day stacked correlation functions as the short‐term average and stack all correlations functions available for the long‐term average. For each frequency band, we analyze segments equivalent to the nearest integer of the longest period and a step between windows of half that period (e.g., a 6 s window and 3 s step for 0.17–0.25 Hz; full parameter table in the supplemental material, available to this article). For each step, the delay time between functions, dt, is calculated in the frequency domain. For each frequency band analyzed, the coda window in which delay times are analyzed is defined by multiples of the longest period, $Tmax$⁠, after Ugalde et al. (2014). Starting from the predicted surface‐wave arrival time (defined as 3 km/s), the minimum bound is $±2Tmax$⁠, and the maximum bound is $±5Tmax$⁠. The minimum window bound is chosen to exclude direct arrivals from waves and also eliminates scatterers near the receiver. The maximum window bound allows multiple wavelengths of all periods to be measured but excludes later times when the coda coherence decreases. The overall delay time, dt/t, is calculated from a weighted least squares regression through the dt points within the coda window and is used to compute dv/v (equation 1). Both the causal and acausal sides of the correlation function are used in the regression, and the line is not forced through the origin.

To further explore how localized velocity changes may be recovered by interferometry methods, we compute synthetic seismograms for multiple velocity model perturbations and compare the seismic velocity variations that are recovered from the MWCS method with the input perturbations. We expect that the modeled velocity change should be accurately recovered in the scenarios in which seismic velocity perturbations are uniform throughout the medium, but underestimated in the localized scenarios, due to the travel‐time delays occurring over part of the path between source and receiver. We compute synthetic seismograms using the frequency–wavenumber integration method implemented in computer programs in seismology (Herrmann, 2013). We simulate a Green’s function point source with a synthetic explosive source that originates at the surface and propagates through a 1D seismic velocity profile representative of the Denver basin region (Table S1). The same source is then propagated through 16 models in which $VP$ and $VS$ parameters have been perturbed by factors of ± 0.1%, 1%, 5%, and 10%, once in the 4–7 km depth interval and once in all layers. The perturbation at the 4–7 km depth mostly overlaps the depths at which Brown et al. (2017) modeled pore pressure increases due to wastewater injection. Modifying velocities from 2 to 4 km also may have been appropriate, but we wanted to avoid changing multiple layers in the model. The second scenario corresponds with a uniform velocity change, which is assumed by the relationship in equation (1). Synthetic seismograms are computed at stations 1–50 km from the source. Ambient noise cross correlations are generally unreliable from stations separated by less than 2–3 times the wavelength of interest (Luo et al., 2015) (⁠$λ=5.9km$ for 0.17 Hz), so interstation distances below 12 km are excluded from the analysis. Although the synthetic seismograms do not incorporate the seismic scattering that gives rise to seismic coda, which is the basis of the seismic noise measurements, they provide insight into the effects on direct seismic waves in a perturbed medium. We measure dv/v from the synthetic seismograms using the same MWCS and dt/tparameters used for the data from the Greeley network (Table S2), and do not target specific phases, although we restrict the analysis to the 0.17–0.25 Hz frequency band. One exception in the synthetic implementation is that we force the regression through the origin, because the correlation lacks an acausal side.

Figure 2 shows the 0.17–0.25 Hz dv/v results averaged across the network compared to seismicity and wastewater injection at high‐rate disposal wells near Greeley. We filter the dv/v time series to signals with periods greater than 3 months, in which injection‐related effects might arise, but no pattern is apparent (Fig. 2, blue line). There is little temporal correspondence between Greeley seismicity and subsurface seismic velocity variations (dv/v). Although approximately 1200 earthquakes were recorded on the Greeley network, only a dozen were greater than a magnitude 2.5, and seismic monitoring began only after the M 3.2 in June 2014—the largest event in the area. Although large earthquakes with extensive damage zones show distinct coseismic velocity drops and postseismic fault healing (e.g., Vidale and Li, 2003; Boschelli et al., 2021), the seismicity in this region is likely too small to cause significant elastic changes after a single event. Similarly, the injection history in and around Greeley does not show any clear correspondence with the overall dv/v time series. The challenge in comparing velocity changes to both the phenomena is that we lack a reference function representing the basin prior to injection or the onset of observed seismicity. In addition, the MWCS method may be ill‐suited for the narrow frequency band used, which limits the number of reliable dt points available for the linear regression used to calculate dt/t.

Many researchers have found that dv/v signals correlate with seasonal variations such as temperature, precipitation, and groundwater. Sens‐Schönfelder and Wegler (2006) found a negative correlation between precipitation and dv/v in Indonesia, with increasing precipitation causing reduced seismic velocities. These observations were for higher frequency correlations for which coda waves directly sampled the saturated medium. Many laboratory and empirical observations indicate that seismic P‐wave velocities increase when a medium becomes saturated. However, Clymer and McEvilly (1981) observed that shallow velocities increased as the water table fell, indicating saturation had reduced S wavespeeds. Thus, the overall effect of ground saturation is generally a reduction in dv/v.

We identify seasonal effects in seismic velocity by applying a low‐pass, zero‐phase filter to the dv/v time series to separate signals with periods of 1 yr or longer (Fig. 2a, green line). Because the lowest frequency of the ambient noise records corresponds to depth sampling of about 5 km, we can constrain the seismic velocity variations to this depth interval, which shows a consistent annual periodic signal from the years 2016 through 2019 with an amplitude of about 0.05% dv/v. Meier et al. (2010) found seasonal velocity variations in southern California and hypothesized that hydrologic loading from groundwater fluctuations could be responsible. To test the observations of Meier et al. (2010), Tsai (2011) modeled the effect of thermoelastic, poroelastic, and direct elastic (i.e., loading) variations on seismic speeds via strain amplitudes in response to groundwater and temperature recordings in Rossmoor, California. Of the three mechanisms, loading correlated best with the phase and amplitude of the observed variations, with groundwater levels having an inverse effect on velocities.

We examine the potential for hydrologic loading and thermoelastic effects to cause the seasonal variations in dv/v for the Greeley, Colorado, network by comparing the phasing of the seismic velocity variations with local water well and temperature measurements, and with the modeled results from Tsai (2011) from southern California. We compare modeled seismic velocity variations from southern California, due to hydrologic loading (Tsai, 2011; Fig. 3a, dashed yellow line), with the seismic velocity variations recorded by the Greeley network. Our observed annual velocity changes (Fig. 3a, green line) slightly precede the phase but match the amplitude of Tsai’s modeled velocity variations (∼0.05 dv/v%) attributed to hydrologic loading. The observed phase offsets could be due to differences in the timing of hydrologic loading effects; however, the exact phase of hydrologic loading in the vicinity of Greeley, Colorado, is difficult to assess, because many shallow wells lack a periodic signal and appear to be influenced by local irrigation practices associated with agriculture. To capture the overall groundwater signal in the region, we analyze 720 groundwater monitoring wells in the South Platte River drainage division (Colorado Department of Water Resources, 2021) and filter the results to time series with a peak periodicity between 0.5 and 1.5 yr, resulting in 125 wells with apparent seasonal variations. To standardize a variety of aquifer systems, depths, and magnitudes of groundwater fluctuations, we normalize recordings by first demeaning and then dividing each time series by its maximum absolute value to provide a sense of the phase of groundwater loading over time. Overall, the well data show a seasonal groundwater peak in the spring (Fig. 3b, gray lines), which indicates that the hydrologic loading signal is roughly consistent in phase with that modeled by Tsai (2011) (Fig. 3b, dashed blue line; he data have been normalized and shifted 15 yr for comparison).

Several authors have attributed annual variations in seismic velocities to seasonal thermoelastic strain. This phenomenon was first related to lateral Global Positioning System displacement signals by Ben‐Zion and Leary (1986), modeling strain fluctuations over tens of kilometers due to surface temperature variations. Tsai (2011) extended this relationship to seismic $VSv$ waves. We collect and plot daily surface temperatures measured at a weather station in Greeley, Colorado, for comparison with the dv/v time series (Fig. 3c; Colorado Department of Water Resources, 2021, station USC00053553). Surface temperatures peak each summer and exhibit phasing that closely resembles our annual dv/v signal from 2016 to 2019. The thermoelastic model of Tsai (2011) indicates that seismic velocity variations may lag behind surface temperature changes by about a month. Given the temporal uncertainty introduced by the 90‐day running average in our measurements, thermoelastic strain via seasonal surface temperature fluctuations could plausibly contribute to the seasonal variation in our results.

It is also possible that these variations are the result of changes in noise source distribution throughout the year due to seasonal changes in ocean activity, though the varying interstation azimuths from station pairs in the network mitigate the effects of directional noise sources on the recovered Green’s functions.

We measure dv/v from the synthetic seismograms for the uniform and depth‐limited seismic velocity perturbations using MSNoise, and plot the measured percent velocity change as a function of distance (Fig. 4). For both the uniform and depth‐limited seismic velocity variations, measured dv/v from the synthetic seismograms is at or near its minimum values at short interstation distances and generally increases with receiver distance. In the uniform model scenarios, in which all layer velocities are increased or decreased proportionally, measured dv/v reaches a peak of about 90% of the imposed value from 44 to 50 km. For the localized scenarios, in which only the velocity of 4–7 km depth profile is modified, measured dv/v generally reaches the maximum of 20% of the imposed value between 30 and 40 km, before decreasing again. Averages across all interstation distances are 67.0% and 11.3% for the uniform and localized scenarios, respectively. Although results are typically averaged over many station pairs and a variety of interstation distances, we would expect the measured values to be more consistent with respect to receiver distance. As mentioned earlier, the narrow filter band used may limit the number of dt points that the MWCS method can calculate, which may explain why the uniform velocity changes could not be fully resolved. In situ velocity variations are rarely known with great accuracy, so producing results with the correct sign and within an order of magnitude is still useful for characterizing processes that induce such changes.

Some considerations in our model set up follow: The synthetic seismograms neither produce coda waves from complex point scattering nor do they reflect direct and scattered waves from 2D and 3D structural effects that would be captured by other methods (e.g., Yang et al., 2022). Incorporation of stochastic velocity perturbations into the seismic velocity model, which is beyond the scope of the 1D modeling performed here, may provide a means to compute a realistic seismic coda that could be directly compared with observations. In addition, our simulation only tests the sensitivity of idealized Green’s functions rather than those derived from daily noise cross correlations. Directly simulating noise sources with variable distributions (i.e., Sager et al., 2022) would allow us to better understand the reliability of our cross correlations. Nevertheless, the modeling approach here helps to illustrate the inherent tradeoffs between the station spacings that are required for waves to intersect deep seismic velocity perturbations and the fraction of the path that accumulates travel‐time delay. Ultimately, retrieving velocity changes within the correct order of magnitude is likely to be beneficial for many monitoring applications, and further work would be needed to interpret the magnitudes of seismic velocity variations that are recovered from various observations.

Our scenarios indicate that broad, vertically localized changes can be detected with some confidence using interferometry techniques. However, the models do not simulate the effects of lateral changes. If wastewater injection does substantially alter seismic properties, these effects are likely concentrated near the disposal site in which pore pressure changes are the highest (e.g., Brown et al., 2017). Although an interstation distance of at least two times the wavelength being analyzed is recommended for reconstructing cross correlations from noise (Luo et al., 2015), our study indicates that dv/v results may be unreliable until spacings of several times the wavelength. At such distances, the cross correlations between stations will have sampled a much larger space than the anomaly of interest. This presents a mismatch between the scale at which potential injection‐related velocity changes may occur and the scale at which they can be reliably measured in our field study.

In the Greeley, Colorado, field study, we observed apparent annual velocity variations in the 0.17–0.25 Hz frequency band with consistent summer peaks on the order of 0.05% dv/v. The phase and amplitude of this signal as well as local groundwater and temperature trends match those observed by Meier et al. (2010) and modeled by Tsai (2011), indicating that our dv/v results primarily detected seasonal changes resulting from hydrologic loading and regional thermoelastic strain. Our results provide further verification of the effect of the environment on seismological measurements, which must be well understood to monitor more subtle processes in the subsurface with ambient noise.

We did not detect long‐term velocity changes or correlation with wastewater injection over the 5 yr of the study. Our modeling shows that velocity variations localized in depth can be detected, provided sufficient seismic velocity variation, but localized seismic velocity changes are inherently underestimated from data recorded at the surface. If these subtle changes were present in the dv/v signal, they might still be obscured by the seasonal changes observed. Further, potential medium changes due to wastewater injection may be even more localized than our model, due to the limited extent of pore pressure increases around injection wells. Finally, like many studies of induced seismicity, our data collection began in response to the onset of felt earthquakes, and, thus, we lack a reference cross‐correlation function that represents the state of the subsurface prior to injection and seismicity. These results indicate that future sites of wastewater disposal would benefit from long‐term monitoring of ambient seismic noise prior to the onset of wastewater injection to compute reference noise correlations for use in noise monitoring efforts and to characterize the environmental effects on seismic velocity that would need to be understood for identification of subtle and small‐amplitude signals. In situ measurements or use of other seismic phases (e.g., scattered waves, Lu and Ben‐Zion, 2022) could provide increased sensitivity to seismic velocity changes at depth.

All seismic data were downloaded through the Incorporated Research Institutions for Seismology (IRIS) webservices (https://service.iris.edu/), including the following seismic network: the XU (Sheehan, 2016). Wastewater injection data were obtained from https://cogcc.state.co.us/data.html (last accessed July 2020) and water well data were obtained from https://cdss.colorado.gov/ (last accessed February 2021). Data from Tsai (2011) were digitized from their figure 4. The codes from computer programs in seismology are available at https://www.eas.slu.edu/eqc/eqccps.html (last accessed November 2021). MSNoise is available at http://www.msnoise.org/ (last accessed May 2020). The supplemental material for this article contains velocity models, Rayleigh wave sensitivity kernels, additional results, and inputs for MSNoise and Computer Programs in Seismology.

The authors acknowledge that there are no conflicts of interest recorded.

This research was funded by the National Science Foundation Award Number 1520846 (Hazards SEES) and U.S. Geological Survey National Earthquake Hazards Reduction Program (NEHRP) program Award Number G18AP00079. Thomas Lecocq, Justin Rubinstein, and an anonymous reviewer provided useful reviews that improved the article. The authors thank Josh Boschelli and Alec Yates for insightful discussions that improved our methods, Kyren Bogolub for leading the 2016 field deployment, Enrique Chon for the Greeley earthquake catalog, and Bob Herrmann for assistance in computing the synthetic seismograms. Any use of trade, firm, or product names is for descriptive purposes only and does not imply endorsement by the U.S. Government.