Aquifer Monitoring Using Ambient Seismic Noise Recorded With Distributed Acoustic Sensing (DAS) Deployed on Dark Fiber

Groundwater is a critical resource for human activities worldwide, and a vital component of many natural ecosystems. However, the state and dynamics of water‐bearing aquifers remain uncertain, mostly due to the paucity of subsurface data at high spatial and temporal resolution. Here, we show that analysis of infrastructure‐generated ambient seismic noise acquired on distributed acoustic sensing (DAS) arrays has potential as a tool to track variations in seismic velocities (dv/v) caused by groundwater level fluctuations. We analyze 5 months of ambient noise acquired along an unused, 23 km‐long telecommunication fiber‐optic cable in the Sacramento Valley, CA, a so‐called “dark fiber." Three array subsections, ∼6 km apart, are processed and the stretching technique is applied to retrieve daily dv/v beneath each location. Near the Sacramento river, dv/v variations in the order of 2%–3% correlate with precipitation events and fluctuations in river stage of ∼1.5 m. In contrast, regions away (2.5 km) from the river do not experience large dv/v variations. These observations reveal short‐scale spatial variability in aquifer dynamics captured by this approach. Dispersion analysis and surface wave inversion of noise gathers reveal that seismic velocity perturbations occur at depths of 10–30 m. Rock physics modeling confirms that observed dv/v are linked to pore pressure changes at these depths, caused by groundwater table fluctuations. Our results suggest that DAS combined with ambient noise interferometry provides a means of tracking aquifer dynamics at high spatial and temporal resolutions at local to regional scales, relevant for effective groundwater resource management.

Within the last few years, seismic methods have emerged as an attractive tool for monitoring aquifer dynamics. In particular, several studies have demonstrated the feasibility of using continuous records of ambient seismic noise for monitoring small changes in subsurface seismic velocities through time that are related to variations of the groundwater table (Clements & Denolle, 2018;Lecocq et al., 2017;Tsai, 2011). Techniques such as coda wave interferometry (Sens-Schönfelder & Wegler, 2006;Snieder et al., 2002) have become a useful tool in detecting small induced changes; given that coda waves are multiply-scattered, they are assumed to have sampled the medium many more times than ballistic waves, hence being more sensitive to very small changes in velocities. These studies commonly use seismic stations sparsely distributed throughout hydrological basins, with separations of several kilometers, which resolves regional changes along the paths that connect pairs of stations. Ambient seismic noise generated by natural sources is commonly exploited, which is characterized by weak signals at low frequencies (< 4 Hz) that have to be stacked for long periods of time (i.e., days to months) yielding monthly or seasonal changes. A few other studies have demonstrated the potential to use infrastructure-generated noise (frequencies in the order of 5-30 Hz) for local monitoring of natural and artificial changes in groundwater levels with spatial resolution in the order of several tens of meters and temporal sampling of days (Fores et al., 2018;Voisin et al., 2016Voisin et al., , 2017. Other techniques, such as seismic interferometry using ballistic waves as opposed to coda waves (Garambois et al., 2019), and single-station methods (Kim & Lekic, 2019) have recently emerged. One of the main objectives of these approaches is exploiting P-wave energy for a better localization and depth resolution of velocity changes caused by aquifer dynamics.
Whereas all these approaches are successful in providing information on aquifer dynamics at intermediate to large scales, there is clear need for a tool that enables efficient monitoring of groundwater dynamics at short spatial scales of a few meters and high temporal resolution of hours to days, while covering large areas for regional characterization. In this study, we show how distributed acoustic sensing (DAS) arrays deployed on unused telecommunication fiber-optic cables can bridge the gap between local and regional studies by continuously monitoring daily variations in groundwater table levels at scales of a few meters for distances of tens of kilometers. DAS systems probe conventional fiber-optic cables buried in the ground with coherent laser pulses and record the phase of the returning light, backscattered at impurities naturally occurring within the fiber. Changes in the optical phase of consecutive backscattered pulses are measured, which are proportional to changes in longitudinal strain along finite sections of the fiber, referred to as the gauge length (Parker et al., 2014). In this way, DAS technology effectively transforms fiber-optic cables into arrays of thousands of single component seismic sensors by measuring strain-rate variations induced by vibrations impacting the cable. DAS enables recording of high-density (∼1 m) seismological data at sampling frequencies that range from the mHz to the kHz, with apertures of 20-30 km, all with a single cable. These capabilities highlight the advantages of using DAS as a seismic sensor in regional, long-term monitoring experiments where obtaining high-resolution subsurface data is critical, an endeavor that would become logistically complicated and prohibitively expensive with traditional sensors. An increasing number of geophysical studies have already successfully applied DAS in active-source imaging in borehole deployments (Daley et al., 2013(Daley et al., , 2016Mateeva et al., 2013) and for recording earthquakes and ambient seismic noise signals in horizontally deployed cables (Dou et Williams et al., 2019). Recent work has demonstrated that DAS can be successfully deployed on existing subsurface telecommunication fiber-optic cables that are currently not being used for data transmission, both onshore and offshore Jousset et al., 2018;Lindsey et al., 2019;Sladen et al., 2019;Williams et al., 2019). These networks components, referred to as "dark" fiber, are widespread and are often available to be leased or purchased by commercial providers. This potential opens up the ability to record continuous seismic data virtually everywhere fiber-optic cables exist, provided that the cables are buried and coupled to the ground and the DAS interrogator is in range and can be powered. The use of DAS on dark fiber networks significantly reduces the cost, time and effort for deployment. In only a few hours, the system can be set up and we can continuously measure the seismic wavefield at tens of thousands of locations.
In this study, coda wave interferometry techniques are applied to 4 months of ambient noise data acquired on a ∼23 km long DAS dark fiber array to the north of the city of Sacramento, in California's Central Valley . Changes in subsurface seismic velocities (dv/v) are estimated at three locations along the array with the goal of investigating not only temporal but also spatial variability within the basin. By exploiting the high spatial density provided by DAS (2 m sampling, 10 m gauge length) and high frequencies contained in infrastructure noise records (∼4-20 Hz), changes in seismic velocities can be spatially localized beneath each subsection of the array under analysis. Correlation between seismic velocity changes and variations in groundwater table dynamics and related changes in pore pressure are investigated by comparing our seismic observations with precipitation and river stage measurements, and validated by a rock physics modeling exercise. Our results suggest that ambient seismic noise interferometry applied to noise records acquired on DAS arrays, potentially deployed on dark fiber networks, is a promising tool for monitoring groundwater table variations at spatiotemporal scales relevant to achieve sustainable management of groundwater resources.

Study Site
The data analyzed in this study was acquired as part of the Fiber-optic Sacramento Seismic Array (FOSSA) experiment . This experiment took place between July 2017 and March 2018 within the southern portion of the Sacramento River Valley (California), which constitutes the northern branch of the Central Valley ( Figure 1a). The study site is situated between the cities of West Sacramento and Woodland, straddling urban and agricultural areas ( Figure 1b).
Geologically, our study area is located on the floodplain of the Sacramento River. Shallow stratigraphy is composed of recent Quaternary alluvial fan, river, flood basin, and stream channel deposits that constitute a mixture of fine sands, gravels, clays and silts which are typically less than 150 ft (∼45 m) thick (Olmsted & Davis, 1961). These deposits contain the shallow aquifer system within the Valley, which is typically unconfined and highly heterogeneous due to the interfingering of coarse-and fine-grained sediments. This alluvium blankets the Sacramento Valley, covering the deeper aquifer system contained in the Pliocene Tehama Formation, slightly more confined as depth increases. These deposits are the result of alluvial fans that formed at the mouth of creeks and streams as a consequence of the uplift of the Coastal ranges to the west. The fans broadened as they extended eastwards, coalescing with adjacent fans. Lithologically, these deposits are mostly composed of pale green, gray, and tan sandstone and siltstones, with lenses of crossbeded pebble and cobble conglomerates. In our study area, which is located approximately along the axis of the Valley, sediments are typically fine-grained and interfinger with deposits resulting from the Sierra Nevada mountains. Based on correlations with borehole logs from natural gas production wells, in this region of the Valley, the Tehama Formation can be as thick as 2,500 ft (∼760 m; Olmsted & Davis, 1961). Similarly to the Quaternary deposits, permeable zones in the Tehama aquifer are also understood to be both discontinuous and heterogeneous (Water Resources Association of Yolo County, 2005).
Regionally, aquifers are recharged by runoff and groundwater from east-facing foothills, by percolation of precipitation, and by infiltration of surface water provided by creeks and streams. Some recharge is likely also derived from applied irrigation water and from the Sacramento River and the associated artificial canals and bypasses constructed by regional and local agencies. However, as a result of the heterogeneity of the alluvium deposits exposed on the surface, infiltration rates within the Valley are highly variable at short length-scales (i.e., within a few kilometers; Water Resources Association of Yolo County, 2005). Groundwater levels are also spatially variable, depending on whether water is stored in the shallow surficial or deep aquifer. Within our study area, groundwater table levels have remained at depths between ∼10 and 0 m below the surface for the last 4 years (California Department of Water Resources, 2019).

DAS Data Acquisition
The DAS data recorded for the FOSSA experiment was acquired along a section of the US Department of Energy (DOE) Energy Sciences Network's (ESnet) Dark Fiber Testbed. This network consists of more than 20,000 km of short-and long-haul single mode fiber-optic cables designed for telecommunications, connecting DOE experimental sites, supercomputing facilities, and associated research networks. Between July 28, 2017 and March 4, 2018, ∼210 TB of continuous DAS data were acquired along a 23.29 km-long transect of this network running between the cities of West Sacramento and Woodland (Figure 1b). Installation notes from the telecommunication company indicate that the cable was largely deployed on conduit buried in soil at depths of 1-1.5 m. Some sections were also placed in shallow horizontal boreholes beneath roads and railway tracks, again in conduit but slightly deeper (3-4 m). Data sets were recorded at a sampling frequency of 500 Hz with a spatial sampling of 2 m, resulting in strain-rate measurements at over 11,000 receiver points. Data was acquired using a commercial DAS unit (Silixa iDAS v.

Data Selection
For this study, ambient seismic noise recorded between October 25, 2017 and March 4, 2018 was utilized. Acquisition parameters of the DAS unit were adjusted on October 25 to improve data quality at long distances. Inspection of the data revealed that waveforms acquired before and after the settings modification were slightly different. These differences would have introduced artificial time lags between waveforms in our time-lapse analysis. Thus, we chose to discard data acquired before October 25 in this study.
With the objective of exploring not only temporal, but also spatial variations in aquifer properties within the study area, three sections along the array were chosen for detailed analysis (subsections A-C, Figures 1b-1e). These profiles are 300 m long, and are located in different settings within the area covered by the array (urban vs. rural). We should note that the distance between the fiber-optic cable and the main channel of the Sacramento River is different for the three sections; while subsections A and B are less than 200 m away from the river, subsection C is at a distance of 2,500 m (Figures 1c-1e).

Ambient Noise Cross-Correlation of DAS Subarrays
Each 1-minute-long DAS ambient noise recording is sequentially processed. The three, 300 m-long array sections of interest are extracted from each file and, subsequently, data for each section is analyzed following well-established ambient noise interferometry approaches Bensen et al., 2007). Data are detrended, demeaned and downsampled to a sampling frequency of 125 Hz (8 ms) after applying an anti-aliasing low-pass filter. Temporal normalization using a running absolute mean with a window of 0.5 s is applied to reduce the effect of earthquakes and other high-amplitude, undesired signals. Next, data are bandpassed between 0.002 and 15 Hz, followed by spectral whitening to balance frequencies in this band. Following these pre-processing steps, cross-correlations are calculated between the southernmost trace of each section, which acts as the virtual source, and all other traces in that section. The results are 1-minute-long virtual shot-gathers for each analyzed section for the duration of the experiment. Due to the effect of cross-correlating coherent noise generated by the DAS instrument, which is present in each trace, the resultant virtual shot-gathers are contaminated by zero-moveout noise that is particularly severe at zero lag-time. This noise is reduced by subtracting the median of all traces for each time sample (Rodríguez Tribaldos et al., 2019). Lastly, all 1 minute-long virtual-shot gathers for each day of the recording period analyzed are stacked using the phase-weighted stacking method of Schimmel and Paulssen (1997) with an exponent of 0.3. For stacking, the causal (positive) and acausal (negative) parts of the cross-correlations are averaged to increase signal-to-noise ratio. As a result of this process, we obtain a virtual shot-gather for each day between October 25, 2017 and March 4, 2018 for each selected section of the DAS array.
The causal (positive) part of a daily cross-correlation stack at each array section is shown in Figures 3a, 3c, and 3e for illustration. At all three sections, the stacked virtual-shot gathers reveal a strong train of surface waves traveling across the sub-arrays. In most stacks, the first 0.5 s are dominated by waves propagating at velocities in the order of 300 m/s. Because of the single component nature of DAS measurements, these arrivals are interpreted as ballistic Rayleigh waves propagating along the direction of the fiber-optic cable (Martin et al., 2018). Secondary arrivals with slightly different velocities immediately follow the main  Rayleigh wave train. These arrivals are interpreted as a mixture of multiply-scattered waves (i.e. coda waves) and higher-order surface wave modes.
In order to investigate temporal coherency of the wave arrivals identified in the stacks, a trace at a specific distance along the sub-array is extracted for each subsection. A large enough distance is needed between virtual source and virtual receiver to allow for the different phases to separate and be distinguishable. A source-receiver separation of 200 m is chosen. Traces are inspected visually and those stacks resulting in noisy traces are removed and not considered for further analysis. Figures 3b, 3d, and 3f show this trace for all daily stacks. In all sections, coherent phase arrivals are present in each daily stack, which confirms the repeatability of the stacked cross-correlated functions in the offset-time domain throughout the experiment.
This repeatability also reveals significant shifts in the arrival times of temporally coherent phases at specific times during the experiment. The most evident time delays are observed for array sections A and B in early January and early March 2018 (Figures 3b and 3d). A hint of these delays is observed for the main Rayleigh wave arrival (i.e. prominent phase before 1 s lag time), but largest time shifts affect coda waves, with the delay slightly increasing with lag time. Consequently, we focus our seismic interferometry analysis on these later wave arrivals (Figure 3). In an attempt to better constrain the depth of the changing medium, this analysis was repeated for cross-correlated data filtered to different frequency bands comprised between 0.002 and 15 Hz. Data filtered between 4 and 15 Hz was chosen, as it reduces the amount of low-frequency arrivals that might not propagate through the array, and shows the larger time lags between consecutive daily stacks.

Coda Wave Interferometry Application
In order to quantify the observed time-shifts, coda wave interferometry techniques are applied to the 200 m offset cross-correlations from each daily stack, filtered between 4 and 15 Hz. We focus our analysis on wave arrivals comprised within a window beginning at 0.8 s and ending at 1.3 s lag time, immediately after the main Rayleigh wave. The stretching technique introduced by Sens-Schönfelder and Wegler (2006) is applied to consecutive pairs of traces throughout the recording period, with the objective of estimating variations in relative seismic velocities (dv/v) between consecutive days. This method assumes that one trace is a stretched version of that same trace for the previous day according to the linear function where τ is time lag, t is time and ϵ is the stretching factor. Assuming that velocity perturbations are homogeneous, ϵ = −dv/v. We choose to measure −dv/v between consecutive traces, rather than comparing the cross-correlation function of the day of interest against the stack of all traces, as it is typically done in coda wave interferometry studies (e.g., Clements & Denolle, 2018;Sens-Schönfelder & Wegler, 2006). The main reason for using this approach is that the relatively large time delay observed between consecutive days at some times throughout the recording period made it impossible to constructively stack all daily cross-correlations without loosing coherency. By stacking all traces, some of the coherent phases disappeared and the reference trace was too different to the daily cross-correlations to obtain any meaningful results using the stretching technique, which relies on the coherency between the two traces that are being compared. Besides, this approach also enables us to evaluate the potential of using DAS for capturing short-term changes that could be masked by stacking all traces, and to apply this technique for real-time monitoring during which changes could be analyzed on a daily basis.
In an iterative approach, stretching factors ranging from −10% to +10% are applied to the later trace of a pair in a moving window with a width of 0.25 s and a step size of 0.02 s. After each stretching operation, the two traces are cross-correlated in these windows. The stretching factor yielding the maximum zero-lag cross-correlation coefficient between the two traces for each analysis window is taken as the optimal dv/v. In our analysis, all estimated dv/v yield cross-correlation coefficients >0.8. Next, dv/v estimates for all windows are considered, and values above and below the 10% quartile are discarded to remove outliers. The median of all retained values is then calculated to yield a single value of dv/v for that day. We chose to use the median, as opposed to the mean, as the representative daily dv/v value because the median is less sensitive to outliers. Uncertainty in dv/v is calculated as the interquartile range of all kept values, as the variability in RODRÍGUEZ TRIBALDOS AND AJO-FRANKLIN 10.1029/2020JB021004 dv/v within each day does not show a perfectly normal distribution, but rather skewed to positive or negative values depending on the day of the recording. This procedure is repeated for each consecutive pair of traces throughout the recording period, resulting in a profile of daily variations in dv/v. Finally, a cumulative median of daily dv/v is calculated for direct comparison of our seismic observations with hydrological data.

dv/v Variations Through Time and Space and Comparison with Hydrological Data
Results from applying the described coda wave interferometry approach to stacks of the three sections under analysis are shown in Figure 4. Subsections A and B reveal a very similar pattern of seismic velocity variations (Figures 4b and 4c). Two significant negative perturbations, on the order of −2% to −3%, are seen in both sections in early January and early March, agreeing with the larger time lags observed in Figures 3b and 3d. In subsection C, however, no significant seismic velocity variations are observed throughout the recording period. These observations suggest spatial variability in seismic velocity perturbations within the study area, which is most likely reflecting variability in aquifer properties, hydraulic boundary conditions, and infiltration rates.
In order to explore this hypothesis, hydrological data contemporaneous to our seismic experiment is investigated for locations near our DAS array. Unfortunately, direct observations of groundwater table elevations at the FOSSA experiment site are sparse in both space and time. Only one groundwater well nearby the fiber-optic cable has been sampled during our recording period (well Sac Bypass, Figure 1b), with only two data points available from December 7, 2017 and February 14, 2018. Despite the lack of direct subsurface measurements, alternative hydrological data connected to groundwater are available for comparison. The top panel of Figure 4 shows measurements of daily precipitation and river stage contemporaneous with our seismic experiment, recorded at environmental stations SMF and BYL, respectively (Figure 1; California Department of Water Resources, 2019). A close comparison reveals that, for those sections of the array located only a few hundred meters away from the river (subsections A and B), both of these data sets are anti-correlated with seismic velocity variations, that is large precipitation events immediately followed by an increase in river stage correspond with a considerable decrease in seismic velocities below these sections. A correlation also exists between the magnitude of change for all three data sets, with a larger dv/v for a larger increase in river stage and larger precipitation event. The temporal correspondence between river stage and changes in seismic velocity perturbations suggest that variations in the river level are directly connected to changes in the subsurface immediately beneath the array. Figure 5 shows contemporaneous measurements of river stage at station BYL and groundwater table elevation at station Sac Bypass, located ∼1,200 m away from subsection A, for a period between January 2016 and January 2019. These observations show that perturbations in groundwater table and river level have agreed in both timing and magnitude for at least the last three years. Thus, we conclude that it is reasonable to use changes in river stage measurements as a proxy for variations in groundwater table elevation for locations near the river (Ha et al., 2008). Previous hydrogeophysical studies have also verified this approach in a near-river aquifer monitoring context (e.g., Baker et al., 2000). Consequently, we interpret seismic velocity perturbations observed below subsections A and B as an indication of fluctuations in the groundwater table at those locations. Moreover, we interpret these changes as variations in shear-wave velocity, as surface-wave coda waves are mostly sensitive to changes in V s .
In subsection C, precipitation and changes in river stage level do not bear any correlation with dv/v, which does not show significant variation through time. These observations are interpreted as either a lack of or a minimal variation in groundwater table levels beneath this section of the array. As previously pointed out, one of the main differences between subsection C and the other two sections of the array is its larger distance to the Sacramento river (2,500 m for subsection C, Figure 1). As part of the development of the Yolo County Integrated Regional Water Management Plan, relative infiltration rates of Quaternary sediments across the county were mapped (Water Resources Association of Yolo County, 2005). According to this study, our fiber-optic array is located in a region of slow to very slow infiltration rates. This observation suggests that recharge of the shallow aquifer from surface water is likely to take longer than just a few days, even for large precipitation events. Hence the increase in groundwater table beneath subsections A and B is likely not a result of infiltration from the surface, but rather a consequence of the increase of the level of the river, which is hydraulically connected with the floodplain. Due to the heterogeneity of the quaternary deposits and consequent discontinuous permeability, aquifer connectivity over thousands of m is difficult, and variations in the level of the Sacramento river in the order of ≤ 2 m do not affect the elevation of the groundwater table beneath subsection C. These observations suggest a strong connectivity between groundwater and the river bed in this basin. Understanding this interaction is extremely important, as significant variations in the groundwater table due to, for example, pumping can have important effects in the river ecosystem.
At even shorter spatial scales, a close look at dv/v variations for subsections A and B reveals a slight difference in the behavior of these velocity perturbations. The magnitude of cumulative dv/v decrease following RODRÍGUEZ TRIBALDOS AND AJO-FRANKLIN 10.1029/2020JB021004 10 of 20 major river stage increase events is in the order of 2.5%-3% in early January and 2% in early March beneath both sections. However, the pattern of recovery of these seismic velocities to the levels shown prior to river stage increase events differs between the two sections. In subsection A, the short term recovery of dv/v agrees with that of river stage levels after the event in early January. River level slowly decreases over 10 days, to stabilize at a level of 2.7 m, ∼0.5 m above the stage previous to the increase event. Seismic velocity perturbations follow a similar behavior, slowly decreasing toward 0% and stabilizing at a value of ∼−1% after 10 days. In the long term, however, river stage levels recover to its mean value of 2 m toward the end of the experiment before the second large increase occurs in early March, whereas seismic velocities remain fairly constant until they decrease again following the river stage increase. In contrast, for subsection B, no short-term recovery is observed in seismic velocities, which remain at a value of −3% after the large perturbation in January to decrease again following the March event. These observations agree with the geological nature of the Quaternary deposits making up the floodplain of the Sacramento river. First, the lack of recovery of the seismic velocities can be interpreted as due to the presence of clay horizons within these deposits. A geotechnical well shown in Ajo-Franklin et al. (2019), located right next to the river bed between subsections A and B, shows layers of clays and silty clays interbedded with fine sands and gravels. Clay horizons would retain part of the water, preventing the full recovery of average seismic velocities. In relation to this effect, the fact that subsection B shows no short-term recovery of dv/v might indicate higher proportions of clays beneath this region. It could also indicate strong hysteresis during drainage/imbibition cycles controlled by the same textural differences.

Depth Localization of Seismic Velocity Changes
Following quantification of seismic velocity variations, we investigate the depth range at which these changes occur. For that purpose, the depth sensitivity of the coda waves used in our analysis is explored by calculating V s sensitivity kernels as a function of depth and frequency. In order to calculate these kernels, information on the frequency content of our data and the shear-wave velocity structure beneath the study area is required. Dispersion analysis followed by surface wave inversion are applied to one of the daily stacks retrieved beneath subsection A. We choose to analyze this section because we considered it to be better constrained, as it is closest to the river stage measurement station (Figure 1). Slant-stacking is applied to the virtual common-shot gather stack from January 8, 2018 to transform the data from the time-offset domain to the frequency-velocity (dispersion) domain. The resultant dispersion spectra, shown in Figure 6a, is characterized by high amplitudes at 4-5 Hz and around 6 Hz. We also calculate the average frequency content of the stack, which also reveals distinct peaks around 4.5 Hz and at 6 Hz ( Figure 6c).
The dispersion curve corresponding to the fundamental mode of the Rayleigh wave is manually extracted from the dispersion spectra and subsequently used as input for the Monte Carlo surface wave inversion algorithm of Maraschini and Foti (2010) Using the retrieved frequency content characteristics and one-dimensional V s structure, V s sensitivity kernels are calculated using the computer programs in seismology (CPS) package (Herrmann & Ammon, 2002). We calculate the sensitivity kernel for the fundamental mode of the Rayleigh wave. Through a series of comprehensive numerical studies, Obermann et al. (2013) show that, at early times, coda waves are dominated by surface waves and are more sensitive to shallow depths. The study shows that, for velocity variations occurring at shallow depths, depth sensitivity is governed by the 1D sensitivity of the fundamental mode of the surface waves. Besides, as mentioned earlier, the axial sensitivity of DAS suggests that the waves recorded by the fiber-optic cable correspond to Rayleigh waves. Thus, the sensitivity kernel calculated here is representative of the coda waves that we have analyzed, and is a good indication of the most likely depth range of occurrence of our dv/v observations. V s kernels as a function of depth are shown for the main frequency components of our stack in Figure 6d. These kernels suggest that the analyzed waves are most sensitive to velocity structure occurring within the top 30 m of the subsurface beneath the array. Peak sensitivities are at 25 m for 4.5 Hz and 18 m for 6 Hz. Following our previous observations, these depths are far below the groundwater table interface, which is estimated to be located at depths of 4-5 m. Thus, we conclude that the dv/v variations observed using coda wave interferometry do not correspond to changes in the location of the groundwater interface itself, but rather to processes that are related to changes in the water level but occur at larger depths.

Beyond Correlation: Physical Links Between Aquifer State and Seismic Properties
A variety of hydrogeologic processes generate seismically observable signatures including direct fluid saturation effects, clay hydration, and pore-pressure perturbations. Before providing a rock physics modeling framework to quantitatively link surficial water table variations to dv/v, we will attempt to rule out several of these processes.
The first of these processes, replacement of air by water near the water table and the capillary fringe, can generate large seismic velocity perturbations for P-waves, and smaller S-wave effects, due to density changes as has been demonstrated in both the laboratory (Knight & Nolen-Hoeksema, 1990) and the field . However, as can be seen in Figure 6, our surface wave observations are dominated by ambient noise in the 4-6 Hz range and are predominantly sensitive to property perturbations between 10 m and 30 m below ground surface, far deeper than the 4.5 m mean groundwater level. Likewise, clay hydration and cohesion controlled phenomenon are mainly relevant to the vadose zone and the capillary fringe, near surface effects that are likely difficult to resolve using our 4-6 Hz measurement band.
The impact of pore pressure perturbations seems a more likely causative mechanism for the observed dv/v changes; an increase in water table height should increase pore pressure in zones with direct fluid communication and decrease effective stress, thus decreasing shear wave velocity. This effect is large at shallow depths and decreases as depth increases, as lithostatic pressure prevails over hydrostatic pressure. Over the relatively shallow depths considered (40 m), pore pressure diffusion times should be small in comparison to the ambient noise integration times (minutes vs. days). Past field hydrogeophysical studies of aquifer state using conventional broadband seismic observations have measured similar effects, mainly a decrease in dv/v correlated with higher water tables or seasonal precipitation cycles (Clements & Denolle, 2018;Lecocq et al., 2017;Tsai, 2011;Voisin et al., 2016).
These past studies have generally assumed a uniform linear stress dependence for V s which seems unlikely, particularly at shallow depths and low effective stress states where a nonlinear effective stress dependence has been firmly established both in laboratory (Prasad, 2002;Zimmer et al., 2007) and field studies . We utilize an established granular media model developed by Walton (1987) to estimate the dependence of V s on effective stress state in the near-surface. The Walton "Smooth" model assumes an uncemented granular media and has been extensively utilized in past studies to provide V s estimates as a function of stress state. We refer the interested reader to Mavko et al. (2020) for the full formulation. Given that the contact with consolidated rocks is well below our zone of investigation, as evidenced by co-linear drilling logs which did not contact such materials , the adoption of a granular media model seems reasonable.
Due to the lack of near-surface property logs, we assumed that soil porosity was a constant value of 38%, typical for unconsolidated sediments, and the coordination number was estimated from the empirical relationship described in García and Medina (2006). Pore pressure was assumed to be hydrostatic and dependent only on depth to water table (hydraulic head) while lithostatic stress was assumed to have a 23 kPa/m gradient. Grain elastic properties were assumed to be those of quartz (Mavko et al., 2020); a grain bulk modulus of 36 GPa, a grain shear modulus of 45 GPa, and a grain density of 2,650 kg/m 3 . To calculate dv/v estimates using this rock physics model, we use the baseline V s estimate provided previously in Figure 6 and used ΔV s estimates based on the model. One small inconsistency is that the Walton model was only used for the ΔV s estimate and does not perfectly fit the model shown in Figure 6b.
With the availability of a rock properties model relating pore pressure to dv/v, we first converted the river stage proxy data discussed previously to water table height. This height was used to calculate pore pressure perturbations, and the parameters listed above were then utilized to estimate dv/v as a function of depth for all times at which measured dv/v values were available. Figure 7 shows the measured river stage values from station BYL in panel (a). Figure 7b compares the measured dv/v from the DAS ambient noise processing to estimates of dv/v from the rock properties model for depths of 18 m (blue) and 25 m (red). These depths correspond to the 6 and 4.5 Hz peaks in surface wave sensitivity shown in Figure 6. As can be seen, the modeled dv/v values derived from water table variations capture the reduction in velocity seen after the major precipitation event which occurred in January 2018 in both sign and magnitude, as well as the late February 2018 event. This suggests that, for reasonable water table variations, the associated pore pressure perturbations should be more detectable using DAS coda wave analysis and surface waves. However, there are several features in the estimate which are not replicated in the coda wave measurements. These features include the November 2017 rainfall event, which may be below our sensitivity levels, and the seasonal river stage reduction seen from late January 2018 to early February. Likewise, the small diurnal variations in river stage are aliased because of the ambient noise stack period. Besides, they may not be large enough to drive water table variations at distance from the river. However, given these caveats, the accuracy in an absolute sense of the measurements is surprising, suggesting that a direct inversion for water table depth or hydrostatic head in confined aquifers is a possibility.

Discussion
Our combined analysis of seismic velocity variations and rock physics modeling indicates that dv/v changes as large as 3% are observed almost immediately after short-term river stage increase episodes, in connection with increasing pore pressures due to an elevation of the groundwater table in the order of ∼1.5 m. These observations are comparable to values reported by a limited number of previous studies monitoring shortterm, local groundwater table fluctuations using ambient seismic noise in the infrastructure frequency band (i.e., 5-30 Hz; Voisin et al., 2016Voisin et al., , 2017. We should note that these dv/v perturbations are considerably larger than those discussed in other hydrogeophysical studies of this type, which mostly focus on monitoring seasonal variations at a regional scale. For example, Sens-Schönfelder and Wegler (2006) observe dv/v variations of up to 1% in a period of 6 months using cross-correlations calculated between stations separated ∼170 m. They attribute this observed dv/v to a groundwater level decrease of ∼25 m inferred from a hydrological model based on precipitation rates. Clements and Denolle (2018) report smaller changes in dv/v on the order of 0.1%-0.15% over several months as a result of seasonal groundwater table perturbations of about 20 m. In a long-term monitoring study focusing on changes occurring over several decades, Lecocq et al. (2017) observe dv/v variations of only ±0.01%. The explanation for this disparity in the magnitudes of the observed seismic velocity variations mostly resides in the frequency content of the ambient seismic noise used for analysis. Sens-Schönfelder and Wegler (2006) high-pass filter their one-day seismic records at 0.5 Hz. Clements and Denolle (2018) analyze ambient seismic noise in the 0.5-2.0 Hz band, whereas Lecocq et al. (2017) explore microseismic noise with frequencies <1 Hz. On the vertical scale, surface waves are sensitive to the average effect of elastic parameters between the surface and a maximum depth governed by their wavelength. Hence, the lower the frequency, the larger the wavelength and the averaging effect. As a result, small changes in subsurface properties will have a very small effect in the average velocities recorded by surface-wave dominated coda waves. Moreover, these investigations estimate seismic velocity changes occurring along the paths between seismic stations that are usually located several kilometers apart. This configuration implies that seismic waves sample a larger volume of the subsurface, which will have a similarly small effect in seismic velocities (Obermann et al., 2013). Lastly, contact theory models similar to the Walton model used in our prior analysis, predict a considerably higher sensitivity to pore pressure variations at shallow depths (at low effective stress states). This would imply a lower predicted hydrogeophysical response at lower seismic frequencies.
In the case of the study by Sens-Schönfelder and Wegler (2006), the combination of shorter station spacing and higher frequency content (>0.5 Hz) with respect to those of Clements and Denolle (2018) is the most likely cause for the larger magnitudes of the observed dv/v values for a similar inferred groundwater table change. Here, we analyze infrastructure generated noise in the 4-15 Hz band, and station separation of 200 m. This enables the detection of smaller perturbations in groundwater table changes which have a larger effect in our recorded seismic velocities. Our observations suggest that, in this type of geological environment, groundwater table variations on the order of 1 m or more are resolvable using coda wave analysis of infrastructure noise recorded using DAS arrays.
An additional concern relating to coda wave interferometry monitoring studies is the ability to recover changes in subsurface properties without the overprint of temporal variations in the properties of noise sources. Prior studies have shown that apparent seismic velocity perturbations as large as 0.05% can be caused by temporal variability in the frequency content of the seismic ambient noise (Zhan et al., 2013). The large magnitude of our recovered dv/v variations and the agreement with auxiliary data and rock physics modeling results suggest that our observations reflect real changes in subsurface properties. Nevertheless, we confirm the reliability of our results by analyzing the stability of the frequency content of our cross-correlated waveforms. As illustrated in Figure 8, day-to-day variability in frequency content is very small for the stacks used for coda wave interferometry, which confirms that our observed velocity variations are due to real changes in subsurface properties, and not spurious measurements introduced by changes in frequency content.
Our study implies that the analysis of seismic velocities in combination with rock physics modeling is a promising approach to infer changes in aquifer state. However, some disagreement remains between estimated dv/v variations and those recorded by coda wave analysis. The most obvious feature is the predicted seasonal recovery occurring as a response of lowering the river level between January and March (Figure 7), which is not represented in our observations. During this time, a decrease in river stage of more than 1 m occurs, which is predicted to cause an increase of 1.5% in dv/v. These fluctuations are in the same range as those observed after the early January increase event, which are visibly captured by our seismic analysis. This observation suggests that the reason for this mismatch is not a lack of sensitivity of the coda wave observations. One possible explanation for this disagreement would be a limited ability to recover long-term dv/v variations due to comparing consecutive traces rather than using the stack of all daily traces as the reference trace. Alternatively, it is worth noting here that, due to a lack of ground-truth geological information at shallow depths at our analysis location, our rock physics model assumes the elastic properties of our model domain to be those of pure quartz, hence implying that sedimentary deposits are composed of 100% sand. However, lithological logs for nearby areas confirm the presence of clay layers, as expected for flood plain deposits. As previously discussed, the retention of some water by these clay layers could delay the water table response and prolong the reduction in seismic velocities visible in our observations. In future studies, a more detailed knowledge of the lithological characteristics of the subsurface and the use of additional ground truth water-level observations can help further investigate the relationship between seismic properties and changes in pore pressure distribution, as well as other physical processes linked to groundwater fluctuations.
Lastly, we should mention that even in cases where monitoring hydrogeophysical state is not useful for direct groundwater management, recovered V s variations can play an important role in providing corrections for deeper monitoring targets. An example would be geothermal systems (e.g., Zeng et al., 2017) or geologic carbon storage (GCS) reservoir monitoring, where perturbations in deeper units might be cloaked by variability in surficial aquifer state. The capacity to correct for near-surface variations provides a path toward recovering true velocity variations at depth. Moreover, the broadband nature of DAS recordings and the ability to continuously acquire data at apertures of tens of kilometers opens up the possibility to apply the interferometric analysis to recordings of a wide variety of frequency and wave types (i.e., body waves as opposed to surface waves) in order to localize seismic velocity variations at a range of depths using a single data set.
In conclusion, our results suggest that seismic interferometry analysis of ambient noise acquired using DAS arrays has significant potential as a tool for monitoring groundwater fluctuations at local to regional scales. DAS technology provides high density of measurements and large aperture, which enable continuous monitoring of spatial variations at high resolution and allows for improved localization of velocity anomalies that can be related to changes in aquifer properties. The possibility of deploying DAS on existing, unused telecommunication fiber-optic cables offers new opportunities for monitoring subsurface processes at minimal effort in complex environments where the deployment of traditional monitoring equipment, such as water-wells, or conventional seismometers could prove challenging and costly. Here, our approach has been applied to three sub-sections of the array that cover different regions within this portion of the Sacramento basin. However, our observations indicate that the DAS array can record ambient noise along 20 km with high signal-to-noise levels, sampling the wavefield at distances as small as 1-2 m. Thus, one could apply this approach to consecutive, overlapping array sub-sections along the entire profile to provide a continuous record of seismic velocity variation and, hence, aquifer property variations along several tens of kilometers with unprecedented spatial and temporal resolution. In this context, monitoring of seismic velocity variations using DAS can provide valuable information on aquifer dynamics at a scale relevant for water agencies to design and implement plans for sustainable management of groundwater resources.
Despite these undeniable advantages, it is worth mentioning that the routine application of dark fiber DAS for long-term monitoring experiments still presents some challenges. The most significant of these barriers is the large data volume generated by these arrays. In the present study, over 11,000 measuring channels at 2 m spacing with a sampling frequency of 500 Hz generated ∼210 TB of raw data in 7 months of recordings.
These unusually large data sizes require innovative approaches in data storage, handling and processing in order to fully exploit the potential of these exceptional data sets. Recent developments in fields such as high performance computing, array processing, and artificial intelligence applied to geophysical data (Bianco et al., 2019;Clements & Denolle, 2020;Dong et al., 2020;Hu et al., 2020;Zhong et al., 2020) promise to alleviate these issues and contribute toward advancing DAS as a revolutionizing tool in subsurface monitoring investigations.

Data Availability Statement
MATLAB was used for data analysis and plotting. Due to the large size of the data set (hundreds of TB), processed data products are made available. Examples of raw DAS data recorded along the array subsections analyzed in the study (shown in Figure 2) and pre-processed data shown in Figure 3 and used to calculate dv/v are available at Rodríguez Tribaldos and Ajo-Franklin (2021). River stage and precipitation data from stations BYL and SMF, respectively, can be accessed at https://cdec.water.ca.gov/. Groundwater level data from station Sac Bypass Shallow can be accessed at: https://water.ca.gov/Programs/Groundwater-Management/ Groundwater-Elevation-Monitoring--CASGEM.