Introduction

One of the worst droughts on record is ravaging the southwestern United States and regions elsewhere in the world. The severity and prolonged duration of this drought bring the water security to a tipping point. To address the water scarcity crisis requires efficient and sustainable management of freshwater resources. A critical component of freshwater is groundwater, which contributes over 60% of the total water supply in dry years in California. The state-of-health of regional groundwater systems has never been more important for agricultural, ecological, industrial, and urban purposes. However, owing to the remoteness and multi-scale, highly heterogeneous structures of water-bearing formations hidden beneath the Earth’s surface, the spatial distribution of aquifers remains poorly understood. Equally limited is our knowledge of the temporal variability of groundwater storage on monthly to decadal timescales. The spatial and temporal constraints of aquifers are not only crucial for the efficient use of groundwater resources but also indispensable for avoiding permanent loss of storage capacity (which may occur when the pore pressure drops below pre-consolidation levels)1,2,3. Observations that characterize effectively the structures and dynamic behaviors of aquifers are thus imperative to sustain groundwater resources.

The traditional observation for tracking groundwater tables is the hydraulic head, which directly measures the water level of specific aquifer layers. In practice, however, the point-scale head measurements are too sparse in time or space to adequately capture variations of the highly nonuniform aquifer systems. Over the past decades, great progress has been made in hydrological investigations by utilizing a variety of non-invasive geophysical methods. Surveys exploiting the electromagnetic properties of the near-surface medium4,5,6,7,8 allow to depict the groundwater systems with impressive detail, but they are expensive to conduct and often limited to static imaging on local scales. Remote sensing techniques have emerged as powerful monitoring tools to inform about the space-time variability of terrestrial water storage. Yet they still have notable limitations. For example, satellite gravimetry1,9 has a low spatiotemporal resolution (hundred-kilometer scale at monthly interval), and GPS and InSAR (Interferometric Aperture Radar) data resolve surface deformation but do not provide direct constraints on structures and processes beneath the surface10,11,12,13,14,15,16,17,18,19.

Seismological observations can serve to add new constraints to the 4D (space-time) fluctuations of groundwater systems. Seismic velocity (v), or, rather, the propagation speed of seismic waves, is determined by elastic moduli and bulk density and can serve as an in-situ measure of the mechanical properties of the underground medium. In particular, relative changes in seismic velocity (Δv/v), which associate with changes in pore pressure, saturation, porosity, and micro-structures20,21,22,23, can be used to infer the volume of stored groundwater. Modeling of seismic wavespeeds, taking into account the poroelastic effect and elastic loading24, has suggested an anti-correlation between Δv/v and groundwater level.

It has been demonstrated in recent years that Δv/v can be measured continuously in time, with high sensitivity and at low cost, by combining two seismic interferometry techniques: (1) interferometry of seismic ambient noise fields, which provides repetitive estimates of Green’s functions of the Earth interior structure25,26,27,28, and (2) interferometry of seismic coda waves, which allows measurements of tiny changes in wavespeed29,30. This protocol has been employed to analyze time-evolving processes in various geological settings31,32,33,34,35,36,37,38, including proof-of-concept studies for monitoring groundwater levels31,39,40,41,42.

In this study, we present a 4D monitoring approach to push such seismological monitoring of groundwater (and other shallow sub-surface) systems to a new level, by not only inferring how Δv/v changes in time (as in previous cases) but also in 3D space43. This advance is made possible by employing recently developed analytical descriptions of seismic coda propagation44 in seismic interferometry. We apply our approach to the drought-stricken, densely populated Los Angeles (LA) coastal region in Southern California (Fig. 1) and measure the space-time evolution of Δv/v during 2000–2020. Our results reveal the spatial patterns of seasonal and long-term fluctuations in groundwater volume at several hundred meters below the metropolitan LA43.

Fig. 1: Hydrogeological settings and station locations in the study area.
figure 1

a Groundwater basins near Los Angeles, California (adjudicated to different local water departments). Primary basins studied in this paper include the San Gabriel Basin, the Los Angeles (LA) Central Basin, the Santa Ana Basin, and the LA West Coast Basin. b The unconfined (forebay) and confined zones in LA Central and Santa Ana basins. The forebay boundary is estimated by the Orange County Water District (OCWD)45 and the Water Replenishment District (WRD) of Southern California46. Confined to semi-confined conditions also exist in the San Gabriel Basin. c Station locations. Yellow triangles denote the broadband, 3-component seismic stations in the CI network operated by the Southern California Earthquake Data Center (SCEDC) available from 2000 to 2020. The blue circle denotes the groundwater monitoring well GGM-1/MP1 operated by OCWD, from which the hydraulic head data are shown in Fig. 2b. The violet diamond denotes the meteorological station ANA at which the precipitation record is shown in Fig. 2a.

Full size image

Results

Seismic interferometry in Coastal Los Angeles Basins

We use records of seismic ambient noise from about 50 broadband seismic stations (Fig. 1c) in the Coastal Los Angeles Basins (CLAB) to calculate Δv/v. CLAB comprises a number of groundwater basins (Fig. 1a). The shallow subsurface is composed of unconsolidated Quaternary-age sediments mainly from marine or alluvial sources. These groundwater basins are mostly bounded by mountains and faults that can act as natural barriers of deep groundwater flow; they are also geographically partitioned by municipal borders (e.g., the boundary between Santa Ana and LA Central basins) and adjudicated to different local water departments.

The northwestern portion of LA Central and Santa Ana basins is characterized as an ‘unconfined zone’ or ‘forebay’45,46 (Fig. 1b), mainly consisting of aquifers made of coarser-grained sediments (sands and gravels) with high permeability. In contrast, the southwestern portion is characterized as a ‘confined zone’45,46, where lenses of aquifers are interbedded with laterally extensive, relatively impermeable lenses of finer-grained sediments (silts and clays) (see Supplementary Note 1). Aquifers in confined to semi-confined conditions also exist in the San Gabriel Valley17. As rainfall is limited to winter and heavy well pumping mostly to dry summers, the groundwater storage in CLAB exhibits strong seasonal variability. In response to the seasonal groundwater fluctuations, confined aquifers show large deformation (leading to up to 60 mm annual uplift or subsidence at the ground surface10,14,15,16,17) due to the high compressibility of clays, whereas the deformation of unconfined aquifers is typically minute due to the low compressibility of sands and gravels.

Variations in seismic wavespeed are inferred from Green’s functions obtained by cross-correlation (at different times) of seismic noise that is recorded continuously at stations in CLAB. We use a broad frequency band of 0.2–2.0 Hz and three sub-bands of 0.2–0.8, 0.4–1.6, and 0.5–2.0 Hz to compute time shifts along the seismic coda, which allows us to probe temporal changes at different depth28,38,47,48. To obtain the spatial distribution of Δv/v, we divide the study area into a 2 km-by-2 km grid (on the latitude-longitude plane) and invert for Δv/v at each grid node for consecutive dates, based on the newly developed coda-wave sensitivity kernels (see, Methods subsection Imaging Δv/v in space).

Temporal fluctuations

Figure 2a shows the Δv/v time series from 2000 to 2020 (blue line) obtained from the 0.2–2.0 Hz band and averaged over all spatial grids in the study area. We plot Δv/v in reverse direction for a more intuitive comparison with groundwater: Upwards corresponds to more groundwater and downwards to less groundwater. Also displayed is the long-term trend of Δv/v (gray line) and the annual cumulative precipitation (light blue bars at the bottom).

Fig. 2: Time series of Δv/v (relative changes in seismic velocity), precipitation, and hydraulic head.
figure 2

a The comparison between Δv/v (bold blue line) and annual cumulative precipitation (light blue bars). Δv/v is measured in 0.2–2.0 Hz and averaged over all spatial grids in the study area. The annual cumulative precipitation is measured at the meteorological station ANA (Fig. 1c). The gray line is obtained by low-pass filtering the blue line to show the long-term trend of Δv/v. The shaded areas (in gray) denote humid periods. b A local-scale comparison between Δv/v (bold blue line) and the hydraulic head (yellow line with circles). Δv/v is averaged over spatial grids within the Santa Ana Basin. Head data are measured at a dedicated groundwater monitoring well GGM-1/MP1 (Fig. 1c).

Full size image

Over the 21-year period, the blue line (that is, the opposite of Δv/v) in Fig. 2a exhibits a predominantly decreasing trend, suggesting an overall decline of groundwater associated with the long-term dearth of precipitation. The decreasing trend is interrupted by brief increasing episodes, corresponding to the years with a few big storms. The blue line also exhibits seasonal fluctuations superimposed on the long-term trends of Δv/v: In each year, Δv/v peaks around January–February and September with annual difference on the order of 0.01%. Featured by their phases and amplitudes, the seasonal variations of Δv/v mainly stem from annual hydrological cycles24. These seasonal and long-term fluctuations will be analyzed separately in the following sections.

For examination on a more local scale, we focus on the Santa Ana Basin. Figure 2b shows Δv/v averaged over spatial grids only within the Santa Ana Basin and the hydraulic head measured at a local monitoring well. This comparison demonstrates that during the time period under study, Δv/v reflects the groundwater level regarding both seasonal fluctuations and long-term trends.

The good agreement between Δv/v and the hydraulic head illustrates the promise of leveraging the large number of existing seismometers in California (and elsewhere worldwide) for tracking groundwater levels. This seismological alternative has three major advantages: First, Δv/v can enhance considerably the spatial and temporal density of the in-situ well measurements, while alleviating the high cost of drilling and maintaining the dedicated monitoring wells. Secondly, Δv/v is affected less by localized heterogeneity and more representative of the average hydrological conditions than the point-scale hydraulic head, because Δv/v integrate the medium properties collectively along the pathways of coda waves. Thirdly, hydraulic head is not always proportional to groundwater storage, and most often the latter is the quantity of practical interests. In this regard, measurements of Δv/v offer a new perspective (on mechanical properties of the aquifer medium), which can be integrated with the head data to constrain the groundwater volume.

Imaging seasonal variations

To investigate the spatial distribution of seasonal changes, we extract the amplitude of seasonal Δv/v (Supplementary Fig. 1) in different frequency bands (0.5–2.0 Hz, 0.4–1.6 Hz, and 0.2–0.8 Hz), which relate to medium properties at different depths (~350, 500, and 700 m, respectively) below the Earth’s surface (see Methods subsection Depth sensitivity).

Figure 3a–c shows that the areas of large seasonal Δv/v (in warm color) migrate southeastwards (from LA Central Basin to Santa Ana Basin) with increasing depth, suggesting that seasonal variations of aquifer storage deepen in that direction. Moreover, in each frequency (depth) range considered, the dominant seasonal Δv/v occur in the confined zones of these basins, with sharp amplitude gradients along the forebay boundary. This indicates that the seasonal variations of hydrologic properties are much more significant in confined aquifers than in unconfined aquifers.

Fig. 3: Seasonal variabilities of Δv/v and surface deformation.
figure 3

a–c Maps showing the seasonal amplitude of Δv/v measured in three decreasing frequency bands corresponding to increasing depths: a 0.5–2.0 Hz, corresponding to about 350 m depth; b 0.4–1.6 Hz, corresponding to about 500 m depth; c 0.2–0.8 Hz, corresponding to about 700 m depth. The depth sensitivities of seismic waves in different frequency bands are demonstrated in Supplementary Fig. 2. d The map of the seasonal amplitude of vertical displacement at the Earth’s surface inferred from InSAR10. Warmer colors indicate stronger seasonal variability, and colder colors weaker seasonal variability. Δv/v maps are derived using records from 2000–2011, and the InSAR map10 during 1992–2011.

Full size image

Seasonal variations can also be inferred from InSAR10 (Fig. 3d), but the two types of data are different and complementary. InSAR measures the deformation at Earth’s surface, which mainly manifests changes in pore pressure and thickness of aquifer layers. In contrast, Δv/v measures the mechanical properties of crustal materials at depth, and thus reflects in-situ variations in pore pressure, saturation, mass density, and micro-structures. Measurements of surface deformation and Δv/v have very different sensitivities to some particular types of subsurface changes. For example, saturation changes in vadose zone may be invisible to InSAR but can lead to considerable changes in Δv/v. Furthermore, Δv/v is less affected by topography, weathering, or vegetation than InSAR.

Despite these differences, the maps of depth-dependent seasonal Δv/v (Fig. 3a–c) and surface deformation (Fig. 3d) show largely consistent patterns of lateral variations, all with dominant changes in the confined zones of LA Central and Santa Ana Basins. The comparison between InSAR and Δv/v maps confirms the expectation that InSAR-derived surface deformation reflects, to some degree, the integration of medium changes across different depths. Besides lateral variations on the surface, Δv/v images further help to characterize the depth and thickness of the aquifer layers, which can be used by decision-makers in water agencies to optimize drilling and pumping strategies and, thus, to promote the efficiency and sustainability of groundwater management.

Imaging cumulative variations

Across the study area, Δv/v indicates an overall reduction of groundwater over the past two decades (Fig. 2a), reflecting a shift to a more arid climate. But in individual basins in CLAB, the long-term behaviors appear different (see Supplementary Fig. 3 for basin-wide Δv/v time series). To assess the spatial variations of long-term changes, we map Δv/v accumulated from 2000 to 2020 in Fig. 4a. It suggests that the San Gabriel and LA Central Basins stores less groundwater now than two decades ago; in the Santa Ana Basin, however, groundwater seems to recover with a slight increase than in the year 2000. In particular, the observations reveal opposite signs of cumulative groundwater change in adjacent basins (San Gabriel and LA Central versus Santa Ana), even though the only boundary between LA Central and Santa Ana basins is a county line, that is, a geopolitical boundary rather than a natural barrier for groundwater flow.

Fig. 4: Long-term Δv/v in different groundwater basins.
figure 4

a Map showing the Δv/v accumulated from 2000 to 2020 (measured in 0.2–0.8 Hz). Colder colors associate with cumulative increases in groundwater and warmer colors with cumulative declines. Note that the boundary between the Santa Ana Basin and the LA Central Basin is a municipal border (i.e., the county line between LA and Orange counties) rather than a natural hydraulic barrier. See Supplementary Fig. 4 for the accumulated Δv/v in 0.5–2.0 Hz (corresponding to shallower depth.) b–d Relationships between two seasonal components of Δv/v (that is, the rainfall-induced decrease in winters and the pumping-induced increase in summers) in b Santa Ana Basin, c LA Central Basin, and d San Gabriel Basin. The rainfall-induced Δv/v are averaged over three-year (backwards) moving windows to downweigh the effect of a few extreme climate events (following the protocol of OCWD45). The red line in b denotes the best-fit line (with a slope close to 1).

Full size image

To explain the distinct cumulative patterns in adjacent basins, we analyze the relationship between the two major seasonal factors that influence the annual hydrologic cycles: the (natural) rainfall that recharges groundwater in winters and the (anthropogenic) well pumping that depletes groundwater in summers. The long-term groundwater change is primarily the year-to-year accumulation of annual differences between these two factors. In Santa Ana Basin (Fig. 4b), the linear relationship between pumping- and rainfall-induced changes suggests that larger volumes of groundwater were pumped out in humid periods and less in arid periods and that in the long run the groundwater drawdown (by pumping) balances the recharge (by precipitation). As a result, total storage has been sustained. In contrast, in the LA Central and San Gabriel Basins (Fig. 4c, d), we find little evidence of a balance between groundwater withdrawal and replenishment. Instead, we observe more years when pumping-induced changes exceed rainfall-induced changes, suggesting that repeated groundwater overdraft leads progressively to excess depletion.

This comparison suggests that the pumping strategy of each municipal water district is the main reason for the distinct Δv/v patterns in individual basins, and that anthropogenic activities, compounding the effect of climate change, play a significant role in shaping the hydrological systems. Effective management of groundwater resources is challenging in water-stressed regions worldwide due to a mixture of environmental, socio-cultural, economic, political, and historical factors. Nevertheless, the Santa Ana Basin case illustrates that putting into place sustainable pumping practices (which can now be verified and quantified via seismic interferometry) can indeed make a difference, setting an encouraging paragon to mitigate the effects of the warming climate.

The basin-wide seismological monitoring presented here can complement the traditional way of assessing groundwater storage (based on streamflow, precipitation, and well measurements), which is subject to uncertainties due to the lack of accurate aquifer models. This seismic approach provides an independent estimate of groundwater budget and safe yield and can thus aid in the policy-making in a sustainable manner.

Discussion

This study presents an unconventional approach to the 4D monitoring of underground hydrologic processes, benefitting from our advances in seismic interferometry techniques. With a pilot application in CLAB, we show that Δv/v recovers the hydraulic head measured in groundwater wells during 2000–2020, illustrating its potential for enhancing the temporal and spatial density of isolated well measurements. Images of Δv/v seasonality agree with surface deformation inferred from InSAR, but also enable the characterization of aquifer behaviors and hydrology at different depths. Based on assessment of long-term Δv/v, we find distinct patterns (decline or slight increase) of groundwater change accumulated in adjacent water districts, due to the effect of anthropogenic pumping practices compounding the climate change. This pilot application shows the promise of leveraging seismometers worldwide to monitor, image, and evaluate underground hydrologic processes. As an in-situ measure of subsurface material property, we anticipate Δv/v to be a powerful 4D geodataset that will bring unique perspectives on constraining near-surface processes and assessing the ever-increasing impact of human activities on shaping the Earth’s shallow subsurface. Integrated with other geodata, the seismological approach presented here can provide a comprehensive understanding of various environmental processes.

Methods

Passive seismic monitoring

The traditional way of passive seismic monitoring provides continuous measurements of Δv/v (relative changes in seismic velocity). This approach incorporates the seismic interferometry in two steps: (1) interferometry of the continuous seismic ambient noise (2) interferometry of seismic coda waves. First, cross-correlating the continuous ambient noise recorded at two seismic receivers allows to estimate the Green’s functions repetitively. Green’s function is the impulse response of the medium at one receiver as if there were a source at the other, and it contains information about structures and elastic properties of the crustal medium between the two receivers25,26. Repeating noise interferometry at different calendar times yield estimates of Green’s functions at consecutive dates. Second, interferometry of coda waveforms gives measurements of Δv/v between different dates. Coda waves are the late arrivals resulted from multiply scattered waves in the reconstructed Green’s functions. By applying cross-spectrum analysis29,47, coda-wave interferometry enables the measurements of tiny perturbations in seismic velocity (Δv/v on the order of 10−4)37,49,50. These two techniques have been successfully applied together to study a wide range of time-varying crustal processes31,32,33,34,35,36,37,38.

In this study, we apply the seismic interferometry using the existing, dense networks of broadband seismic stations in Los Angeles Metropolitan area (Fig. 1c) from the Southern California Earthquake Data Center (SCEDC). We make use of the continuous records of seismic ambient noise at all three components from about 50 stations during 2000–2020. The noise records are pre-processed by removal of instrumental response, whitening over 0.08–8.0 Hz in spectral domain and 1-bit normalization in time domain. We compute daily cross-correlations of pre-processed noise between all station pairs within 50 km and stack over 20 days. Then we apply interferometry between coda-waveforms from the cross-correlations at consecutive times (with 5-day step).

Imaging Δv/v in space

In the traditional framework of Δv/v passive monitoring described above, the coda-wave interferometry is based on the assumption of homogeneous medium perturbation29,30. In other words, previous Δv/v studies only consider the temporal variations, whereas in space the inhomogeneous distribution of Δv/v are typically ignored and only averaged changes are estimated. However, in most monitoring applications, the spatial distributions of the perturbations are crucial for understanding the physical origin and mechanisms of the dynamic processes.

In this study, we go one step further to image the inhomogeneous Δv/v distribution in space, by advancing the recent efforts to solve inverse problems51,52,53,54,55 based on coda-wave sensitivity kernels44,56,57,58,59. The sensitivity kernels (Supplementary Fig. 5) describe the propagation of coda waves in a statistical sense by delineating the likelihood of the coda waves’ travel path. Here we make use of the newly developed sensitivity kernels based on radiative transfer equation44,59. These new kernels provide more accurate descriptions of sensitivity taking into account the actual anisotropy of the scattered wave fields (compared to previous kernels derived under diffusion approximation56,57). Our measurements concern the early coda, which is primarily scattered surface waves53,54,59 and are thus particularly important for Δv/v imaging at shallow depth (e.g., the typical depth of groundwater aquifers).

The coda wave sensitivity kernels are then used to connect the spatial distribution of Δv/v and travel-time shifts (Δt) between coda waveforms. In Eq. 1, the right-hand side is the spatial integration of Δv/v at different spatial locations that are weighted by coda sensitivity kernel K, and the left-hand side is the Δt measured at the corresponding lapse time. Note that here the sensitivity kernel is derived for two-dimensional infinite space44 and the integration is conducted over the latitude-longitude plane.

$$\Delta t\left(t\right)=\int\limits_{S}\frac{\Delta v}{v}\big(\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}\big)\cdot {{{{{\rm{K}}}}}}\big(\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup},t{{{{{\rm{;}}}}}}\mathop{{{{{\boldsymbol{{{{{\mathscr{S}}}}}}}}}}}\limits^{\rightharpoonup},\mathop{{{{{\boldsymbol{{{{{\mathscr{R}}}}}}}}}}}\limits^{\rightharpoonup}\big)\cdot {dS}\big(\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}\big)$$
(1)

Next, an inverse problem (Eq. 2) is set up and resolved to give the spatial distribution of Δv/v. The model parameters we invert for in vector m include Δv/v at different spatial locations (Eq. 3); the observations, vector d, are the Δt measurements at different lapse times between different station pairs (Eq. 4); and the design matrix G is constructed by coda-wave sensitivity kernels at the corresponding travel-times for all station pairs (Eq. 5). Supplementary Fig. 5 shows the sum of all kernels for all station pairs used in this study.

$${{{{{\bf{Gm}}}}}}={{{{{\bf{d}}}}}}$$
(2)
$${{{{{\bf{m}}}}}}=\left({\begin{array}{c}\frac{\Delta v}{v}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{1}\big)\\ \frac{\Delta v}{v}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{2}\big)\\ \frac{\Delta v}{v}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{3}\big)\\ \vdots \\ \frac{\Delta v}{v}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{{{{{{\rm{n}}}}}}}\big)\end{array}}\right)$$
(3)
$${{{{{\bf{d}}}}}}=\left(\begin{array}{c}\Delta {t}_{11}\\ \Delta {t}_{12}\\ \vdots \\ \Delta {t}_{16}\\ \Delta {t}_{21}\\ \vdots \\ \Delta {t}_{26}\\ \Delta {t}_{31}\\ \vdots \\ \Delta {t}_{p6}\end{array}\right)$$
(4)
$${{{{{\bf{G}}}}}}=\left(\begin{array}{ccc}\begin{array}{c}\begin{array}{ccc}{{{{{{\rm{{ K}}}}}}}}_{11}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{1}\big) & {{{{{{\rm{{ K}}}}}}}}_{11}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{2}\big) & {{{{{{\rm{{ K}}}}}}}}_{11}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{3}\big)\\ {{{{{{\rm{{ K}}}}}}}}_{12}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{1}\big) & {{{{{{\rm{{ K}}}}}}}}_{12}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{2}\big) & {{{{{{\rm{{ K}}}}}}}}_{12}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{3}\big)\\ & \vdots & \end{array}\\ \begin{array}{ccc}{{{{{{\rm{{ K}}}}}}}}_{16}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{1}\big) & {{{{{{\rm{{ K}}}}}}}}_{16}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{2}\big) & {{{{{{\rm{{ K}}}}}}}}_{16}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{3}\big)\\ {{{{{{\rm{{ K}}}}}}}}_{21}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{1}\big) & {{{{{{\rm{{ K}}}}}}}}_{21}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{2}\big) & {{{{{{\rm{{ K}}}}}}}}_{21}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{3}\big)\end{array}\\ \vdots \end{array} & \begin{array}{c}\begin{array}{c}\cdots \\ \cdots \\ \ddots \end{array}\\ \begin{array}{c}\cdots \\ \cdots \end{array}\\ \cdots \end{array} & \begin{array}{c}\begin{array}{c}{{{{{{\rm{{ K}}}}}}}}_{11}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{{{{{{\rm{n}}}}}}}\big)\\ {{{{{{\rm{{ K}}}}}}}}_{12}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{{{{{{\rm{n}}}}}}}\big)\\ \vdots \end{array}\\ \begin{array}{c}{{{{{{\rm{{K}}}}}}}}_{16}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{{{{{{\rm{n}}}}}}}\big)\\ {{{{{{\rm{{K}}}}}}}}_{21}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{{{{{{\rm{n}}}}}}}\big)\end{array}\\ \vdots \end{array}\\ \begin{array}{c}\begin{array}{ccc}{{{{{{\rm{{ K}}}}}}}}_{26}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{1}\big) & {{{{{{\rm{{ K}}}}}}}}_{26}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{2}\big) & {{{{{{\rm{{ K}}}}}}}}_{26}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{3}\big)\\ {{{{{{\rm{{ K}}}}}}}}_{31}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{1}\big) & {{{{{{\rm{{ K}}}}}}}}_{31}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{2}\big) & {{{{{{\rm{{ K}}}}}}}}_{31}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{3}\big)\end{array}\\ \vdots \end{array} & \begin{array}{c}\begin{array}{c}\cdots \\ \cdots \end{array}\\ \ddots \end{array} & \begin{array}{c}\begin{array}{c}{{{{{{\rm{{ K}}}}}}}}_{26}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{{{{{{\rm{n}}}}}}}\big)\\ {{{{{{\rm{{ K}}}}}}}}_{31}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{{{{{{\rm{n}}}}}}}\big)\end{array}\\ \vdots \end{array}\\ \begin{array}{ccc}{{{{{{\rm{{ K}}}}}}}}_{{{{{{\rm{p}}}}}}6}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{1}\big) & {{{{{{\rm{{ K}}}}}}}}_{{{{{{\rm{p}}}}}}6}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{2}\big) & {{{{{{\rm{{ K}}}}}}}}_{{{{{{\rm{p}}}}}}6}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{3}\big)\end{array} & \cdots & {{{{{{\rm{{ K}}}}}}}}_{{{{{{\rm{p}}}}}}{{{{{\bf{6}}}}}}}\big({\mathop{{{{{{\bf{r}}}}}}}\limits^{\rightharpoonup}}_{{{{{{\rm{n}}}}}}}\big)\end{array}\right)$$
(5)

Here we solve this inverse problem using a Bayesian approach60 that maximizes the posterior probability assuming Gaussian distributions. In Eq. 6, CD is the covariance matrix determined by the cross-correlation coefficient of the two small windows of coda waveforms for Δt calculation, and CM is the smoothing matrix that decays exponentially with distance.

$$\widetilde{{{{{{\bf{m}}}}}}}={\left({{{{{{\bf{G}}}}}}}^{{\prime} }{{{{{{\bf{C}}}}}}}_{{{{{{\rm{D}}}}}}}^{-1}{{{{{\bf{G}}}}}}+{{{{{{\bf{C}}}}}}}_{{{{{{\rm{M}}}}}}}^{-1}\right)}^{-1}{{{{{{\bf{G}}}}}}}^{{\prime} }{{{{{{\bf{C}}}}}}}_{{{{{{\rm{D}}}}}}}^{-1}{{{{{\bf{d}}}}}}$$
(6)

The workflow of our data processing is summarized in Supplementary Fig. 6.

In this study, we divide the study area into 2 km-by-2 km grids (on horizontal plane) and calculate the coda-wave sensitivity kernels using a homogeneous model with scattering mean free path of 100 km. We then invert for the Δv/v distribution on each grid. For each station pair, we compute Δt at 6 lapse-time via coda-wave interferometry. Using Δt calculated for all station pairs, the spatial inversions of Δv/v are performed repetitively on consecutive dates (with 5-day steps) from 2000 to 2020. In the main text, we average the time series of Δv/v over all grids in the study area (Fig. 2a) for an overall evaluation of temporal behavior; we also average over grids only within Santa Ana Basin (Fig. 2b) for a local-scale analysis.

Depth sensitivity

To characterize the near-surface aquifers, in this study we use the early coda (with arrival time between 15 to 50 s), mainly consisting of scattered surface-wave content52,53,58 that are most sensitive to shallow depth. As commonly used in surface wave studies, the seismic waves in different frequencies propagate and thus probe medium properties at different depth. Here we derive the depth sensitivity of Δv/v measured in different frequency bands using kernels of Rayleigh waves61 based on a 1D average velocity model62 in Los Angeles basin (Supplementary Fig. 2a). As shown in Supplementary Fig. 2b, the three frequency bands, 0.5–2.0, 0.4–1.6, and 0.2–0.8 Hz, used for Δv/v measurements are most sensitive to medium perturbations at around 350, 500, and 700 m beneath the surface.