Monitoring terrestrial water storage, drought and seasonal changes in central Oklahoma with ambient seismic noise

Signiﬁcant imbalances in terrestrial water storage (TWS) and severe drought have been observed around the world as a consequence of climate changes. Improving our ability to monitor TWS and drought is critical for water-resource management and water-deﬁcit estimation. We use continuous seismic ambient noise to monitor temporal evolution of near-surface seismic velocity, dv/v, in central Oklahoma from 2013 to 2022. The derived dv/v is found to be negatively correlated with gravitational measurements and groundwater depths, showing the impact of groundwater storage on seismic velocities. Seasonal cycling of dv/v follows atmospheric temperature changes with a phase shift, which can be explained by thermo-elastic strain in the uppermost crust and sedimentary cover. The occurrences of droughts appear simultaneously with the local peaks of dv/v, demonstrating the sensitivity of near-surface seismic velocities to droughts. The results illustrate the potential of using seismic data

manuscript submitted to Geophysical Research Letters Abstract Significant imbalances in terrestrial water storage (TWS) and severe drought have been observed around the world as a consequence of climate changes.Improving our ability to monitor TWS and drought is critical for water-resource management and waterdeficit estimation.We use continuous seismic ambient noise to monitor temporal evolution of near-surface seismic velocity, dv/v, in central Oklahoma from 2013 to 2022.The derived dv/v is found to be negatively correlated with gravitational measurements and groundwater depths, showing the impact of groundwater storage on seismic velocities.
Seasonal cycling of dv/v follows atmospheric temperature changes with a phase shift, which can be explained by thermo-elastic strain in the uppermost crust and sedimentary cover.The occurrences of droughts appear simultaneously with the local peaks of dv/v, demonstrating the sensitivity of near-surface seismic velocities to droughts.The results illustrate the potential of using seismic data for monitoring TWS and drought at regional to local scales.

Plain Language Summary
Terrestrial water storage (TWS) is fundamental to the well-being of inhabitants on Earth.However, current approaches to measure TWS variations have limited temporal or spatial resolution.In this study, we use near-surface seismic velocity variations, dv/v, derived from continuous seismic recordings to monitor changes of TWS in central Oklahoma.A negative correlation between the long-term trend of dv/v with gravity measurements reflects the impact of groundwater recharge/discharge on near-surface seismic velocity.In addition, a seasonal cycling of dv/v has similar periodicity to recordings of air temperature, which can be explained by thermo-elastic strain at the subsurface.Comparisons between dv/v and drought index further show the possibility of using near-surface seismic velocity as a proxy for monitoring severe drought for local communities.Considering the high temporal sampling and flexible spatial deployment, seismometers may be used to monitor subsurface water distributions.This can be useful for sustainable water management and reliable water-deficit estimation.
The groundwater and surface water are dominant prerequisites for agricultural irrigation, and have important social-and economical-impacts to modern society (Rodell et al., 2009;Famiglietti et al., 2011;Scanlon et al., 2012).Severe climate changes have also led to frequent occurrences of drought and flooding.In addition, increasing demands on water resources for economical and social developments have imposed tremendous stresses on TWS (Feng et al., 2013;Long et al., 2013).Hence, it is critical to accurately and timely monitor the change of TWS in order to maintain sustainable water-resources management and rational water-inadequacy estimation (Alsdorf & Lettenmaier, 2003).
Direct measurements of TWS, such as installing gauges in wells, can accurately assess the levels of aquifers.However, due to the expense of well drilling and instrument maintenance, a limited number of wells typically lead to TWS estimates with a low spatial resolution (Alsdorf & Lettenmaier, 2003).In contrast, remote sensing techniques, such as Global Positioning System (GPS) and Interferometric Synthetic Aperture Radar (InSAR), allow measuring the deformation of the Earth's surface, which can then be used to infer groundwater variations and underground fluid migration (Bawden et al., 2001;Argus et al., 2005).In spite of the high temporal resolution of GPS measurements (K.H. Ji & Herring, 2012), its point-sampling characteristics fail to accurately map lateral heterogeneity of TWS.On the other hand, radar and InSAR can be used to delineate surface deformation with a much higher spatial resolution (Watson et al., 2002;Lanari et al., 2004), but they typically suffer from comparatively low temporal resolution, which depends on the orbital frequency of satellites.Since 2003, the Gravity Recovery and Climate Experiment (GRACE), a joint mission supported by the National Aeronautics and Space Administration (NASA) and the German Aerospace Center, provides detailed information on gravity variations on the Earth's surface by measuring the position changes of twin satellites (Tapley, Bettadpur, Watkins, & Reigber, 2004).These temporal mass changes near the Earth's surface primarily result from near-surface water circulation and migration.As a state-of-the-art technique, GRACE measurements have been widely used to assess TWS variations from regional to global scales (Rodell et al., 2007;Tian et al., 2017), and also enable estimating ice sheet and glacier melting in Greenland and Antarctica (Velicogna & Wahr, 2006a, 2006b;Harig & Simons, 2012).However, its limited temporal sampling (one month) and spatial resolution (around 300−400 km) (Tapley, Bettadpur, Watkins, & Reigber, 2004) cannot satisfy the urgent requirement of managing sustainable water-resource at regional to local scales.
There have been many applications for time-lapse seismic velocity analysis of water injection (Landrø, 2001;Nakata et al., 2022), reservoir monitoring (Rickett & Lumley, 2001;Angerer et al., 2002), and carbon capture and storage (D.Lumley, 2010;Zhu et al., 2019).Because of data availability and computational cost, these classical approaches to monitor temporal changes of the seismic velocities is limited by the occurrence of earthquakes in small regions, which commonly leads to sparse temporal sampling (D.E. Lumley, 2001;Kamei & Lumley, 2017).Over the past decades, continuously recorded seismic ambient noise has been widely used to delineate spatial and temporal variations of seismic velocities within the Earth's crust and upper mantle.It alleviates the restriction of using ballistic wave propagation between earthquakes and seismometers (Campillo & Paul, 2003;Shapiro et al., 2005;Lin et al., 2008;Yao & Van Der Hilst, 2009;Ritzwoller et al., 2011), thus enabling us to image Earth's structure for tectonically inactive regions.
Since 2008, central Oklahoma has experienced a significant increase in seismicity, some of which led to infrastructure destruction, such as the 2011 Mw 5.7 Prague earthquake (Sun & Hartzell, 2014;Sumy et al., 2017) and the 2016 Mw 5.8 Pawnee earthquake (Barbour et al., 2017;Chen et al., 2017).In response to these events, the United States Geological Survey (USGS) and Oklahoma Geological Survey (OGS) have deployed dense seismic networks to monitor earthquake activity and examine their relation with industry operations, such as saltwater injection during unconventional shale gas production (Keranen et al., 2013;McGarr, 2014;Yeck et al., 2017).These large amounts of continuous seismic recordings provide an important opportunity to investigate environmental changes via seismic ambient noise, and in particular to monitor variations of TWS and droughts in central Oklahoma.Two major bedrock aquifers in central Oklahoma are the Garber-Wellington (GW) in the south and the Vamoosa-Ada (VA) in the east (Figure 1A).Several minor reservoirs, such as the Enid Isolated Terrace (EIT), Cimarron River (Ci), and North Canadian River (NC) also contribute to the hydrological complexity in the study region (Osborn & Hardy, 1999).As an agricultural state, irrigation in Oklahoma highly depends on the supply of groundwater.Detected by the U.S. Drought Monitor (Svoboda et al., 2002), central Oklahoma has suffered from severe droughts in the last decade (Kuwayama et al., 2019).The U.S. is the world's third-largest wheat exporter, and severe droughts in these two major wheat-producing states, Kansas and Oklahoma, led to 7-8% reduction of wheat in 2021 in comparison to their five-year average (Hegarty, 2022).Improving the capability to better monitor water storage in Oklahoma can contribute to sustainable water resources management and steady global food supplys.

Data and Methods
Our study area ranges from 97.90 • W to 96.30 • W in longitude and 35.20 • N to 36.50 • N in latitude (Figure 1A).Nine-year-long (from 2013 to 2022) continuous seismic recordings for all available seismometers are collected as the dataset for this study.Due to the limited number of permanent seismometers deployed in Oklahoma, we also use tempo-rary stations with duration longer than two years to compensate for the data availability.In total, 54 seismometers, mainly from networks NX, OK, GS and O2, are used in Figure 1A.We set 15 and 100 km as the minimum and maximum for the inter-station distance when grouping station pairs and making time-lapse measurements.
Only vertical components of the broadband seismic recordings are used in the analysis.Due to the imbalanced sampling (data only observed at the Earth's surface) and high susceptibility of shallow structures to failure, large shallow velocity changes can be erroneously mapped to smaller variations at greater depths (Juarez & Ben-Zion, 2020).
Here, we estimate the depth sensitivity of our measurements based on fundamental mode Rayleigh waves dispersion with a 1-D velocity profile provided by the OGS (Figure S2 in the Supporting Information).In order to monitor groundwater storage near the Earth's surface, we filter the continuous recordings from 0.1 to 1.0 Hz, allowing us to map timelapse velocity changes down to 1.0 km depth.Our data preprocessing procedures include re-sampling, bandpass-filtering, spectral whitening, instrument response deconvolution, and one-bit normalization (Bensen et al., 2007).The continuous recordings, after preprocessing, are cut into one-hour-long segments.For each station pair, the data segments with the same timestamp are cross-correlated and stacked to construct daily cross-correlation functions (CCFs).Here, we choose a 60-day stacking window with 30-day overlap, which gives us the best performance and stability.The reference CCF for each station pair is the stacked result of all available daily CCFs, while for all daily and refernce CCFs, the positive and negative branches are stacked to boost the signal-to-noise ratio.Once we gather the daily and reference CCFs, their relative time shift, dt, can be computed by using the Moving Window Cross Spectrum method (MWCS) (Clarke et al., 2011), during which the surface and coda waves are isolated by applying dynamic time windows with apparent velocities of 3.0 and 2.0 km/s, respectively (Figure S3 in Supprting Information).With the assumption of homogeneous velocity perturbation, the relative time delay (dt/t) between the daily and reference CCFs should be negatively proportional to the velocity variation (dv/v), i.e., dv/v=-dt/t (Poupinet et al., 1984).More detailed information on the post-processing procedure of ambient noise CCFs can be found in Section S2 of the Supporting Information.
The used dv/v measurements have acceptable uncertainties with respect to frequencies and azimuthal angles.With the same dataset and workflow, four different frequency bands (0.1-1.0 Hz, 0.5-2.0Hz, 0.2-0.8Hz, 0.4-1.6Hz) give us similar dv/v patterns, es-pecially for recovering seasonal variations and long-term trends (Figure S6 in the Supporting Information).The CCFs may also be affected by the heterogeneous distributions of noise sources, such as periodic ocean tides (Ardhuin et al., 2011) and wind speeds (Young et al., 1994).Here, we group the station pairs and re-calculate dv/v according to their azimuthal angles and the similarity among four azimuth groups (0 − 90 • , 90 − 180 • , 180 − 270 • , 270 − 360 • ) suggests that the measured dv/v mainly reflect near-surface seismic velocity variations in the study area, rather than heterogeneous noise source distributions (Figure S7 in the Supporting Information).More details about the uncertainty estimation can be found in Section S3 of the Supporting Information.

Seasonal changes and long-term trends
We compute the mean and median of dv/v among all station pairs in order to illustrate the averaged temporal evolution of near-surface seismic velocity in central Oklahoma (Figure 1B).The raw dv/v time series are smoothed by a Gaussian filter with standard deviation σ=6 days, in order to eliminate unrealistic high-frequency fluctuations.The uncertainty of dv/v is comparatively small, within 5% for most measurements (blue shade in Figure 1B).Some temporary stations were not working, resulting in relatively larger uncertainty from 2016 to 2018.Since 2014, dv/v reduces gradually until the middle of 2017 and then rebounds to its peak value (+0.03%) around the end of 2018 (Figure 1B).After that, it declines further to -0.02% until 2020, and then fluctuates in a relatively low-value zone.Besides this long-term fluctuation, we also observe clear seasonal cycling in dv/v measurements.Statistically stacking the intra-annual pattern of each year, we find that dv/v declines annually to a trough during the summer time (April to June), whereas its peak value commonly appears in the winter season from October to December (Figure 1D).Both seasonal and long-term changes can be clearly identified in the time-frequency analysis of the dv/v measurements by using moving window Fourier transform (Figure 1C).A continuous response, centering around one year period, reflects the seasonal cycling, while another strong response existing at periods greater than two years represents the long-term trend in our measurements.

Comparison with GRACE observations
Figure 2 compares the dv/v results with the GRACE measurements (Syed et al., 2008;Landerer & Swenson, 2012;Richard Peltier et al., 2018), expressed in terms of water equivalent thickness in centimeters.The gap in the GRACE measurements around 2018 comes from an observational gap between the GRACE and GRACE-FO missions, with the latter launched in May 2018.The sampling rates for both GRACE and GRACE-FO missions are one month.To compare signals, both dv/v and GRACE data are interpolated and filtered into the same frequency band.We further quantify the similarity between GRACE and dv/v results by applying a moving window cross-correlation, with a 800-day-long sliding window, between these two time-series (Figure 2B).
Overall, the GRACE data is negatively correlated with our dv/v results, and the absolute cross-correlation coefficient in Figure 2B is greater than 0.7 with almost zero time lag, suggesting that gravity perturbations and seismic velocity variations reflect similar environmental changes in central Oklahoma.Furthermore, we use a least-square regression to determine a linear relation between dv/v and GRACE results (∆h), which gives dv/v = −1.68×10−3 ∆h−9.61×10−3 , with the final data misfit as 0.00848.The confidential ellipse with 2σ (black dashed ellipse in Figure 2C) covers most data points, suggesting a good linear fitting between these two independent datasets.In comparison to the annual statistical stacking of dv/v shown in Figure 1D, we also observe seasonal cycling in the GRACE data, but with an opposite intra-annual pattern (Figure 2D).The strong annual similarity between these two independent time series again suggests that they may reflect similar physical processes, such as variations of TWS near the Earth's surface (Tapley, Bettadpur, Ries, et al., 2004;Rodell et al., 2007;Famiglietti et al., 2011).
As a state-of-the-art tool for monitoring TWS (Rodell et al., 2018), the limited spatial resolution (300 to 400 km) and temporal sampling (one month) of the GRACE measurements cannot satisfy current requirements of ground and surface water storage management at regional to local scales.In contrast, seismic ambient noise measurements have the considerably higher temporal and spatial resolution, depending on the specific deployment configuration, and can be used to monitor the intra-seasonal persistence of TWS deficits and surpluses in a relatively small basin with near real-time fashion and low costs.
Furthermore, since more than 13,000 permanent seismometers have already been deployed around the world and can be openly accessed through the Data Management Center of the Incorporated Research Institutions for Seismology (IRIS-DMC), continuous seismic noise recordings provide a compelling supplement to the GRACE measurements for monitoring water balance on the Earth's surface.

Correlation with groundwater measurements
As a low-latitude inland state, surface water storage and snowfall are negligible in Oklahoma (Swenson et al., 2008), so groundwater storage and soil moisture are the major contributors to the TWS in Oklahoma.Negative correlations between groundwater levels and dv/v have been observed in California (Clements & Denolle, 2018;Mao et al., 2022;Qin et al., 2022) and Germany (Lecocq et al., 2017).Also, the spatiotemporal variations of seismic velocity, after projecting into 2-D maps by using coda wave sensitivity kernels (Mao et al., 2022), are coherent with the discharge and recharge of aquifers in California.The sensitivity of near-surface seismic velocity to groundwater storage can be directly examined by comparing dv/v with groundwater levels.
Deployed by the OGS, the gauge in the Spencer well measures groundwater depths for monitoring the status of the Garber-Wellington (GW) aquifer (Mashburn et al., 2014), which can be obtained from the Oklahoma Water Resource Board (OWRB).Considering lateral variations of seismic velocities, a local dv/v is measured from three seismometers (OK.CHOK, OK.SMO, OK.SWND) surrounding the Spencer well by using the workflow described in Section 2. Starting at 15.5 ft in 2016 (Figure 3A), a monotonic increase of the groundwater depth, indicating a discharge of the GW aquifer, is well correlated with the dynamic increase of seismic velocity around the Spencer well.In contrast, the groundwater level gradually reduces from 17.0 ft in 2019 to 15.0 ft in 2022, representing the recharge of GW aquifer, while the seismic velocity correspondingly decreases by about 0.06 % during the same period (Figure 3A).Similar to the gravitational variations shown in Figure 2A, this negative correlation between groundwater levels and long-term trend of local dv/v further illustrates the sensitivity of near-surface seismic velocities to groundwater recharge/discharge, and demonstrates the capability of using seismic data to monitor TWS with higher temporal and spatial resolution.The seasonal variations around the long-term changes are modeled in terms of thermo-elastic strain (Figure 3B-D) and discussed in Section 4.

Comparison with precipitation and drought index
Next, we compare the dv/v measurements with recordings from precipitation and drought monitoring.Figure 4A shows the history of drought index in central Oklahoma, collected from the USDM (Noel et al., 2020), which classifies the drought condition in five levels, D0 to D4, representing abnormal, moderate, severe, extreme, and exceptional droughts, respectively.From the USDM recordings for central Oklahoma, up to 80% of areas suffered from different levels of drought (D1-D4) from January 2014 to April 2015.
Another two severe droughts appeared from November 2016 to May 2017 and December 2017 to October 2018, when 60% of areas in central Oklahoma were under moderate drought (D1) and up to 20% of areas were exceptionally dry (D4).
The drought index is negatively correlated with the precipitation data (Figure 4B) collected from the Oklahoma Climatological Survey (Boone et al., 2012).For instance, no drought is observed in May 2015 and April 2019 when there were high precipitation volumes (greater than +10 inches precipitation anomaly).In contrast, due to comparatively less precipitation in April 2014 (-10 inches precipitation anomaly), 40% of areas in central Oklahoma were in D4-level drought.Comparably, the dv/v results decrease since February 2014 and reach a minimum in March 2016 when the drought in central Oklahoma is less severe, while they increases during the next severe drought until January 2018 when 90% of areas were in D2-level drought.Because of the frequent and heavy rains in 2019, dv/v fluctuates in a relatively low-value zone where only abnormal droughts were detected in central Oklahoma.It is interesting to note that in central Oklahoma, almost every major drought season coincides with a rapid increase of dv/v.These correlations among the drought index, precipitation, and dv/v again illustrate that nearsurface seismic velocities are influenced by climate changes.

Potential impact of pore pressure on seismic velocity changes
Seismic velocities within the Earth's uppermost crust may be affected by a variety of factors, including environmental changes (Hillers et al., 2015;Lecocq et al., 2017;Clements & Denolle, 2018;Mao et al., 2022), earthquakes (Peng & Ben-Zion, 2006;Brenguier, Campillo, et al., 2008;Bonilla et al., 2019;Qiu et al., 2020), and volcanic activities (Sens-Schönfelder & Wegler, 2006;Brenguier, Shapiro, et al., 2008;Brenguier et al., 2014Brenguier et al., , 2016)).Multiple factors can lead to changes in effective confining pressure or dynamic stresses, and eventually result in seismic velocity perturbations with different magnitudes and characteristic time scales.For instance, tides, temperature changes, snow loading, or sea level changes typically produce 10 −5 to 10 −3 velocity changes with time scales ranging from hours to years (Yamamura et al., 2003;Taira et al., 2018), whereas tectonic/volcanic activities might also lead to 0.1% to more than 10% velocity perturbations (Niu et al., 2008;Brenguier et al., 2014;Pei et al., 2019).Here, it is important to note that the absolute amplitude of derived dv/v depends on the temporal sampling rate used in the analysis (Bonilla et al., 2019).Longer time windows reflect averages over larger time scales and spatial extent, and result in smaller amplitudes than local changes that may be resolved by very short time windows.
Seismic velocity is known to be sensitive to decreasing/increasing effective confining pressure, which affect opening/closure of the microcracks and/or pore space during recharge/discharge of water storage (Birch, 1960;Simmons, 1964;Nur & Simmons, 1969).
As we can monitor in-situ seismic velocity variations, a corresponding sensitivity allows estimating changes of effective confining pressure that is relative to water storage in Oklahoma.Here we attempt to calculate the potential stress sensitivity based on the linear regression between dv/v and GRACE measurements (Figure 2C).Since the GRACE measurement is expressed as equivalent water thickness in centimeters (∆h), the slope value connecting gravity and dv/v (−1.68×10 −3 cm −1 ) can be transfered to the potential stress sensitivity of dv/v by ∆P = ρg∆h.Converting all variables into standard units, the estimated stress sensitivity is −1.72×10 −7 Pa −1 , which is in the same magnitude as the in-situ measurements of previous studies, 10 −7 Pa −1 (Yamamura et al., 2003;Silver et al., 2007).This linear slope also allows estimating, in future studies, the changes of gravity for areas that are too small to be sampled by the GRACE mission.In addition, it provides a possibility to estimate regional stress perturbation from measurements for relative seismic velocity changes.
The clear correlations between our dv/v results with GRACE (Figure 2) and groundwater well measurements (Figure 3A) suggests that when the pore space of sedimentary rocks is filled with water during aquifer recharge, the increasing pore-pressure, and related decreasing effective confining pressure and rock rigidity, generates the reduction of near-surface seismic velocity (Dong & Lu, 2016;Qin et al., 2022).In contrast, the increasing confining pressure and rigidity during drought periods lead to increases in seis-mic velocity.This helps to explain the correspondence between severe droughts with local peaks in the derived dv/v results shown in Figure 4A.

Potential impact of rock density changes on near-surface seismic velocity
Besides pore pressure changes, bulk density changes may also lead to seismic velocity perturbations.It has been concluded that the bulk density ρ, comparing with shear modulus G, may play a more significant role on shear velocity, with respect to water saturation (Nur & Simmons, 1969;M. Wyllie et al., 1962;Mavko & Jizba, 1991).It has recently been validated in a laboratory experiment (Li et al., 2018) that increasing water saturation, from 0 to about 100%, leads to a larger decay in shear velocity than shear modulus solely.We, therefore, discuss the potential impact of density changes on seismic velocity to interpret our results.
To quantitatively investigate the influence of bulk density, we introduce a two-layer conceptual model, with a dry upper layer and a water-saturated lower layer, to represent the aquifer contact.Taking the thicknesses of the dry upper layer and the water table as L 1 and L 2 , respectively, the apparent velocity of this two-layer model can be computed by using the following Voigt-Reuss-Hill approximation (M.R. J. Wyllie et al., 1956;Mavko et al., 2020), where V 1 and V 2 stand for the velocities of the upper and lower layers.We assume a homogeneous and isotropic medium in each layer, V s = G/ρ (Mavko et al., 2020), and a constant shear modulus G that is unchanged by water saturation for simplicity (Nur & Simmons, 1969;M. Wyllie et al., 1962).With the definition of dry and saturated rocks (Gassmann, 1951), the apparent shear velocity of this two-layer model can be re-written as: where ρ dry and ρ sat , controlled by the porosity ϕ, denote the averaged densities of dry and fully saturated rocks, and ρ 0 and ρ w represent the densities of host rock and pure water, respectively.
To represent the circumstance in Oklahoma, we use G = 24.0GP a, ϕ = 0.3, and ρ 0 = 2.65 g/cm 3 (D.C. Wyllie & Mah, 2004) to approximate the sediment, with ρ w = 1.00 g/cm 3 .The total thickness of the two-layer model is L = L 1 +L 2 = 100 m, which is equivalent to the average thickness of the GW aquifer (Mashburn et al., 2014).When the groundwater level L 1 changes from 5 to 6 m representing the discharge of the GW aqufier, the apparent shear velocity V s,app reduces by 2.9 m/s with dv s /v s = −0.03%,which is comparable to the measured dv/v (±0.04%) from seismic ambient noise recordings (Figure 1B).In spite of the simplification of the model, this synthetic velocity perturbation analysis provides a negative correlation between time-lapse seismic velocity changes and groundwater variations through the change of apparent rock density in water reservoirs.
We provide two possible mechanisms, pore-pressure and bulk density, for interpreting the observed negative correlation between groundwater storage and long-term seismic velocity variations.However, at the current stage, we cannot further distinguish these two mechanisms, or assess the non-uniqueness in modeling and inversion, because of limited knowledge of rock properties in Oklahoma, such as averaged porosity around the water reservior, elastic moduli of dry rocks, approximated water saturation in dry and wet seasons, confining pressure around aquifers, etc, which require more detailed investigations in the future.

Thermo-elastic effects on the seasonal changes of near-surface seismic velocities
Annual atmospheric temperature variations have been used to explain seasonal changes of dv/v in many regions, including California (Meier et al., 2010;Hillers et al., 2015;Qiu et al., 2020;Clements & Denolle, 2023), Kyushu Island (Q.Wang et al., 2017), northern Chile (Richter et al., 2014), and Mars (Qin et al., 2023).The nonlinear relation between air temperature and seismic velocity changes can be modeled in terms of thermoelastic strain at the subsurface (Berger, 1975;Ben-Zion & Leary, 1986;Tsai, 2011).This process has been investigated in laboratory experiments by heating and cooling rock sam-ples (Snieder et al., 2002).Here, we use the following experssions of thermo-elastic strain from Ben-Zion and Leary (1986) and Tsai (2011) to explain the relation between air temperature and seismic velocity in central Oklahoma, Here T 0 and ω denote the mean value and frequency of the temperature record, α th and κ are the linear thermo-expansion coefficient and thermo-diffusivity, respectively, λ and µ are two Lamé parameters, and ν stands for the Poisson's ratio.The values of these model parameters for the study region are collected from previous studies (Deming & Borel, 1995;S. Ji et al., 2010;Zhai et al., 2019), and are listed in Figure 3C.The parameter m in Equation 3 represents the second Murnaghan elastic constant, while y b represents the thickness of an unconsolidated cover layer.To simulate velocity variations induced by temperature changes, we determine y b and m by searching for the minimum misfit between observed dv/v and simulated dv/v from Equation 3.This grid-search (Figure 3D) gives values of y b = 2.89 m and m = −1.78× 10 7 GP a, which are reasonable approximations for central Oklahoma.As shown in Figure 3B, the time series of the simulated (green) and measured (blue) dv/v have good correspondence in both amplitudes and phases with a cross-correlation coefficient of 0.843.In addition, the averaged temporal shift between the measured dv/v and temperature fluctuations (red) is around 63 days.The similarity between the simulated and observed seasonal variations of dv/v suggests that air temperature changes are another important environmental factor affecting changes of near-surface seismic velocity.A mismatch between the observed and thermo-elastic simulated dv/v is observed from 2016 to 2019, where the time shift ranges from 54 days in 2016 to 12 days in 2019.Potential causes of these local deviations include changes in the effective layer thickness y b resulting from soil moisture variations, and larger uncertainty due to reduced data coverage associated with fewer temporary seismometers during this period (Section 2).

Figure 1 .
Figure 1.Relative seismic velocity changes measured by ambient noise cross-correlation analysis in central Oklahoma.Available seismometers from different networks are shown as triangles in panel A. Thin white lines denote fault traces mapped at the Earth's surface (Marsh & Holland, 2016), while thick blue lines represent the boundaries of the Garbar-Wellington aquifer (GW), the Vamoosa-Ada aquifer (VA), the North Canadian River (NC), and the Cimarron River (Ci) (https://www.owrb.ok.gov/maps).Panel B shows the mean (blue) and median (red) values of the estimated velocity variations, while the blue shades represent the standard error of the measurements.Panel C illustrates the time frequency analysis of the measured dv/v by using short time Fourier analysis.Panel D shows the statistical stacking of the annual pattern of dv/v.

Figure 2 .
Figure 2. Comparison between relative seismic velocity change (dv/v) with GRACE measurements.Panel A compares dv/v (red) with GRACE measurements (blue) expressed in termes of equivalent water thickness in centimeters.Panel B illustrates the local cross-correlation map between dv/v and GRACE observations.The gap in Panels A and B comes from the survey gap between GRACE and GRACE-FO missions.Panel C shows the relation between dv/v and GRACE results (h) through a linear regression.Yellow bars represent the standard errors of dv/v and GRACE measurements.Panel D shows the statistical stacking of the annual pattern of the GRACE measurements.

Figure 3 .
Figure 3. Analysis of potential causes of long-term and seasonal cycling of measured seismic velocity variations.The black line in panel A represents the historical measurement of groundwater depths for the Spencer well (purple dot in Figure 1A), while the red line is the local dv/v computed from three seismometers in the vicinity of the well (OK.CHOK, OK.SMO, OK.SWND).Panel B compares simulated velocity variations (green dashed line) from a thermoelastic modeling (Ben-Zion & Leary, 1986; Tsai, 2011) with the historical air temperature (red) recorded in central Oklahoma.Panel C gives model parameters used in the thermo-elastic calculation, in which the incompetent layer thickness y b and the Murnaghan constant m are determined by a grid search shown in panel D.

Figure 4 .
Figure 4. Comparison among dv/v, drought index and precipitation records in central Oklahoma.Panel A compares measured dv/v (black) with drought index from the U.S. Drought Monitor, expressed as the area percentage under different drought categories (droughtmonitor.unl.edu).The severity of drought increases from D0 to D4. Panel B presents the precipitation records in central Oklahoma, collected from the Oklahoma Climatological Survey (http://climate.ok.gov).
Taking advantage of recently deployed seismometers in central Oklahoma, we estimate relative seismic velocity variations (dv/v) using continuous ambient seismic noise recordings.The negative correlation between dv/v and GRACE measurements, as well as groundwater levels, can be explained by changes of near-surface seismic velocities due to changes of pore pressure and/or bulk density induced by water saturation.The time delay between the seasonal cycling of dv/v and air temperature recordings can be explained in terms of thermo-elastic strain at the subsurface.The simultaneous occurrences of severe droughts and local peaks of dv/v in central Oklahoma illustrate the sensitivity of near-surface seismic velocities to meteorological droughts.Considering the high temporal sampling rates and flexible spatial resolution of seismic recordings, using analyses of the type performed in this paper can improve the ability to monitor terrestrial water distribution, which is critical for sustainable water-resource management and accurate water-deficit estimation.The historical precipitation recordings for central Oklahoma are accessed from the Oklahoma Climatological Board (https://www.mesonet.org/index.php/weather/monthlyrainfall table).The historical air temperature recordings in Oklahoma are downloaded from Weather Underground (www.wunderground.com/history).
physical Consortium at the University of Texas at Dallas.The raw continuous seismic recordings are collected and pre-processed by using Obspy.We use the MSNoise package to measure temporal evolution of near-surface seismic velocity from ambient noise recordings.The groundwater data are collected from the Oklahoma Geological Survey and the Oklahoma Geological Water Resource (https://www.owrb.ok.gov/maps/PMG/ owrbdata GW.html).Drought Mitigation Center at the University of Nebraska-Lincoln provided the historical drought index for central Oklahoma (https://droughtmonitor.unl.edu/DmData/DataDownload/DSCI.aspx).