Compaction and subsidence of the Groningen gas field in the Netherlands

The Groningen gas field in the Netherlands is Europe’s largest gas field. It has been produced since 1963 and production is expected to continue until 2080. The pressure decline in the field causes compaction in the reservoir which is observed as subsidence at the surface. Measured subsidence is characterized by a delay at the start of production. As linear compaction models cannot explain this behavior, alternative compaction models (e.g. Rate Type Compaction Model and Time Decay model) have been investigated that may explain the measured subsidence. Although the compaction models considered in this study give a good match to this delay, their forecasts are significantly different. Future measurements of subsidence in this area will indicate which type of compaction model is preferred. This will lead to better forecasts of subsidence in future. The pattern of overand underestimation of the subsidence is similar for the compaction models investigated and tested. The pattern can be explained by differences in modeled porosity and aquifer activity illustrating the improvement of subsurface knowledge on the reservoir using subsidence measurements.


Introduction
The Groningen gas field, situated in the north of the Netherlands, is the largest European gas field and has been in production since 1963.Pressures have declined from about 350 bar to, on average, 100 bar in 2012.Pressure reduction in the reservoir causes compaction which is observed at the surface as subsidence.Predictions of surface subsidence at the end of field life, in 2080, reach 60 cm.Due to the proximity to sea level surface subsidence is an important issue and is therefore closely monitored.The surface subsidence measurements show a delay of subsidence at the start of production.Conventional linear elastic compaction models cannot explain this behaviour.Other compaction models, such as Rate Type Compaction models (de Waal, 1986) and Time Decay model (Mossop, 2012) can describe this behaviour.We have applied these compaction models to the Groningen reservoir in order to estimate and predict subsidence.The resulting pattern of fit between the modelled and measured subsidence provides more insight into the reservoir properties.

Laboratory compaction experiments
In the laboratory compaction is measured on core samples of reservoir rock.The rocks are subjected to different loading rates and deformation is measured.The amount of compaction depends on the rate of applied stress on the sample.A faster rate gives a stiffer response of the sample (less compaction).For a constant rate, ongoing compaction (creep) is observed.The result of a typical laboratory compaction experiment is shown in Fig. 1.Initially, the sample is loaded at a rate of 2300 bar h −1 , this is subsequently lowered to 62, 6.2 and 0.62 bar h −1 .The stress-strain path belonging to the imposed loading rate is followed until the loading rate is changed.After a change in loading rate, a clear direct strain response is visible.Following this direct response a more gradual response is seen until the stress-strain path corresponding to the new loading rate is reached.
The dependence of compaction on the rate of applied stress offers a possible explanation for the observed subsidence delay above depleting reservoirs and ongoing subsidence for fields after abandonment.
Published by Copernicus Publications on behalf of the International Association of Hydrological Sciences.Waal et al., 2015).The lines of constant stress-strain rate are called isotachs (green) and are lines of constant loading rate.

RTiCM model
The Rate Type Compaction Model (RTCM) has been proposed by de Waal (1986) to explain rate type dependent behavior for sandstones.The original RTCM model is not able to simulate the continuous transition from one loading rate to another.Also the transition from a constant loading rate to a creep phase cannot be simulated.Therefore in Pruiksma et al. (2015), the RTCM model has been rewritten in isotach formulation, now called RTiCM.The isotach concept was developed in geotechnics (den Haan, 1996).Isotachs are lines of constant loading rate in the stress-strain diagram.Every loading rate corresponds to a unique isotach which is reached eventually if a sample is kept subjected to this loading rate.The isotach concept unifies compaction behavior for changes in loading rate with the compaction behavior for creep phases (de Waal, 2015).
The RTiCM model can be summarized in three equations: Here σ is the vertical effective stress in the sample, and ε is the vertical strain.This strain is composed of a vertical direct strain ε d , which is linear elastic and a vertical creep strain ε s .The model has three input parameters (c m,d , c m,ref , b) and two state parameters (σ ref , σref ).The constant ε 0 is defined by Equation ( 1) describes the direct response, while Eqs.( 2) and (3) represent the rate dependent behavior.

Fit to laboratory experiment
The vertical effective stress as a function of time σ (t) is known from the load on the sample.The vertical strain as a function of time ε(t) can then be computed as the output of the model.An explicit Euler forward scheme is used for the solution (for details see Pruiksma et al., 2015).

Time Decay Model
The Time-Decay model (Mossop et al., 2012) equation reads: with ε the strain, c m the compaction coefficient, P (t) the pressure depletion as a function of time and τ a decay constant.Changes in pressure depletion lead to changes in strain with a certain delay time.Contrary to the RTiCM, the Time Decay lacks a direct response and is not based on laboratory experiments.

Application to field case: Groningen
Both the RTiCM and the Time Decay model have been applied to model the compaction and, from there calculate the subsidence of the Groningen gas field.Compaction is determined by the volume of gas bearing rock, the pressure depletion and the compaction coefficient which in turn are dependent on the rock type.For the Rotliegend sandstone which forms the reservoir rock in Groningen, laboratory experiments have shown a polynomial dependency on porosity (van Thienen-Visser et al., 2015).Both compaction models have been constrained using subsidence measurements over the whole of the Groningen field.

Subsidence data
The subsidence data consists of optical levelling and InSAR measurements.The optical levelling campaigns started in 1964 and were then limited to the central and southern part of the Groningen field.In subsequent campaigns in 1972 and 1987 the density and coverage of the optical levelling network was improved to cover the entire field.From 1987 to

Reservoir description
The Groningen gas field is a Rotliegend sandstone reservoir at a depth between 2600 and 3200 m.The thickness of the reservoir gradually ranges from ∼ 100 m in the southsoutheast to ∼ 300 m in the north-northwest.As the thickness changes, the reservoir lithology changes from mainly sandstone to more claystone.From east to west the thickness is relatively uniform (Mijnlieff and Geluk, 2011).The seal is formed by Zechstein salt covering the whole of the field and strongly varying in thickness over the field (Fig. 2).
A static geological reservoir model (NAM, 2013)  for detailed information on the porosity, reservoir thickness and depth.Porosity in the wells was determined from petrophysical analyses of well logs.The porosity of the reservoir varies from 0.12 to 0.22 leading to a compaction coefficient range between 5 × 10 −6 to 2 × 10 −5 bar −1 .A reservoir dynamic model (NAM, 2013) was used to match the historical gas flow and pressures in the reservoir.For both the Time Decay and the RTiCM model compaction was modeled using the dynamic reservoir model.

Approach
The above mentioned compaction models have been used to calculate subsidence at the benchmark locations at the surface.The translation from compaction to subsidence is based on a 1-D model (van Opstal, 1974).In this method, the non-compacting area around the reservoir (side-, overand underburden) is assumed to be homogeneous and elastic.In the case of Groningen, the overburden is highly nonhomogeneous as illustrated in the thickness of the Zechstein (Fig. 2).Additionally a rigid basement at a depth of 5 km is assumed, which determines the shape of the subsidence bowl (van Opstal, 1974) at the edges of the reservoir.The difference between measured and modeled subsidence is used to fit the resulting unknown parameters: τ for Time Decay and To fit the unknown parameters an ensemble of model realizations is generated that spans a large range for the parameters (α c m between 0.4 and 1.0, τ between 0 and 12 years, c m,d c m,ref between 0.4 and 0.8 and σref between 10 −10 and 10 −3 bar yr −1 ).The differences between the modeled and measured subsidence are input in a Bayesian "Red Flag" procedure (Nepveu et al., 2010).The Red Flag procedure gives the highest probability and the lowest χ 2 value to the parameter combination that fits the measured subsidence best.χ 2 is defined as where sig i is the standard deviation of the ith subsidence measurement.The results of the Red Flag procedure are shown for the top three models for the RTiCM model in Table 1 and for the Time Decay model in Table 2.The χ 2 value is smaller for the RTiCM model and the parameter set is more distinct for the RTiCM model, which is indicated by the high probability (99.9 %) of the first parameter set and low probabilities (order ×10 −4 ) for the second and third parameter set.

Results
Figure 3 shows the fit to one benchmark in the center of the Groningen field for both the RTiCM compaction model and the Time Decay model.
Figure 4 shows the compaction in the reservoir in the end of 2011 calculated with the RTiCM model.Superimposed on the compaction are the optical levelling benchmarks and the fit of the modeled subsidence to the measured subsidence.
Both compaction models fit the data, within one standard deviation.A further distinction between these models will only be visible, if significant changes in production rate will occur.Here it is noted that early 2014 the production in the central part of the field has been strongly reduced in view of mitigating the induced seismicity; the effect on compaction of that measure is still to be analyzed, when sufficient data become available.

Discussion
The compaction models (RTiCM and Time Decay) both fit the delay character of the observed subsidence in the first 10 years after the start of the gas production (Fig. 3).They underpredict the maximum subsidence in the center of the subsidence bowl by 2-3 cm for the RTiCM model and 5-6 cm for the Time Decay model at the end of 2011.Spatially both compaction models show the same pattern of overesti-mation and underestimation (Fig. 4).An overestimation of the subsidence occurs in the eastern part and in the northwestern part of the field.An underestimation exists in the southwestern part of the field.The differences between the compaction models are in the amplitude of maximum compaction (RTiCM larger than Time Decay) and the shape of the subsidence bowl at the edges of the field.The RTiCM model predicts a slightly steeper subsidence bowl than the Time Decay model.
As is clear from Fig. 4, relatively large misfits (up to 8 cm) occur over the field.In the van Opstal (1974) method a depth of a rigid basement is assumed, which governs the shape of proc-iahs.net/372/367/2015/Proc.IAHS, 372, 367-373, 2015 the subsidence bowl.A variation of rigid basement would most likely follow the depth of the reservoir which varies from south (2600 m) to north (2900 m) and therefore cannot lead to the observed mismatch patterns.Another assumption in the van Opstal (1974) method is a mechanically homogenous over-and underburden.The Zechstein salt layer which forms the seal of the gas reservoir varies greatly over small distances, which makes the overburden not homogeneous (Fig. 2).The thickness map of the Zechstein, however, does not correlate with the pattern of misfit.Also the presence of shallow soft soils such as peat and clay does not correspond to the misfit pattern (Wassing and Dost, 2012), indicating that shallow causes cannot explain the discrepancies.We conclude that the misfit pattern has to be related to the reservoir itself.Compaction in the reservoir depends on the volume of gas bearing rock, the compaction coefficient and the pressure depletion.As the compaction coefficient has a polynomial dependency on porosity a slight change in porosity will induce a large change in compaction.The pattern of the porosity (TNO, 2013) shows a high in the same area as the overestimation in the eastern part of the field.Figure 5, left, shows the pattern of misfit if a 0.85 reduction factor is applied on the porosity in a circle with a radius of 4 km in the area of high porosity.The fit to the measured subsidence is now within 2 cm, which is within two standard deviations and thus acceptable.For the other misfit areas no clear correlation with porosity is found.In both the northwestern as well as the southwestern areas an aquifer adds to the uncertainty of the reservoir model and thus to uncertainties of the pressures in those parts of the field.In the southwestern part of the field a less active aquifer (i.e. more depletion) would cause a larger modeled subsidence.In Fig. 5, right, the compaction modeled with a reservoir model consisting of a less active aquifer is shown.The misfit in the area improves considerably (van Thienen-Visser et al., 2014).The anomaly to the northwest is more difficult to explain: a more active aquifer leads to higher pressures and thus less modeled subsidence.The aquifer in the northwest is already modeled as a large active aquifer and an even larger aquifer would not fit the pressures and the water measurements in the northwest of the field.This area of the field is characterized by a couple of faults with large throws.The porosity as well as the other reservoir parameters are modelled continuously over these faults, which is not realistic due to the large throws of these faults.As a hypothesis, the mismatch in the northwestern area could possibly be explained by the presence of gas-blow-free water level, inhibiting the aquifer activity.

Conclusions
The delay character of subsidence measured over the Groningen field can be fitted using a Rate Type Compaction Model.The RTiCM model gives the smallest differences between measured and modeled subsidence.The pattern of over-and underestimation of the subsidence is similar for both compaction models investigated and tested here.The pattern can be explained by differences in modeled porosity and aquifer activity illustrating the improvement of subsurface knowledge on the reservoir using subsidence measurements.

Figure 1 .
Figure 1.A compaction experiment (black) with the fit to the compaction experiment (red) with the RTiCM compaction model (de Waal et al., 2015).The lines of constant stress-strain rate are called isotachs (green) and are lines of constant loading rate.

Figure 2 .
Figure 2. Thickness (m) of the Zechstein salt layer.The difference between calculated and modeled subsidence (m) is indicated at the benchmark locations (label: subsidence diff).A red color indicates that the measured subsidence is larger than the modeled subsidence (from: TNO, 2013).
ref , σref for the RTiCM model, as well as a multiplier (α c m ) on the porosity compaction coefficient relation for both compaction models.Additionally, the RTiCM model uses the laboratory experiments for the b value, and to define a range of c m,d c m,ref values used in the Red Flag procedure (Nepveu et al., 2010).

Figure 3 .
Figure 3. Different compaction models (RTiCM, linear isotach and Time Decay) and their fit to a levelling benchmark in the center of the field (from: TNO, 2013).

Figure 4 .
Figure 4. Compaction in the Groningen reservoir at January 2012 calculated with the RTiCM model (from TNO, 2014).The difference between calculated and modeled subsidence is indicated at the benchmark locations (label: Subsidence diff).A red color indicates that the measured subsidence is larger than the modeled subsidence.

Figure 5 .
Figure 5. (a) Pattern of misfit (m) in the eastern part of the field, with an adaptation for the porosity of 0.85 applied in the green circle with a radius of 4 km.(b) Pattern of misfit (m) using a more active aquifer in the southwestern part of the field.

Table 1 .
Red Flag output for the RTiCM model.

Table 2 .
Red Flag output for the Time Decay model.