Subsidence Induced by Gas Extraction: A Data Assimilation Framework to Constrain the Driving Rock Compaction Process at Depth

Surface movement can be induced by many human subsurface activities: production of natural gas, geothermal heat extraction, ground water extraction, phreatic groundwater level lowering, storage of natural gas and CO₂. In this manuscript, we focus on subsidence caused by gas production. While geological interpretations, seismic campaigns and flow modeling often provide a relatively rich pre-existing knowledge, understanding of the driving mechanisms for production-induced subsidence is still poor and forecasts are often very uncertain. This is related to the multiple poorly constrained models that translate gas production to ground surface displacements. Currently, a biased constraint of these models is inferred by arbitrarily pre-selecting a subset of those. Here, we have devised and deployed an integrated approach of the entire chain of models from the flow simulations to the ground surface displacements which, for the first time, accounts for all our pre-existing knowledge in terms of processes and uncertainties attached to them. More specifically for the transfer between reservoir depletion to compacting volume at depth, four reservoir-rock compaction models are a-priori considered, ranging from linear elastic model to nonlinear time-dependent viscous-type model. After assimilation of the geodetic observations (i.e., the ground-surface displacements) with ensemble-smoother algorithms, we demonstrate that even when all the a-priori known complexities were present in all steps of the modelling chain, the model parameter uncertainties of each model could be reduced. Interestingly we demonstrate that one can discriminate which reservoir-rock compaction model driving subsidence is activated at depth. This identification of the activated compaction model at depth is crucial to build confidence in our subsidence forecasts. The predictive power of the integrated approach is demonstrated with an ensemble of synthetic but complex reservoir flow simulations mimicking all the characteristics and uncertainties representative for real gas fields in the north of the Netherlands.

Introduction

Various types of energy-related anthropogenic subsurface activities, such as reservoir gas production, geothermal heat extraction, ground water extraction, storage of natural gas or CO₂, can lead to surface movement (Zoback 2007). The relatively rich pre-existing knowledge of the underground around anthropogenic subsurface activities (e.g., through geological interpretations of seismic campaigns, borehole analysis, and reservoir flow modelling, e.g., see NAM 2016 for the Groningen gas field, the Netherlands) could make one think that the driving physical processes of anthropogenically-induced subsidence are well constrained. This is unfortunately not the case and the contribution of multiple driving physical processes is still highly debated (e.g., Mossop and Segall 1999; van Thienen-Visser et al., 2015a; Fokker et al., 2016).

In this contribution we focus on one instance of anthropogenic subsurface activities which can potentially lead to subsidence: reservoir gas production. During the extraction of natural gas from a gas field, the reservoir pressure decreases, leading to compaction. This reduction in volume at reservoir depth may induce surface subsidence (Doornhof 1992), with consequences for the environment and for human activities (e.g., van Thienen-Visser et al., 2015b; Simeoni et al., 2017). The consequences of reservoir production and its changes call for subsidence forecasting approaches. However, the link between reservoir production and subsidence is non-trivial. Some gas fields in the Netherlands and elsewhere in the world have shown a non-linearity between pressure depletion and subsidence, or even a delay between the start of production and the onset of subsidence and a continuation of subsidence even after production had stopped (Hettema et al., 2002; van Thienen-Visser et al., 2015a). This non-linear relationship between reservoir production and subsidence has been explained by the potential delay between reservoir depletion and compaction and multiple reservoir-compaction models have been already developed (e.g., Mossop 2012; De Waal 1986; NAM 2015; Pijnenburg et al., 2018, 2019). An efficient procedure to discriminate which compaction model is the most likely to explain the subsidence observations is a crucial step to gain confidence in our subsidence predictions. Defining and describing this procedure is the objective of the current contribution.

Subsidence caused by gas extraction is the result of multiple inter-connected physical processes: from the flow of fluids within the reservoir to the ground surface lowering and in between the reservoir rock compaction. As in many scientific realms, like weather forecasting and hydrology, each of these processes is attached to uncertainties which must be honored. In this regard a probabilistic ensemble-based approach (e.g., Reggiani and Weerts 2008; Jaynes 2003, Evensen 2003; Emerick and Reynolds 2013a; Emerick and Reynolds 2013b) can be appropriate and fruitful. An “ensemble” here implies that we build multiple subsidence realizations, based on the possible choices of processes and subsurface parameters. Later, deploying an ensemble-smoother algorithm, this prior ensemble is confronted with the subsidence observations to refine the predictions and to identify the most likely driving processes; and more specifically in our case to discriminate which compaction model is the most likely to explain the subsidence observations.

Ensemble-based inversion procedures of surface subsidence have already been developed in the past (see e.g., Fokker et al., 2012; Baù et al., 2015; Fokker et al., 2016; Gazzola et al., 2021). However, these pioneering works were designed for specific applications, where the full spectrum of uncertainties from the reservoir flow to the modelling of subsidence was not accounted for. A-priori omitting additional complexities of potentially activated processes, the “apparent” constraint of the subset of uncertain models was thus biased. For example (Fokker et al., 2016) assumed a linear compaction model and only varied the reservoir compaction coefficient and the subsurface elastic moduli. Instead in the present contribution our objective is to introduce an integrated approach, coined Ensemble-based Subsidence Interpretation and Prediction (ESIP), which honors all the pre-existing known complexities in terms of processes and subsurface parameters.

For evaluation purposes, our integrated approach is applied on numerical flow simulations aiming to mimic the characteristics of real gas fields in the north of Netherlands. First, we describe our fit-for-purpose exercise and the challenges attached to it. Second, we describe the integrated approach from each forward model to the inversion procedure. Third, the results of the integrated approach applied on our tailored exercise are presented and discussed. We demonstrate that while accounting for all the pre-existing known complexities, our integrated approach enables to identify the correct compaction process activated at depth.

Description of the Exercise and Challenges

In the exercise we designed we use the reservoir flow model of a typical gas field in the north of the Netherlands. The complexities of such a field are reproduced with a synthetic reservoir model, and multiple reservoir simulations are performed to mimic the prevalent uncertainty. One of the realizations is selected as the synthetic truth model and used for creating synthetic subsidence data, while the others are employed to show how the knowledge about the reservoir compaction and the subsidence forecasting can be improved with these data.

It is often the case that over-simplified synthetic reservoir models are used for testing purposes. This is the case of the simplified disk-shaped reservoir with spatially uniform pressure drop presented by Baù et al. (2015). Instead here our synthetic reservoir model preserves all the complexities of a real gas field (see Section 3.1) in terms of geometry and spatio-temporal pressure drop history.

Making use of model-based synthetic subsidence data presents pros and cons. The main inconvenient is that synthetic data can be missing additional complexities that are unmapped in our model-based simplified version of the reality. However, it is important to stress that for our exercise these additional complexities were accounted for by an appropriate evaluation of the uncertainties of the observations formalized with the full data covariance matrix (see Section 3.6). The main advantage of using synthetic data is that in our model-based realm, one can virtually control everything and thus it makes it ideal for testing purposes. More specifically in contrast to real data, we know upfront which processes and parameters are at the origin of the synthetic data. One can thus test if our integrated approach can effectively recover these processes and parameters.

Our exercise aims to estimate and forecast ground-surface displacements from 2008 to 2015 based on synthetic data from 1986 to 2007. We ran reservoir simulations for multiple scenarios from 1986 to 2015. Because the synthetic data correspond to one member of the prior ensemble, we have data and models from 1986 to 2015. Our predictions from 1986 to 2015, based on model conditioning from 1986 to 2007, can therefore be directly checked against our synthetic measurements from 1986 to 2015.

Besides testing the predictive and forecasting capabilities of our integrated approach, the prime scientific challenge that we are tackling with our exercise is: does our integrated approach combining forward models (and their uncertainties) with ensemble-smoother algorithms help to discriminate which type of compaction model best explains the synthetic data? Four types of compactions models, ranging from linear elastic model to nonlinear time-dependent viscous-type model, were considered: the linear, the bilinear, the time decay, and the rate type model (see Section 3.3). We selected the bilinear compaction model to generate the synthetic data and thus the goal of the exercise is to test if our integrated approach enables to retrieve this information.

Description of the Integrated Approach

This section details the integrated approach that we designed to constrain the compaction process at depth and consequently to improve forecasting capability of subsidence. The goal of our ensemble-based integrated approach is first to perform a global uncertainty analysis where all the uncertain parameters of each process of the entire chain of models are vary simultaneously. The succession of processes to be modelled is: 1) flow of fluids within the reservoir, 2) compaction at depth of the reservoir caused by pressure depletion, 3) ground surface displacement caused by the reservoir compaction and mechanical response of the surrounding medium around it. For each of these processes, multiple parameters can be involved. The a-priori range and distribution of each of these model parameters is pre-defined by our expert knowledge and the pre-existing literature for our case-study: typical gas fields in the north of the Netherlands. The distribution of each model parameter is uniform, and we use the Latin hypercube sampling method to randomly select parameter values from the overall multidimensional distribution.

Ensemble of Pressure-Depletion Scenarios

A synthetic gas field was built with Shell’s propriety reservoir simulator MoReS, mimicking the complexities of a typical real field in the north of the Netherlands. These complexities (see Figure 1) are: 1) a fault potentially compartmentalizing the reservoir with an uncertain sealing capacity (that is an uncertain fault transmissibility mapped into a fault permeability multiplier with an uniform a-priori distribution between ${10}^{-4}$ and $1$), 2) an uncertain reservoir permeability including the possibility of a flow barrier (with an uniform a-priori distribution between ${10}^{-4} mD$ and $1 mD$) and a high permeability streak (with an uniform a-priori distribution between $15 mD$ and $1000 mD$), 3) an uncertain amount of residual gas in the aquifer (with an uniform a-priori distribution between $0$ and $0.2$), and 4) irregular reservoir boundaries and grid size. The variability in the permeability of the flow barrier results in the variability of the aquifer activity. We built an ensemble of pressure-depletion scenarios to map our a-priori belief in terms of reservoir geometry, geology, flow properties, and their uncertainty. Figures 2, 3 display the first and the last of the 76 members of the prior ensemble of flow simulations. Each member is simulating 30 years of field production; the total volumes of gas and water produced are representative of real fields in the north of the Netherlands. The depletion level of a reservoir-connected aquifer is regularly uncertain in real scenarios. This is also the case for our synthetic scenario, where the prior ensemble is intentionally mapping uncertain pre-existing knowledge on the degree of aquifer depletion. The pressure depletion of the aquifer is clearly variable between members of the prior ensemble (see Figures 2, 3).

FIGURE 1

FIGURE 1. Cross-sections of the reservoir and connected aquifer for one member of the ensemble of flow simulations. Top: initial water saturation (fraction of the pore space occupied by water) of the model indicating the progressive interface between reservoir (west side) and aquifer (east side). Middle/Bottom: horizontal/vertical permeability-structures (in mD).

FIGURE 2

FIGURE 2. Cross-sections of the pore pressure distribution in the reservoir and connected aquifer for two members of the ensemble of flow simulations. Top: pore pressure pre-production [1986]. Bottom: pore pressure at the end of production [2015].

FIGURE 3

FIGURE 3. Top-views of the pore pressure distribution in the reservoir and connected aquifer for two members of the ensemble of flow simulations. Top: pore pressure pre-production [1986]. Bottom: pore pressure at the end of production [2015]. The blue rectangle of zero pore pressure, on the south-side of the model, corresponds to an intentionally cropped-volume (non-active zone) intended to mimic the complex geometry of real field cases. Still for the sake of capturing structural complexities of real field cases, east-west and north-south corridors of relatively smaller grid-block volumes are also visible.

Upscaling

The 3-D pressure fields (Figures 1–3) are upscaled to 2-D maps in two steps: the vertical averaging of the pressure in the reservoir layers and the horizontal averaging over a coarser grid. This averaging process reduces the computational runtime and is justified by the fact that a distribution of compaction over an upscaled cell in a thin and deep reservoir has a limited effect on the surface subsidence (Geertsma, 1973).

The upscaling needs to be performed such that the amount of compaction for the upscaled 2-D grid and the initial layered 3-D grid is identical. Therefore, the vertical and horizontal averaging of the 3-D pressure field need to be weighted, to take into account the variability in compaction between each individual grid block of the reservoir model. This variability in compaction is controlled by differences in pressures, volumes, and porosities between each grid block. The vertically averaged pressures $P_{av}$, volumes $V_{tot}$, porosities ${\varnothing }_{av}$ and grid block center positions $X_{av}$ and $Y_{av}$ for one vertical column of the 3-D grid can be expressed as:

$P_{av}=\frac{{{\displaystyle \sum }}_{k=1}^n{C}_{mk}\left({\varnothing }_k\right).P_k.V_k}{C_m\left({\varnothing }_{av}\right).V_{tot}}(1)$

$V_{tot}={\displaystyle \sum }_{k=1}^n{V}_k(2)$

${\varnothing }_{av}={\displaystyle \sum }_{k=1}^n{\varnothing }_k.\frac{V_k}{V_{tot}}(3)$

$X_{av}={\displaystyle \sum }_{k=1}^n{X}_k.\frac{V_k}{V_{tot}}(4)$

$Y_{av}={\displaystyle \sum }_{k=1}^n{Y}_k.\frac{V_k}{V_{tot}}(5)$

In Equation 1 $C_{mk}\left({\varnothing }_k\right)$ corresponds to the compaction coefficient of one grid block as a function of its porosity, pressure, and time. This specific relationship is constrained by uniaxial experiments (see Section 3.3). The volume $C_{mk}\left({\varnothing }_k\right)$ considered for each individual grid block is the net volume, that is only including reservoir parts. After this first step of vertical upscaling we have a 2-D irregular fine grid of $P_{av}$, ${\varnothing }_{av}$, $V_{tot}$ at upscaled positions ($X_{av}$, $Y_{av}$).

The second step of the upscaling consists in the horizontal averaging of pressures, volumes, and porosities of the 2-D fine grid $(X_{av},Y_{av})$ falling within the limits of a 2-D regular coarse grid covering the area of the reservoir model. The regular coarse grid size can be flexible and is typically between 500 and 1000 m. In order to maintain the same amount of compaction before and after upscaling, the same averaging procedure needs to be followed as for the vertical averaging in Equations 1–5. After this second step of horizontal upscaling we end-up with a 2-D coarse irregular grid with new $(x_{av},y_{av})$ positions, pressures, volumes, and porosities. Figure 4 presents the upscaled 2-D pressure fields for the first and last member of the prior ensemble of reservoir flow simulations.

FIGURE 4

FIGURE 4. Upscaled pore pressure distributions for two members of the ensemble of flow simulations. The dark rectangles indicate the surface locations of the geodetic benchmarks; their uneven spatial distribution is mimicking real scenarios of the north of the Netherlands.

From Pressure Depletion to Compaction

Translating the pressure depletion between two epochs to volume compaction is non-trivial. Models for rock compaction based on laboratory measurements and/or field measurements are still highly debated (De Waal, 1986; Hettema et al., 2002; van Thienen-Visser et al., 2015a; NAM, 2015; Pruiksma et al., 2015; Spiers et al., 2017; van Eijs and van der Wal, 2017; Pijnenburg et al., 2018, 2019; Smith et al., 2019). For our exercise, four types of compaction models, ranging from linear elastic model to nonlinear time-dependent viscous-type model, were considered: the linear, the bilinear, the time decay, and the rate type model.

The linear elastic compaction model assumes a linear relationship between pressure depletion $dP$ and compaction $V_{comp}$:

$V_{comp}\left(x,y,t\right)=C_m\left(x,y\right) V\left(x,y\right) dP(x,y,t)(6)$

where $V$ is the grid block net volume. Only one material parameter, the compaction coefficient of the reservoir rock $C_m$ is needed for the linear compaction model. This last model parameter can be chosen as a single value or as a function of porosity $\phi$. A polynomial regression derived from uniaxial laboratory experiments on rock samples representative of the gas reservoirs of the north of the Netherlands (NAM, 2017) yielded a third-order polynomial:

$C_m\left[{10}^{-5}ba{r}^{-1}\right]=267.3 {\phi }^3-68.72 {\phi }^2+9.85 \phi +0.21(7)$

The three other compaction models (bilinear, time-decay, rate type) intend to capture non-linear time-dependent viscous-type and/or pressure-dependent behavior. They have been formulated in response to geomechanical and geodetic measurements pointing out nonlinearities and delays in compaction and subsidence relative to the production of gas fields (De Waal, 1986; Hettema et al., 2002; van Thienen-Visser et al., 2015a; NAM, 2015; Pruiksma et al., 2015; Spiers et al., 2017; van Eijs and van der Wal, 2017). Recent observations also indicate continuing subsidence even when production has stopped.

In the bilinear compaction model (NAM 2017), the compaction coefficient is assumed to be pressure-dependent. It combines two linear relationships between pressure depletion and reservoir rock compaction as:

$V_{com{p}_{pre}}\left(x,y,t\right)=C_{m_{pre}}\left(x,y\right) V\left(x,y\right) \left(P_0\left(x,y\right)-P\left(x,y,t\right)\right)(8)$

$V_{com{p}_{post}}\left(x,y,t\right)=C_{m_{pre}}\left(x,y\right) V\left(x,y\right) \left(P_0\left(x,y\right)-P_{trans}\right)+C_{m_{post}}\left(x,y\right) V\left(x,y\right) \left( P_{trans}-P\left(x,y,t\right)\right)(9)$

where $P_0$ and $P_{trans}$ respectively define the initial pressure before the start of production and the transition pressure. The first relationship, Equation 8, should fit the stiff behavior at early stages using a low $C_{m_{pre}}$ value from the onset of pressure depletion up to the transition pressure $P_{trans}$. The second relationship, Equation 9 addresses the later, weaker behavior using a high value for $C_{m_{post}}$ for higher pressures. Two material parameters $C_{m_{pre}}$, $C_{m_{post}}$ and the pressure $P_{trans}$ are required to compute the bilinear compaction.

The time-decay model was inspired by the ubiquitous diffusive behavior of many physical systems pushed into non-equilibrium high-energy states, which slowly decay to their low-energy equilibrium states again (Mossop 2012). Such a diffusion-type process can be modeled using a convolution of a linear relationship between pressure depletion and reservoir rock compaction with an exponential time decay function:

$V_{comp}\left(x,y,t\right)=C_m\left(x,y,t\right) V\left(x,y,t\right) dP\left(x,y,t\right){\ast }_t\frac{1}{\tau }\exp \left[-\frac{t}{\tau }\right](10)$

where ${\ast }_t$ is the convolution operator with respect to time. To compute the time decay compaction, the material parameter $C_m$ and the characteristic time decay constant $\tau$ are required.

The rate type compaction model was derived from laboratory measurements designed for capturing the loading-rate dependency in sandstones (De Waal 1986). The onset and arrest of production can be seen as changes in loading rate of the reservoir rocks. At the change in loading rate, a first direct elastic strain response ${\epsilon }_d$ is modeled, followed by a more gradual creep strain response ${\epsilon }_s$. We followed (Pruiksma et al., 2015) with a formulation based on the linear isotach compaction. The rate type isotach compaction was implemented as an explicit Euler finite-difference scheme keeping a constant time step $\Delta t$ (Pruiksma et al., 2015). To calculate the compaction of one grid block grid $(x, y)$ the applied numerical scheme can be divided into five steps:

1) From the current effective vertical stress $\sigma \left(t\right)$ and strain $\epsilon \left(t\right)$, calculate the creep strain rate as:

${\dot{\epsilon }}_s\left(t\right)=\left(\frac{\epsilon \left(t\right)-{\epsilon }_0}{\sigma \left(t\right)}-C_{m_d}\right){\dot{\sigma }}_r{\left(\frac{\epsilon \left(t\right)-{\epsilon }_0}{\sigma \left(t\right) C_{m_{ref}}}\right)}^{-1/b}(11)$

The vertical stress is derived from the reservoir depth and the mean density ${\rho }_{mean}$ of the subsurface up to the reservoir top $z_r$, and the effective stress is its difference with the pressure:

$\sigma \left(t\right)= {\rho }_{mean} g z_r-P\left(t\right)(12)$

At the onset of pressure depletion/production $t_0$, the direct elastic strain ${\epsilon }_d\left(t_0\right)$ and creep strain ${\epsilon }_s$ ($t_0)$ are both taken zero, and thus the total strain $\epsilon \left(t_0\right)$ is set to zero.

The reference total strain is expressed as:

${\epsilon }_0=-C_{m_{ref}} {\sigma }_r(13)$

with the reference vertical effective stress ${\sigma }_r=\sigma \left(t_0\right)$.

Three material parameters ($C_{m_{ref}}$, $C_{m_d}$, and $b$) and one state parameter (${\dot{\sigma }}_r$) are needed to compute the rate type compaction. The parameter $b$ is an empirical laboratory-derived constant. The material parameters $C_{m_{ref}}$ and $C_{m_d}$ are respectively the reference and direct compaction coefficients, where $C_{m_{ref}}$ is the compaction coefficient corresponding to the pre-depletion loading rate, and thus by definition relatively high. Parameter $C_{m_d}$ is dedicated to map out the direct effect at the change of loading rate. In the scenario of the change of loading rate due to the onset of pressure depletion, $C_{m_d}$ is expected to be relatively low in order to mimic the stiff response of the reservoir rocks. The state parameter ${\dot{\sigma }}_r$ represents the reference vertical effective stress rate at the start of the reservoir depletion.

2) The second step of the Euler scheme consists in calculating the increase in creep strain as:

$\mathrm{\Delta }{\epsilon }_s={\dot{\epsilon }}_s\left(t\right) \mathrm{\Delta }t(14)$

and update the creep strain as:

${\epsilon }_s\left(t_{+1}\right)\to {\epsilon }_s\left(t\right)+\mathrm{\Delta }{\epsilon }_s(15)$

3) The time is updated as $t_{+1}\to t+\mathrm{\Delta }t$

4) Following a linear stress-strain relationship one can calculate the direct elastic strain as:

${\epsilon }_d(t+\mathrm{\Delta }t)=C_{m_d} (\sigma \left(t+\mathrm{\Delta }t\right)-{\sigma }_r)(16)$

5) Finally One can Calculate the Total Cumulative Strain as

$\epsilon \left(t+\mathrm{\Delta }t\right)={\epsilon }_s(t+\mathrm{\Delta }t)+{\epsilon }_d(t+\mathrm{\Delta }t)(17)$

And the total cumulative compaction as:

$V_{comp}\left(t+\mathrm{\Delta }t\right)=-V_0 \epsilon \left(t+\mathrm{\Delta }t\right)(18)$

with $V_0$ the grid block net volume, assumed constant over time. After this last fifth step the workflow returns to the first step for the next time step.

From Compaction to Ground-Surface Displacements

In order to propagate the 2-D compaction fields to surface subsidence we employ influence functions, also called Green functions. They are rotationally symmetric surface displacement profiles for a unit volume of compaction (a nucleus of volumetric strain) and they are used in conjunction with the 2-D compaction fields, in order to calculate the total surface displacements at the desired geodetic benchmarks locations. We used the semi-analytical approach developed by (Fokker & Orlic 2006), which is based on the nucleus of strain concept (Mindlin & Cheng 1950) and which extends Geertsma’s analytical solutions (Geertsma 1973; Van Opstal 1974) in the sense that it can handle a layered elastic subsurface. The influence functions give the vertical and horizontal surface displacements for a single “nucleus” (a unit center of compression) located at reservoir depth, and for given elastic parameters. The required input parameters are the reservoir depth, the depths of the layer interfaces, Young’s modulus and Poisson’s ratio of each layer. For the present study, an elastic two-layer subsurface model has been considered. The interface between the layers was put at −3800 m depth, 80 m below the reservoir (−3720 m). Uncertainties of the model parameters (Young’s modulus and Poisson’s) of both layers have been mapped out in the prior ensembles (see Section 3.5) to cover the broad range of possibilities from a homogeneous elastic subsurface (no elastic differences between both layers) to a subsurface with a stiff basement (that is with a strong contrast of elasticity between the top and bottom layer).

Prior Predictions

For each type of compaction model and based on the prior information, a prior model vector ensemble of $N_e$ vectors $M_0=\left(m_1,m_2, \dots m_{Ne}\right)$ is created. Following the Latin hypercube sampling procedure, these vectors are generated by stochastically selecting values in the prior distributions of the driving input parameters. In this sense the available prior knowledge is mapped to the ensuing ground-surface displacement, by means of the geomechanical forward models for both the compaction and the influence functions to translate compaction at depth to ground-surface displacements. For each member of the 76 prior flow simulations, 10 sets of parameters for the geomechanical forward models (compaction model + influence function) are selected with the Latin hypercube sampling approach. 760 members of prior ground surface displacements are thus initially computed for each compaction model. For each prior ensemble, the length $N_e$ of the prior model vector ensemble $M_0$ is thus initially 760. Note here that the prior model parameters used to generate the ensemble of reservoir flow simulations are not included in $M_0$ and thus will not be updated during the data assimilation procedure (detailed in Section 3.7).

The prior distributions of all the model parameters used for the compaction models and influence function are presented respectively in Tables 1, 2. For some model parameters, stochastic factors are employed to generate the prior ensembles. The prior ensembles of some compaction coefficients ($C_m$, $C_{m_{pre}}$, and $C_{m_d}$) are generated by multiplying the initial compaction coefficient value $C_m$ (as defined by Equation 7) with their respective stochastic factors ($C_m fac$, $C_{m_{pre}} fac$, and $C_{m_d} fac$). The early stage low compaction coefficient $C_{m_{pre}}$ of the bilinear compaction model, is then multiplied by the stochastic factor $C_{m_{post}} fac$ to obtain the prior ensemble of relatively higher later stage compaction coefficient $C_{m_{post}}$. The prior ensemble of pre-depletion relatively higher compaction coefficient $C_{m_{ref}}$ of the rate type compaction model, is generated by dividing $C_{m_d}$ with the stochastic factor $C_{m_{ref}} fac$.

TABLE 1

TABLE 1. Bounded uniform prior distributions for each model parameters used to build each type of compaction models.

TABLE 2

TABLE 2. Bounded uniform prior distributions for each model parameters used to build the semi-analytical influence functions for an elastic two-layers cake subsurface model. The subscript “r” indicates the elastic properties of the top layer (from the ground surface to depth −3800 m). The subscript “b” indicates the elastic properties of the bottom layer (from −3800 m to an infinite depth).

From our a-priori knowledge one can generate an ensemble of surface displacements predictions for each component of surface displacements (horizontal displacement in West-East direction u₁ an North-South direction u₂, and vertical displacement u₃) at the specific measurement locations, that is the geodetic benchmarks locations. From these ensembles of surface displacements one can generate a prior ensemble of 90 double differences for the training period [1986–2007]. A double difference indicates a difference in time of a level difference in space (van Leijen et al., 2017), which circumvents the effect of the choice of references in time and space. In other words, one double difference is the difference of ground surface displacements (either vertical or horizontal displacements) between two benchmarks locations and between two epochs/times of geodetic campaigns. The combinations of modeled 90 double differences (in terms of benchmark locations and time of the geodetic campaigns) correspond to those of the measured (synthetic in our case) 90 double differences (see Section 3.6). These 90 double differences include all the benchmark locations presented in Figure 4.

The geomechanical forward model (compaction model + influence function) is indicated by the function $\mathit{ }G$ working on each vector of prior model parameters. For the sake of clarity, the ensembles $G{M}_0$ for each type of compaction model, will from now on only refer to the ensemble of $N_e$ prior double difference $dd$ predictions:

$G{M}_0=d{d}^{prior}=\left(d{d}_1^{prior},d{d}_2^{prior}, \dots d{d}_{Ne}^{prior}\right)(19)$

For each type of compaction model one can define a mean and a covariance matrix of the prior double differences predictions. The mean over the $N_e$ members is defined as:

$<G{M}_0>= \mu \left[d{d}^{prior}\right]=\frac{1}{N_e}{\displaystyle \sum _{i=1}^{Ne}d{d}_i^{prior}}(20)$

The covariance over the $N_e$ members between the jth location and the kth location is defined by

$c_{jk}^{prior}=\frac{1}{N_e-1} {\displaystyle \sum _{i=1}^{N_e}\left(d{d}_{ij}^{prior}-\mu \left[d{d}_j^{prior}\right]\right)}\mathrm{·}\left(d{d}_{ik}^{prior}-\mu \left[d{d}_k^{prior}\right]\right)(21)$

which can be written in matrix notation as:

$C^{prior}=\left(G{M}_0^{\prime }G{M}_0^{\prime T}\right)/(N_e-1) (22)$

And where

$G{M}_0^{\prime }=G{M}_0-<G{M}_0> (23)$

The mean of the priors $<G{M}_0>$ can be seen as the best prior estimates for the various model types. The ensembles can be regarded as the implementation of the prior knowledge, and, assuming a production scenario from time T₁ to T₂ (that is from 2008 to 2015), they can be used to generate the best prior estimates for the various types of compaction model from T₀ to T₂ (from 1986 to 2015). For a real scenario, the times T₁ and T₂ correspond respectively to the present day and a future geodetic campaign, that is the forecasting period.

Geodetic Observations

Two types of ground-surface displacement measurements are considered for our exercise: optical levelling and GPS. Their spatial distribution with regard to the reservoir mimics real cases of the north of the Netherlands (Figure 4). The levelling measurements give vertical displacements only (subsidence or heave, indicated by u₃); GPS also measures the horizontal displacements (West-East component u₁, and North-South u₂). The additional assimilation of the horizontal displacements is particularly interesting for constraining the shape of the subsidence bowl. Our integrated approach allows flexible incorporation of other measurement types like InSAR, if necessary.

The selected 90 double differences for the training period [1986–2007], and including both optical levelling and GPS, were chosen from one member (which was not used anymore in the further analysis) of the prior ensemble of the bilinear compaction model.

Along with the simulated observations, the uncertainties of the observations due to measurement noise, and temporal and spatio-temporal idealization noise were considered, formalized with the full data covariance matrix according to (van Leijen et al., 2017). A detailed description of the method used to simulate the noise of the measurements and the full covariance matrix is beyond the scope of the present publication. We simply note here that the full data-covariance matrix, describing the uncertainty/variance of each measurement, but also the temporal and spatio-temporal correlations between the measurements, even if often neglected (e.g., Fokker et al., 2012), is essential for a proper conditioning of the model parameters. Indeed, the data covariance matrix is key to properly weight the relative importance of each measuring point during the conditioning step (see later Section 3.7).

Figure 5 presents the full data covariance matrix of the 90 double differences for the training period [1986–2007]. The temporally correlated idealization noise, mostly giving positive off-diagonal numbers, describes the effect of benchmark instability. Most negative off-diagonal numbers are related to the spatio-temporal components reflecting additional surface motion not associated with the signal of interest, e.g., compaction of shallow soft soil layers. Additional complexities, such as additional ongoing processes which are not mapped in our forward modelling strategy, are thus accounted for in our integrated approach by the data covariance matrix.

FIGURE 5

FIGURE 5. Full covariance matrix of the 90 double differences (in mm²).

Conditioning of the Models With the Data

The next step consists in confronting or conditioning the prior estimates up to T₁ with the geodetic data acquired up to T₁ (that is from 1986 to 2007 [1986–2007]).

The ensemble-smoother algorithm consists in an inversion scheme, for which the goal is to maximize an objective function J (or minimize–log [J]) of the form (Menke 1989; Tarantola 2005):

$J\left(m\right)=\exp \left[-\frac{1}{2}\left({\left(G\left(m\right)-dd\right)}^T C_{dd}^{-1} \left(G\left(m\right)-dd\right)+{\left(m-m_0\right)}^T C_m^{-1} \left(m-m_0\right)\right)\right] (24)$

where $m$ and $G\left(m\right)$ are respectively the “optimized” (posterior) vector of model parameters and double differences predictions (that is $d{d}^{post}$) from time T₀ to T₁ [1986–2007].

In other words, the objective function is integrated in an inversion scheme seeking the solution for the vector $m$ of model parameters that optimizes the match with the data and with the prior information $m_0$. This way, the ensemble-smoother conditioning step updates both models and predictions.

The optimal “least-square” solution of Equation 24 for one particular realization and assuming a linear inverse problem is given by Tarantola (2005):

$\widehat{m}=m_0+C_m{G}^T{\left(G{C}_m{G}^T+C_{dd}\right)}^{-1}\left(dd-G{m}_0\right)(25)$

Converting this expression into an ensemble expression, we can translate this equation to a forward model that does not have to be a matrix multiplication–i.e., using the equation $G{M}_0^{\text{'}}$ (Equation 23) as the difference of $G{M}_0$ and its mean. The estimate then is replaced by the new ensemble and Equation 25 becomes:

$\widehat{M}=M_0+C_M{G}^T{\left(G{C}_M{G}^T+C_{dd}\right)}^{-1}\left(DD-G{M}_0\right)(26)$

which can be rewritten as:

$\widehat{M}=M_0+\frac{M_0^{\prime }{M}_0^{\prime T}}{N_e-1} G^T{\left\{G\frac{M_0^{\prime }{M}_0^{\prime T}}{N_e-1}{G}^T+C_{dd}\right\}}^{-1}\mathrm{·}\left(DD-G{M}_0\right)=M_0+M_0^{\prime }{M}_0^{\prime T}{G}^T{\left\{G{M}_0^{\prime }{M}_0^{\prime T}{G}^T+\left(N_e-1\right)C_{dd}\right\}}^{-1}·\left(DD-G{M}_0\right)=M_0+M_0^{\prime }{\left[G{M}_0^{\prime }\right]}^T{\left\{G{M}_0^{\prime }{\left[G{M}_0^{\prime }\right]}^T+(N_e-1) C_{dd}^{-1}\right\}}^{-1}·\left(DD-G{M}_0\right)(27)$

In Equation 27 $M_0$ is the prior ensemble of model parameters; $G{M}_0$ represents the result of the non-linear geomechanical forward model working on all the members of the prior ensemble, that is the ensemble of prior double differences predictions. The expression $C_M=M_0^{\text{'}}{M}_0^{\text{'}T}/(N_e-1)$ corresponds to the model covariance matrix which includes the known and belief bandwidths of the model parameters. Primes indicate anomalies with respect to the ensemble mean as $M_0^{\text{'}}=M_0-<M_0>$. Finally $DD=\left(dd+{\epsilon }_1, dd+{\epsilon }_2,\dots dd+{\epsilon }_{N_e}\right)$ corresponds to an ensemble of double differences data realizations created adding to the vector data $dd$ random noise vectors $\epsilon$ corresponding to the uncertainty range of the measurements. This procedure ensures a posterior error covariance that is consistent with the theory. The background of this was discussed by Burgers et al. (1998).

After the data assimilation step has been performed, that is computing Equation 27, posterior model parameters included in $\widehat{M}$ are used to re-run the compaction models and influence functions to generate the posterior ensemble of ground-surface displacements. This implementation is called the Ensemble Smoother with one Single Step of data assimilation (ES-SS).

The non-linearity in the geomechanical forward model can be more appropriately handled by applying multiple steps of data assimilation, so called Ensemble Smoother with Multiple Data Assimilation (ES-MDA) (Emerick and Reynolds 2013a). In ES-MDA, the output ensemble of the smoother is used as input for the next update. The same data are used for all assimilation steps; to compensate for the effect of the multiple application, the data covariances used in the update steps are increased with respect to its actual value. For each step, a data covariance multiplication factor ${\alpha }_i$ must be chosen, so that ${\displaystyle \sum }_{i=1}^{N_a}\frac{1}{{\alpha }_i}=1$ (with $N_a$ the number of assimilation steps). Following Emerick and Reynolds (2013a), we apply four steps with decreasing ${\alpha }_i$ with increasing $i$, giving progressively larger weights to later update steps.

ES-MDA updates model parameters in a more “progressive” way, relative to the rather “aggressive” single step of the ES-SS. Indeed, even if at the end of the four steps of the ES-MDA the actual data covariance is honored, during each step of the ES-MDA the update can be seen as relatively “gentle” because of the data covariance inflation. Therefore, one can expect a higher ability of the ES-MDA for avoiding ensemble collapse. The downside of ES-MDA is that it is computationally slower; in our case, 4 times slower as we pick four steps of assimilation.

Performance Assessment

The performance of the ensemble-smoother algorithms (ES-SS and ES-MDA), in updating the prior ensembles towards posterior ensembles consistent with the data, can quantitatively be evaluated by different metrics. First we define a metric based on the Absolute Error ($AE$) (Baù et al., 2015; Gazzola et al., 2021), which is possible since the synthetic data were generated with known parameters:

$AE\left(dd\right)=\frac{1}{N_e·N_{dd}}{\displaystyle \sum _{j=1}^{N_e}{\displaystyle \sum _{i=1}^{N_{dd}}\left|d{d}_{i,j}-d{d}_{i,true}\right|}}(28)$

$AE\left(m\right)=\frac{1}{N_e}{\displaystyle \sum _{j=1}^{N_e}\left|m_j-m_{true}\right|}(29)$

A second metric, the Average Ensemble Spread ($AES$) can be defined independently from the true values: the variation of the values with regard to their average (Baù et al., 2015; Gazzola et al., 2021):

$AES\left(dd\right)=\frac{1}{N_e·N_{dd}}{\displaystyle \sum _{j=1}^{N_e}{\displaystyle \sum _{i=1}^{N_{dd}}\left|d{d}_{i,j}-{\overline{dd}}_i\right|}}(30)$

$AES\left(m\right)=\frac{1}{N_e}{\displaystyle \sum _{j=1}^{N_e}\left|m_j-\overline{m}\right|}(31)$

Both metrics can be defined independently for the data $\left(dd\right)$ and for the model parameters $\left(m\right)$. For the comparison of data with model results, the average is determined. For the model parameters we calculate the value for each parameter separately, because we want to know the behavior of the separate parameters, and their absolute values differ considerably. Low values of $AE$ correspond to a model close to the truth. Low values of $AES$ indicate that the spread of the ensemble and thus the model uncertainties are small.

To evaluate the update of the ensemble-smoother algorithms (ES-SS and ES-MDA) we evaluate the variation between the prior and the posterior metrics with respect to the prior metrics: $\frac{A{E}_{prior}-A{E}_{posterior}}{A{E}_{prior}}$ and $\frac{AE{S}_{prior}-AE{S}_{posterior}}{AE{S}_{prior}}$. The closer their value is to unity, the better the estimate is (for $AE$) or the stronger the constraining power of the smoother, in the sense of reducing the uncertainty of the estimate (for $AES$). In Sections 4, 5, we will only use the normalized values, but we will still call them $AE$ and $AES$ for convenience.

A third class of quality metrics is formed by a chi-square $\left({\chi }^2\right)$ test. It evaluates the difference between model results and data with account of their covariances. We have defined two formulations. One formulation $\left(\frac{{\chi }_c^2}{N_{dd}}\right)$ in which the covariance is formed by the combination of covariances of both 1) the posterior model ensemble and 2) the data (Equation 32). One formulation $\left(\frac{{\chi }^2 }{N_{dd}}\right)$ in which only the covariance of the data (Equation 33) is used:

$\frac{{\chi }_c^2}{N_{dd}}=\frac{1}{N_{dd}}{(dd-\mu \left[d{d}^{post}\right])}^T.{\left[C_{dd}+C_{GM}\right]}^{-1}.\left(dd-\mu \left[d{d}^{post}\right]\right)(32)$

$\frac{{\chi }^2 }{N_{dd}}=\frac{1}{N_{dd}}{(dd-\mu \left[d{d}^{post}\right])}^T.C_{dd}^{-1}.\left(dd-\mu \left[d{d}^{post}\right]\right)(33)$

The chi-square test (Equations 32 and 33) judge if the model posterior predictions and the data are consistent. Loosely speaking, the average squared mismatch should be of the order of the covariance. A chi-square test close to unity means that the model ensemble is matching the data and that the quality of the match is in agreement with the error covariance of the data.

Results

This Section 4 presents all the results obtained with four prior ensembles, corresponding to each type of compaction model, and one identical member of the bilinear ensemble as data. Supplementary Appendix A.1 presents the results obtained with another set of four different prior ensembles and bilinear data. Similar conclusions can be drawn for the results presented in Section 4 and Supplementary Appendix A.1.

Model Parameters

Figures 6, 7 display the prior and posterior ensembles of model parameters for both the ES-SS and ES-MDA. The posterior ensembles are the ones inferred after calibration during the training period [1986–2007]. The less aggressive multi-step update of the ES-MDA helps to “better constrain” most of the model parameters. “Better constrain” means that the reduction of the parameters uncertainties is more effective with the ES-MDA than with the ES-SS. For the model parameters of the influence function, used to propagate the reservoir compaction in terms of ground-surface displacement, only the ratio between the Young’s modulus of the top layer ($E_r$ from the ground surface to -3800 m depth, that is including the reservoir) and the Young’s modulus of the bottom layer ($E_b$ that is the basement) is updated significantly during the data assimilation procedure. Interestingly this update of the Young’s modulus is very similar for all types of compaction models.

FIGURE 6

FIGURE 6. Prior and posterior (following the ES-SS algorithm and for the training period [1986–2007]) ensembles of model parameters.

FIGURE 7

FIGURE 7. Prior and posterior (following the ES-MDA algorithm and for the training period [1986–2007]) ensembles of model parameters.

For the bilinear compaction model, where we know the true model parameters used to generate the synthetic data, the visual-inspection-based conclusion about the parameter distributions is confirmed by the metrics $AE$ and $AES$. The metrics show more constraining power for the ES-MDA than for the ES-SS (Figure 8): the posterior ensembles of model parameters obtained after the fourth assimilation step of the ES-MDA are less scattered and more centered around the true set of model parameters. Figure 8 also shows that the Young’s modulus of the basement ($E_b$) is more constrained by the data assimilation than the elasticity of the top layer ($E_r$). The conditioning of the ratio of Young’s modulus ($E_r$/ $E_b$) towards the true, as visualized in Figures 6, 7, is thus dominated by the update of the Young’s modulus of the basement ($E_b$).

FIGURE 8

FIGURE 8. $AE$ and $AES$ performance metrics (over the training period [1986–2007]) for the model parameters of the bilinear ensembles.

Ground-Surface Displacements

Figures 9–12 display the prior and posterior ensembles of ground-surface displacements for both the ES-SS and ES-MDA. The results are shown for two specific geodetic benchmark locations, one on top of the reservoir (location 18, Figure 4) and one on top of the aquifer (location 8, Figure 4). The results for other geodetic benchmark locations are presented at Supplementary Appendix A.1. For both ensemble-smoother algorithms (ES-SS and ES-MDA), a qualitative assessment of the curves indicates that the spread of each posterior is smaller than the spread of their prior. The vertical ground-surface displacements of the posterior ensembles are mostly “well aligned” with the data in the sense that the spread of the posterior ensembles encompasses the data. Only for the linear compaction model, and for the combination rate-type compaction model with ES-SS, the alignment is suboptimal. With the exception of the rate-type compaction model with the ES-SS, the horizontal ground-surface displacements (West-East component u₁), even if of low magnitudes, are relatively “well constrained” by the posterior ensembles for both the reservoir benchmark and aquifer benchmark.

FIGURE 9

FIGURE 9. Prior and posterior ground-surface displacements for the linear ensembles. Thin blue curves: prior ensembles. Thin red curves: posterior ensembles. Blue thick curves: mean of the prior ensembles. Red thick curves: mean of the posterior ensembles. Dark thick curves: data. Top row: with the ES-SS algorithm. Bottom row: with the ES-MDA algorithm and with the posterior obtained after four steps of data assimilation. On top of each graph the benchmark location (see Figure 4) and the component of the ground-surface displacement is indicated. The vertical dashed gray line corresponds to the end of the training period.

FIGURE 10

FIGURE 10. Prior and posterior ground-surface displacements for the bilinear ensembles. Thin blue curves: prior ensembles. Thin red curves: posterior ensembles. Blue thick curves: mean of the prior ensembles. Red thick curves: mean of the posterior ensembles. Dark thick curves: data. Top row: with the ES-SS algorithm. Bottom row: with the ES-MDA algorithm and with the posterior obtained after four steps of data assimilation. On top of each graph the benchmark location (see Figure 4) and the component of the ground-surface displacement is indicated. The vertical dashed gray line corresponds to the end of the training period.

FIGURE 11

FIGURE 11. Prior and posterior ground-surface displacements for the time decay ensembles. Thin blue curves: prior ensembles. Thin red curves: posterior ensembles. Blue thick curves: mean of the prior ensembles. Red thick curves: mean of the posterior ensembles. Dark thick curves: data. Top row: with the ES-SS algorithm. Bottom row: with the ES-MDA algorithm and with the posterior obtained after four steps of data assimilation. On top of each graph the benchmark location (see Figure 4) and the component of the ground-surface displacement is indicated. The vertical dashed gray line corresponds to the end of the training period.

FIGURE 12

FIGURE 12. Prior and posterior ground-surface displacements for the rate type ensembles. Thin blue curves: prior ensembles. Thin red curves: posterior ensembles. Blue thick curves: mean of the prior ensembles. Red thick curves: mean of the posterior ensembles. Dark thick curves: data. Top row: with the ES-SS algorithm. Bottom row: with the ES-MDA algorithm and with the posterior obtained after four steps of data assimilation. On top of each graph the benchmark location (see Figure 4) and the component of the ground-surface displacement is indicated. The vertical dashed gray line corresponds to the end of the training period.

The posterior ground-surface displacements at the benchmark location on top of the aquifer (benchmark location 8, Figure 4) are in better agreement with the data than the movements calculated with the prior ensemble. Actually, the mean posteriors moved very close to the data during the assimilation procedure. This is remarkable given 1) the high uncertainty of the aquifer depletion and 2) the low magnitudes of the ground-surface displacements on top of the aquifer compared to those on top of the reservoir. Interestingly, for the benchmark location on top of the aquifer, the relative bandwidths of the posterior (with the ES-MDA after four steps of assimilation) horizontal ground-surface displacements (defined as the ratio between the standard deviation and the mean of the ensemble) are smaller than those of the posterior vertical displacements. As an example, for the bilinear compaction model, the relative bandwidth of the posterior horizontal ground-surface displacements on top of the aquifer is 0.57 whereas for the posterior vertical ground-surface displacements it is 1.01. This points out the importance of the use of the horizontal ground displacements to constrain the shape of the subsidence bowl on the sides of the reservoir and ultimately to constrain the aquifer activity.

For each compaction model, both metrics $AE$ and $AES$, are better for the ES-MDA than for the ES-SS (Figure 13). For the bilinear, time decay, and rate type compaction model, both metrics have reduced considerably after the fourth step of assimilation with the ES-MDA. This demonstrates the good fit of the mean posterior ensembles and the reduction of the initial large prior bandwidths of ground-surface displacements after application of the ES-MDA. Even if the time decay compaction model is slightly under-performing (indicated by a lower $AE$), the performance differences as evaluated by the $AE$ and $AES$, between the posterior ensembles obtained with the ES-MDA for the bilinear, time decay, and rate type compaction are not striking. However, the $\frac{{\chi }_c^2}{N_{dd}}$ and $\frac{{\chi }^2}{N_{dd}}$ metrics indicate that the posterior ensembles of the bilinear compaction model is the one best explaining–that is with the values closest to unity–the synthetic data (Figure 13).

FIGURE 13

FIGURE 13. $AE$, $AES$, $\frac{{\chi }_c^2}{N_{dd}}$ and $\frac{{\chi }^2}{N_{dd}}$ performance metrics for the ensembles of double differences. The horizontal dark dashed lines indicate an optimum posterior ensemble ($\frac{{\chi }_c^2}{N_{dd}}=1$ and $\frac{{\chi }^2}{N_{dd}}=1$).

The suboptimal update obtained with the ES-SS for the ground-surface displacements of the rate-type compaction model (Figure 12) is confirmed by the lower $AE$ and $AES$ than for the bilinear and time decay ensembles (Figure 13). For the rate-type compaction model, the performance differences between the ES-SS and ES-MDA are larger than for the other compaction models. It seems that for a complex compaction model (with a higher number of model parameters) as the rate type, the need for a less aggressive assimilation procedure with multiple steps is more pronounced. The suboptimal alignments of the posterior ground-surface displacements of the linear compaction model (Figure 9) are confirmed by their poor values for $AE$, $\frac{{\chi }_c^2}{N_{dd}}$ and $\frac{{\chi }^2}{N_{dd}}$ (i.e. $\frac{{\chi }_c^2}{N_{dd}}\gg 1$ and $\frac{{\chi }^2}{N_{dd}}\gg 1$) compared to those derived for the other compaction models (Figure 13).

An important question concerns the reliability of our subsidence forecasts. Figure 13 presents the quality metrics both for the training period [1986–2007] and for the full period [1986–2015] including the forecasting period [2008–2015], for which no measurements were used in the assimilation. We see that the metric for the full period is not significantly different from the metric for the training period. Especially the $\frac{{\chi }_c^2}{N_{dd}}$ and $\frac{{\chi }^2}{N_{dd}}$ for the bilinear model do not increase with the extension of the evaluation period. This indicates that the forecast of the bilinear model is trustworthy.

Discussion

The main objective of our exercise was to test if our integrated approach could effectively identify the bilinear compaction model as the best compaction model to explain the synthetic data (themselves generated with a bilinear compaction model). Our analysis has demonstrated that with the combination of the $AE$ and $AES$ metrics it was possible to identify the bilinear and rate-type as best performing compaction models. The additional $\frac{{\chi }_c^2}{N_{dd}}$ and $\frac{{\chi }^2}{N_{dd}}$ metrics allowed to clearly identify the bilinear compaction model as better performing than the rate-type model. Our integrated approach is thus indeed capable to identify the main driving compaction model. This is especially remarkable as our exercise incorporated the full spectrum of uncertainties from the reservoir flow to the modelling of subsidence. To our best knowledge, there are no previous studies that have honored this full spectrum. Even if large uncertainties in the ensembles of pressure-depletion scenarios and influence functions were effectively “blurring” the spatio-temporal signature of the reservoir compaction, the integrated approach successfully recognized the bilinear compaction as best model to explain the synthetic data.

It is important to stress here that the pressure-depletion realizations were not subjected to the assimilation procedure. Even if the uncertainties in the pressure depletion were kept constant throughout the assimilation procedure, the model parameters of the bilinear compaction model were well constrained and aligned with the true-synthetic data. Both footprints in the subsidence signal, of 1) the pressure depletion and 2) the bilinear compaction, were thus effectively disentangled by our integrated approach.

In addition, our analysis has demonstrated that the update of the model parameters of the influence function is independent of the type of compaction model. Indeed, for each type of compaction model, the posterior ensembles of the Young’s modulus ratio (between the two subsurface layers) similarly converged towards the data. This last result indicates that again our integrated approach was capable to disentangle the unique footprint of the elastic layering on the subsidence signal from the footprints of the pressure depletion and compaction.

The main limitation of the current workflow is that the ensemble of reservoir flow simulations is unchanged upon assimilation. Ideally one should update the controlling reservoir parameters and make a new reservoir simulation for every ensemble member. Introducing potential future implementations, it is important to utilize the versatility of applications offered by ensemble-smoother algorithms. Our integrated approach was developed for the production of a natural gas field, but it can easily be adapted for other settings of exploitation of subsurface resources causing ground-surface movement. Typical examples are the injection of CO₂ (Vasco et al., 2010), underground gas storage (Teatini et al., 2011; Castelleto et al., 2013), steam injection (Khakim et al., 2012), geothermal systems (Mossop & Segall 1999; Allis et al., 2009; Vasco et al., 2013), groundwater extraction (Galloway and Burbey 2011; Zhu et al., 2015), phreatic groundwater level management (Fokker et al., 2019) and salt mining (Fokker and Visser 2014). We have recently deployed the ES-SS to identify the driving mechanisms of induced seismicity at the Groningen gas field on a real dataset (Candela et al., 2021). This specific scenario demonstrates that the ensemble-smoother algorithms can be also applied to a discrete dataset, such as an earthquake catalog. The flexibility of the implementation of the ensemble-smoother also offers the possibility to assimilate multiple types of data. As an example, jointly assimilating subsidence geodetic observations and seismicity data, our understanding of the driving mechanisms of both induced subsidence and seismicity might be improved. An active research area in the Netherlands is the disentangling of the relative contribution between deeply seated sources of subsidence (as natural gas extraction) and shallow seated sources of subsidence (Fokker et al., 2019). These sources often interfere. Incorporating in ESIP the additional forward models for shallow seated sources of subsidence, such as soil compaction and oxidation caused by lowering of the phreatic level, might help to achieve this disentangling (Candela et al., 2020).

Conclusion

We have presented an integrated approach for subsidence evaluation, coined ESIP (Ensemble-Based Subsidence Interpretation and Prediction). The key elements of ESIP are:

- it combines fast physics-based forward models with ensemble-smoother algorithms for the data assimilation;

- it covers all the known a-priori complexities in terms of processes and subsurface parameters, and their uncertainties;

- the error covariance matrix of the subsidence data enables to properly weight the contribution of each measuring point during the conditioning;

The predictive power of ESIP has been demonstrated for a synthetic gas reservoir-aquifer system incorporating all the characteristics of real cases in the north of Netherlands. The key conclusions of the presented exercise are:

- the model parameter uncertainties of each fast forward model are reduced with the assimilation of subsidence data;

- our hindcast and forecast are improved with the assimilation of subsidence data;

- the bilinear compaction process gives the best hindcast and forecast with well-defined error bounds in agreement with the data;

- the bilinear compaction process is thus successfully identified as the driving compaction process at depth.

The identification of the driving compaction process is especially important for the subsidence forecasts. In the case of a reservoir production shutdown, the subsidence forecast obtained with each reservoir compaction model would be drastically different. As an example, for an identical pressure diffusion scenario after shutdown, the bilinear compaction model would result in a more severe subsidence than for the linear compaction model.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Author Contributions

TC: Conceptualization, Methodology, Software, Writing-Original draft preparation. AC: Software. EP: Visualization, Writing-Original draft preparation. MP: Writing-Original draft preparation. DH: Writing-Original draft preparation. KK: Writing-Original draft preparation. PF: Supervision, Conceptualization, Methodology, Writing-Original draft preparation.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

The data generated by the integrated approach are available on request from the corresponding author.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/feart.2022.713273/full#supplementary-material

References

Allis, R., Bromley, C., and Currie, S. (2009). Update on Subsidence at the Wairakei-Tauhara Geothermal System, New Zealand. Geothermics 38, 169–180. doi:10.1016/j.geothermics.2008.12.006

Baù, D., Ferronato, M., Gambolati, G., Teatini, P., and Alzraiee, A. (2015). Ensemble Smoothing of Land Subsidence Measurements for Reservoir Geomechanical Characterization. Int. J. Numer. Anal. Meth. Geomech. 39 (2), 207–228. doi:10.1002/nag.2309

Burgers, G., Jan van Leeuwen, P., and Evensen, G. (1998). Analysis Scheme in the Ensemble Kalman Filter. Mon. Wea. Rev. 126 (6), 1719–1724. doi:10.1175/1520-0493(1998)126<1719:asitek>2.0.co;2

Candela, T., Koster, K., Stafleu, J., Visser, W., and Fokker, P. (2020). Towards Regionally Forecasting Shallow Subsidence in the Netherlands. Proc. Iahs, 382, 427–431,. doi:10.5194/piahs-382-427-2020

Candela, T., Pluymaekers, M., Ampuero, J. P., van Wees, J. D., Buijze, L., Wassing, B., et al. (2021). Controls on the Spatio-Temporal Patterns of Induced Seismicity in Groningen Constrained by Physics-Based Modeling with Ensemble-Smoother Data Assimilation. Rev. Geophys. J. Int.

Castelletto, N., Ferronato, M., Gambolati, G., Janna, C., Marzorati, D., and Teatini, P. (2013). Can Natural Fluid Pore Pressure Be Safely Exceeded in Storing Gas Underground? Eng. Geology. 153, 35–44. doi:10.1016/j.enggeo.2012.11.008

De Waal, J. A. (1986). On the Rate Type Compaction Behaviour of sandstone Reservoir rock.Doornhof, D. (1992). Surface Subsidence in the Netherlands-the Groningen Gas-Field. Geologie en Mijnbouw 71 (2), 119–130.

Doornhof, D. (1992). Surface Subsidence in The Netherlands: The Groningen Gas Field. Geologie en Mijnbouw 71, 119–130.

Emerick, A. A., and Reynolds, A. C. (2013a). Ensemble Smoother with Multiple Data Assimilation. Comput. Geosciences 55, 3–15. doi:10.1016/j.cageo.2012.03.011

Emerick, A. A., and Reynolds, A. C. (2013b). Investigation of the Sampling Performance of Ensemble-Based Methods with a Simple Reservoir Model. Comput. Geosci. 17 (2), 325–350. doi:10.1007/s10596-012-9333-z

Evensen, G. (2003). The Ensemble Kalman Filter: Theoretical Formulation and Practical Implementation. Ocean Dyn. 53 (4), 343–367. doi:10.1007/s10236-003-0036-9

Fokker, P. A., Gunnink, J. L., Koster, K., and de Lange, G. (2019). Disentangling and Parameterizing Shallow Sources of Subsidence: Application to a Reclaimed Coastal Area, Flevoland, the Netherlands. J. Geophys. Res. Earth Surf. 124. doi:10.1029/2018JF004975

Fokker, P. A., and Orlic, B. (2006). Semi-analytic Modelling of Subsidence. Math. Geol. 38 (5), 565–589. doi:10.1007/s11004-006-9034-z

Fokker, P. A., and Visser, J. (2014). Estimating the Distribution of Salt Cavern Squeeze Using Subsidence Measurements. In 48th US Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association.

Fokker, P. A., Visser, K., Peters, E., Kunakbayeva, G., and Muntendam-Bos, A. G. (2012). Inversion of Surface Subsidence Data to Quantify Reservoir Compartmentalization: A Field Study. J. Pet. Sci. Eng. 96-97, 10–21. doi:10.1016/j.petrol.2012.06.032

Fokker, P. A., Wassing, B. B. T., van Leijen, F. J., Nieuwland, D. A., and Nieuwland, D. A. (2016). Application of an Ensemble Smoother with Multiple Data Assimilation to the Bergermeer Gas Field, Using PS-InSAR. Geomechanics Energ. Environ. 5, 16–28. doi:10.1016/j.gete.2015.11.003

Galloway, D. L., and Burbey, T. J. (2011). Review: Regional Land Subsidence Accompanying Groundwater Extraction. Hydrogeol J. 19, 1459–1486. doi:10.1007/s10040-011-0775-5

Gazzola, L., Ferronato, M., Frigo, M., Teatini, P., and Zoccarato, C. (2021). “Ensemble Smoother with Multiple Data Assimilation Analysis,” in Challenges and Innovations in Geomechanics. IACMAG 2021. Lecture Notes in Civil Engineering. Editors (Springer, Cham), 125. doi:10.1007/978-3-030-64514-4_93

Geertsma, J. (1973). Land Subsidence above Compacting Oil and Gas Reservoirs. J. Pet. Tech. 25 (06), 734–744. doi:10.2118/3730-pa

Hettema, M., Papamichos, E., and Schutjens, P. (2002). Subsidence Delay: Field Observations and Analysis. Oil Gas Sci. Tech. - Rev. IFP 57 (5), 443–458. doi:10.2516/ogst:2002029

Jaynes, E. T. (2003). Probability Theory: The Logic of Science. Cambridge University Press.

Khakim, M. Y. N., Tsuji, T., and Matsuoka, T. (2012). Geomechanical Modeling for InSAR-Derived Surface Deformation at Steam-Injection Oil Sand fields. J. Pet. Sci. Eng. 96-97, 152–161. doi:10.1016/j.petrol.2012.08.003

Menke, W. (1989). International Geophysics Series. Geophys. Data Anal. Discrete inverse Theor. 45.

Mindlin, R. D., and Cheng, D. H. (1950). Nuclei of Strain in the Semi‐Infinite Solid. J. Appl. Phys. 21 (9), 926–930. doi:10.1063/1.1699785

Mossop, A. (2012). An Explanation for Anomalous Time Dependent Subsidence. In 46th US Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association.

Mossop, A., and Segall, P. (1999). Volume Strain within the Geysers Geothermal Field. J. Geophys. Res. 104, 29113–29131. doi:10.1029/1999jb900284

NAM (2017). Ensemble Based Subsidence Application to the Ameland Gas Field – Long Term Subsidence Study Part Two (LTS II). Available at https://nam-feitenencijfers.data-app.nl/download/rapport/2a2da56c-faee-4453-a8a7-c64c43b418a6?open=true.

NAM (2016). Technical Addendum to the Winningsplan Groningen 2016; Production, Subsidence, Induced Earthquakes and Seismic Hazard and Risk Assessment in the Groningen Field. Available at http://www.nam.nl/algemeen/mediatheek-en-downloads/winningsplan-2016.html.

NAM (2015). Wadden Sea Long Term Subsidence Studies – Overview Report. Available at: https://nam-feitenencijfers.data-app.nl/download/rapport/2ca6c8d8-c0d4-4c10-8460-672f93b4cdaa?open=true.

Pijnenburg, R. P. J., Hangx, S. J. T., and Spiers, C. J. (2019). Inelastic Deformation of the Slochteren sandstone: Stress-Strain Relations and Implications for Induced Seismicity in the Groningen Gas Field. J. Geophys. Res. Solid Earth, 124(5), 5254-5282.

Pijnenburg, R. P. J., Verberne, B. A., Hangx, S. J. T., and Spiers, C. J. (2018). Deformation Behavior of Sandstones from the Seismogenic Groningen Gas Field: Role of Inelastic versus Elastic Mechanisms. J. Geophys. Res. Solid Earth 123(7), 5532–5558. doi:10.1029/2018jb015673

Pruiksma, J. P., Breunese, J. N., van Thienen-Visser, K., and de Waal, J. A. (2015). Isotach Formulation of the Rate Type Compaction Model for sandstone. Int. J. Rock Mech. Mining Sci. 78, 127–132. doi:10.1016/j.ijrmms.2015.06.002

Reggiani, P., and Weerts, A. H. (2008). A Bayesian Approach to Decision-Making under Uncertainty: An Application to Real-Time Forecasting in the River Rhine. J. Hydrol. 356 (1-2), 56–69. doi:10.1016/j.jhydrol.2008.03.027

Simeoni, U., Tessari, U., Corbau, C., Tosatto, O., Polo, P., and Teatini, P. (2017). Impact of Land Subsidence Due to Residual Gas Production on Surficial Infrastructures: The Dosso Degli Angeli Field Study (Ravenna, Northern Italy). Eng. Geology. 229, 1–12. doi:10.1016/j.enggeo.2017.09.008

Smith, J. D., Avouac, J.-P., White, R. S., Copley, A., Gualandi, A., and Bourne, S. (2019). Reconciling the Long-Term Relationship between Reservoir Pore Pressure Depletion and Compaction in the Groningen Region. J. Geophys. Res. Solid Earth 124, 6165–6178. doi:10.1029/2018JB016801

Spiers, C. J., Hangx, S. J. T., and Niemeijer, A. R. (2017). New Approaches in Experimental Research on Rock and Fault Behaviour in the Groningen Gas Field. Neth. J. Geosciences 96 (5), s55–s69. doi:10.1017/njg.2017.32

Tarantola, A. (2005). Inverse Problem Theory and Methods for Model Parameter Estimation. SIAM 89. doi:10.1137/1.9780898717921

Teatini, P., Castelletto, N., Ferronato, M., Gambolati, G., Janna, C., Cairo, E., et al. (2011). Geomechanical Response to Seasonal Gas Storage in Depleted Reservoirs: A Case Study in the Po River basin, Italy. J. Geophys. Res. 116, F02002. doi:10.1029/2010JF001793

Van Eijs, R. M. H. E., and Van der Wal, O. (2017). Field-wide Reservoir Compressibility Estimation through Inversion of Subsidence Data above the Groningen Gas Field. Neth. J. Geosciences 96 (5), S117–S129. doi:10.1017/njg.2017.30

van Leijen, F., Samiei-Esfahany, S., van der Marel, H., and Hanssen, R. (2017). Technical Report. Delft, Netherlands: Delft University of Technology.Uniformization of Geodetic Data for Deformation Analysis; Contribution to the Research Project: Second Phase of the Long-Term Subsidence Study in the Wadden Sea Region (LTS2)

Van Opstal, G. H. C. (1974). The Effect of Base-Rock Rigidity on Subsidence Due to Reservoir Compaction. Proc. 3rd Congr. Int. Soc. Rock Mech. 2, 1102–1111.

van Thienen-Visser, K., Breunese, J. N., and Muntendam-Bos, A. G. (2015b). Subsidence Due to Gas Production in the Wadden Sea: How to Ensure No Harm Will Be Done to Nature. 49th US Rock Mechanics/Geomechanics Symp.

van Thienen-Visser, K., Pruiksma, J. P., and Breunese, J. N. (2015a). Compaction and Subsidence of the Groningen Gas Field in the Netherlands. Proc. IAHS 372, 367–373. doi:10.5194/piahs-372-367-2015

Vasco, D. W., Rucci, A., Ferreti, A., Novali, F., Bissell, R. C., Ringrose, P. S., et al. (2010). Satellite-based Measurements of Surface Deformation Reveal Fluid Flow Associated with the Geological Storage of Carbon Dioxide. Geophys. Res. Lett. 37, L03303. doi:10.1029/2009gl041544

Vasco, D. W., Rutqvist, J., Ferreti, A., Rucci, A., Bellotti, F., Dobson, P., et al. (2013). Monitoring Deformation at the Geysers Geothermal Field, California Using C-Band and X-Band Interferometric Synthetic Aperture Radar. Geophys. Res. Lett. doi:10.1002/grl.50314

Zhu, L., Gong, H., Li, X., Wang, P., Chen, B., Dai, Z., et al. (2015). Land Subsidence Due to Groundwater Withdrawal in the Northern Beijing plain, China. Eng. Geology. 193, 243–255. doi:10.1016/j.enggeo.2015.04.020

Zoback, M. D. (2007). Reservoir Geomechanics. Cambridge Press. doi:10.1017/cbo9780511586477