Physics-Based Relationship for Pore Pressure and Vertical Stress Monitoring Using Seismic Velocity Variations

Abstract

Previous studies examining the relationship between the groundwater table and seismic velocities have been guided by empirical relationships only. Here, we develop a physics-based model relating fluctuations in groundwater table and pore pressure with seismic velocity variations through changes in effective stress. This model justifies the use of seismic velocity variations for monitoring of the pore pressure. Using a subset of the Groningen seismic network, near-surface velocity changes are estimated over a four-year period, using passive image interferometry. The same velocity changes are predicted by applying the newly derived theory to pressure-head recordings. It is demonstrated that the theory provides a close match of the observed seismic velocity changes.

1. Introduction

Seismic waves contain information about the subsurface; for instance, subsurface seismic properties such as shear modulus and density can be derived from observations of wave propagation. In Earth sciences, seismic data are therefore an important source of information. Relevant physical and chemical information can be found in seismic properties and especially in their variations. We distinguish two types of variations: spatial and time-lapse ones. Spatial variations in seismic velocity can indicate layer boundaries, faults or more subtle subsurface heterogeneity, e.g., [1]. Time-lapse variations in seismic velocity can be caused by changes in stress, temperature or composition.

A technique to retrieve time-lapse velocity variations is passive image interferometry [2]. This technique consists of two steps. First, the Green’s function is estimated using cross-correlations of seismic noise recorded at two receivers. Secondly, the cross-correlation result using, e.g., a day’s worth of data, is compared with a reference Green’s function to obtain the relative change in seismic velocity. This technique has been applied to find volcanic precursors [3], to monitor stress changes [4], ensure the safety of civil structures [5], and to monitor groundwater tables [6].

A groundwater table is the depth in the subsurface where the soil or rock is fully saturated with water. It is the border between the unsaturated upper part (vadose zone) and the lower phreatic zone. Over the seasons, the soil moisture in the vadose zone varies due to precipitation and evaporation. At the same time, the groundwater table fluctuates due to drainage and inundation. Changes in both zones affect the loading and pore pressure of the deeper subsurface. Locally, the pore pressure at depth can be measured with a piezometric well. For measurement over a larger region and depth range, recent developments in seismic methods are promising.

In many monitoring studies, the seismic velocity change is empirically linked to other observations, for which the velocity change is used as a measurement. For instance, Clements and Denolle [7] fitted seismic velocity changes linearly to changes in the groundwater table. This is also the case for other groundwater and pore pressure monitoring studies [2,6,7,8,9,10,11,12]. An empirical relationship between seismic velocity and the groundwater table or pore pressure can be useful, but without theory it cannot provide understanding on how seismic velocity changes are caused.

The Groningen region in the Netherlands has been a test bed for monitoring studies, both due to presence of a large seismic network [13] and pronounced changes in the subsurface due to gas extraction [14] and a thick layer of unconsolidated sediments [15]. The large azimuthal coverage of the network was used to test different quality parameters for passive image interferometry [16]. Moreover, seismic interferometry was applied to a string of geophones in the reservoir to measure compaction [17]. In another study, a dense surface network of stations was used to detect velocity changes in refracted waves over a one-month period [18]. Using the same dense array of 417 stations, a novel implementation of passive image interferometry was tested [19]. Following heavy rain, they found a velocity reduction that propagates downward with time. This they explained with effective-pressure diffusion.

This study provides a physics-based model connecting fluctuations in the pore pressure and vertical compressional stress to seismic velocity variations through changes in effective stress. We combine the relationship between induced stress and velocity change, derived from first principles by Tromp and Trampert [20], with wave propagation theory, basic hydrology and geomechanics to relate velocity change directly to fluctuations in the pore pressure and vertical compressional stress. Using data collected in Groningen, the Netherlands, we validate this relationship, and show that it enables us to monitor near surface (i.e., top 500 m) changes in pore pressure using passive image interferometry.

2. Theory

We derive a physics-based relationship between shear-wave velocity change and fluctuations in the groundwater table and pore pressure. First, we rewrite the existing relationship between seismic velocity and induced stress by Tromp and Trampert [20]. Then, we use basic hydrology and geomechanics to model the induced stress. Last, we study the relative contribution of different terms in the final equation. A similar derivation for compressional-wave velocity change can be found in Appendix A.

2.1. Velocity Change Due to Induced Stress

Relationships between seismic velocity and induced stress were derived from first principles by eq. 38 (Tromp and Trampert [20]). They showed that for an induced stress ${\mathbf{T}}^0$ (defined as negative for compression), written in terms of an induced pressure $p^0=-\frac{1}{3}tr\left({\mathbf{T}}^0\right)$ and a symmetric trace-free induced deviatoric stress ${\tau }^0={\mathbf{T}}^0-\frac{1}{3}tr\left({\mathbf{T}}^0\right)\mathbf{I}$, the shear-wave velocity can be written as

$$\begin{array}{cc}\hfill \rho {\tilde{\beta }}^2=(\mu +{\mu }^{\prime }{p}^0)& +\frac{1}{2}(1-{\mu }^{\prime })\widehat{\mathbf{k}}·{\tau }^0·\widehat{\mathbf{k}}\hfill \\ & -\frac{1}{2}(1+{\mu }^{\prime })\widehat{\mathbf{a}}·{\tau }^0·\widehat{\mathbf{a}},\hfill \end{array}$$

(1)

with mass density $\rho$, shear-wave velocity $\beta =\sqrt{\mu /\rho }$ without induced stress, its perturbation $d\beta =\tilde{\beta }-\beta$ due to induced stress, shear modulus $\mu$, its pressure derivative ${\mu }^{\prime }=\partial \mu /\partial p$, direction of propagation $\widehat{\mathbf{k}}$, and direction of motion $\widehat{\mathbf{a}}$. We rewrite this equation to find relative velocity change $d\beta /\beta$ by subtracting $\mu =\rho {\beta }^2$ from both sides of the equation

$$\begin{array}{cc}\hfill \rho ({(\beta +d\beta )}^2-{\beta }^2)={\mu }^{\prime }{p}^0& +\frac{1}{2}(1-{\mu }^{\prime })\widehat{\mathbf{k}}·{\tau }^0·\widehat{\mathbf{k}}\hfill \\ & -\frac{1}{2}(1+{\mu }^{\prime })\widehat{\mathbf{a}}·{\tau }^0·\widehat{\mathbf{a}},\hfill \end{array}$$

(2)

and dividing by $\mu$:

$$\begin{array}{cc}\hfill \frac{d\beta }{\beta }(2+\frac{d\beta }{\beta })=\frac{{\mu }^{\prime }{p}^0}{\mu }& +\frac{1-{\mu }^{\prime }}{2\mu }\widehat{\mathbf{k}}·{\tau }^0·\widehat{\mathbf{k}}\hfill \\ & -\frac{1+{\mu }^{\prime }}{2\mu }\widehat{\mathbf{a}}·{\tau }^0·\widehat{\mathbf{a}}\hfill \end{array}$$

(3)

Since velocity changes are usually small $\left(\left|\frac{d\beta }{\beta }\right|\ll 2\right)$, this relationship can be approximated as

$$\frac{d\beta }{\beta }=\frac{{\mu }^{\prime }{p}^0}{2\mu }+\frac{1-{\mu }^{\prime }}{4\mu }\widehat{\mathbf{k}}·{\tau }^0·\widehat{\mathbf{k}}-\frac{1+{\mu }^{\prime }}{4\mu }\widehat{\mathbf{a}}·{\tau }^0·\widehat{\mathbf{a}}.$$

(4)

2.2. Velocity Change Due to Surface Load and Pore Pressure

Wherever pore pressure plays a role, poroelastic theory states that deformation is proportional to the effective stress, defined by $\tilde{\mathbf{T}}=\mathbf{T}+{\alpha }_{B}u$ p. 32 [21], where ${\alpha }_B$ represents the Biot constant and u is the pore pressure. The solid framework carries the part $\tilde{\mathbf{T}}$ of the total external stress $\mathbf{T}$, while the remaining part, ${\alpha }_{B}u$, is carried by the fluid. An increase in pore pressure essentially reduces the forces on the contact surfaces between the grains (i.e., $\tilde{\mathbf{T}}$ becomes less negative, since (effective) stress is defined as negative for compression). The remaining pore pressure, $(1-{\alpha }_B)u$, is counteracted by internal stresses in the solid. In unconsolidated or weak material, ${\alpha }_B$ is close to 1. We therefore need to formulate Equations (1)–(4) in terms of the effective stresses and the pore pressures that affect them.

We examine changes in effective stress that are induced by (near) surface processes such as fluctuation of the groundwater table or the atmospheric pressure. Such processes cause almost instant changes in the stress, while pore pressure diffusion is highly dependent on the permeability of the layers and viscosity of water. Therefore, we assume induced stress ${\mathbf{T}}^0\left(t\right)$ from surface processes to be independent of depth, while we allow induced pore pressure $u^0(z,t)$ to vary with depth. In a one-dimensional model, this automatically satisfies the zero divergence condition $\nabla ·{\mathbf{T}}^0=0$. For the loading effect, we only consider changes in vertical stress $T_{33}^0\left(t\right)$.

We write the induced effective stress as

$${\tilde{\mathbf{T}}}^0={\mathbf{T}}^0+u^0\mathbf{I}=\left(\begin{array}{ccc}{u}^0& 0& 0\\ 0& u^0& 0\\ 0& 0& u^0+T_{33}^0\end{array}\right),$$

(5)

with corresponding induced effective pressure

$${\tilde{p}}^0=-\frac{1}{3}tr\left({\tilde{\mathbf{T}}}^0\right)=-u^0-\frac{1}{3}{T}_{33}^0,$$

(6)

and induced effective deviatoric stress

$${\tilde{\tau }}^0={\tilde{\mathbf{T}}}^0-\frac{1}{3}tr({\tilde{\mathbf{T}}}^0)\mathbf{I}=-\frac{1}{3}{T}_{33}^0\left(\begin{array}{ccc}1& 0& 0\\ 0& 1& 0\\ 0& 0& -2\end{array}\right).$$

(7)

We combine Equations (4), (6) and (7) to write the velocity change as

$$\begin{array}{cc}\hfill \frac{d\beta }{\beta }& =\frac{{\mu }^{\prime }}{2\mu }(-u^0-\frac{1}{3}{T}_{33}^0)\hfill \\ & -\frac{1}{3}\frac{1-{\mu }^{\prime }}{4\mu }{T}_{33}^0({\widehat{k}}_1^2+{\widehat{k}}_2^2-2{\widehat{k}}_3^2)\hfill \\ & +\frac{1}{3}\frac{1+{\mu }^{\prime }}{4\mu }{T}_{33}^0({\widehat{a}}_1^2+{\widehat{a}}_2^2-2{\widehat{a}}_3^2)\hfill \end{array}$$

(8)

We distinguish three types of shear waves: shear waves with vertical propagation and horizontal motion (${\widehat{k}}_1={\widehat{k}}_2=0$; ${\widehat{k}}_3=1$; ${\widehat{a}}_1^2+{\widehat{a}}_2^2=1$; ${\widehat{a}}_3=0$), shear waves with horizontal propagation and horizontal motion (${\widehat{k}}_1^2+{\widehat{k}}_2^2=1$; ${\widehat{k}}_3=0$; ${\widehat{a}}_1^2+{\widehat{a}}_2^2=1$; ${\widehat{a}}_3=0$), and shear waves with horizontal propagation and vertical motion (${\widehat{k}}_1^2+{\widehat{k}}_2^2=1$; ${\widehat{k}}_3=0$; ${\widehat{a}}_1={\widehat{a}}_2=0$; ${\widehat{a}}_3=1$). Respectively, this simplifies Equation (8) to

$${\left(\frac{d\beta }{\beta }\right)}^{Vertical}=-\frac{{\mu }^{\prime }}{2\mu }{u}^0-\frac{{\mu }^{\prime }-1}{4\mu }{T}_{33}^0,$$

(9)

$${\left(\frac{d\beta }{\beta }\right)}_{SH}^{Horizontal}=-\frac{{\mu }^{\prime }}{2\mu }{u}^0,$$

(10)

and

$${\left(\frac{d\beta }{\beta }\right)}_{SV}^{Horizontal}=-\frac{{\mu }^{\prime }}{2\mu }{u}^0-\frac{{\mu }^{\prime }+1}{4\mu }{T}_{33}^0.$$

(11)

3. Model Validation

To validate the relationships derived in the previous section, we model surface-wave velocity changes based on pressure-head measurements and compare their results to independent measurements obtained with passive image interferometry. From pressure-head measurements, recorded by piezometer B08C0952 [22] in Groningen, the Netherlands (Figure 1), we estimated changes in pore pressure and vertical stress, which are used as inputs to Equations (9)–(11) to model shear-wave velocity changes. Corresponding surface-wave velocity changes were then obtained through surface-wave dispersion modeling. To obtain independent measurements of surface-wave phase-velocity changes, we applied passive image interferometry to seismic noise. To this end, we used seismic data from borehole geophones [23] at 200 m depth, located within a 10 km radius from the piezometer (Figure 1).

3.1. Static Model

The static parameters in Equations (9)–(11) are $\mu$ and ${\mu }^{\prime }$. They are estimated from detailed subsurface models available in the Groningen area. As a starting point, we used a local shear-wave velocity and density model from Kruiver et al. [24] and a compressional-wave velocity model from Romijn [25] at $[X_{RD},Y_{RD}]=[261.96\mathrm{km}, 582.86\mathrm{km}]$, visualized in Figure 2a–c. Note that the chosen model location below receiver G424 very close to the piezometer is a point location, whereas the used surface waves sample the entire region.

Using compressional-wave (a), shear-wave (b) and density models (c) we computed the bulk modulus (d), the shear modulus (e), and the confining pressure (f). A smoothed derivative of the shear modulus with respect to the confining pressure ${\mu }^{\prime }$ is shown in Figure 2g. Note that $\mu$ increases with depth, whereas ${\mu }^{\prime }$ decreases with depth. ${\mu }^{\prime }$ has values above 50 for most of the unconsolidated-sediment depth range (upper 800 m) and much smaller values for the underlying chalk rock.

3.2. Stress Model

The borehole piezometer shown in Figure 1 [22] has registered pressure heads (i.e., the levels to which groundwater rises in a frictionless tube due to pore pressure) at depths $z_i=7.3$, $27.3$, $105.3$, $132.3$, and $170.8$ m. Figure 3a shows the change in pressure head $dh(z_i,t)$ with respect to the average pressure head between 1 January 2017 and 31 December 2019. If there were a high permeability between the measurement depths $z_i$, all measurements $dh(z_i,t)$ should be identical. Figure 3a shows this is not the case for depths $z<105.3$ m. This corresponds with the presence of low-permeability clay layers in this depth range.

From a linear interpolation of the pressure head, we obtained changes in pore pressure $u^0(z,t)={\rho }_{w}gdh(z,t)$ as a function of depth and time. $dh(z_i,t)$ values at the deepest three measurement levels are nearly identical. For that reason, we extrapolated $dh(z,t)$ below depths of 170.8 m as a constant with depth down to 840 m and set $dh$ to zero for $z>840$ m, from where the sediments become consolidated. The resulting $u^0(z,t)$ field is shown in Figure 3b for parameters ${\rho }_w=1000$ kg/m³ and $g=9.8$ m/s².

To estimate the order of magnitude of the induced vertical compressional stress, we assumed that changes in pressure head at the shallowest level are similar to changes in the groundwater table. Hence, we used $T_{33}^0\left(t\right)\approx -\varphi {\rho }_{w}gdh(z=7.3$ m$,t)$. Figure 3c shows the result for porosity $\varphi \sim 0.25$.

3.3. Shear-Wave Velocity Change

To compare the effects of groundwater table loading $T_{33}^0$ and induced pore pressure $u^0$ on the seismic velocity, we studied the orders of magnitude that can be predicted using Equations (9)–(11). Typical values for stress and pore pressure can be found at 1 January 2018, when the induced compressional stress was estimated at $T_{33}^0\approx -1000$ Pa, and the average pore pressure was $u^0\approx 2000$ Pa. Using the information of Figure 2e,g (${\mu }^{\prime }\approx 80$ and $\mu \approx 5\times {10}^8$ Pa) in Equations (9) and (11), we found, for the groundwater loading, a relative increase in shear-wave velocity $d\beta /\beta =-\frac{{\mu }^{\prime }\mp 1}{4\mu }{T}_{33}^0\approx 0.004\%$, and we found a relative decrease in shear-wave velocity for the induced pore pressure $d\beta /\beta =-\frac{{\mu }^{\prime }}{2\mu }{u}^0\approx -0.02\%$. Since the effect of vertical compressional stress is considerably smaller than the pore pressure effect in our Groningen setting, and we have no accurate measurements of $T_{33}^0$, we chose to neglect the effect of the vertical compressional stress. Hence, we simplify Equations (9)–(11) to

$$\left(\frac{d\beta }{\beta }\right)\approx -\frac{{\mu }^{\prime }}{2\mu }{u}^0.$$

(12)

In accordance with Equation (12), we used the information from Figure 2e,g and Figure 3b to construct shear-wave velocity changes as a function of time and depth. The result in Figure 4 forms the basis in the forward modeling of surface-wave velocity changes. The effects of the pore pressure $u^0$ are mainly responsible for the time-lapse variations, while the shear modulus $\mu$ and its pressure derivative ${\mu }^{\prime }$ regulate the amplification of shear-wave velocity changes as a function of depth. For instance, at the interface at 50 m depth, the amplitude of velocity change rapidly decreases due to an increase in the shear modulus (Figure 2e), and a decrease in the pressure derivative of the shear modulus (Figure 2g).

3.4. Surface-Wave Dispersion Forward Modeling

Given our physics-based relationship regarding shear-wave velocity change as a function of depth (simplified as Equation (12)), and our observations of surface-wave phase velocity changes as a function of frequency (Section 3.5), we can relate them through surface-wave dispersion modeling. Assuming that the lateral variations in our region of interest are small, we can use a one-dimensional average static model of the subsurface and use the adjoint technique of Hawkins [26] to compute sensitivity kernels which relate partial derivatives of the change in surface-wave phase velocities (Love and Rayleigh) to the small stress-induced changes in shear-wave velocities. These partial derivatives can be used to estimate the effect on relative changes in surface-wave phase velocities using

$$\frac{dv}{v}\left({\omega }_i\right)=\frac{1}{v\left({\omega }_i\right)}\sum _j\frac{\partial v}{\partial {\beta }_j}\left({\omega }_i\right){\beta }_j\frac{d{\beta }_j}{{\beta }_j},$$

(13)

where $\partial v/\partial {\beta }_j\left({\omega }_i\right)$ is the partial derivative of the surface-wave phase velocity with respect to ${\beta }_j$ at frequency ${\omega }_i$, $d{\beta }_j$ is the small perturbation to ${\beta }_j$, and $v\left({\omega }_i\right)$ is the actual surface-wave phase velocity. While this approach is approximate and discounts the impact of anelasticity and lateral heterogeneity, it is sufficient to demonstrate that the stress-induced changes in shear-wave phase velocities give the correct changes in observed surface-wave velocities.

For the static models introduced in Section 3.1, we computed the sensitivity kernels of fundamental Rayleigh and Love waves for density, compressional-wave velocity and shear-wave velocity. Figure 5 shows the kernels at a frequency of 1 Hz. This figure shows that, in the Groningen setting, fundamental mode Love waves are more sensitive to shallower changes than Rayleigh waves, which peak around 200 m depth.

Figure 6 shows several eigen modes for Love and Rayleigh waves at 1 Hz, indicating that there should be reasonable amounts of Rayleigh-wave energy recorded at a depth of 200 m (depth of the borehole geophones we are going to use below), but little fundamental Love-wave energy. Overtones of both Love and Rayleigh will be recorded if excited by the seismic noise; however, our analysis suggests that the fundamental mode is dominant, at least in ZZ and TT cross-coherence’s (Figure 7), as shown in Figure 8.

Using the adjoint method [26], we modeled surface-wave velocity changes corresponding to the shear-wave velocity changes shown in Figure 4. Figure 9 shows the relative velocity changes for the example frequency of 1 Hz as a function of date (a) and for the example date 31 Aug 2019 as a function of frequency (b). We show relative velocity change of Love waves (black), Rayleigh waves (red), and its Voigt average ${(dv/v)}_{Voigt}=\frac{2}{3}{(dv/v)}_{Rayleigh}+\frac{1}{3}{(dv/v)}_{Love}$ (blue). Love waves show a larger induced velocity change than Rayleigh waves. This can be understood from their higher sensitivity to shallow depths (Figure 5), where the largest shear-wave velocity changes occur (Figure 4). For low frequencies, we used the Voigt average of Rayleigh and Love waves, but higher frequencies (Figure 8) show less energy for Love waves than for Rayleigh waves. Therefore, we used Rayleigh wave velocity change for frequencies $f\ge 1$ Hz, which has a much lower amplitude than the Voigt average.

3.5. Passive Image Interferometry

In an approach similar to Fokker and Ruigrok [16], building upon the work of Sens-Schönfelder and Wegler [2], we retrieved relative changes in the seismic velocity using passive image interferometry, consisting of two main processes. First, the time-varying Green’s functions between two seismic receivers were estimated. Second, these estimates were used to determine the relative changes in arrival time, corresponding to the relative changes in velocity.

To estimate a daily Green’s function (i.e., a one-day average of a Green’s function), we computed the cross-coherence, i.e., the spectrally normalized cross-correlation, [28] of ambient seismic noise, recorded by seismic receivers at $x_A$ and $x_B$:

$$\widehat{H}(x_B,x_A,\omega )=\frac{\widehat{u}(x_B,\omega ){\widehat{u}}^{*}(x_A,\omega )}{|\widehat{u}(x_B,\omega )\left|\right|\widehat{u}(x_A,\omega )|}.$$

(14)

The frequency domain is indicated by a hat and the star denotes a complex conjugation. We stacked cross-coherences calculated from 50% overlapping time windows of 20 min durations, where the first time window ranges from 0:00 to 0:20 UTC, the second from 0:10 to 0:30 UTC, etc. We repeated this procedure for the data of every day between 1 January 2016 and 1 January 2020. The cross-coherences were computed for all component combinations. To speed up computation, we applied rotations to radial (R) and transverse (T) components after cross-correlation. Applying rotation after correlation gives a small yet acceptable error (Appendix B). Figure 7 shows the average cross-coherence between 1 January 2017 and 31 December 2019 as a function of distance between receivers. A different time period was chosen for the reference cross-coherence, since in 2016 quite a lot of data are missing. In this way, the reference cross-coherence contains equal contributions from all seasons, while available data from 2016 can still be compared to a well estimated Green’s function.

We determined velocity changes using the stretching method in the time domain [29]. Relative velocity changes $dv/v=ϵ$ were found at the maximum correlation coefficient $CC\left(ϵ\right)$ between lapse (daily) cross-coherence $H_{lapse}$, stretched in time with factor $(1-ϵ)$, and reference cross-coherence $H_{ref}$,

$$CC\left(ϵ\right)=\frac{{\int }_{t_1}^{t_2}{H}_{\mathrm{lapse}}\left[t(1-ϵ)\right]H_{\mathrm{ref}}\left[t\right]dt}{\sqrt{{\int }_{t_1}^{t_2}{\left(H_{\mathrm{lapse}}\right)}^2\left[t(1-ϵ)\right]dt}\sqrt{{\int }_{t_1}^{t_2}{\left(H_{\mathrm{ref}}\right)}^2\left[t\right]dt}}.$$

(15)

The obtained velocity changes are relative to the reference cross-coherence, here defined as the average cross-coherence between 1 January 2017 and 31 December 2019 (Figure 7). We decided to only use the coda of the cross-coherence (red areas in Figure 7), because the coda is less dependent on changes in illumination, and velocity change causes larger changes in traveltime for late arrivals. This results in more accurate and stable measurements of the velocity change. As time windows (integration boundaries in Equation (15)) for the cross-coherence, we therefore used $(x/300+5)$ s $<\left|t\right|\le 100$ s, where the first term is the distance x between the two receivers divided by the minimum expected propagation velocity. An additional 5 s was added to exclude the direct field with certainty. We applied bandpass filters to the cross-coherences, varying both the center frequency and the frequency span to obtain surface wave velocity changes as a function of frequency.

A disadvantage of the stretching method is spurious velocity changes when coda-wave amplitudes vary [30]. The use of a spectral normalization of the cross-correlation, the cross-coherence, limits amplitude variations, but when the ambient noise amplitudes are weaker than the normalization water level used, the amplitudes of the cross-coherence change, resulting in spurious arrivals. For this reason, we analyzed the velocity change distributions of 78 receiver pairs (Figure 1). Spurious arrivals can easily be spotted by an inconsistent distribution of velocity change. Figure 10 shows an example of an inconsistent distribution of velocity change. These spurious velocity changes can especially be observed during summers and at low frequencies ($f<0.2$ Hz). During these seismically quiet periods, the 4.5 Hz geophones record a significant proportion of instrumental noise, which results in a poor estimate of Green’s function.

Employing coda waves measured on the horizontal components, we cannot distinguish between Rayleigh, Love and body waves, since the coda consists of multiply-scattered waves. To show that surface waves dominate the Groningen noise, we applied multichannel analysis of surface waves [27] to the reference cross-coherences, measured on the radial, vertical and transverse components. Figure 8 shows a strong presence of both Rayleigh and Love waves. The surface-wave dispersion can be seen up until about 1 Hz. The dispersion cannot be discerned at higher frequencies, probably due to the near-surface heterogeneity over the area sampled with the 78 receiver pairs (Figure 1). Since surface waves dominate Green’s function estimates (the cross-coherences), we treat the average velocity change in the coda for lower frequencies as the Voigt average of the Rayleigh and Love wave phase velocity changes, and for higher frequencies as Rayleigh wave phase velocity changes. Figure 11 shows that for the coda of the cross-coherence every horizontal component configuration (i.e., RR, RT, TR, TT) leads to similar velocity changes. This indicates that arrivals of both Love and Rayleigh waves can indeed be measured on the coda of all horizontal components. An average over all horizontal component configurations will lead to a cleaner distribution of velocity change. For the cross-coherence of vertical motion at both receivers, Rayleigh waves dominate. Hence, we treat velocity change obtained with coda waves recorded on the vertical component as phase velocity changes in the Rayleigh wave.

3.6. Model Validation

Relative seismic phase-velocity variations were retrieved for all combinations of the seismic receivers indicated in Figure 1 for frequency bandwidths of [0.3 0.4], [0.7 0.8], [1.0 1.2], and [1.3 1.6] Hz. Figure 12 presents the mean velocity change (black) and its distribution (colored background) of 78 station combinations of one vertical (ZZ; left) and four horizontal-component configurations (i.e., the average velocity change of components RR, RT, TR, TT; right). The observations of surface-wave velocity change (black) match the results of the independent forward model (purple and red) quite well, both in shape and magnitude. For the vertical component and the higher frequencies of the horizontal components, the Rayleigh-wave velocity-change model is shown. For the lower frequencies of the horizontal components, a Voigt average was taken over the Rayleigh and Love velocity-change models (Figure 9).

4. Discussion

To retrieve velocity variations of surface waves, we used the stretching method [29]. This method is based on the assumption that velocity change is homogeneous. This assumption is valid for velocity changes due to fluctuations in the pore pressure, because the causes of these fluctuations, i.e., precipitation and evaporation, are similar for the whole area. Local velocity changes, for which the assumption is not valid, cannot be retrieved using this method.

In a similar study, Fokker and Ruigrok [16] retrieved velocity variations in a region 15 km northwest of the area sampled in this study. Compared to their results, our velocity variations are three times smaller. This discrepancy can likely be explained by differences in depth sensitivity for fundamental modes and the overtones of Rayleigh waves, and the dominant amplitude of motion at 50 and 200 m depths. At a depth of 200 m, used in the present study, motion of the fundamental mode dominate over motion of the overtones (Figure 8). At a depth of 50 m, used by Fokker and Ruigrok [16], motion of the first overtone dominates over the motion of the fundamental mode. In Mordret et al. [19], it is shown that, in the Groningen setting, the first overtone of the Rayleigh wave has a higher sensitivity to velocity changes than the fundamental mode.

Note that, for some frequencies, we observed the first Rayleigh overtone on the RR component (Figure 8), while we used the fundamental mode in the modeling (Figure 9 and Figure 12). The match between model and observation will likely improve when overtones are included.

In our model for velocity change, we excluded the Love wave contribution for frequencies $f\ge 1$ Hz due to the low Love-wave energy (Figure 8) and its relatively small amplitude at 200 m depth (Figure 6). If the Voigt average (Figure 9 blue) was taken at these higher frequencies, the amplitude of the velocity change model would be much higher. The Voigt averaged model for velocity change would overestimate the velocity change observations at higher frequencies (Figure 12f,h).

In our stress model (Section 3.2), we interpolated the change in pressure head between the measurement depths, we extrapolated $dh$ as uniform changes between 170.8 m and 840 m depth, and we assumed no changes below the consolidation interface at 840 m depth. We expected the permeability to decrease at this interface, limiting the changes in pore pressure. However, pore pressure changes at depths below 840 m depth are not relevant for the purpose of this study, since the waves we study have no sensitivity to changes at depths below 500 m (Figure 5).

We observed that some studies found a positive correlation between groundwater fluctuations and velocity changes [6,11], while others found an anti-correlation [2,7,8,31,32]. This can be explained by the different mechanisms presented in Section 2. The contributions of $u^0$ and $T_{33}^0$ in Equations (9) and (11) have opposite effects: an increase in pore pressure results in a decrease in seismic velocity, while an increase in compressional stress ($T_{33}$ becomes more negative) results in an increase in seismic velocity. The balance between the mechanisms depends, on the one hand, on the permeability determining the size of $u^0$, relative to $T_{33}^0$. On the other hand, it depends on the presence of Rayleigh and Love waves, and their responses to changes in S-wave velocities. Both negative and positive correlations with, respectively, pore pressure and surface weight were found by Wang et al. [10], who modeled the pore pressure from precipitation measurements and used measurements of snow depth for surface loading.

This article argues that shear-wave velocity changes are caused by fluctuations of the effective stress through changes in the shear modulus. There is, however, a second mechanism to couple velocity change to induced effective stress. Effective pressure leads to compaction, affecting the density and hence the shear-wave velocity. This effect can be quantified using the density derivative of the shear-wave velocity, $d\beta /\beta =-\frac{1}{2}d\rho /\rho$, and the definition of the Bulk modulus, $d\rho /\rho =p^0/\kappa$. A rise in effective pressure ${\tilde{p}}^0=-u^0-\frac{1}{3}{T}_{33}^0$ would induce a relative change in density $d\rho /\rho =-(u^0+\frac{1}{3}{T}_{33}^0)/\kappa$, resulting in an increase in shear-wave velocity $d\beta /\beta =(u^0+\frac{1}{3}{T}_{33}^0)/\left(2\kappa \right)$. For a typical bulk modulus in the Groningen subsurface, $\kappa \sim 5·{10}^9$ Pa (Figure 2d), this would result in velocity changes that are two to three orders of magnitude smaller than the ones observed. Therefore, we neglected the mechanism of compaction.

This study focuses on seismic velocity changes due to changes in the pore pressure. Here, we address sources for velocity changes that could not have been measured due to our choices for surface waves, frequency bandwidths, temporal resolution and methodology.

In some studies, atmospheric temperature variations are correlated to seismic velocity change, e.g., [33]. The temperature oscillations in the subsurface, due to yearly atmospheric temperature variations, can mathematically be described as a highly damped wave propagating downwards. Most of the effect is located above 20 m. For depths below 20 m, the temperature changes for quartz, with a thermal diffusivity ${\alpha }_q=1.4$ mm²/s, would be limited to 0.1 °C. The surface waves, with frequencies we used, have no sensitivity to changes at depths shallower than 20 m (Figure 5). Furthermore, surface waves are sensitive to changes over a large range of depths, while the wavelength of the thermal “wave” is approximately 25 m, resulting in temperature contributions of a range of seasons. This would not be detectable using surface waves at frequencies lower than 2 Hz.

Velocity changes caused by the Lunar tides were also left out, because we have only one velocity measurement per day, whereas Lunar tides cause two oscillations of velocity change per day. Velocity changes induced by Groningen earthquakes (maximum magnitude 3.4, at 3 km depth) would not be detectable due to the relatively small induced stress in the shallow subsurface and the heterogeneous nature of the velocity changes. Changes in water saturation are not relevant as a potential source either, since these changes can only be observed above the groundwater table, which around Groningen is at a depth of approximately 1 m.

The derived relationship between shear-wave velocity and induced effective stress is validated in the context of surface waves, traveling through the shallow subsurface of the Earth. It provides a new understanding of time-lapse variations in the subsurface. We postulate that it can be directly applied to monitoring studies using shear waves or compressional waves (Appendix A) and may provide new insights in monitoring civil structures.

5. Conclusions

In this study, we developed a theory relating seismic velocity change to fluctuations in the pore pressure and vertical stress. By combining a relationship between seismic velocity and induced stress, derived from first principles [20], with basic hydrology and geomechanics, we derived a physics-based relationship for seismic velocity changes as a function of induced pore pressure, vertical compressional stress and elastic parameters. To validate this relationship, we modeled seismic surface-wave velocity changes, based on measurements of the pressure head, using the newly derived relationship and the adjoint method [26]. Surface-wave velocity changes were independently retrieved by applying passive image interferometry [2] to seismic noise measured in the subsurface of Groningen, the Netherlands. The close match between model and observation shows the validity of the derived theory, and justifies the use of seismic velocity variations for pore pressure monitoring.