Reconciling the Long‐Term Relationship Between Reservoir Pore Pressure Depletion and Compaction in the Groningen Region

The Groningen gas reservoir, situated in the northeast of the Netherlands is western Europe's largest gas reservoir. Due to gas production measurable subsidence and seismicity has been detected across this region, attributed to the deformations induced by reservoir pore pressure depletion. We investigate the surface displacement history using a principal component analysis‐based inversion method to combine a diverse set of optical leveling, interferometric synthetic aperture radar, and Global Positioning System data to better constrain reservoir compaction and subsidence history. The generated compaction model is then used in combination with prior pressure depletion models to determine a reservoir uniaxial compressibility. The best fitting model of uniaxial compressibility is time independent but spatially variable. The absence of evidence for any significant time delay between changes in depletion and compaction rates supports an instantaneous poroelastic reservoir response. The absence of nonlinear yielding at the largest compaction strains suggests that anelastic deformations are a minor part of reservoir compaction.


Introduction
The Groningen gas field, situated in the northeast of the Netherlands (Figure 1), has been in production since 1963, with 70% of the estimated 2,800 billion cubic meters (bcm) of gas now extracted (e.g., Bourne et al., 2014). Prior to gas extraction the region was considered naturally aseismic with no recorded historical events. However, since the 1990s small-magnitude earthquakes have been detected, with these shallow events causing nonstructural building damage and public concern (Dost et al., 2012). The last decade has shown an increase in both the mean frequency and the magnitude of earthquakes. On 17 January 2014 gas production of the Groningen field was limited to 42 bcm per year, in response to the increasing frequency and magnitude of seismicity across the region. This reduction was achieved through reduced extraction from the northern gas extraction sites, in the Loppersum area ( Figure 1). In 2015, the reservoir extraction rate dropped further to 27.5 bcm per year (van Thienen-Visser & Fokker, 2017).
The seismicity is thought to be due to strain induced by the decrease of the bulk reservoir volume Dempsey & Suckale, 2017;Nederlandse Aardolie Maatschappij, NAM, 2016), but the details of the underlying mechanisms remain poorly understood (Spiers et al., 2017). Progressive fluid extraction from this reservoir has lead to a documented decrease in reservoir pore fluid pressure by 25 MPa as of 2017 (Bourne & Oates, 2017). The resulting increase in the effective vertical stress must have driven compaction of the reservoir, due to the combination of linear poroelasticity (Wang, 2000) or to anelastic compaction (Schneider et al., 1996), leading to surface deformation ( Figure 1).
An improved understanding of the geomechanical properties of the reservoir would help reconcile the compaction response to changing pressure depletion. For a purely poroelastic medium, the elastic deformation of the porous material is coupled to fluid flow, as described by Biot's theory of a porous elastic medium saturated with a compressible fluid. In that case, strain at any point within the reservoir is linearly proportional to the local pressure change. In contrast, in the case of an anelastic response of the reservoir, strain at any point within the reservoir would depend on the entire time evolution of the reservoir pressure depletion and would result in a nonlinear relationship between compaction and pressure change. This formulation would be sensitive to the entire time evolution of the reservoir pressure depletion, with a inferred time decay response to reservoir pressure depletion. One way to assess the reservoir properties is to measure surface deformation and model how it relates to pressure variations in the reservoir. The objective of this study is to measure the evolution of reservoir compaction in response to reservoir pore pressures decreasing over the last 50 years, taking advantage of the wealth of data now available, which are documented in the surface deformation.

Geological Setting
The Groningen field is situated in the Groningen structural high, formed during the late Kimmerian tectonic phase. The subsurface is well known from various geophysical investigations (seismic reflection and seismic refraction), borehole core samples, and borehole logging ( Figure 2). The reservoir levels are located within the Upper Rotliegend Group, at depths between 2.6 and 3.2 km, with thickness varying from 100 to 300 m. The Upper Rotliegend comprises the Slochteren sandstone and Ten Boer claystone formation. The gas has migrated from its source in the highly faulted underlying Pennsylvanian Carboniferous limestone containing coal layers (Figure 2a). The overlying thick and impermeable evaporite and anhydrite layers of the Permian Zechstein group provide the seal for the reservoir. This caprock formation comprises a laterally discontinuous 50-m-thick basal anhydrite and 0.2-to 1-km-thick evaporite deposit with an internal disconnected anhydrite layer of 50-m thickness. Seismic reflection, seismic refraction, sonic log, and well core samples have been combined to produce a 3-D elastic model of the Groningen region (NAM, 2013), with knowledge of the reservoir depth and thickness determined from this model. A cross section of the P wave velocity model and a cartoon representation of the separate lithologies are shown in Figures 2a and 2b.

Geodetic Data
Surface deformation monitoring in the Groningen region commenced in 1964 with optical leveling, providing measurements of relative height differences between reference benchmarks, with repeated measurements taken every 1-8 years. Since 1992, surface deformation monitoring has been expanded to incorporate The colors correspond to the separate data types, with blue representing optical leveling, green Global Positioning System and red interferometric synthetic aperture radar. interferometric synthetic aperture radar (InSAR) measurements, measurements along line of sight (LOS). These data have been processed by Skygeo (http://www.skygeo.com/) and Ketelaar et al. (2006), for the selection of coherent persistent scatter points; they have modeled and removed atmospheric signals. (Supporting information section S1.2 discusses this further.) To expand further the temporal resolution of surface deformation, continuous Global Positioning System (cGPS) stations were deployed by NAM (2013) starting in March 2013, allowing us to measure a three-dimensional displacement vector time series. This network was gradually expanded to include 13 stations by April 2014. The temporal span and sampling of each of these techniques is shown in Figure 3.

Optical Leveling
Optical leveling is a surveying technique that provides measurements of height differences between reference benchmarks. We use the data provided by NAM for the period spanning 1964 to 2013, in its raw format (without any processing such as outlier rejections; NAM, 2016). The array geometry of the network benchmarks changes significantly across the period. The initial campaign in 1964 was restricted to sparse coverage of the central and southern portion of the Groningen field, with further densification of the array in 1972 and 1987 covering the entire region (supporting information Figure S1).
We chose a reference benchmark external to the zone of induced subsidence. The data were spatially differenced relative to a benchmark at Rijksdriehoekscoördinaten coordinate RD-Easting = 249080 and RD-Northing = 554740, about 10 km south of the gas field. This benchmark is external to the deforming region, and its measured elevation is stable for the study period. Error analysis of the optical leveling data is still poorly constrained, with no independent measure of error given between each benchmark. The information regarding the acquisition pathways is not available, which limits our ability to determine the data uncertainties. We therefore assumed a standard uncertainty on each benchmark of 1 cm, the minimum expected uncertainty for the data acquisition between benchmarks (NAM, 2013).

PS-InSAR
Persistent scatter interferometric synthetic aperture radar (PS-InSAR, Ferretti et al., 2001) is a technique that constrains motion between subsequent synthetic aperture radar (SAR) acquisitions of high coherence, persistent scatterers. This technique measures the component of motion along the LOS of the satellite, expanding on the only vertical deformation of optical leveling.
The Groningen region has been monitored by radar acquisitions since the first acquisition of European Remote Sensing 1 (ERS-1) in 1992. The first study using such data used satellite acquisitions of ERS-1 and Envisat (Ketelaar et al., 2006). These measurements were validated through comparison with optical leveling acquisitions (NAM, 2016). Since 2009, the InSAR data sets have been augmented substantially by the acquisition of RADARSAT-2 and TerraSAR-X data. The available SAR data have been processed systematically by SkyGeo (https://skygeo.com). Here we use all the InSAR data available except for the sparse ERS measurements prior to 2004, due to their poor temporal coverage and coherence.
Since the PS-InSAR data represent surface displacements for different persistent scatterers for each radar track, some common datum must be chosen to combine the individual InSAR tracks. In addition, the large persistent scatterer (PS) numbers in this region mean that processing is computationally expensive, without any additional information gained from the large data set due to the smooth pattern of ground deformation. A quadtree decomposition, a method that subdivides the signal depending on the variance of the signal, is often adopted. However, this method is independent for each satellite track and does not use a common grid between individual tracks. Instead, for an initial outlier rejection method each PS-InSAR data set is decomposed onto a common regularized grid at 300-m spacing, with LOS displacements determined as the median value and grid cell weighting determined by the global variance of the InSAR measurements scaled by the number of persistent scatterer points for each grid cell. This gridding effect on the RADARSAT-2 can be seen in supporting information Figure S4, where supporting information Figure S4a represents the persistent scatterer deformation across the period 2009-2017, supporting information Figure S4b represents the gridded deformation across the period 2009-2017, and Figure S4c represents the gridded normalized 10.1029/2018JB016801 weighting for each of the grid cells. This method creates a deformation signal consistent with that of the raw data.
We used the PS-inSAR results from SkyGeo, which provided filtered time series of LOS measurements for each of the observed scatterers corrected for tropospheric effects using two different atmospheric models. We consider the time series obtained with these two corrections as well as the uncorrected time series. The two atmospheric corrections supplied are a normal atmospheric model based on linear, polynomial, and sinusoidal fitting of the InSAR observations removing any residuals and an atmospheric estimation made using a 90-day Gaussian filter.

cGPS
Deployment of a cGPS network started with one station in March 2013 and acquisition from two prior stations, with expansion to 13 stations by April 2014. This cGPS expands on the PS-InSAR data by recording a three-dimensional deformation signal. The temporal acquisition of each of the stations is shown in Figure 3. The deployment fixing type can be separated into building-attached antennae (0647, DZYL, EEMS, FROO, OVER, TENP, TJUC, VEEN, ZAND, and ZDVN) and ground-fixed antennae (STED, USQU, and ZEER). The different types of antenna location could be important in analyzing the deformation signal, as discussed below, since the ground-fixed stations are coupled better to the subsurface and could be less affected by seasonal temperature variation than are GPS measurements on building structures.

Compaction Modeling
Surface deformation in the Groningen region is clearly located in the region of gas extraction ( Figure 1). Given the spatial distribution of the subsidence, we assume that surface deformation results from reservoir compaction. Deformation causes a surface displacement with a minimum lateral spatial wavelength similar to the depth to the reservoir (Segall et al., 1994;Segall, 2010). Consequently, the subsurface acts as a low-pass filter restricting our resolution when studying reservoir compaction to a length scale comparable to the reservoir depth. Shallow surface deformation should be manifest as short-wavelength features superimposed on the regional signal . Possible near-surface deformation (e.g., soil mobility and salt flow effects; Marketos et al., 2015) are expected to be small and are neglected in this work. In addition, a region to the south of the reservoir that comprises the Veendam dissolution salt mine is removed from the data set to minimize the effect from this unrelated shallow surface deformation.
Earlier investigations of reservoir compaction have coupled geodetic inversion with dynamic reservoir simulation models Fokker & Van Thienen-Visser, 2015;NAM, 2016;and van Thienen-Visser & Fokker, 2017). NAM (2016) compaction model is calculated using a combination of reservoir thickness, spatial pressure variation, and porosity derived reservoir compressibility and calibrated to match the history of gas production from each production well and the limited fluid influx from observation wells (see NAM, 2016, for more details). Earlier studies suggest that the misfit between inverted compaction models and reservoir simulation models show a nonlinear relationship, hypothesized to be due to a complex compaction behavior of the reservoir. This nonlinear relationship is then replicated using either a time decay model (NAM, 2013) or a rate-dependent compaction model (RTiCM, van Thienen-Visser et al., 2015). These dynamic reservoir simulation models are based on sparse geological borehole core data constraining the spatial reservoir compressibility and hydrological properties. For these reasons, in this study we do not implement a coupled approach using a dynamic reservoir simulation model but instead invert surface subsidence measurements directly to estimate reservoir compaction.
In reservoir dynamics, provided that the reservoir lateral variations are greater than the reservoir height, the deformation caused by volume depletion is predominantly in the vertical direction. The total reduction of the reservoir height can therefore be expressed due only to the vertical strain in the reservoir z . The expression of the reservoir height change at a instantaneous point can be expressed as follows: where H is the reservoir height, ΔH is the vertical compaction of the reservoir (referred to as m throughout the remainder of the article), and C m the uniaxial compressibility, the vertical response of the reservoir due 10.1029/2018JB016801 to pressure depletion. Averaged over the thickness of the reservoir (H) this formulation can be reduced to As a starting point we assume a linear poroelastic reservoir embedded in a homogeneous elastic half-space. Surface displacements, d, can be related to the compaction, m, within the reservoir (divided into cuboids) by the linear equation where ⃗ d is the input surface deformation data vector, ⃗ m the compaction within the separate reservoir cuboids, and G is the d × m Green's Function matrix that relates the reservoir compaction to the surface displacements, using the formulation of Kuvshinov (2017) for a strain, also adopted in previous study (Bierman et al., 2015). The coefficients of G depend on the geometry (the position of the point at the surface relative to the source of strain at depth) and the elastic properties of the half-space. The reservoir is discretized as a series of cuboids of 500 m × 500 m in horizontal extent and with a vertical extent the same as the local reservoir thickness. This problem is overdetermined with >9,000 observation points per epoch related to 3,942 modeled parameters.
We describe the inversion procedure used to estimate the compaction ( ⃗ m) at each epoch based on the geodetic data in the next section. Once ⃗ m is estimated, it can be compared with NAM pressure depletion models in combination with modeled compressibility models from porosity derived measurements.

Principal Component Analysis Inversion Modeling
The time evolution of a deformation source can be derived by inverting the displacement data available at each epoch. Such an approach is computationally expensive for large data sets and would not impose any coherent time evolution of the source, since it would yield independent models at each epoch. We use a principal component analysis inversion method (PCAIM, Kositsky & Avouac, 2010) to overcome these limitations by decomposing the data into a set of linearly uncorrelated principal components, each associated with its own spatial and temporal function. The spatial function for each of these principal components can then be inverted separately for compaction and combined with the time functions and principal component amplitudes to determine the compaction history. The computational gain of the method comes from the fact that the original data set can generally be reconstructed within uncertainties with a small number of principal components. In addition, the method makes use of a weighted low rank approximation of the variance-covariance matrix (Srebro & Jaakkola, 2003), which allows weighting of the data according to their uncertainties and proper handling of missing data (Kositsky & Avouac, 2010).
The appropriate number of principal components is chosen based on the chi-square statistics, determined by the output of the PCA and the observed surface deformation. The number of components is selected such that the misfit between the model and the data is of the same magnitude as the measurement uncertainties (i.e., where 2 red 1.0; Kositsky & Avouac, 2010). The reduced chi-square is where N is the number of data values, P is the number of parameters, X obs is the observations, X PCA is the PCA modeled displacements, and obs is the 1-sigma uncertainty of the data and the sum is over the observations.
Regularization is imposed as the Laplacian smoothing on the reservoir compaction, with overall weighting defined by the smoothing parameter 2 . Since the surface subsidence data are irregularly distributed, areas with dense data coverage should have lower smoothing on regional compaction. Following the approach suggested by Tarantola (2005), and detailed in Ader et al. (2012), we can weight the imposed smoothing with a resolution matrix that takes the form where G is the Green's function matrix, C d is the data covariance matrix, and fj is the Laplacian matrix.
The best fitting compaction model depends on the value of the smoothing parameter, with higher values leading to greater smoothing. The smoothing parameter is then selected by the trade-off between solution roughness (the average second-order finite difference sum of each fault patch) and the surface deformation misfit function (see supporting information section S3.2).
Combining the methods outlined above, the governing inversion follows the general form where ⃗ W is the vector of grid cell weighting for each geodetic data set; G is the Green's function matrix that relates the reservoir compaction to the LOS displacements using the nucleus of strain formulation of Geertsma (1973); ∇ 2 is the finite difference approximation of the Laplacian operator, which acts to smooth the compaction distribution, the importance of which is governed by the value of the scalar smoothing factor 2 ; R is the resolution matrix of the compaction cuboids for each geodetic data set; ⃗ m is the vector of compaction for difference cuboids; and ⃗ d is the input LOS deformation data vector for the geodetic data set.

Leveling Results 1964-2013
The PCA was conducted for the entire leveling array across the region. The spatial and temporal deformation signals for a one-component and two-component PCA analysis are shown in supporting information Figures S2 and S3, respectively. The variances of reduced chi-square values suggest that only one component is statistically significant (supporting information Table S1). The additional second component primarily accounts for nonspatially correlated large displacements at a few leveling points. These displacements are due to local effects and are not representative of the regional deformation. We therefore use only the first component that accounts for the majority of the original data variance (as shown in supporting information Table S1).
The temporal evolution of the optical leveling data shows a long-term subsidence since the start of gas extraction ( Figure 4). Initially, the subsidence rate was slow between 1964 and 1974, with an increase in the subsidence rate from 1974 to 2013. This surface deformation pattern is consistent with earlier studies (van Thienen-Visser & Breunese, 2015; Bierman et al., 2015;NAM, 2016), showing the center of subsidence in the northern section of the reservoir. Comparison between the maximum subsidence and the mean reservoir pressure depletion (Figure 4) shows good agreement across the period of 1964-2013. The short-time scale deviations from this common trend between 1990-1992 and 1997-2000 are possibly artifacts due to significant changes in the survey geometry between the acquisitions (see supporting information Figures S2 and S3).

InSAR Results 2004-2016 5.2.1. Envisat
The descending track was processed using a PCA during the period 2004-2010. A three-component analysis is shown in supporting information Figure S8. The first component displays a long-term monotonic deformation across the period of interest. The second principal component displays a similar spatial pattern to the first principal component, but a temporally varying deformation signal, although this deformation signal is typically close to the noise in measurements. We combined these data with the leveling data to refine the model for the 2004-2016 period, as described above, with additional studies showing agreement between Envisat and optical leveling data sets (NAM, 2013).

RADARSAT-2
The RADARSAT-2 data are the longest-running InSAR data set across the Groningen region, spanning late 2009 to 2017. This is the only geodetic data set that has a significant acquisition before and after the change in gas extraction rate of 2014. This data set was decomposed into its principal components for the period spanning 2009-2017 using the method described above. To investigate the effects of the atmospheric corrections, we applied a PCA to the results obtained with the two different atmospheric models used by SkyGeo and to the SAR measurements (supporting information Figures S5-S7). All displayed nearly identical spatial component in the first principal component with a clear signal restricted to the reservoir extent area. The time functions associated with the first components show for each a monotonic trend. For the atmospheric-corrected data sets the second and third components showed significantly different deformation trends and spatial distributions that correlate only poorly with the extent of the reservoir and a seasonal effect in the first component. For the raw data set, uncorrelated for atmospheric effects, the first component is monotonic with the higher components not spatially correlated with the Groningen region. The deformation signals associated with the second and third components are ∼1 mm across the period of interest, which is within the source of error expected from the atmospheric effects and is not significant in the reduced chi-square analysis (supporting information Table S1). Due to the similarity of the first component in the atmospherically corrected and uncorrected data sets, and the PCA analysis only statistically significant to the first component for the atmospherically modeled data, we used only the first component derived from the results with no atmospheric correction ( Figure 5).

TerraSAR-X
TerraSAR-X has been acquired across the region in three viewing geometries, one ascending and two descending tracks. Acquisition of the combined TerraSAR-X data spans 2011-2017, with each track displaying differing initial starting time (TerraSAR-Asc, 25 June 2013, TerraSAR-Dsc1, 10 June 2013, TerraSAR-Dsc2, 14 December 2011. Acquisitions for each of these tracks do not completely cover the region of the Groningen gas reservoir, but the acquisition geometry allows for a combined coverage of the reservoir, with some regions of overlap so that a combined satellite track inversion can reduce ambiguity in a single LOS measurement. The spatial PCA for each of these satellite tracks is shown in supporting information Figure S9 and the temporal component in supporting information Figure S10. Each track shows statistical significance for a single principal component, with all tracks displaying a long-term monotonic deformation signal, in agreement with the RADARSAT-2 data.

InSAR and GPS Results After 2014
Since 2014 the Groningen region has been densely monitored by both InSAR and GPS, with data acquired from four satellite tracks (RadarSAT-2 Descending, TerraSARX-Ascending, TerraSARX-Descending 1, and TerraSARX-Descending 2) and 13 cGPS stations ( Figure 6). To reduce the orbital and atmospheric effects in the InSAR data and the poor spatial coverage of the GPS network, we implemented a combined inversion technique incorporating all the InSAR tracks and GPS stations.
The PS-InSAR tracks were decomposed into their principal components, with each track displaying only a single statistically significant component, and all tracks showing a long-term monotonic temporal signal ( Figure 6, and more information in supporting information section S1.2). The cGPS stations are decomposed into three statistically significant principal components, which show a long-term monotonic deformation and two seasonal signals offset by half a year (supporting information Figure S12). The seasonal components in the cGPS analysis show no clear spatial correlation with the reservoir region and are unlikely to be related to gas extraction, so are not discussed further in this study (for further information see supporting information section S1.3). A common monotonic trend is consistent for all temporal components of all the satellite tracks and cGPS, as shown in Figure 6b. Since all satellite tracks now share the same common temporal evolutions, the InSAR and cGPS can be combined into a joint inversion of the spatial pattern of subsidence (Figure 6a), with the smoothing factor determined as 2 = 0.8 (supporting information Figure S17). This spatial-compaction component was then recombined with the common temporal trend, the averaged temporal trend, to determine compaction variation across the period 2014-2017.
Comparison of the RADARSAT-2 and InSAR-GPS inversion models displays different spatial patterns of compaction, with the InSAR-GPS showing more localized centers and greater compaction in the southern portion of the reservoir (supporting information Figure S20). However, assessment of the RADARSAT-2 data alone is biased by having only a single acquisition geometry and fewer PS points in grid cells, giving biased deformation signals over the same areas. This combined method gives a better evaluation of the deformation in the southern portion of the reservoir, an area which shows poor spatial sampling from RADARSAT-2 data alone.

Compaction Since Start of Gas Extraction
In order to determine the compaction since the start of gas extraction, we combined all the geodetic inversions, with optical leveling spanning 1964-2004, Envisat 2004-2009, RADARSAT-2 2009-2014, and InSAR-GPS 2014-2017. To minimize the model roughness needed to fit the data, we compared the data fit obtained with different values of the regularization parameters ( ), with a Pearson correlation coefficient determined between the compaction models for the overlapping period. Optical leveling and RADARSAT-2 data sets showed correlation at a significance of 15%, with low correlation expected due to differences arising in the optical leveling acquisition or LOS bias in the separate inversion. A significance level of 15% was used with the minimal smoothing parameters selected for each of the geodetic inversions (see supporting information Figure S21). These geodetic data sets are then combined to determine the long-term compaction since the start of gas extraction.

Reservoir Compressibility
Current modeling procedures of NAM (2013) implement a spatial reservoir uniaxial compressibility determined by defining an empirical relationship between well porosity and uniaxial compressibility measurements (Figure 7a). These values are then extended away from the well using petrological logging, sonic log, and seismic data. This procedure is hindered by possible errors in the evaluation of the porosity away from the control points and also by the possibility that uniaxial compressibility could vary due to factors not represented by the porosity. Prior studies have accommodated these porosity derived compaction estimates, comparing with the surface subsidence history to determine possible anelastic deformation of the reservoir due to depletion. The Rate Type isopach Compaction Model (RTiCM;  contains a direct elastic response and secular strain (creep), with the work of van  suggesting that a RTiCM model with a secular strain time decay at 3-8 years, fits the optical leveling derived surface subsidence better than an elastic response only model. Such a procedure is limited to a single geodetic data set, with proposed time decay equal to the repeat acquisition of the geodetic data set. We reinvestigate the possibility of a time-delayed response by determining the uniaxial compressibility through the combination of the NAM (2013) pressure depletion model and multigeodetic derived compaction model and investigating the spatial and temporal variability of these uniaxial compressibility values.
Pressure depletion models are generated using the Shell "MoRES" flow models (NAM, 2013), determining the time evolution of pressure depletion at a finite-element mesh of minimum edge length of 50 m. This procedure incorporates extraction rates, pressure gauge measurements, flow gauge measurements, and tracer timing measurements in history matching to determine the optimal transmissivity and diffusivity across the reservoir. To allow direct comparison with the geodetically derived compaction model, the pressure depletion is determined for the same regular grid structure and is smoothed spatially by the same function used in the geodetic inversion. The relationship of the geodetically derived compaction is compared to the reservoir modeled pressure diffusion thickness and is shown in Figure 7b. Lines represented in black are the comparisons at well locations and in gray for locations within the pressure grid. Reservoir uniaxial compressibility can be determined by formulating the gradient for the relationship between pressure depletion thickness and geodetic derived compaction, governed by the relationship shown in equation (2). The variation in these uniaxial compressibility values for all strain volumes is demonstrated in the normalized density map of Figure 7c. Figure 7c shows no discernible variation of uniaxial compressibility with pressure depletion, demonstrating that the reservoir is an invariant of uniaxial compressibility with pressure depletion. The values for the invariant compressibility were determined from the mean uniaxial compressibility for each strain volume, with values displayed in Figure 8a. A PCA is applied to the uniaxial compressibility variations for each strain volume through time. The analysis finds a statistical significance to a single principal component, with this single component explaining 97% of the data variance. The temporal component of the principal component is used to determine the percentage variation of the mean value of uniaxial compressibility through time and is represented in Figure 8b. Figure 8b demonstrates scattered variation with no discernible trend from the mean uniaxial compressibility value, with values within 1% of the mean value. The lack of variation of uniaxial compressibility with pressure depletion demonstrates that the reservoir has a dominant poroelastic response across the length scales of smoothing for geodetic derived compaction models; secular strain deformation could be present at very local length scales, although would not be recoverable in surface displacement.
The modeled compaction from reservoir simulations is then determined using the modeled pressure-invariant uniaxial compressibility, as shown in Figures 9e-9h. Comparisons are then made between the reservoir compaction models and geodetic compaction models in Figure 9, with the variability between the two models representing possible errors from the assumed pressure-invariant uniaxial compressibility. The compaction model determined using the geodetic derived uniaxial compressibility measurements shows good spatial agreement across all acquisition periods, with minor disagreements attributed to spatial smoothing. The epoch of 2014-2017 (Figures 9d and 9h) represents a period of reduced production in the northern reservoir, in the Loppersum area. This variation in production is observed in both geodetic and reservoir derived compaction models, demonstrating a poroelastic response to changes in production rates, with no discernible delay between production variation and surface subsidence. Such a direct response to changes in gas extraction disagrees with the prior work of van  that attributed the lack of surface subsidence change to a time delay in rock, compaction compared to my suggested pressure diffusion behavior.

Discussion and Conclusion
This study demonstrates the usage of a PCAIM to determine reservoir compaction from a diversity of data sets; optical leveling, PS-InSAR and cGPS. With this approach we were able to generate a model of compaction of the Groningen reservoir since the start of gas extraction, reducing the bias from single geodetic acquisition. We applied this procedure using the InSAR data sets preprocessed for the removal of atmospheric bias. Atmospheric effects actually account for a significant percentage of the original data variance, and it is possible that our results are effected by the possible bias in the atmospheric modeling. These sources of bias are alleviated by the fact that we combine the InSAR data with the leveling and GPS measurements that have different and smaller sensitivity to atmospheric effects. It would be instructive in future studies to use an independent component analysis, as in principle it would not require atmospheric correction a priori. The method should allow separating the geodetic signal and atmospheric effects provided they follow independent time evolution. This should reduce the LOS acquisition error and could allow the detection of additional features missed in the current processing procedure. Such a method is promising for the reservoir monitoring as it could help extract subannual geodetic signals better than with the PCA-based technique used in this study.
This compaction model is then used in conjunction with a reservoir depletion model (NAM, 2013) to investigate the reservoir uniaxial compressibility. We identify no temporal variability in the reservoir uniaxial compressibility, a result that contrasts with some previous studies (e.g., van Thienen-Visser et al., 2015).
We instead determine a pressure invariant but spatially variable reservoir compressibility and demonstrate that this model is sufficient to fit within uncertainties of the various measurements of surface deformation assembled in this study. We demonstrate the production variation imposed in 2014 to reduce compaction in the Loppersum area are observed in both reservoir and geodetically derived compaction models, with small spatial variations attributed to imposed smoothing and detection levels of geodetic data sets. The presence of a dominant pressure-invariant uniaxial compressibility, with the lack of a detectable anelastic effect on reservoir uniaxial compressibility, demonstrates that the reservoir is either represented by a dominant poroelastic response or the material has not yet reached a yielding point where anelastic deformation is more dominant than poroelastic. If the yielding point has not yet been reached, further production induced reservoir strain would supply a greater stress change than an anelastic medium. Therefore, this higher stress could be released in higher-magnitude induced seismicity. In addition, the dominance of a poroelastic response allows the use of modulations in reservoir production, and its relationship with surface subsidence, to better constrain reservoir properties (e.g., diffusivity). These calculations could allow a better constraint of reservoir compaction and its relationship to induced seismicity and help mitigate future seismic hazard.