4D Physics‐Based Pore Pressure Monitoring Using Passive Image Interferometry

This study introduces a technique for four‐dimensional pore pressure monitoring using passive image interferometry. Surface‐wave velocity changes as a function of frequency are directly linked to depth variations of pore pressure changes through sensitivity kernels. We demonstrate that these kernels can be used to invert time‐lapse seismic velocity changes, retrieved with passive image interferometry, for hydrological pore pressure variations as a function of time, depth, and region. This new approach is applied in the Groningen region of the Netherlands. We show good recovery of pore pressure variations in the upper 200 m of the subsurface from passive seismic velocity observations. This depth range is primarily limited by the reliable frequency range of the seismic data.

Recently, Fokker et al. (2021) provided a physical model for pore pressure monitoring using surface-wave phase-velocity changes. Building on the theory of Tromp and Trampert (2018), they showed that pore pressure changes induce shear-wave velocity variations through changes in effective stress. Using surface-wave dispersion modeling (Hawkins, 2018), they showed that pore pressure changes explain the measured phase-velocity changes both in phase and amplitude.
In the current study, we demonstrate that measured velocity changes can be inverted for pore pressure variations as a function of time and space. We introduce pore pressure sensitivity kernels for surface-wave phase-velocity changes, and compute velocity variations by applying passive image interferometry to seismic ambient noise measurements in Groningen, The Netherlands. An inversion of these velocity changes results in models of pore pressure variation as a function of time, depth and region. Different regions of Groningen show a different temporal behavior that coincide with the jurisdictions of two independent water boards.

Groningen Setting, Data and Models
The Groningen region in the Netherlands has been studied extensively in the context of induced seismicity (e.g., Bourne et al., 2018;Hettema et al., 2017;Nepveu et al., 2016;Trampert et al., 2022) and subsidence (e.g., Van der Wal & Van Eijs, 2016;Van Thienen-Visser & Fokker, 2017;Van Thienen-Visser et al., 2015). The installation of a large dense network of borehole geophones  enabled intensive research activity. Seismic measurements on multiple depth levels were used to estimate shallow 1D velocity and attenuation profiles (Hofman et al., 2017;Ruigrok et al., 2022) and to estimate soil amplifications (Van Ginkel et al., 2019), while the large azimuthal coverage of the network was used to test different quality assessment parameters for passive image interferometry (Fokker & Ruigrok, 2019). The great amount of geological and geophysical models, provided by previous studies, and the presence of the large seismic network make Groningen an ideal region to test our approach of physics-based pore pressure monitoring.
The Groningen region can be divided into water board Noorderzijlvest in the northwest and water board Hunze en Aa's in the southeast. The borders between different water boards are shown in Figure 1 in light blue. Different water boards in the Netherlands can have different policies regarding groundwater management, and thus the pore pressure variations may be region dependent. In the southeastern region, at the location shown in Figure 1 as the blue dot, a deep borehole piezometer (Dinoloket, 2022) takes direct continuous measurements of the pore pressure at multiple depth levels up to 170 m. Shallow direct measurements of pore pressure variation can be found throughout the whole region (Grondwatertools, 2022).  (Dinoloket, 2022). Different regions are indicated by circles. The color coding is used in Figure 4, Figures S2, S4-S7 in Supporting Information S1 to distinguish regional results. The outline of the Netherlands and the Groningen gas field are shown as black and red lines, while the borders between different water boards are shown in light blue.
Hydrologically, we can classify the shallow subsurface in the Groningen area roughly into three layers ( Figure S1 in Supporting Information S1). In the first 25 m we find an unconfined aquifer. Pore pressure variations within this layer are a direct result of changes in the groundwater table. From 25 m to roughly 75-100 m, we find an aquitard, spanning the entire region with only sparse openings. Due to the low permeability of this clay layer, pore pressure diffusion cannot fully penetrate this layer and hence we do not expect large seasonal pore pressure variations. A confined aquifer can be found below the clay from 75-100 m to 200-300 m depth. The pore pressure in this layer is determined by the groundwater table at the recharge locations. Therefore, the spatial pore pressure variability is expected to be small within this layer.
From the seismic network in Groningen  we use data from the 4.5 Hz borehole geophones at 200 m depth at the locations shown in Figure 1 by the black triangles. We chose the deepest geophones from the borehole network, because they register the highest power of coherent noise from distant sources, compared to the power of incoherent noise from close sources. Each colored circle indicates a subregion that we inves tigate. For each subregion we gather shear-wave velocity and density models from Kruiver et al. (2017) and a compressional-wave velocity model from Romijn (2017). From these models we compute all elastic parameters needed in this study ( Figure S2 in Supporting Information S1).
The models for compressional-wave velocity, shear-wave velocity and density (Figures S2a-S2c in Supporting Information S1) allow us to compute the bulk modulus, the shear modulus and the confining pressure (Figures S2d-S2f in Supporting Information S1). The pressure derivative of the shear modulus, needed for the sensitivity kernel, can be computed by a pointwise derivative of the shear modulus with respect to the confining pressure. At layer interfaces, however, the shear modulus can change abruptly due to a change in material from one layer to another. This will result in an unrealistic estimate for its pressure derivative. A smoothing operation with a robust weighing function and positivity constraint removes outliers that occur at such a layer intersection. Figure S2g in Supporting Information S1 shows our model for the pressure derivative of the shear modulus dμ/dp at the center of the corresponding region.

Passive Image Interferometry
To compute seismic velocity changes we apply passive image interferometry (Sens-Schönfelder & Wegler, 2006) to seismic ambient noise measured in Groningen, The Netherlands. This method consists of two processes. First, the Green's function between two seismic receivers is estimated using cross-correlations of ambient seismic noise. Second, time-lapse variations in arrival times are identified, corresponding to velocity variations.
To estimate the Green's function for one lapse period, we compute the cross-coherence of seismic noise, recorded by seismic receivers at locations x A and x B . The cross-coherence represents the spectrally normalized cross-correlation, and can be computed in the frequency domain (Wapenaar, Slob, et al., 2010): where u is ground velocity. The frequency domain is indicated by a hat and the star denotes a complex conjugation. We stack cross-coherences calculated from 50% overlapping time windows of 20 min duration, where the first time window ranges from 0:00 to 0:20 UTC, the second from 0:10 to 0:30 UTC, etc., for a lapse period of 21 days. We repeat this procedure for lapse periods between 01 January 2017 and 01 January 2020. The cross-coherences are computed for vertical components. Figure S3 in Supporting Information S1 shows an example of cross-coherences in the time domain as a function of date, for receiver combination G014-G104 in the orange region ( Figure 1) and frequency range [1.3 1.6] Hz.
We then determine velocity changes using the stretching method in the time domain (Lobkis & Weaver, 2003). Relative velocity changes dv/v = ϵ are found at the maximum correlation coefficient CC(ϵ) between lapse cross-coherence H lapse , stretched in time with factor (1 − ϵ), and reference cross-coherence H ref , The reference cross-coherence is defined as the 3-year average from 01 January 2017 0:00 UTC to 01 January 2020 0:00 UTC, hence the retrieved velocity change is relative to the average within this period.
The coda of the cross-correlation is more likely to contain stable parts of the Green's function, because this only requires a stable background noise structure (Hadziioannou et al., 2009), while direct waves also require well-illuminated Fresnel zones (Wapenaar, Draganov, et al., 2010). For this reason, we omit all arrivals of direct waves, and choose our time windows (integration boundaries in Equation 2) for the cross-coherence as τ < |t| < 2τ, where τ = (x/v low + 5) s. v low is the fundamental-mode Rayleigh wave phase velocity in the model of Figures S2a-S2c in Supporting Information S1. An additional 5 s is added to exclude the direct Rayleigh waves with more certainty. This narrow window excludes most body waves in the coda and should mainly leave closely scattered surface waves.
We filter the cross-coherences with a bandpass filter before we estimate the velocity change for the chosen frequency range. To obtain velocity variations as a function of frequency range, we repeat this process for multiple frequency ranges. We compute an average velocity change for the regions indicated by the circles in Figure 1, using all receivers pairs within the indicated circles. This also allows us to compute the standard deviation of the sampling distribution of velocity change ∕ = ∕ √ , as an indication of the measurement uncertainty on the one hand, and the intrinsic variability over a region on the other hand.
We use the coda of the cross-coherence evaluated for the vertical components to estimate velocity changes. Likely, the velocity changes are caused by fundamental-mode Rayleigh waves, but contributions from higher modes, Love and body waves cannot a-priori be excluded. We repeat the approach of Fokker et al. (2021) to find what type of waves is the main contributor to the observed velocity change by making a forward calculation for the region containing the piezometer. Figure S4 in Supporting Information S1 shows velocity changes for five frequency ranges, retrieved using passive image interferometry (purple), and fundamental-mode phase-velocity changes for Rayleigh (red dashed) and Love (blue dashed) waves, modeled from the pore pressure variations measured by Dinoloket (2022). The velocity variations closely resemble fundamental-mode Rayleigh-wave velocity changes. Therefore, we treat the velocity changes measured on the vertical components as fundamental-mode Rayleigh-wave phase-velocity changes. We tried the same modeling with a Voigt average of Love and Rayleigh (Fokker et al., 2021), but this degraded the fit to the piezometer data.

Pore Pressure Sensitivity Kernels
To connect Rayleigh-wave phase-velocity change to pore pressure variation, we combine the physics-based relationship derived by Fokker et al. (2021) with shear-wave sensitivity kernels to construct pore pressure sensitivity kernels. Building on Tromp and Trampert (2018), Fokker et al. (2021) derived that a change in pore pressure u 0 via effective stress induces shear-wave velocity change with shear-wave velocity β, shear modulus μ, and pressure derivative of the shear modulus μ′ = dμ/dp. A positive change in pore pressure thus results in a negative change in shear-wave velocity.
Changes in the shear-wave velocity directly induce Rayleigh-wave phase-velocity changes with Rayleigh-wave phase velocity v, and shear-wave sensitivity kernel K β . We can now substitute Equations 3 in 4, resulting in represents the pore pressure sensitivity kernel for Rayleigh-wave phase velocity.
Shear-wave sensitivity kernels for Rayleigh-wave phase velocity can be calculated using the adjoint method (Hawkins, 2018) together with one-dimensional models for compressional-wave velocity v p , shear-wave velocity v s , and density ρ. Figure 2a shows the shear-wave sensitivity kernel for the region centered at receiver G424 (purple region in Figure 1), constructed from the elastic model shown in Figures S2a-S2c in Supporting Information S1. The fraction −μ′/2μ shown in Figure 2b is calculated using the shear modulus and its pressure derivative (Figures S2e and S2g in Supporting Information S1). In accordance with Equation 6, we multiply Figures 2a and 2b to obtain the pore pressure sensitivity kernel shown in Figure 2c.

Inversion for Pore Pressure Variation
To invert surface-wave velocity change for pore pressure variation as a function of depth and time, we need to discretize the linear relation described by Equation 5. We expand pore pressure change u 0 as where function S j (z) is chosen to be a cubic natural spline function, and m j (t k ) its coefficients at time t k , which is the center of the 21 day lapse period (Section 3). We then rewrite Equation 5 as For each lapse time t k , this can be written as a linear forward problem, represents the data, the forward operator, and m j (t k ) the model coefficients of the pore pressure change. Model coefficients m j (t k ) can be retrieved using the explicit least-squares formulation (Tarantola, 2005), with data covariance C d and prior model covariance C m . Based on the pressure head measurements in the southeastern region we expect a variance in pore pressure of 10 6 Pa 2 , hence we choose the model covariance as C m = 10 6 I, where I represents the identity matrix. Since we are interested in the mean velocity change dv/v(ω i , t k ) per region, we define the data covariance as the variance in the set of cross-coherences per region (see Figure 3a, error bars). We note that this variance can reflect the cross-coherence variability per region and/or direct observational uncertainty. We therefore use The resolution R(t k ) of the inverted model representation ̃( ) can be obtained by substituting the data d in Equation 12 for the forward operator G, and the posterior model covariance can be found bỹ After inversion for model representation m j (t k ), we repeat the process for all lapse times t k , and compute our final model for pore pressure variation using Equation 7. Figure 3 shows the steps in the inversion scheme for the region centered at receiver G424 (purple region in Figure 1). Velocity changes retrieved using passive image interferometry form the data of this inversion (Figure 3a, error bars; two example frequency ranges). We use velocity variations of multiple frequency ranges with varying center frequency and frequency span (Figure 3b), and we define 10 spline functions S j (Figure 3c). Following Equation 11, we construct forward operator G ij (Figure 3d). Figure 3e shows pore pressure variations as retrieved using Equations 7 and 12, and Figure 3f shows the posterior model covariance as computed using Equation 15. The uncertainty of the retrieved model can then be computed using the square root of the diagonal of the posterior model covariance. Pore pressure changes smaller than this uncertainty are colored gray in Figure 3e. The resolution matrix is computed using Equation 14 (Figure 3g), indicating that we only have sufficient resolution to confidently infer the model coefficients corresponding to the first six splines. Therefore, pore pressure variations can only be retrieved at depths smaller than about 200 m. The resolution matrix shows that deeper pore pressure models have contributions from splines 2 and 6-10, and are thus smeared out over a large depth range. To show how well the pore pressure model explains the velocity variations, we use Equation 9, the forward operator G, and the inferred pore pressure model ̃ to predict the data. Figure 3a (solid lines) shows the result.
We construct a four-dimensional pore pressure model by repeating the inversion procedure for all regions shown in Figure 1. We compute velocity changes (Figure 4a shows five example frequencies) and construct pore pressure sensitivity kernels based on the elastic parameters shown in Supporting Information S1 ( Figure S2). The inversion leads to pore pressure models as a function of time, depth and region. Figure 4b shows in purple the inferred model in the region centered at receiver G424 for five depths, compared to the independent direct measurements of pore pressure variation in black (Figure 1, blue point; Dinoloket, 2022). The four-dimensional model of pore pressure variations is illustrated in Figure 4c, where for five depth levels and seven dates the pore pressure is shown in a colored map view. Detailed comparisons between pore pressure models and comparisons with shallow independent piezometric measurements are shown in Supporting Information S1 (Figures S5 and  S6). The comparison of shallow pore pressure models in the northwest and the southeast shows significant spatial variations, while lateral variations of deeper pore pressure models could not be classified as significant. The shallow pore pressure models also compare well in phase and amplitude to the direct independent measurements of pore pressure change. The relative misfit between velocity change measured using passive image interferometry and predicted based on the inferred pore pressure model is shown in Supporting Information S1 ( Figure S7), indicating that measured velocity variations between 0.7 and 1.8 Hz are well explained by our pore pressure model. In the lower frequency ranges, that is, larger depths, the model does not explain the data, in agreement with the information displayed on the posterior covariance and resolution matrix.

Hydrologic Interpretation
The inferred pore pressure models reveal the characteristics of the hydrologic classification (Section 2, Figure S1 in Supporting Information S1).
Within the confined aquifer, pore pressure models compare well to the direct measurement in the southeast (Figure 4b) and models for the different regions are very similar to each other ( Figures S5d-S5f in Supporting Information S1). The seasonal trends show lower pore pressures during summers and higher pore pressures during winters. The source for pore pressure change in this lower layer is due to locations where the clay layer is absent or very thin and pore pressure diffusion can reach this aquifer. Therefore, the pore pressure in this aquifer represents groundwater fluctuations at the recharge locations.
Within the aquitard, we observe small pore pressure variations that show neither a clear seasonal pattern, nor consistency over the different regions. Within this layer we expect much smaller pore pressure variations, because the hydraulic conductivity in the order of 1 mm per day is too low for pore pressure diffusion to reach the core of this layer. In the inversion process, pore pressure variations must therefore have leaked from depths corresponding to neighboring splines. The resolution in Figure 3g shows that this is possible.
Within the unconfined aquifer, pore pressure variations are a direct result of the changing groundwater table.
Changes in the groundwater table are very site dependent, since their sources (i.e., precipitation, topography, groundwater extraction, and groundwater management) can vary from region to region. Interestingly, there is a significant ( Figure S5 in Supporting Information S1) difference in amplitude between shallow pore pressure variations in the southeast (purple and blue areas) and the northwest (red and orange areas). Independent shallow piezometric measurements of the pore pressure (Grondwatertools, 2022) show for this aquifer an amplitude increase in seasonal variations from the southeast to the northwest. The amplitude differences between the regions coincide with the jurisdictions of two different water boards that may have different policies for groundwater management. The mismatch between shallow pore pressure models and the direct measurements shown in Figure 4b can potentially be explained by local topography or the presence of clay, since the direct measurements are taken at a point location, while the models represent an average over a lateral area of 250 km 2 . The spatial variability shown by other pore pressure measurements from this region (purple area in Figure S6 in Supporting Information S1) supports this hypothesis. Other shallow pore pressure measurements (Grondwatertools, 2022) show closer agreement with the shallow models ( Figure S6 in Supporting Information S1).

Discussion
In this study we obtained seismic velocity changes using the stretching method (Lobkis & Weaver, 2003). However, Zhan et al. (2013) showed that varying amplitudes in the noise can lead to spurious velocity changes. This is what we observe at frequency ranges containing the frequencies of 0.63 or 1.24 Hz, which are eigenfrequencies of nearby wind turbines (Van der Vleut, 2019). With varying wind direction, the swinging direction of the wind-turbine masts changes and therefore the directions, into which Rayleigh and Love waves are excited, will change. This causes substantial amplitude variations and hence spurious velocity changes. For this reason we excluded all frequency ranges containing these eigenfrequencies.
The advantage of the stretching method mostly lays in the ability to detect weak velocity changes using low signal-to-noise ratios. However, it makes use of the assumption of homogeneous velocity change. Using this method we can therefore only retrieve an average velocity change over a relatively large region. Alternatively, one could estimate velocity change using the moving window cross-spectral method (Clarke et al., 2011;James et al., 2017), dynamic time warping (Mikesell et al., 2015), or the wavelet method (Mao et al., 2020). These methods can be used for a higher-resolution spatial inversion of velocity change, taking into account the sensitivities of different wave types at different arrival times and frequencies (James et al., 2019;Mao et al., 2022;Margerin et al., 2016;Obermann et al., 2013).
By using the coda of the cross-correlations of vertical components close after the arrival of the fundamental-mode Rayleigh wave, we excluded most Love-wave energy. If the ratio of Love to Rayleigh energy were known in the Groningen area, one would be able to add velocity change measured on the horizontal components (i.e., RR, RT, TR, TT). The pore pressure sensitivity kernels for Rayleigh and Love would need to be averaged accordingly. A Voigt average between Rayleigh and Love as used by Fokker et al. (2021) would be too rough an approximation for pore pressure inversion, since the ratio of Love to Rayleigh energy varies as a function of frequency (Juretzek & Hadziioannou, 2016).
Velocity changes are linked to pore pressure variations through pore pressure sensitivity kernels. To compute these kernels for Rayleigh-wave velocity change, we determined pressure derivatives of the shear modulus by a point-wise comparison between the shear modulus and the confining pressure. While this is a reliable method to determine the pressure derivative within a layer of one material, at interfaces this can lead to spurious values. A smoothing operation with a weighing function can remove such outliers at the cost of resolution. Alternatively, one could conduct a lab experiment to determine the pressure derivative of the shear modulus as a function of depth and hence maintain a better vertical resolution.
There are unexplained low-frequency data ( Figure S7 in Supporting Information S1). For frequencies below 0.5 Hz we are pushing the 4.5 Hz geophones to their limits. With much instrumental noise at these frequencies, the retrieved velocity variations are of low quality. However, for the inversion part the quality of the low-frequency velocity variations does not really matter, since the resolution shows that the pore pressure models below 200 m cannot be interpreted anyway.
In this study we showed that the velocity variations between 0.7 and 1.8 Hz can be attributed to pore pressure changes. While in Groningen pore pressure change is the main source for velocity variation, other sources also need to be addressed. Locally, earthquakes can cause subsurface damage, resulting in a velocity drop (e.g., Brenguier et al., 2008;Wegler et al., 2009). However, this local effect has only been reported for much larger earthquakes than the ones observed in the Groningen area. Also temperature variations can induce seismic velocity changes (e.g., Colombero et al., 2018;Richter et al., 2014). Seasonal temperature variations by thermal diffusion through quartz, however, are naturally restricted to 0.1°C for depths below 20 m, and thermal energy storage systems only induce local temperature changes that cannot be resolved with our spatial resolution. Moisture variations within the vadose zone cause changes in density that can affect surface-wave velocities (e.g., Knight et al., 1998). In Groningen, however, the groundwater table can be found at approximately 1 m depth, which leaves a very small vadose zone and therefore a limited sensitivity to changes therein. For these reasons, we do not expect that other mechanisms should notably affect the seismic velocity, and therewith the pore pressure models at depths below 20 m.
Within the inversion procedure for depth variations of pore pressure, we used well-defined data and model covariances, enabling the use of the explicit Bayesian formulation. When data or model covariances are not available, it is still possible to carry out a damped least squares inversion. One can search for an optimum weight for the residual norm minimization and the solution norm minimization. Additionally, one could use the correlation coefficient CC max (ω, t) (Equation 2) as proxy for the quality of the retrieved velocity changes, since Fokker and Ruigrok (2019) showed that the standard deviation of retrieved velocity changes σ(ω i , t k ) correlates strongly with 1 − CC max (ω i , t k ). Therefore, this can be used as an alternative to the data covariance presented in this study (Equation 13).

Conclusions
This study introduces a new technique for pore pressure monitoring using passive image interferometry. We derived that pore pressure sensitivity kernels can be used to link surface-wave velocity change as function of frequency directly to pore pressure change as function of depth. In Groningen, The Netherlands, most sensitivity to pore pressure changes lays in the very shallow subsurface (i.e., top 200 m), much shallower than the sensitivity to shear-wave velocity change. We showed that pore pressure sensitivity kernels can be used to invert surfacewave velocity changes for pore pressure variations as a function of depth, resulting in four-dimensional pore pressure models, agreeing with independent measurements of pore pressure variation and showing hydrological features.