Seismic wave speed monitoring approach

Groundwater levels have been shown to affect seismic wave speeds, enabling a wide range of monitoring applications (e.g. Andajani et al., Reference Andajani, Tsuji, Snieder and Ikeda2020; Clements & Denolle, Reference Clements and Denolle2018; Mao et al., Reference Mao, Lecointre, van der Hilst and Campillo2022; Roumelioti et al., Reference Roumelioti, Hollender and Gueguen2020; Sens-Schönfelder & Wegler, Reference Sens-Schönfelder and Wegler2006, Yang et al., Reference Yang, Wang, Yuan and Ge2018). Also, groundwater levels in the Groningen region show this relationship: groundwater levels anti-correlate with seismic wave speeds. Figure 1 exemplifies this anti-correlation using observations of Fokker et al. (Reference Fokker, Ruigrok, Hawkins and Trampert2023) and Grondwatertools (2024). Most studies, however, only link groundwater levels through empirical correlations.

Figure 1. Comparison of seismic wave speeds (a) and pressure head data (b). The changes in seismic wave speeds were computed by Fokker et al. (Reference Fokker, Ruigrok, Hawkins and Trampert2023, fig. 4a, red curve at frequency range [1.4–1.6] Hz) from seismic ambient noise measured in the Groningen subsurface (Dost et al., Reference Dost, Ruigrok and Spetzler2017; KNMI, 1993), while the pressure head data were obtained from Grondwatertools (2024, well-id: B03D0016, filter: 001).

Fokker et al. (Reference Fokker, Ruigrok, Hawkins and Trampert2021) first derived the physics describing how groundwater levels affect seismic wave speeds. Two mechanisms are considered here. Groundwater can both exert a pore pressure and a weight load on the system. Pore water pressure reduces the effective pressure, effectively lowering seismic wave speeds, while the weight load increases the effective pressure and hence seismic wave speeds. Using these theoretical relationships and knowledge of the elastic properties of a material, we can model seismic wave speed variations as a result of groundwater dynamics. Figure 2 gives an example of such a modelling exercise: changes in pressure head from the Groningen subsurface in the Netherlands (Grondwatertools, 2024) were used in combination with elastic parameters from this region to predict changes in seismic wave speed (Fig. 2, purple, red). We note that this forward model has not been fitted but follows from physics-based modelling.

Figure 2. Models and observations of seismic wave speed variations (dv/v) shown by Fokker et al. (Reference Fokker, Ruigrok, Hawkins and Trampert2021, fig. 12e). The observations were obtained using passive image interferometry (Sens-Schönfelder & Wegler, Reference Sens-Schönfelder and Wegler2006) on seismic ambient noise measurements (Dost et al., Reference Dost, Ruigrok and Spetzler2017; KNMI, 1993) between frequencies of 1.0 and 1.2 Hz. The background colours indicate the probability of such a value, while the black curve shows the mean observation. Models of wave speed changes are based on pressure head measurements (Grondwatertools, 2024) and shown for frequencies of 1.0 (purple) and 1.2 Hz (red).

Seismic wave speeds can independently be measured using interferometric methods, for instance, passive image interferometry (Sens-Schönfelder & Wegler, Reference Sens-Schönfelder and Wegler2006). By comparing seismic ambient noise measurements at different locations, we can retrieve information about the seismic propagation speed in the subsurface. By repeating this measurement, we find time-lapse changes in seismic wave speeds. Figure 2 shows the comparison between modelled and observed seismic wave speed changes in the Groningen region. The seismic wave speeds tend to slightly increase in summer and decrease in winter, corresponding to the prediction based on pressure head data.

Understanding the physics involved, we can leverage seismic wave speed data to infer pore water pressure variations as a function of time, depth and region. Fokker et al. (Reference Fokker, Ruigrok, Hawkins and Trampert2023) showed this possibility for the upper 200 m of the subsurface. Figure 3 shows their result in a four-dimensional visualisation. Interestingly, this model reveals hydrological characteristics, corresponding roughly to the present lithology.

Figure 3. Pore water pressure inference (top) as obtained by Fokker et al. (Reference Fokker, Ruigrok, Hawkins and Trampert2023, fig. 4c) from seismic wave speed measurements using a physics-based inversion scheme. The horizontal and vertical axes show time and depth, while each subfigure shows a map view of pore water pressure changes for seven different subregions. The bottom figure shows maps indicating these seven different subregions in Groningen and the Zegveld area in the Netherlands, our area of interest. The colour coding corresponds with Figs. 6 and 7.

Hydrologically, the shallow subsurface in the Groningen area can be roughly divided into three layers. Within the first 25 m, there exists an unconfined aquifer where changes in pore water pressure directly correlate with fluctuations in the groundwater table. The likelihood of finding clay layers, characterised by limited permeability, is largest within the depth range from 25 down to 75 m. Within such aquitards, pressure diffusion is naturally restricted, resulting in minimal seasonal variations in pore water pressure. Below the clay layers, down to a depth of 200−300 m, lies a confined aquifer. Pore water pressure in this layer is influenced by the groundwater table at recharge points, leading to relatively uniform spatial variability in pore water pressure within this layer. A hydrogeological model of this region is presented in Section 5, Fig. 6.

We observe some of these characteristics in the pore water pressure inference in Fig. 3: relatively high spatial variability in the shallow unconfined aquifer, smaller values for the aquitard depth range and almost uniform pore water pressure variations in the confined aquifer.

Most monitoring studies employing seismic wave speed variations focus on the saturated subsurface. As the unsaturated subsurface facilitates water exchange between surface and subsurface water (Vereecken et al., Reference Vereecken, Huisman, Bogena, Vanderborght, Vrugt and Hopmans2008), it determines the well-being of ecosystems (Rodriguez-Iturbe et al., Reference Rodriguez-Iturbe, D’Odorico, Laio, Ridolfi and Tamea2007) and the recharging of groundwater reserves (Dobriyal et al., Reference Dobriyal, Qureshi, Badola and Hussain2012), hence it may be interesting to study saturation levels using geophysical methods. We anticipate two monitoring approaches using seismic wave speed data: either using seismic waves that are sensitive to pressure changes in the saturated subsurface or using seismic waves that are sensitive directly to changes within the unsaturated layers.

As the presence of water in the unsaturated subsurface increases the total mass, the shallow layers exert a weight load on the full system. It has been shown empirically (Wang et al., Reference Wang, Brenguier, Campillo, Lecointre, Takeda and Aoki2017) and theoretically (Fokker et al., Reference Fokker, Ruigrok, Hawkins and Trampert2021) that a weight load can increase the seismic wave speeds. We propose a saturation monitoring approach using seismic waves that are sensitive to such a loading mechanism. The presence of water in the shallow subsurface affects the pressure at all depths; hence, lower-frequency surface waves, which are generally more sensitive to changes in the deeper subsurface, are good candidates for this approach. We expect this to be feasible and valuable, particularly in areas where groundwater pressure heads respond slowly to precipitation and evaporation as is the case in De Veluwe, The Netherlands (Zaadnoordijk et al., Reference Zaadnoordijk, Bus, Lourens and Berendrecht2019).

It has also been shown that water seepage through the vadose zone directly affects in-situ seismic wave speeds (e.g. Blazevic et al., Reference Blazevic, Bodet, Pasquet, Linde, Jougnot and Longuevergne2020). To assess where and how changes in the vadose zone can be monitored using seismic wave speed variations, all properties affecting seismic wave speeds need to be considered, including capillary stress, adsorptive stress, atmospheric pressure, water temperature, density and phase transitions (Linneman et al., Reference Linneman, Strickland and Mangel2021; Mordret et al., Reference Mordret, Gradon and Brenguier2022). A physics-based analysis can then determine which mechanisms are dominant in the vadose zone, and which properties can potentially be inferred from seismic wave speed measurements. Although laboratory experiments get closer to understanding various aspects of wave speed changes in the unsaturated subsurface (e.g. Smits et al., Reference Smits, Vasconcelos, Willingshofer and Beekman2024; Gaubert et al., Reference Gaubert, Bordes, Garambois, Brito and Voisin2022), the intricate interplay of all physical processes needs further study.

Electrical conductivity monitoring approach

Carpentier et al. (Reference Carpentier, Meekes, Frumau, Verberne, Candela and Koster2024) developed a monitoring approach for SMC and groundwater levels using electrical conductivity data. They established an empirical relationship between electrical conductivity and SMC (Fig. 4) and designed an electromagnetic induction (EMI) based monitoring approach for SMC and groundwater levels. TNO and SoilMasters employed a mobile electromagnetic mapping system to recover soil moisture and groundwater levels in a managed peatland near the city of Gouda. These innovative measurements were further used to spatially predict subsidence and CO₂ emissions. The results matched closely with precipitation and drought patterns, highlighting the validity and potential of the approach.

Figure 4. Empirical relationship between electrical conductivity and volumetric water content, also known as soil moisture content, for various soil types (Carpentier et al., Reference Carpentier, Meekes, Frumau, Verberne, Candela and Koster2024, fig. 6, after Brevik et al., Reference Brevik, Fenton and Lazari2006). Soils represented in this figure are Clarion (fine-loamy, mixed, superactive, mesic Typic Hapludolls), Nicollet (fine-loamy, mixed, superactive, mesic Aquic Hapludolls), Knoke (fine, smectitic (calcareous), mesic Vertic Endoaquolls) and Canisteo (fine-loamy, mixed (calcareous), superactive, mesic Typic Endoaquolls).

They also explain how ‘A cross-domain technology was introduced that uses an electromagnetic induction mapping system, supported by InSAR, GPS, in-situ probes, CO₂ flux data, and improved prior shallow geological models. Shallow EMI is a high-resolution version of this for ultra-shallow application (Altdorff et al., Reference Altdorff, Bechtold, Van der Kruk, Vereecken and Huisman2016). InSAR derived displacements are estimates of the terrain elevation with millimetre resolution and measures subsidence in time-lapse measurements. The combination of land-based and airborne EMI with InSAR enables the surveying of huge patches of land in short time’.

Carpentier et al. (Reference Carpentier, Meekes, Frumau, Verberne, Candela and Koster2024) conducted three time-lapse EMI field experiments on the Zegveld peat observatory site (in March, June and September 2021) and one field test in Cabauw (in March 2021) to ascertain the correlation between SMC and groundwater level. In the end, too much interference by electromagnetic devices that are permanently present at the Cabauw test field for other monitoring studies resulted in the abandonment of this location after the March measurement campaign. The study was continued at the Zegveld pilot site.

’The EMI field data was subsequently confronted with monitoring data that was previously collected for different studies (Boonman et al., Reference Boonman, Hefting, van Huissteden, van den Berg, van Huissteden, Erkens and van der Velde2022): InSAR (Sentinel-1), probe-based in-situ SMC (Royal Dutch Meteorological Institute (KNMI) at Cabauw, National Research Programme on Greenhouse Gas Dynamics in Peatlands and Organic soils (NOBV) at Zegveld), and CO₂ flux (KNMI at Cabauw, NOBV at Zegveld). These datasets were used in data assimilation and calibration procedures to predict SMC.

In-situ located SMC point measurements after calibration correlated well to the average electrical conductivities at those locations during the three timelapses. A solid empirical relation could be established between SMC and electrical conductivity (Fig. 4), allowing for spatial prediction of SMC on the Zegveld field into plausible maps (Fig. 5). Figure 5 shows that at the Zegveld site both the SMC and the electrical conductivity follow the behavior of precipitation and the groundwater level. The usual pattern arose that a short- or long-term precipitation event would take place, then with short delay the groundwater level would rise and after more delay the SMC and the electrical conductivity would rise. For a drought the opposite effect occurred with the same delays: first a period without precipitation, then the groundwater level drops, and subsequently the SMC and the electrical conductivity drop. This is a helpful insight that enables us to not only spatially predict the SMC from electrical conductivity, but also make assessments of how the SMC will react after a wet and dry period’.

Figure 5. Soil moisture content (SMC) time-lapse maps by Carpentier et al. (Reference Carpentier, Meekes, Frumau, Verberne, Candela and Koster2024, fig. 7) from Zegveld in March, June and September 2021 at 0.3 m depth (three panels right) derived from the observed variations in electrical conductivity at 0.3m (3 panels left). X- and Y axis represent GPS coordinates in the Dutch RD coordinate system. The colorbar units in the left three panels represent electrical conductivity in milliSiemens and in the right three panels SMC in percentage.

Strengths of both monitoring techniques

Both groundwater levels and soil moisture content have been retrieved and mapped using geophysical methods. We highlighted two independent monitoring approaches, each having its own advantages and limitations. Other workers in the field have attempted before to combine seismic-based and EM-based measurements for obtaining hydrogeological parameters, for example, Garambois et al. (Reference Garambois, Sénéchal and Hervé Perroud2002) and Vereecken et al. (Reference Vereecken, Huisman, Pachepsky, Montzka, van der Kruk, Bogena, Weihermüller, Herbst, Martinez and Vanderborght2014). They achieved partial success, mainly because the seismic and EM methods were operated independently and the results were combined.

In this paper, we intensify the integration of seismic and EM methods by having both methods constrain each other and therefore reducing the non-uniqueness of the solution. The strength of the seismic approach mostly lies in the ability to obtain passive observations using ambient noise measurements on existing seismic stations. Since all of the Netherlands is relatively close to the North Sea, seismic ambient noise originating from microseisms is always present. This makes the method highly repeatable with minimal financial and logistic effort. The approach however requires a relatively long sensor deployment, in the order of a year, and the physics-based approach requires substantial knowledge of elastic parameters. The obtained depth resolution is rather low, in the order of 10 m near the surface to 50 at 200 m, and can not be enhanced much further at greater depth, because the sensitivity to changes decreases rapidly with depth. However, the approach allows for a wide spatial coverage, making it ideal to compute average changes over a large area. Lateral resolution and vertical resolution in the very shallow subsurface can be enhanced by using denser sensor spacing, provided that higher-frequency noise sources are persistently present in the area of interest. We expect that also the unsaturated subsurface can be monitored using seismic wave speed variations.

The EMI approach allows for much higher resolutions, on the order of 1 m, dependent on the frequency and depth, while the acquisition can be done relatively quickly and easily. One of the limiting factors originates from the distance between the transmitter and the material of interest. Airborne EMI is most sensitive to the shallowest few meters of the subsurface, while the sensitivity of ground-based EMI extends to 50 m depth. However, the area that can be covered on a single acquisition day is significantly larger for the airborne approach. This trade-off makes airborne surveys particularly useful for monitoring the unsaturated subsurface with a large spatial coverage. It is important to keep in mind that other electromagnetic devices can cause interference, making this approach unfeasible in urban environments.

Illustration of an integrated electromagnetic and seismic approach

By using information from both electromagnetic and seismic measurement techniques, we can construct a more complete picture of groundwater dynamics. As the EMI method excels in mapping clay content, this can be used as prior information for the seismic approach: pressure head variations need only to be inferred in aquifers, while variations in aquitards can be set to zero a priori. The soil moisture content observed by EMI can potentially be used to calibrate the seismic approach, so less information is required. Combining these measurement techniques enriches the available information and enables large-scale, higher-resolution inferences.

The electromagnetic and seismic data used in the previous sections stem from two different locations and cannot be combined. At this moment, several EM profiles have been measured in the Groningen area but are not available yet in the form of layered models. To still illustrate the value of combining these methods, however, we constructed layer models of sand and clay content as they could have been measured using electromagnetic surveys. We then used the layer models as prior information to the seismic inversion approach. The hydrogeological model has been collected from TNO–GDN (2024, REGIS-II) in the Groningen region.

The layering model was constructed as follows. Hydrogeological models were collected from TNO–GDN (2024) for multiple locations within the regions of interest. For every formation, we computed the average depth and thickness. We stacked all formations to obtain an average hydrogeological model, indicating the layer type, sand, clay or complex, as a function of depth. All stacked sand-sand or clay-clay interfaces were joined to one layer, after which layers with a thickness of less than 10% of their depths were discarded. The results are shown in Fig. 6, where the colours correspond to the regions in Fig. 3.

Figure 6. Layering model derived from TNO–GDN (2024) as it could have been measured using electromagnetic methods.

Having knowledge of the subsurface layering, we can use this for model discretization and prior information. As derived in Fokker et al. (Reference Fokker, Ruigrok, Hawkins and Trampert2023), surface wave velocity c directly depends on the pore water pressure u:

(1) $${dc \over c}\left(\omega \right)=\int _{0}^{\infty }K_{u}\left(\omega, z\right) du\left(z\right) dz,$$

where K _u represents the pore pressure sensitivity kernel as a function of frequency ω and depth z. Assuming uniform changes in pore water pressure within a layer, we can use Fig. 6 to discretize Equation 1 to

(2) $${\boldsymbol{d}}\left(t_{k}\right)={\boldsymbol{Gm}}\left(t_{k}\right),$$

where $d_{i}\left(t_{k}\right)={dc \over c}\left(\omega _{i},t_{k}\right)$ represents surface-wave phase-velocity change at frequency range ω _i and time range t _k , the forward operator is given by $G_{ij}=\int _{0}^{\infty }K_{u}(\omega _{i},z) B_{j}(z) dz$ with boxcar function B _j (z), yielding 1 for depth z within layer j and 0 otherwise. Model m _j (t _k ) represents uniform pore water pressure change within layer j at time range t _k and is mathematically coupled to pore pressure changes as $du(z,t_{k})=\sum _{j}B_{j}(z)m_{j}(t_{k})$ . To invert for pore water pressure changes within each layer, we employ the explicit least-squares formulation (Tarantola, Reference Tarantola2005),

(3) $$\tilde{\boldsymbol{m}}\left(t_{k}\right)=\left({\boldsymbol{G}^{T}}{\boldsymbol{C}_{\boldsymbol{d}}^{-1}}\left(t_{k}\right)\boldsymbol{G}+\boldsymbol{C}_{\boldsymbol{m}}^{-1}\right)^{-1}\boldsymbol{G}^{T}\boldsymbol{C}_{\boldsymbol{d}}^{-1}\left(t_{k}\right)\boldsymbol{d}\left(t_{k}\right),$$

where C _d is the data covariance, given by the squared standard deviation of velocity change on the diagonal, and C _m represents the prior model covariance, indicating the expected variance per layer. Based on pressure head measurements in the region, we expect seasonal pore water pressure changes with a variance of 10⁶ Pa² (Grondwatertools, 2024; Zaadnoordijk et al., Reference Zaadnoordijk, Bus, Lourens and Berendrecht2019). Within clay layers, however, pressure diffusion is rather slow, hence we expect much smaller changes. We therefore construct the prior model covariance on the diagonal as 10⁶ Pa² for the sand and complex layers, and 10⁴ Pa² for clay layers.

Figure 7 shows the pore pressure variation models inferred with Equation 3 for the hydrogeological model and parametrization shown in Fig. 6. The coloured curves indicate the posterior model, while the black curves show independent piezometric measurements (Grondwatertools, 2024), within the depth range of the shown model and located within the coloured areas in Fig. 3. Overall, there is a good match between the inferred and measured pore pressure changes. The matches in the deeper layers up to 200 m depth show the value of this approach. However, some layers show a mismatch between inference and independent observation. This can possibly be explained by horizontal variability, while the approach works best with little to no horizontal variability. This leaves room for improvement, possibly by reducing the horizontal scale.

Figure 7. Pore water pressure models inferred from surface-wave phase-velocity variations using prior hydrogeological knowledge. The colour coding corresponds to the locations shown in Fig. 3 and the model parametrization in Fig. 6. The coloured curves represent the inferred pore water pressure models, while the black curves show independent piezometric measurements within the same region and depth range.

The inferred pore pressure models show variations within sand and complex layers, while only small changes are present in the clay layers. The resolvability of pore pressure changes in sand layers has slightly increased by the explained approach but at the cost of low resolvability in clay layers. We deem this reasonable, as pore water pressure variations are expected to be small to non-existent in clay layers.

We note that the sensitivity of seismic wave speeds to changes in pore water pressure still decreases with depth, limiting the resolvability of seismic surface waves. This implies that relatively thinner layers cannot be resolved, especially in the deeper subsurface.

Last, we need to note that the presented hydrogeological model is currently inspired by electromagnetic data. We hope to use real EM data in the near future.