CC-BY 4.0
1. Introduction
1.1. Broad overview
The Critical Zone (CZ) is Earth’s near-surface region where freshwater, air, soil, rock, and living organisms interact, extending from the vegetation canopy at the surface down to freely circulating groundwater stored in aquifers (Giardino and Houser, 2015). Monitoring soil moisture in the CZ and understanding its dynamics is critical to improving agricultural management, monitoring droughts, floods, and groundwater recharge, and monitoring the stability of slope failures, and thus has implications for food security, climate systems, and surface hazards (floods, landslides, and earthquake-induced liquefaction). The spatial heterogeneity of soils’ properties (e.g., hydraulic conductivity) and their potential for evapotranspiration, geological structure, meteorological forcing (e.g., precipitation, wind, and radiation), and vegetation control the spatial heterogeneity of soil moisture (e.g., Chaney et al., 2015). Conventional soil moisture monitoring is through (1) in situ soil moisture probes (Su et al., 2014) that have local sensitivity or (2) remote or airborne data that use electromagnetic wave signals (amplitudes and phase) to estimate soil moisture and have regional or global sensitivity (Babaeian et al., 2019). As stated by Flinchum, Holbrook, et al. (2022), “a major challenge in CZ science is to obtain measurements that have both high spatial resolution and cover large areas”.
Passive ambient noise monitoring refers to the methodology of using repeated cross-correlation of the ambient seismic field, recorded at a single or a pair of seismic stations, to estimate the changes in seismic properties in the subsurface. We focus this review on the scattered waves in the coda of coherent ambient seismic waves, which are typically surface waves in the context of soils and sediments. “Ambient noise monitoring” is transforming environmental seismology as it provides new passive methods to monitor the mechanical response of the shallow subsurface (e.g., Larose et al., 2015; Le Breton et al., 2021; Wang and Yao, 2020). Temporal variations in near-surface seismic properties stem from a range of interacting factors whose relative importance is frequently site-dependent, underscoring the spatial heterogeneity of soils. This review sets out to build a unified framework for monitoring hydromechanical behavior in the critical zone with ambient-noise seismology, examining the method’s strengths and limitations.
Relative changes in shallow seismic velocities, dv/v hereafter, have been observed to vary between 0.1% and 50% under many tectonic and volcanic factors (e.g., Wegler and Sens-Schönfelder, 2007; Brenguier, Campillo, et al., 2008; Brenguier, Shapiro, et al., 2008; Obermann, Froment, et al., 2014; Yukutake et al., 2016; Feng, Huang and Wu, 2020; Cabrera-Pérez et al., 2023; Donaldson et al., 2019; Okubo et al., 2024; Esfahani et al., 2024) and solid Earth tides (e.g., Sens-Schönfelder and Eulenfeld, 2019; Mao, Campillo, et al., 2019; Takano and Nishida, 2023; Kramer, Y. Lu and Bokelmann, 2023). In addition to the rise in analyses monitoring temporal changes in the subsurface induced by earthquakes or volcanic activity, there are numerous environmental factors that change seismic properties of the subsurface: varying snow load (Q.-Y. Wang, Brenguier, et al., 2017; Donaldson et al., 2019), rainfall load and diffusion (Q.-Y. Wang, Brenguier, et al., 2017; Andajani et al., 2020; Feng, Huang, Hsu, et al., 2021; Clements and Denolle, 2023), atmospheric pressure fluctuations (Gradon et al., 2021; Kramer, Y. Lu and Bokelmann, 2023; Kramer, Y. Lu, Q.-Y. Wang, et al., 2025), river and lake levels (e.g., Rodríguez Tribaldos and Ajo-Franklin, 2021), thawing and freezing of soils (e.g., James, Knox, Abbott, Panning, et al., 2019). Figure 1 represents the vast range of time and spatial scales at which the hydrological phenomena have been reported to carry seismic signatures in the literature reviewed in this manuscript.
Range of spatial and time scales of the hydrological processes reported in the literature. Left: lateral spatial extent ranging from centimeters to 1000 km. Right: depth scale from centimeters to 100 m. Solid boxes represent the range of time and spatial scales where seismic signatures correlate with specific hydrological processes. Long-term climate effects have been reported at all spatial scales and depths.
The core principle of ambient noise monitoring is that repeated measurements of the seismic wavefield, measured through the cross-correlation of the ambient seismic field between channels and stations, capture the seismic signatures of changes in effective stress or pore pressure. Coda waves extracted from “noise correlation functions” (NCFs) are commonly used in ambient noise monitoring, as they are less affected by variations in noise sources than direct waves (e.g., Hadziioannou, Larose, Coutant, et al., 2009; Colombi et al., 2014). Recently, using the direct, or ballistic, waves of the NCF, typically either P-waves or surface waves (e.g., Garambois et al., 2019; Takano, Brenguier, et al., 2019; Mordret et al., 2020; Brenguier, Courbis, et al., 2020; Sheng et al., 2024), or specific scattered phases such as P–S (e.g., Y. Lu and Ben-Zion, 2021), has helped improve the spatial localization of velocity changes. However, the use of ballistic waves for monitoring requires some knowledge about stability and location of repeated anthropogenic sources (e.g., trains, Higueret et al., 2024). In addition, coda waves are more sensitive than ballistic waves to small changes in dv/v, as they sample the medium many more times due to multiple scattering (Snieder, 2002; Colombi et al., 2014). With the assumption that the scatterers themselves do not change positions, the perturbation of the seismic velocity averaged along the coda wave path is linearly proportional to the accumulated phase shift: dv/v = −dt/t, where t is the phase lag and dt is the phase shift in the coda (Poupinet et al., 1984). This review focuses on the observations made with coda waves but will highlight the potential for ballistic wave analysis in relevant cases. Beyond resource monitoring, many of the most damaging natural hazards, such as landslides, earthquake-induced liquefaction, and nonlinear site amplification, are also governed by hydromechanical processes in the Critical Zone. These hazards arise from the same physical mechanisms discussed throughout this review: changes in effective stress, pore pressure, saturation, and soil structure. In particular, the response of near-surface soils to strong shaking is modulated by water content and stress history, linking seismic hazard to hydrological variability. By framing these hazards as manifestations of CZ dynamics, ambient field seismology offers a unifying observational approach: the same seismic observables (e.g., dv/v, attenuation (Q), and the Horizontal-to-Vertical Sprectral Ratio (HVSR)) used to track groundwater or vadose (unsaturated) processes can also be used to detect precursors to slope failures, assess site amplification conditions, and evaluate liquefaction susceptibility. This review aims to bridge these application areas and demonstrate how environmental seismology contributes to an integrated understanding of both hydrogeophysical processes and geohazards in the CZ.
1.2. Seismic structure of the near surface
The seismic structure of the shallow subsurface, especially in the critical zone, is marked by a surface low-velocity zone with shear wavespeeds (VS) ranging from 130–1000 m/s spanning unsaturated and saturated soil, densities ranging from 1.7–2 g/cm3, and P wavespeed (VP) of 350–2500 m/s (e.g., Befus et al., 2011; Boore, 2016; Nagashima and Kawase, 2021). The shallowest layers are often highly weathered and have a greater sensitivity to hydrological processes at and above the water table depth than more consolidated sediments (Anderson et al., 2013; Okay and Özacar, 2023). Scattering attenuation in sediments, generated by higher porosity and crack density in shallow sediments, is often equal to or greater than intrinsic attenuation, which is dispersed as heat (Jemberie and Langston, 2005; Eulenfeld and Wegler, 2016; Y.-P. Lin and Jordan, 2023). Combined intrinsic anelastic and scattering effects may yield strong attenuation in shallow sediments, with total quality factors (Q) as low as 2–10 for S and P waves in shallow soils (Malagnini, 1996).
2. Relating seismic properties to water content
Relating seismic measurements into hydrological processes requires understanding three topics: (1) hydrological data and, when lacking measurements, how seismologists infer it; (2) rock physics models to estimate elastic and anelastic properties with varying water content; (3) seismic wave spatial and depth sensitivities depending on the wave types.
2.1. Hydrological models
Water content, either in the form of water saturation or pore pressure, is a key element in the hydromechanical response of the critical zone. Due to the logistical challenges of obtaining high-resolution in situ measurements of water content and the resulting sparsity of data coverage, numerical and process-based models are often used to understand hydrological processes within the critical zone. The water budget includes sources (precipitation, snow melt, and runoff from uphill) and sinks (drainage, runoff, evapotranspiration).
The main approach that seismologists have used to estimate pore pressure in the fully saturated regime is through coupling infiltration from rainfall or snow melt with simplified outflow models based on either baseflow (Sens-Schönfelder and Wegler, 2006; Feng, Huang, Hsu, et al., 2021) or the poroelastic response to surface water input (E. Roeloffs, 1996; E. A. Roeloffs, 1998; Talwani et al., 2007; Tsai, 2011; Andajani et al., 2020; Q.-Y. Wang, Brenguier, et al., 2017). For regions with a significant soil moisture reservoir, the antecedent precipitation index has been used to represent the capacity of the vadose zone (Illien, Andermann, et al., 2021). Clements and Denolle (2023) found that baseflow and poroelastic-based models of pore pressure, along with the cumulative deviation from the moving mean of precipitation (Smail et al., 2019), equally worked as long-term proxies of groundwater level when compared to hydrologically-induced dv/v.
Hydrological models that account for the full water budget with inputs and outputs (runoff, downward drainage, and evapotranspiration) are more conducive to modeling vadose zone dynamics. Lecocq et al. (2017) modeled water storage as a budget between input (precipitation) and output (river runoff and evapotranspiration) over a time period of three decades. Z. Shen et al. (2024) modeled water storage as a mass balance between precipitation (input), evapotranspiration modulated by surface moisture measurements and temperature (output), and deep water drainage (i.e. groundwater downward draining as an output). Additionally, groundwater levels can be inferred from river (Rodríguez Tribaldos and Ajo-Franklin, 2021) or lake (Clements and Denolle, 2023) levels, as surface body water levels are often hydraulically connected to the groundwater table.
2.2. Hydromechanical models of Earth materials
Theoretical, laboratory, and field studies have demonstrated that both VP and VS vary with water saturation in rocks (Nur and Simmons, 1969; O’Connell and Budiansky, 1974; Kuster and Toksöz, 1974). In partially saturated rocks, seismic velocity perturbations caused by changes in water content are controlled mainly by the density effect. Increasing bulk density through water saturation, achieved by replacing air with water, results in a gradual reduction in both VP and VS. When the rock achieves nearly fully saturated conditions (80%–90%), the non-compressibility of the fluid in the pores dominates the density effect, and VP experiences a sharp increase. Commonly, VS continues to decrease with saturation. In a few laboratory cases, however, studies have shown that a small increase in VS can also be observed when saturation is close to 100%, which may be due to the effect of surface tension in small pores causing an increase in the shear rigidity of the rock (Knight and Nolen-Hoeksema, 1990). This increase in VP and VS at full saturation will mostly be observed near the water table and the capillary fringe. Beyond this general behavior, studies have shown that other hydrological processes can influence the behavior of seismic velocities in the unsaturated zone (e.g., see below discussion of Section 5.1). For example, strong matric suction effects can lower the effective stress and cause a consequent increase in VP, which competes against the VP reduction caused by the density effect (Linneman et al., 2021).
Below the water table, under fully saturated conditions, pore structure controls the relative change in seismic velocity. For low-porosity rocks with crack-shaped pores, VP increases with full saturation, whereas VS remains unchanged with saturation (Nur and Simmons, 1969). For high porosity rocks with circular pores, VP increases with saturation whereas VS decreases with saturation (Toksöz et al., 1976). At shallow depths, where confining pressure is low, increases in pore pressure thus increase VP and decrease VS. At the near-surface, Rayleigh waves are far more sensitive to variations in VS than VP (Song et al., 1989). Thus, numerous observations have been made of a decrease in Rayleigh wave speed (VS) with increasing groundwater levels or pore pressure. The constitutive relations for nonlinear elasticity show that relative change in elastic wavespeeds, dv/v, is proportional to the dilatational strain 𝜖kk (Murnaghan, 1937; Ostrovsky and Johnson, 2001), and observations suggest that:
| \begin {equation} \mathrm {d}v/v = -\beta \epsilon _{kk}, \end {equation} | (1) |
| \begin {equation} \mathrm {d}v/v = - \beta s_{sk} \,\mathrm {d}h. \end {equation} | (2) |
| \begin {equation} \mathrm {d}v/v= \mathrm {cov}(\mathrm {d}v/v,p)/\mathrm {var}(p) , \end {equation} | (3) |
In the partially saturated regime, above the water table, partial saturation effects on seismic wavespeeds are less trivial. The first-order effect is the change in density that occurs by replacing the air in the pores with water:
| \begin {equation} \rho = (1 - \phi ) \rho _s + \phi ( S_w \rho _w + (1- S_w)\rho _a ), \end {equation} | (4) |
Generic seismic profile of a critical zone from Solazzi et al. (2021). Left plot: saturation increases from the surface until the water table (left). Middle plot: elastic physical properties (density 𝜌, P wave velocity VP, S wave velocity VS) as a function of depth. Right plot: Surface wave dispersion as a function of frequency with (blue lines) and without (red dashed lines) capillary forces included in the model. Solazzi et al. (2021) showed that when the VP/VS ratio, which is usually a metric for water content, is about constant in the partially saturated zone, surface wave velocities are sensitive to saturation due to capillary suction effects.
It should be noted that VS’s sensitivity to water level greatly depends on the lithology: it is strong for clays but small for sand (Dong and N. Lu, 2016). Such a difference may explain why most seismological studies exhibit spatially variable sensitivity.
Changes in seismic attenuation, Q, in the partially saturated zone are less obvious. P-wave attenuation is greatly sensitive to partial saturation due to the contrast in the compressibility of mixed fluids (air vs. water) (e.g., Chapman et al., 2021; Amalokwu et al., 2014), whereas S waves are insensitive to compressible fluids. Numerous studies have correlated values of Qp/Qs⩽1 with the presence of air/gas in a partially saturated medium (e.g., Hudson et al., 2023). Nevertheless, theoretical and empirical (Winkler and Nur, 1982; Diewald et al., 2024) studies have shown that S-wave or S-wave-dominated seismic attenuation varies with saturation. Seismic attenuation within partially saturated porous materials is mainly caused by the relative movement between the solid matrix and the pore fluid. This movement generates shear stresses in the fluid, leading to energy loss through viscous friction (Mavko and Nur, 1979).
2.3. Wave types and their various spatial and depth scales
Direct ballistic waves can be used to probe large changes in the subsurface, as phase shifts in the direct waves may be precise enough to estimate the seismic velocity changes. By using direct waves, the localization of the velocity perturbation is achieved simply through time-lapse tomography. Using ambient noise seismology, studies have focused on time-lapse shear wave tomography using surface wave dispersion curves (e.g., Dou et al., 2017; Mordret et al., 2020; Cheng et al., 2022; Sun et al., 2025). The typical frequencies used in these analyses range from 2 to 15 Hz, providing sensitivity to the shallow seismic structure from a few meters to approximately 40–50 m in depth. The spatial resolution of these measurements is typically on the order of tens of meters. Ballistic wave-based tomograpy is limited by the density of the seismic channels and the quality of the high-frequency signals. Besides surface waves, several studies have leveraged distinct seismic phases, such as refracted P waves (Brenguier, Courbis, et al., 2020) or converted phases (Y. Lu and Ben-Zion, 2021), when possible. One key element of using ballistic waves is that their spatial and depth sensitivity depends on the source-receiver geometry and the Earth’s structure (ibid.). Because time-lapse tomography, which infers absolute seismic velocities in space (laterally and in depth) and time, rock physics models can extract temporal and spatial variations in hydro-thermo-mechanical properties of the critical zone (e.g., Sun et al., 2025).
When distinct seismic phases are undistinguishable or when the changes are subtle, seismologists use coda waves instead to image changes in the subsurface. Because scattered waves accumulate travel-time perturbations due to changes in seismic velocities, phase shifts are increased in late coda waves. For homogeneous velocity perturbations, such delays accumulate linearly regardless of the wave type involved (e.g., Obermann, Planès, Larose, et al., 2013).
The depth sensitivity of coda waves is most often assumed to be that of fundamental mode Rayleigh wavespeed, given a perturbation in VS as a function of depth and frequency. Sensitivity kernels, which link observed wavespeed perturbations (dv/v) to true VS perturbations between seismic stations often rely on solutions from the radiative transfer model for short lag-time measurements (Pacheco and Snieder, 2005; Planès et al., 2014). A good rule of thumb is that the penetration depth of the Rayleigh wave is around one horizontal wavelength (Barajas, Margerin, et al., 2022). When seismic velocity and density are relatively constant with depth, scattered body waves dominate the late coda of NCFs, though sensitivity kernels may couple surface and body wave sensitivity at depths below the critical zone (Obermann, Planès, Larose, et al., 2013; Obermann, Planès, Hadziioannou, et al., 2016; Obermann, Planès, Larose, et al., 2019; Margerin et al., 2019; Barajas, Margerin, et al., 2022).
3. Review of sensing technologies
We now focus on the technologies that have served the passive seismological community to investigate critical zone hydrology.
3.1. Seismological observation networks
3.1.1. Permanent seismic networks
The most commonly used sensors for ambient noise monitoring of the critical zone are broadband seismometers, typically deployed directly at the surface. Broadband seismometers have been deployed since the 1960s, although continuous digital recording began in the late 1990s. By 2002, most seismic networks had begun continuous recording of their entire broadband networks, which fueled the use of the ambient seismic field for seismological research (e.g., Campillo and Paul, 2003; Shapiro et al., 2005). The deployment of broadband sensors serves two purposes. First, the global detection and monitoring of explosions and earthquakes drive the deployment of these sophisticated stations worldwide, particularly on oceanic islands. With this purpose in mind, broadband sensors typically record at high gain to retrieve weak seismic waves from distant earthquakes and explosions. Second, regional and national seismic networks, with the specific goal of monitoring earthquake sources, have favored the deployment of broadband seismometers on bedrock and near fault zones. Because seismic waves amplify in soft sediments, network engineers have chosen to avoid these sites and focus on the station’s sensitivity to fault zone waves. These locations are not ideal for subsurface hydrology, as the groundwater reservoirs and soils are predominantly located in basins. However, the advantages of broadband stations are their long-standing, multi-decadal, continuous records. Because their sensitivity spans a broad range of frequencies, seismologists have repurposed broadband seismometers to probe the near-surface in the 0.1–10 Hz frequency range, covering both powerful microseismic ambient fields and stationary anthropogenic seismic noise (e.g., Hasselmann, 1963; Díaz, 2016).
Accelerometers are also used for permanent seismic monitoring and typically record at lower digital gain to retrieve high fidelity of strong ground accelerations and higher frequencies (above 0.1 Hz). Regional seismic networks (e.g., the Southern California Seismic Network SCSN, the Northern California Earthquake Data Center NCEDC, and the Pacific Northwest Seismic Network PNSN) comprise several types of seismometers, ranging from broadband seismometers to accelerometers, for both monitoring microseismicity and early warning of strong earthquakes. Regional and global networks of seismometers (e.g., GSN, ORFEUS, Geoscope) archive their data on publicly accessible data servers, either on-premise servers or in cloud storage (e.g., E. Yu et al., 2021). Most of these installations are made at the surface only.
Several observatories pair surface with borehole instrumentation enabling top and bottom analysis of the near-surface environment. After the 1995 Hyogo-ken Nanbu Earthquake, NIED (National Research Institute for Earth Science and Disaster Prevention) constructed an uphole/downhole observation network, KiK-net (Kiban Kyoshin network), with approximately 700 stations (Aoi et al., 2004). Each KiK-net station has a borehole of 100 m or more at depth, and accelerometers have been installed both on the ground surface and at the bottom of the boreholes. Such pairs of sensors enable direct measurement of near-surface dynamics or physical properties between the two sensors. The KiK-net accelerometers in Japan have been unrivaled by other downhole networks in the volume of available acceleration seismograms, given the quality, homogeneity of the information describing subsurface layers, and the rate of seismicity. In addition to national observatories, another model for focused surface-to-borehole observations exists in the geotechnical community, referred to as “liquefaction arrays”. In the United States, such arrays include the Wildlife Liquefaction Array (Holzer and Youd, 2007) and the “Seattle Liquefaction Array” (Steidl and Hegarity, 2017), which feature piezoelectric pore pressure sensors, along with several shallow borehole seismometers, to continuously sample water pressure and earth vibration.
Near-surface structure information of the site on which these seismic observations are made is sometimes available. For instance, the velocity profiles and geological information, as well as the observed records, are widely accessible on the NIED website. C. Zhu, G. Weatherill, et al. (2021) created an open-source site database that contains site characterization parameters directly derived from available velocity profiles, including average wave velocities, bedrock depths, and velocity contrast. Additionally, the site database includes topographic and geological proxies inferred from regional models or maps. Geotechnical arrays also typically conduct surveys of the seismic and hydrological structure (e.g., Youd et al., 2004).
3.1.2. Temporary seismic arrays
Most seismic imaging of the near-surface entails a small aperture array of a few 100 m with seismic channel spacing of a few meters. When analyzing deeper and larger groundwater reservoirs, such as aquifers in sedimentary basins, the analysis aperture may range from 10 km to 100 km. Given these spatial considerations, seismologists use direct or refracted P waves (e.g., Flinchum, Holbrook, et al., 2022), direct surface waves (e.g., Sobolevskaia et al., 2024; Sun et al., 2025), and coda waves dominated by multi-scattered S and surface waves (e.g., Z. Shen et al., 2024) to explore the spatial and temporal variations in the CZ seismic properties.
For temporary seismic deployments, low-cost geophones, called “nodes” such as the FairField 3C or the SmartSolo nodes, have also been used for short-term experiments to image and monitor water table, precipitation, soil moisture, or seasonal frozen ground freeze–thaw cycles using seismic methods (Oakley et al., 2021; Hua et al., 2023; Gochenour et al., 2024; H. Liu et al., 2025; Behm et al., 2019; Hua et al., 2023; Grobbe et al., 2021). Nodes are easy to deploy in “large-N” arrays, fast to install in rapid response to earthquakes or sudden environmental events; their disadvantage is a short battery life of a month that limits long-term monitoring.
A novel technology transforming the field of seismology, particularly in near-surface seismic imaging and monitoring, is distributed acoustic sensing (DAS). DAS technology transforms traditional fiber-optic cables into dense arrays of sensors, enabling the measurement of distributed strain or strain-rate changes along their length (e.g., Hartog, 2017). Current DAS technology can continuously record the seismic wavefield at a spatial resolution of a few meters along cable lengths that can exceed 100 km. Several studies have also demonstrated that DAS can record strain-rate signals over a broad range of frequencies, ranging from 0.01 to hundreds of Hertz (N. J. Lindsey et al., 2019; Paitz et al., 2020), though near-source, high ground motion observations are often “clipped” due to DAS’s limited sensitivity at high strain rates (C.-R. Lin et al., 2025). The primary advantage of DAS lies in its capacity for continuous, real-time monitoring across extensive areas or structures sampling the wavefield over many closely spaced (∼3–10 m) virtual channels using either pre-existing telecommunication fiber-optic cables or short-term fiber deployments. DAS measures strain rate along the axial direction of the cable, making it especially suitable for ambient noise studies, in which Rayleigh waves traveling along the axis of the cable can be analyzed for subsurface imaging and monitoring. DAS arrays produce substantial amounts of raw strain data (∼TBs), which makes storing, sharing, and analyzing the data computationally difficult.
3.2. Hydrological observations networks
Seismologists have attempted to explain the hydrological signature by observing time-dependent seismic wavespeeds. Occasionally, in-situ measurements are possible near or at seismic stations, such as groundwater wells (e.g., C.-Y. Lin et al., 2024) or soil moisture probes (e.g., Illien, Andermann, et al., 2021; Oakley et al., 2021; Sobolevskaia et al., 2024). Typically, the distance between the seismometer and the in-situ hydrological sensor is up to a few kilometers, and the hydraulic conductivity between the two sites is either unknown or highly uncertain. Piezoelectric probes are frequently used to infer water levels and have been employed in several local studies (Voisin, Guzmán, et al., 2017; Gaubert-Bastide et al., 2022).
Networks of groundwater wells are increasingly made available to the public (Goodall et al., 2008), such as the National Ground Water Monitoring Networks maintained in individual countries; however, data rights still pose challenges for data sharing (Condon et al., 2021). Overall, groundwater wells are either pumped for agricultural, domestic, or industrial uses or used to monitor groundwater levels (Margat and Van der Gun, 2013) passively. In the past, one challenge seismologists have faced in accessing groundwater well data to compare with their seismic proxies has been the decentralization of groundwater well observation networks and the non-uniformity in the temporal sampling of the wells.
Determining the soil moisture (SM) content of the CZ is crucial for effective resource management and hazard mitigation during extreme droughts and floods. There are several methods for inferring or measuring soil moisture. Babaeian et al. (2019) provides an extensive review of these approaches, which are categorized as remotely sensed (using satellite or airborne sensors), proximally sensed (e.g., GPS reflectometry), or terrestrially sensed (e.g., time-domain reflectrometry). These methods are only sensitive to the upper 50 cm, limiting analysis to the upper part of the critical zone and lacking sensitivity for deeper processes such as the weathered bedrock (2–30 m). Additionally, ground-based sensors are particularly local, while remotely sensed SM measurements typically have a footprint of 10s of kilometers: empirical downscaling strategies are often required to bridge the scales.
When lacking in situ or remotely sensed measurements, seismologists have turned to inferring water content from a hydrological model (see Section 2.1). For example, Lecocq et al. (2017) used precipitation data together with a runoff, discharge, and evapotranspiration model to estimate groundwater storage.
4. Review of methodologies in seismology
4.1. Active near-surface imaging with body waves
Near-surface geophysics (NSG) typically refers to temporary experiments that involve deploying sensors, transmitting a source, and inverting data to determine water content and structural properties. NSG is a mature field of geophysics, with a vast set of methodologies and dedicated curriculum (Butler, 2005; Everett, 2013). Non-invasive methods for characterizing critical zone water content include electrical and electromagnetic surveys (Siemon et al., 2009), magnetic susceptibility and resonance methods (Hertrich, 2008), ground penetrating radar (Klotzsche et al., 2018), gravity survey (Gehman et al., 2009; Syed et al., 2008), and seismic, on which this review will expand.
Seismic methods for inferring near-surface water content typically use active seismic surveys and short-period seismometers (geophones) (Befus et al., 2011; Flinchum, Holbrook, et al., 2022). Usually, travel time measurements of the direct and refracted P waves are used to generate reflection and refraction images of the near-surface (Befus et al., 2011). The spatial scales of these surveys vary from 10 to 100s of meters, depending on the spacing of the geophones (e.g., Uecker et al., 2023; Flinchum, Grana, et al., 2024). Sophisticated methods, such as full-waveform inversion, which were initially developed to image regional-scale subsurface stratigraphy, are now being applied to the critical zone (X. Liu et al., 2022).
4.2. Time-lapse surface wave tomography
Surface waves are dispersive, and their dispersion relation depends on the vertical shear wave profile. Using either active or passive (e.g., ambient noise cross correlations) sources, surface wave imaging in near-surface geophysics has significantly advanced the monitoring of water content changes in the critical zone. By analyzing surface wave dispersion, these methods utilize the frequency-dependent velocity of Rayleigh waves to infer subsurface elastic properties, which are sensitive to variations in water content. Dispersion inversion techniques, such as multichannel analysis of surface waves (MASW) or the frequency-wavenumber method (Park et al., 1999), allow the creation of shear wave velocity profiles, and when repeated over time, they correlate with changes in soil moisture and porosity (Socco et al., 2010; Cheng et al., 2022; Sobolevskaia et al., 2024; Sun et al., 2025) or shallow aquifer structure (Gochenour et al., 2024; Eppinger et al., 2024).
4.3. Changes in seismic velocity dv/v
The relative change in seismic velocity, dv/v, derived from ambient noise cross-correlation is sensitive to near-surface changes in saturation, pore pressure, and overall effective stress. dv/v monitoring of hydrological changes in the critical zone spans from point measurements using single station cross-channel correlations (Sens-Schönfelder and Wegler, 2006; Donaldson et al., 2019; Kim and Lekic, 2019; Lindner et al., 2021; Feng, Huang, Hsu, et al., 2021; Clements and Denolle, 2023) to kilometer-scale variations using inter-station cross-correlations (Lecocq et al., 2017; Q.-Y. Wang, Brenguier, et al., 2017; Clements and Denolle, 2018; Mao, Lecointre, et al., 2022; Mao, Ellsworth, et al., 2025). We illustrate the clear evidence that dv/v correlates with water content in Figure 3.
Temporal evolution of dv/v relative to water storage seen for (a) soil moisture (green and red lines) by Illien, Andermann, et al. (2021), (b) groundwater (blue line) by many, including Clements and Denolle (2018), and (c) in regional-scale karstic groundwater storage (gray line) by Almagro Vidal, et al. (2021).
The temporal resolution of the dv/v transient depends on the time scale of the process under study, which spans from seconds for earthquake shaking (Illien, Turowski, et al., 2025) and liquefaction (Viens, Bonilla, et al., 2022) to decades for climatic processes (Lecocq et al., 2017; Clements and Denolle, 2018; Mao, Lecointre, et al., 2022; Mao, Ellsworth, et al., 2025), as shown in Figure 1. For long-term monitoring studies and at regional scales, NCFs are averaged in 1–30 day moving windows to increase signal-to-noise ratio (SNR) (Seats et al., 2012), while other studies with short aperture may report sub-hourly temporal resolution in dv/v (Xie et al., 2023; Rodríguez Tribaldos and Ajo-Franklin, 2021). The non-stationarity of the seismic sources and their coherence is the main limiter of temporal resolution and robustness of dv/v measurements. Typically, short aperture arrays or single-station measurements can yield high temporal resolution (e.g., minutes to hours) at high frequency signals (e.g., Xie et al., 2023) while greater array aperture and lower frequency signals tend to become stable over longer time periods (e.g., days to months) (e.g., Mao, Ellsworth, et al., 2025; Clements and Denolle, 2023). Improving the quality and temporal resolution of cross-correlation is an area of research (A. M. Baig et al., 2009; Hadziioannou, Larose, A. Baig, et al., 2011; Viens and Van Houtte, 2019; Yang et al., 2023; He et al., 2024). Retrieving dynamics at time scales shorter than minutes, ambient field seismology relies on the stationarity of the ambient field. It tends to be limited to single-channel analysis or near-sensor processes.
Coda waves multiply scatter in the subsurface and thus have increased, though still weak, sensitivity to perturbations in the medium (Colombi et al., 2014). Given that the sensitivity is weak, observed dv/v values typically under-predict the perturbation at depth (Obermann, Planès, Larose, et al., 2013; Obermann, Planès, Hadziioannou, et al., 2016; Barajas, Margerin, et al., 2022; Yuan et al., 2021). Furthermore, dv/v is measured as decreasing with lag time in the coda (e.g., Gassenmeier et al., 2014; Takano, Brenguier, et al., 2019), which is also expected from the low sensitivity of body waves to perturbations (Obermann, Planès, Larose, et al., 2013; Barajas, Margerin, et al., 2022).
dv/v can be measured in both the frequency and time domains. The moving window cross-spectral (MWCS) technique developed by Poupinet et al. (1984) and the stretching technique developed by Lobkis and Weaver (2003) and first used in an environmental context by Sens-Schönfelder and Wegler (2006) are the most commonly used methods in the frequency and time domains, respectively. The stretching technique works well at low SNR levels (Hadziioannou, Larose, Coutant, et al., 2009), whereas the MWCS is computationally faster to calculate (Clarke et al., 2011). Recently Okubo et al. (2024) found that the stretching and MWCS methods may yield different results depending on which scattered wave dominates in the coda. Additionally, Mikesell et al. (2015) developed a method based on dynamic time warping to measure local phase shift measurements in the coda, Mao, Mordret, et al. (2019) developed a wavelet-based method to measure phase shifts, and Yuan et al. (2021) suggested combining the wavelet and stretching techniques to improve the frequency-dependence of the measurements relative to all previously stated methods. Overall, studies have attempted to measure dv/v over a range of seismic frequencies, which is desirable in the attempt to estimate the depth of the changes.
4.4. Changes in seismic attenuation
Attenuation is highly sensitive to water content, but is always comprised of two effects, intrinsic attenuation and scattering attenuation, which are often difficult to disentangle in observations. A first method for measuring attenuation as related to fluids estimates the amplitude decay as a function of distance in the frequency domain through direct single wave measurements (e.g., Hudson et al., 2023) or using the exponential decay in the spectral autocorrelation method (e.g., Aki, 1957; Prieto et al., 2009; Nakahara, 2012; Lawrence et al., 2013). Another class of methods proposed by Malagnini et al. (2022) utilizes a generalized inverse technique for regression among pairs of earthquakes and stations. A second method for attenuation measurements focuses on the waveform envelope shape for which the radiative transfer theory (Sato et al., 2012; Mayor et al., 2014) is frequently used to model the dissipation of seismic energy and incorporate explicit functional forms for intrinsic and scattering attenuation, which may be necessary to extract the seismic signature of fluids (e.g., Hirose, Nakahara, et al., 2019; Hirose, Ueda, et al., 2022; Hirose, Q.-Y. Wang, et al., 2023; van Laaten and Wegler, 2024). Overall, measurements are made with absolute values of attenuation and evaluated over time. Recently, Pinzon-Rincon et al. (2025) observed 5% variations of seismic attenuation due to changes in shallow water storage in the first tens of meters of the subsurface (12–23 Hz).
4.5. Changes in ground motion site response
The response of the near-surface to incoming earthquake ground motions is referred to as “site response” or “site effects”, and they often include linear effects of wave propagation in the shallow soil layers and the damage and nonlinear deformation of these materials when ground motion amplitudes are sufficiently high. Site effects arise due to the presence of impedance (the product of wave velocity and density) contrasts and subsurface topographies. Site effects also depend on the attenuation of waves in the near-surface environments. Modeling site effects, therefore, requires in-depth knowledge of the site’s velocity structure and the various factors controlling attenuation.
Site response can be precisely measured using instrument-based methods, such as spectral ratio (SR) techniques (Borcherdt, 1970), the generalized inversion technique (GIT, Andrews, 1986), and residual analysis (Abrahamson and Youngs, 1992), within ground motion models (GMMs). SR methods require a reference-target recording pair, such as a borehole station or a station on subsurface reference rock. In contrast, GIT and GMM residual analysis depend on a regional seismic network that includes the site of interest and has recorded multiple past earthquakes with comprehensive azimuthal and distance coverage. SR methods have also been utilized for measuring the environmental impacts on slope stability (e.g., Borgeat et al., 2024).
Where direct measurement or modeling is not possible, amplification can be predicted by generic prediction models linking amplifications to various predictor variables. Predictive variables can be site parameters from in situ geotechnical/geophysical measurements, e.g., shear wave velocity in the shallowest 30 m (VS30), sediment thickness, and site resonant frequency (e.g., Hassani and Atkinson, 2017; Kwak and Seyhan, 2020) or indirect data from existing regional models or maps, such as topographic slope, terrain and geology (e.g., Wald and Allen, 2007; G. A. Weatherill et al., 2020). These indirect methods have recently shown good performance at predicting site response due to the emergence of quality databases and the associated use of machine learning methods (C. Zhu, G. Weatherill, et al., 2021; C. Zhu, Cotton, Kawase and Nakano, 2023; C. Zhu, Cotton, Kawase and Bradley, 2023).
Another metric of site response is the horizontal-to-vertical spectral ratio (HVSR). HVSR identifies fundamental resonance frequencies of soil layers, aiding in seismic microzonation, site amplification studies, and subsurface characterization (Nakamura, 1989). It is particularly valuable for evaluating the thickness and stiffness of sediments, seismic hazard assessment, and earthquake-resistant design (see review by Molnar et al., 2018). The influence of environmental variations on site responses is just beginning to be explored (e.g., Vassallo et al., 2022; Kula et al., 2018; Colombero et al., 2018). Händel et al. (2025) have also identified seasonal variations of high-frequency ground motion characteristics in northeastern Hokkaido and Honshu.
5. Seismic monitoring of hydrological resources and hazards
5.1. Vadose-zone dynamics
Seismological methods offer valuable insights into the partially saturated vadose zone, where complex water dynamics, including infiltration, evapotranspiration, drainage, and runoff, interact with the mechanical properties of near-surface soils and rocks. These processes have been modeled using a range of hydrological and poroelastic frameworks, and increasingly constrained by seismic observations of velocity changes (dv/v) that respond to water saturation and stress conditions.
Several studies have used seismological data to model hydrological processes. For instance, Fores et al. (2018) modeled ±0.2% dv/v variations in the 6–8 Hz frequency band by integrating a water budget, including rainfall, evapotranspiration, and drainage, with gravity data and Biot-Gassmann-based petrophysical models, effectively tracking changes in water saturation in a karstic setting. Similarly, Illien, Andermann, et al. (2021) developed a two-layer hydrological model incorporating an antecedent precipitation index and rainfall thresholds to explain ±2% dv/v variations in the 4–8 Hz frequency range driven by transpiration and runoff. Their model effectively captured responses to heavy monsoonal precipitation in Nepal. In permafrost settings, Sobolevskaia et al. (2024) and Sun et al. (2025) analyzed dv/v in terms of thermoelastic, poroelastic, and saturation-driven processes caused by thawing, precipitation, and air temperature variations.
Vadose zone seismology has been enabled by both long-term sparse networks and short-term dense deployments. The Crépieux–Charmy water exploitation site in France exemplifies a long-term test bed where ambient noise seismology has been combined with hydrological observations using permanent and temporary densified arrays, as well as piezo-electric sensors (Voisin, Guzmán, et al., 2017; Garambois et al., 2019; Gaubert-Bastide et al., 2022). The Susquehanna Shale Hills Critical Zone Observatory in the United States has supported year-round seismic monitoring to detect shallow responses to meteorological and hydrological forcing (Oakley et al., 2021). More recently, dense distributed acoustic sensing (DAS) deployments have enabled higher spatial resolution. For example, several studies used DAS to detect sub-seasonal dv/v decreases following precipitation and long-term increases associated with drying, with a depth sensitivity down to 150 m in semi-arid eastern California and in permafrost regions of Alaska (Z. Shen et al., 2024; Sobolevskaia et al., 2024; Sun et al., 2025).
Interpreting vadose zone dynamics requires careful attention to the differential sensitivity of seismic waves to fluid and matrix properties. S waves are particularly responsive to changes in effective stress and soil stiffness. Garambois et al. (2019) found an approximately linear relationship of a 1% drop in S-wave dv/v and 1 m increase in water table height at frequencies above 5 Hz, although with deviations from that linear relation with a hysteresis observed during wetting and drying (e.g., Garambois et al., 2019; Gaubert-Bastide et al., 2022; J. Shen and T. Zhu, 2025) can be explained with capillary effects in the partially saturated zone (Dong and N. Lu, 2016; Solazzi et al., 2021). In contrast, P-wave velocities are more sensitive to changes in fluid compressibility; in the same study, increases in P-wave velocity were attributed to greater saturation and reduced pore compressibility. Shallow monitoring studies, such as Oakley et al. (2021) and Sun et al. (2025), further highlighted that near-surface P-wave velocities may also respond to thermoelastic effects caused by temperature fluctuations and moisture changes in the top few meters. DAS technologies, which primarily sense strain rates related to mixed wave modes, offer a promising avenue to observe these joint sensitivities over depth and time (Z. Shen et al., 2024; Sobolevskaia et al., 2024; Sun et al., 2025).
5.2. Groundwater resource management
Groundwater monitoring with dv/v is a developing technique for water resource management at the aquifer scale. There is a strong anti-correlation between dv/v and relative changes in groundwater levels, i.e., when groundwater levels decrease, dv/v increases and vice-versa. Clements and Denolle (2018) showed that an increase of up to 16 m of groundwater level over a ∼20 × 40 km area during the winter of 2004/2005 was commensurate with a 0.15% decrease in dv/v in the 0.5–2.0 Hz frequency range (sensitivity to ∼500 m depth). Spatial changes in dv/v and groundwater level were the largest in the deepest part of the basin. Mao, Lecointre, et al. (2022) compared dv/v and groundwater levels over the ∼40 × 80 km greater Los Angeles coastal aquifers in the 0.5–2.0 Hz frequency band (sensitivity to ∼350 m depth). They measured a long-term increase of up to 0.2% dv/v in the central Los Angeles (LA) basin associated with groundwater pumping from 2000–2020. Whereas the Santa Ana basin with active groundwater replenishment had a 0.2% decrease in dv/v. Almagro Vidal, et al. (2021) measured dv/v in the 0.7–0.9 Hz frequency range (depth sensitivity to ∼200 m) over an area of ∼60 × 60 km in Northern Italy. They found that 0.05% dv/v variations were correlated to 1 × 108 m3 changes in total groundwater storage in a karstic aquifer. At a much smaller spatial scale (∼500 × 500 m area) in the Crépieux–Charmy field site, Gaubert-Bastide et al. (2022) measured a 3% decrease in dv/v in the 2–5 Hz frequency range corresponding to a 3 m increase in groundwater level. Figure 4 shows the progress of increasing the spatial resolution of groundwater monitoring with dv/v for water management.
Using seismological proxies to groundwater levels. Left: Clements and Denolle (2018) first developed a smooth map of groundwater level change in the San Gabriel basin, CA. Spatial extent: ∼20 × 60 km. Middle: Mao, Lecointre, et al. (2022) mapped dv/v at 2 × 2 km with 2D sensitivity kernels in the greater Los Angeles, CA area. Spatial extent: ∼40 × 80 km. Right: Gaubert-Bastide et al. (2022) mapped dv/v at 20 × 20 m in the Crépieux–Charmy field site, Lyon, France. Spatial extent: ∼500 × 500 m.
Beyond measuring relative changes in groundwater levels, seismological studies have attempted to infer groundwater dynamics. For example, Delouche and Stehly (2023) found that deep karstic and confined aquifers had a rapid response to seasonal rain (30–60 days) and a reduction in velocities. Mao, Ellsworth, et al. (2025) found that intense rainfall events during the 2023 winter storms in California replenished shallow aquifers but not the deep aquifers that are the main contributors to Calfornia’s water storage.
5.3. Landslide monitoring
Landslide dynamics are profoundly affected by subsurface hydrology (Greco et al., 2023). Rapid changes in the stress field or the water storage, e.g., due to massive rainfall events or earthquake shaking, are often causes for slope failures (Dai et al., 2002). Because hydrological data alone do not always correlate with landslide movements (Schulz et al., 2018), additional information about the rheology and transient processes, such as heterogeneous soil swelling in shear zones (ibid.) or rheological and hydrological heterogeneity (Tacher et al., 2005), is required. The pioneering work of Mainsant et al. (2012) demonstrated that pore-pressure-driven dv/v decreased 7% in the 10–12 Hz frequency range in the four days precedding a catastrophic landslide. Since then, ambient noise monitoring has become a popular way to monitor continuously and forecast landslides: Le Breton et al. (2021) recently provided a comprehensive review of how ambient noise monitoring has been applied to and benefited landslide research, and in particular, the advantages and challenges in making the dv/v measurements robust at short time scales for monitoring landslides. More recently, researchers have monitored the water content in landslides, tailing dams, and unstable rock slopes with dv/v, especially in response to rainfall events (Liu et al., 2024; Ouellet et al., 2022; Borgeat et al., 2025). Marc et al. (2021) showed that the landslide susceptibility rate in the epicentral zones of large earthquakes followed a similar recovery as dv/v in the year following strong ground motion. Collins et al. (2024) used HVSR to estimate the thickness of the slow-moving Barry Arm landslide. Overall, ambient noise methods provide in-situ, high temporal resolution monitoring of the coupled mechanical and hydrological processes that drive landslides. These ground failures, whether triggered by seismic or hydrological loading, should be viewed as dynamic responses of the critical zone, where stress redistribution and fluid–solid coupling govern stability.
5.4. Seismic hazards
Site effects can significantly impact seismic hazards by amplifying seismic waves and altering soil behavior in the shallow subsurface (i.e., the upper ∼50 m). Under weak ground motions, shallow soil layers typically behave in a linear elastic manner, with amplification primarily dependent on the incoming seismic wave strength. During strong shaking, soft soils may undergo substantial dynamic strains and exhibit non-linear responses (Ishihara, 1996). This non-linearity often results in a relative de-amplification of high-frequency ground motions and a shift of energy toward lower frequencies (Beresnev and Wen, 1996). Numerous studies have highlighted the potential of seismic interferometry to monitor the non-linear behavior of shallow sediments from single seismometers (Bonilla and Ben-Zion, 2020; Hobiger et al., 2014; Hobiger et al., 2016; Viens, Denolle, et al., 2018; Esfahani et al., 2024), surface-borehole arrays (Bonilla, Guéguen, et al., 2019; Nakata and Snieder, 2011; Qin et al., 2020), and DAS arrays (Viens, Bonilla, et al., 2022). dv/v can decrease as much as 50–60% in the few seconds following strong motion at frequencies above 1 Hz (Bonilla, Guéguen, et al., 2019). In extreme cases, phenomena such as liquefaction (Nakata and Snieder, 2011), dynamic compaction (Viens and Delbridge, 2024), and significant variations in soil stiffness and strength (Chandra et al., 2015) have been observed. Recent studies have also highlighted the impact of environmental changes, such as rainfall and freeze-thaw cycles, on the subsurface response to ground motions (Roumelioti et al., 2020; Händel et al., 2025). Importantly, seismic hazard itself, via liquefaction and ground failure, is a hydromechanical outcome of CZ processes, reinforcing the need for an integrated perspective.
5.5. Permafrost
Permafrost, which underlies vast areas of the Northern Hemisphere and other cold regions, is warming at an increasing rate due to climate change. The degradation of permafrost poses significant risks, including the release of greenhouse gases such as methane and carbon dioxide, which further exacerbate global warming, as well as threats to ecosystems, hydrology, and infrastructure (Biskaborn et al., 2019). Conventional monitoring systems include temperature measurements through the Global Terrestrial Network for Permafrost and conventional active geophysics surveys, such as resistivity and ground-penetrating radar (Kneisel et al., 2008). Environmental seismology has demonstrated great potential in monitoring permafrost phenomena, mainly through the dv/v and HVSR methods. dv/v monitoring constrains the timing of seasonal freeze and thaw in the active layer, which comprises the upper layer of soil. Thawing of the active layer in early summer induces a 5–10% drop in dv/v at frequencies of 1–30 Hz, whereas refreezing in the late fall increases dv/v a similar amount (James, Knox, Abbott and Screaton, 2017; James, Knox, Abbott, Panning, et al., 2019; Guillemot, Helmstetter, et al., 2020; R. Steinmann, Hadziioannou, et al., 2021; James, Minsley, et al., 2021). This 5–10% dv/v drop is much less than ∼75% difference in velocity between the unfrozen active layer and the frozen permafrost body (Kneisel et al., 2008), which suggests that most dv/v studies are sensitive to depths deeper than the active layer. Seasonal fluctuations in the thickness of the active layer also lead to a sharp velocity contrast at the active layer-permafrost boundary, which can be monitored with the HVSR on land (Kula et al., 2018; Weber et al., 2018; Köhler and Weidle, 2019) and in submarine environments (Overduin et al., 2015; Angelopoulos et al., 2024). In addition to seasonal variations in dv/v, Albaric et al. (2021) and Lindner et al. (2021) recovered long-term decreases in dv/v in permafrost over 3 and 15 years, respectively, suggesting an expansion and warming of the active layer. DAS arrays deployed in permafrost allow for detailed inversion of permafrost structure using tomographic methods (Ajo-Franklin et al., 2017) and dv/v monitoring of both the timing and spatial extent of the freeze and thaw cycle (Cheng et al., 2022; Sobolevskaia et al., 2024). Recently, Steinmann, Seydoux, et al. (2022) employed an artificial intelligence-based method to identify the pattern of ground freezing from continuous ambient seismic noise in an urban environment, and Sun et al. (2025) utilized a rock-physics model to disentangle the relative contributions of thawing, precipitation, thermoelastic, and poreelastic stresses, each of which may dominate at specific depth ranges and depending on the water table depth (C. Yu et al., 2025).
5.6. Snow load and melt
The effects of snow load and melt on dv/v can be multi-faceted. Guillemot, van Herwijnen, et al. (2021) measured a 2–3% decrease in dv/v between 10–25 Hz after dry snowfall. They interpreted snowfall as adding new mass onto the surface without increasing the rigidity, thereby decreasing seismic velocity. In contrast, vertical loads induced by the weight of snow (and or pre-snow precipitation) have been interpreted with a positive increase in dv/v at frequencies below 10 Hz (Cannata et al., 2017; Makus et al., 2023). Heavy, wet snow has been interpreted as (1) increasing water content in the pores as correlated with snow-water equivalent thickness (Hotovec-Ellis et al., 2014) or (2) increasing pore pressure in fluid-filled cracks (Silver et al., 2007; Taira and Brenguier, 2016).
6. Summary of seismological evidence for ground-water hydrology
Our review highlights the ubiquitous correlation between groundwater levels, dv/v, and other seismic markers. To achieve such a correlation, studies may have had to disentangle the hydrological signature in dv/v from other contributing factors (e.g., tectonic, thermoelastic) using modeled and/or observed data. From the numerous studies mentioned in this review, we collected the values of dv/v, their corresponding dH groundwater level changes, the type of seismic network configuration, and the frequencies (i.e., depth) involved. Figure 5 shows how dv/v has been reported to vary against groundwater levels. The studies are compiled based on comparisons with aquifer or water-table levels. We categorize them based on the experimental setting: inter-station cross-correlations in regional seismic networks, small aperture arrays, or single-station measurements.
Change in seismic velocities as a function of groundwater level change. Values are shown from individual studies that reported water levels, or pore pressure, which we converted to water level change using a hydrostatic pressure p =𝜌 g dH, where 𝜌 = 1000 kg/m3 and g = 9.81 m/s2 the gravity constant. Squares represent values obtained from inter-station correlations from regional seismic networks with station spacing on the order of 10 km. Triangles are measurements made with inter-station correlations of small arrays of apertures of about 1 km. Circle markers are measurements obtained from single-station methods. Gray triangle is a digitally extracted data from Garambois et al. (2019). The color bar represents the middle of the frequency band used in each study. The dashed red line is simply showing dv/v = −GWL without any fitting.
Overall, the simple relation dv/v(%) = −GWL(m) is strikingly accurate. The deviation from this relation is explained by the following: (1) a choice of inter-station correlations with large separation, and (2) the choice of a low-frequency band of analysis. Both settings yield a lower sensitivity of dv/v to perturbations at depth. The remaining variability is likely due to differences in lithology, as studies range from karstic carbonate environments to alluvial fans and large sedimentary basins. The sensitivity of shear waves to sediments is drastically different, especially between sand and clay (Dong and N. Lu, 2016).
7. Future directions and research priorities
7.1. Bridging data gaps with enhanced sensor deployment
Seismology is well-positioned to fill a hydrological data gap between remotely sensed observations and ground-based hydrological sensors, which have a spatial sensitivity of 10–100 m and a depth sensitivity of 1 cm to 100 m. Despite the recent progress in Earth observation, data gaps persist in the spatial distribution of seismic sensors at regional scales. Future studies should prioritize the co-located deployment of seismic and hydrological sensors to improve ground truth data and validate models. Multi-sensor data fusion will likely be a desirable mode of analysis (see Figure 6), such as supplementing critical zone observatories (e.g., CZO and OZCAR Brantley et al., 2017; Gaillardet et al., 2018) with added geophysical monitoring capabilities (e.g., seismic and geodetic networks). This approach will enhance our understanding of subsurface processes by directly correlating seismic observations with hydrological measurements. Emerging technologies, such as distributed acoustic sensing (DAS), could play a pivotal role in addressing this gap, enabling cost-effective, high-density sensor networks.
Critical Zone processes are probed by geophysical measurements through ground base sensors sensitive to the subsurface (seismometers, distributed acoustic sensing) and airborne or satellite measurements that probe the surface or very near surface properties. Art made by Elena Hartley.
Additionally, it should be noted that seismological data tends to be publicly accessible, especially through open-access repositories such as the Earthscope Data Services and the European Integrated Data Archive (EIDA). There is also a growing availability of hydrological data sets, such as those provided by the United States Geological Survey, the Hydro-Climatic Data Network (Lins, 2012), the National Water Information System (U.S. Geological Survey, 2020), NASA’s Earth-Observing Giovanni system, and Copernicus Climate Change Service. Open seismic and hydrological data thus allow for both timely and long-term critical zone science.
7.2. Developing unified physical models
A critical challenge is the lack of a unified theoretical framework for integrating seismic observations with hydrological and geomorphological processes. Bridging the disciplines of continuum mechanics, granular physics, and soil science could lead to comprehensive models that capture the coupled effects of environmental influences. For instance, some studies have focused on nonlinear elasticity as a linear superposition of dilatational strains (e.g., Donaldson et al., 2019) or extensional strains (Okubo et al., 2024), or directly modeling dv/v terms (e.g., Sun et al., 2025), while others have focused on induced stresses (e.g., Fokker et al., 2021), and others use petrophysical models to estimate P and S wavespeed with water content (Solazzi et al., 2021; Z. Shen et al., 2024; Sobolevskaia et al., 2024).
An additional challenge is to separate elastic and poroelastic effects: seismic wavespeeds decrease in unconfined aquifers with increasing saturated but increase in confined aquifers due to the increased hydrological load (Delouche and Stehly, 2023). Furthermore, the coupling of thermoelastic stresses with hydrological processes is likely to play a role, given the large thermoelastic stresses in the critical zone (Sens-Schönfelder and Eulenfeld, 2019; Oakley et al., 2021; Sun et al., 2025). Collaborative efforts among seismologists, soil scientists, and hydrologists will be key to developing a unified physical model.
7.3. Standardizing data processing
Drawing unified models from the diverse literature on ambient noise monitoring for hydrological work is challenged by the individualized data processing scheme. Processing strategies for seismic data are currently study-specific, limiting the generalization and reproducibility of findings. Each study has a given sensor configuration, seismic frequency of interest (e.g., depth sensitivity), and lapse time in the coda (e.g., different wave types and spatial sensitivities for early vs. late coda), yielding a large variability in the proportionality between dv/v and groundwater levels (see Figure 5). Efforts to standardize preprocessing workflows, combined with the adoption of automated feature extraction and machine learning techniques (e.g., recent examples to predict river discharge or soil saturation directly from ambient noise Abi Nader et al. (2023), Cunha Teixeira et al. (2025)), can accelerate the identification of hydrological signatures in seismic data. These advancements will also facilitate the development of transferable methodologies for broader applications in CZ science.
8. Conclusion
This review highlights the emerging role of ambient field seismology in bridging hydrological and geophysical sciences for monitoring the Critical Zone (CZ). Advances in ambient noise interferometry, coda wave analysis, and time-lapse tomography have demonstrated that near-surface seismic velocity and attenuation variations are intimately linked to changes in water saturation, pore pressure, and effective stress.
The range of applications covered in this review, from groundwater monitoring, vadose zone dynamics, slope stability, to seismic hazard assessment, may appear distinct, yet they reflect shared physical mechanisms: water-driven changes in subsurface stress and mechanical properties. Recognizing this, we advocate for a unified conceptual framework in which hydrological processes and geohazards are understood as coupled systems responding to climatatic and tectonic forcings.
Seismic velocity changes do not always behave uniformly with increasing water content: while pore pressure increases in saturated media typically reduce velocity, partial saturation and capillary effects in the vadose zone may increase soil stiffness and seismic speeds. These contrasting effects underscore the importance of site-specific rock physics, soil structure, and frequency-dependent wave sensitivity. The interpretation of seismic observables, whether dv/v or Q, must account for these regimes.
Across applications, the magnitude and timescale of seismic responses vary widely from sub-percent dv/v changes over hundreds of meters and tens of years in deep aquifers to rapid, multi-percent dv/v changes in shallow landslides. While the range of aforementioned spatial and temporal scales underscore the versatility of ambient noise monitoring, there is a need for standardized sensitivity metrics and frameworks to compare observations across hydrogeologic settings.
Fully realizing the potential of ambient field seismology will require transdisciplinary collaboration. Integrating seismic observations with hydrological modeling, soil physics, geotechnical engineering, remote sensing, and climate science is essential to build interpretable, scalable, and actionable monitoring systems. Beyond collaboration, this integration demands co-development of research questions, shared instrumentation strategies, and open data standards across traditionally siloed disciplines.
To translate these insights into practice, we outlined in Section 7 a set of specific future directions: from multi-sensor data fusion and physically informed modeling to standardization and reproducibility efforts. Together, these strategies can help realize the promise of passive seismic methods for understanding and forecasting critical zone dynamics in a changing world.
Acknowledgments
This research was partially funded by the David and Lucile Packard Foundation (MD). This manuscript has been co-authored by LV with number LA-UR-25-20095 by Triad National Security, LLC, under contract with the U.S. Department of Energy (no. 89233218CNA000001). LV was partially supported by the Los Alamos National Laboratory Center for Space and Earth Science Chick Keller Fellowship. TC was partially supported by a USGS Mendenhall fellowship. The authors thank the broader community that has contributed to the discussion with MD about this review, namely Kuan-Fu Feng, Shujuan Mao, Stéphane Garambois, and Christophe Voisin. The authors thank two anonymous reviewers and Stephanie James for reviews that improved the manuscript. Much of the open-access data presented in the publications are hosted by the Earthscope Consortium Data Services through the National Science Foundation’s Seismological Facility for the Advancement of Geoscience (SAGE) Award under Cooperative Agreement EAR-1724509. In European literature, EIDA is the European Integrated Data Archive infrastructure within ORFEUS to provide access to seismic waveform data in European archives.
Declaration of interests
The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.
Appendix
We summarize the values of dv/v and groundwater changes published in the literature in Table 1. The readers will find Python scripts to reproduce Figure 5 on https://github.com/Denolle-Lab/dvv-hydro-CompGeo-review-2024.
Literature review of studies that have measured dv/v and compared with groundwater levels
| Region | Wave type | Network | Freq (Hz) | ΔGWL (m) | Hydro estimate method | dv/v (%) | Reference |
|---|---|---|---|---|---|---|---|
| Southern CA, USA | Coda | Array | 0.2–2 | 16 | Well | 0.10 | Mao, Lecointre, et al. (2022) |
| CA, USA | Coda | Single | 4–8 | 6 | Rainfall | 1.15 | Clements and Denolle (2023) |
| San Gabriel, CA | Coda | Array | 2–5 | 25 | Well | 0.15 | Clements and Denolle (2018) |
| New Zealand | Coda | Small array | 6–8 | 2 | Piezo, seismometers | 2.00 | Voisin, Garambois, et al. (2016) |
| Lyon, France | Coda | Small Array | 2–5 | 2 | Piezo | 2.50 | Gaubert-Bastide et al. (2022) |
| Lyon, France | Coda | Small array | 3–20 | 3 | Piezo | 2.50 | Voisin, Guzmán, et al. (2017) |
| Ketzin, Germany | Coda | Array | 1.4–3 | 0.6 | Well | 0.60 | Gassenmeier et al. (2014) |
| Texas, USA | Coda | Single | 0.01–8 | 25 | Well | 2.00 | Kim and Lekic (2019) |
| Germany | Coda | Array | 0.1–0.8 | 0.1 | Well, water balance | 0.01 | Lecocq et al. (2017) (*) |
| Malta | Coda | Single | 3–10 | 0.5 | Well, rainfall | 1.00 | Laudi et al. (2023) |
| CA, USA | Coda | Small array | 4–15 | 1.5 | River gage | 2.50 | Rodríguez Tribaldos and Ajo-Franklin (2021) |
| Australia | Coda | Small array | 20 | 0.3 | 2.00 | Olivier et al. (2017) | |
| Kyushu, Japan | Coda | Array | 0.1–0.39 | 0.01 | Rainfall | 0.02 | Q.-Y. Wang, Brenguier, et al. (2017) (*) |
| OK, USA | Coda | Array | 0.1–1 | 0.15 | GRACE | 0.02 | Zhang et al. (2023) (***) |
| Greece | Late coda | Single | 0.3–1 | 2–5 | Wells | 0.03 | Delouche and Stehly (2023) |
| Italy | Coda | Array | 0.1–0.9 | 100 | Geodetic inversion | 0.10 | Almagro Vidal, et al. (2021)(**) |
| Germany | Late coda | Small array | 1.5–3 | 0.4 | Well | 0.20 | Gassenmeier et al. (2014) |
| Lyon, FR | Surface waves | Small array | 3–40 | 2.5 | Piezo and wells | 1.50 | Garambois et al. (2019) |
| Mojave, CA | Coda | Array | 0.1–2 | 12 | Well | 0.10 | Tsai (2011) and Meier et al. (2010) |
| Taiwan | Coda | Single | 0.1–0.9 | 7 | Well | 0.40 | Feng, Huang, Hsu, et al. (2021) |
Regions, waves used, array configuration, sensor configuration, frequency band of measurements, groundwater level change ΔGWL (m), dv/v(%) values reported, and references. The “rainfall” method refers to diffusion, either through poroelastic or baseflow models, of rainfall data. “single” refers to single station measurements, “small array” refers to arrays of geophones with less than 1 km in aperture, and “array” refers to regional seismic networks with inter-station spacing of about 10 km. (*) Converted from the hydraulic pressure dH = p/g𝜌, where g = 9.81, 𝜌 = 1000 kg/m3. (**) GWL is computed for 1000 km2 basin. (***) Zhang et al. (2023) finds the following relation: dv/v = −1.68 × 10−3ΔGWL − 9.61 × 10−3.