Regulating Subsidence and its uncertainty in the Dutch Wadden Sea

At the start of gas production its effects on land subsidence are not certain. There are uncertainties in mechanisms, models and parameters. Examples are non-linear deformation of reservoir rock, fault transmissibility, behaviour of overlaying salt and aquifer activity. Looking back at historical cases in the Netherlands, a factor two or three difference between initial prediction and final outcome is quite common. As the Dutch regulator, SSM is tasked with assuring proper management by operators of the risks associated with land subsidence from natural gas production in The Netherlands. Large initial uncertainties can only be tolerated if operators can demonstrate that timely actions can still be taken when predefined subsidence limits are at risk of being exceeded now or in the future. The applied regulatory approach is illustrated by the case history of gas production induced subsidence in the Dutch Wadden Sea area. This environmentally highly sensitive UNESCO World Heritage Site is a natural gas province. Extensive legal, technical and organisational frameworks are in place to prevent damage to its natural values. Initial uncertainties in the predicted subsidence (rate) were later exacerbated by the detection of strong non-linear effects in the observed subsidence behaviour, leading to new concerns. It was realised that – depending on the underlying physical cause(s) – there will be a different impact on future subsidence. To assure proper management of the additional uncertainty by the operator, several improvements in the regulatory approach have been implemented. Possible underlying mechanisms had to be studied in depth and improved data analysis techniques were requested to narrow down uncertainties as time progresses. The approach involves intensified field monitoring, scenario's covering the full range of uncertainties and a particle filter approach to handle uncertainties in predictions and measurements. Spatial-temporal double differences, production data and the full covariance matrix are used to confront scenario predictions against measurements and to assess their relative probability. The regulator is actively involved in assuring this integrated control loop of predictions, monitoring, updating, mitigation measures and the closing of knowledge gaps. The regulator involvement is supported in the Mining law and by appropriate conditions in the production plan assent. With the approach it can be confidently assured that subsidence (rate) will remain within the allowed range.

1 Introduction

The Dutch regulator State Supervision of Mines (SSM) is, amongst other things, tasked with the supervision of the risk associated with land subsidence due to natural gas production. With half the country near sea level and some of the production activities near the unique natural areas of the Wadden Sea tidal flats, land subsidence from gas production is not an issue that is taken lightly in the Netherlands. Experience over the last 50 years demonstrates that accurate prediction of induced land subsidence at the early stages of a gas production project is difficult. A 100 %–200 % deviation in outcome compared to what was considered the most likely scenario at the start of production is not unusual (de Waal et al., 2012, 2015b). From the perspective of the regulator, initial predictions with such large uncertainties can only be accepted if the operator can demonstrate – at any moment during the production period – that timely actions can still be taken when predefined subsidence limits are at risk of being exceeded now or in the future. A scenario-based approach with uncertainties reducing as time progresses and a continued demonstration that control can be retained if reality starts to deviate outside the expected range are crucial elements of the regulatory  approach. To achieve this SSM has stimulated the development and application of advanced monitoring and data analysis tools by operators as integral part of the regulatory framework (de Waal et al., 2012, 2015b). The need and realisation of these further advances in the regulatory approach are illustrated in this paper by the case history of gas production induced subsidence in the environmentally sensitive Wadden Sea area in the Netherlands.

2 The Wadden Sea case history

The Wadden Sea (Fig. 1) is a large temperate coastal wetland system behind a chain of coastal barrier islands. It is one of the world's most important wetlands and it is on the UNESCO world heritage list, based on its rich diversity, unique morphodynamical features and its wildlife values. It is one of the Netherlands most notable nature conservation areas protected under the European Birds and Habitats Directives. Gas production from fields underneath the Wadden Sea is nevertheless permitted, but only under very strict conditions.

Figure 1 The Dutch Wadden Sea is part of a semi-contiguous UNESCO world heritage site that spans the coastline of The Netherlands, Germany and Denmark. (© UNESCO; source: http://whc.unesco.org, last access: 28 March 2020).

3 Existing regulatory framework

The key question underpinning the applied approach for the Wadden Sea has always been: how much subsidence is acceptable and at which rate? And from the regulator perspective: how can it be reliably assured that (future) subsidence will stay within these limits? To address the issue the concept of “effective subsidence capacity” was developed (van Herk et al., 2010). It is the maximum human-induced subsidence that the affected area can robustly sustain. To determine this “effective subsidence capacity”, the maximum volumetric rate of relative sea level rise that can be accommodated in the long term, without damaging the natural values of the Wadden Sea area, was established first. The volume of sediment that can be transported and deposited by nature into the tidal basin where the subsidence occurs ultimately determines the “limit of acceptable average subsidence rate”. The capability of the tidal basins to “capture” sediment is the overall rate-determining step (Fig. 2). Effective subsidence capacity is then the maximum average subsidence rate available for human activities. It is obtained by subtracting the subsidence volume rate “consumed” by natural relative subsidence in the area (sea-level rise plus natural shallow compaction) from the total long-term acceptable subsidence volume rate limit.

Figure 2 When subsidence and sea level rise are balanced by (extra) sedimentation, the intertidal flats are maintained, and the natural values are likely preserved.

In the operational procedure for mining companies, six-years-average expectation values of subsidence rates are used to calculate the maximum allowable production rates. This is done under the provision that production will be reduced or halted if the expected or actual subsidence rate (natural + man induced) cannot be guaranteed not to exceed the limit of the acceptable subsidence rate now or in the future. The approach is known as the “Hand on the Tap” method. Monitoring and management schemes ensure that predicted and actual subsidence rates stay within the limit of acceptable subsidence rate and that no damage is caused to the protected nature. For further details see (van Herk et al., 2010) and (de Waal et al., 2012).

4 Challenges from observed subsidence delay

Initial uncertainties in the predicted subsidence (rate) of the Ameland field underneath the Wadden Sea were exacerbated by the observation of strongly non-linear (delayed) subsidence during later stages of the production period, see e.g. (Houtenbos, 2015). These were not accounted for in the original subsidence scenarios. Initially the operator assumed that the unexpected behaviour resulted from bi-linear elastic behaviour of the reservoir rock. Later, a continued more or less constant subsidence rate at much reduced reservoir depletion rates was observed, and it became clear that this is too limited an approach (NAM, 2015).

Figure 3 Prognosis for the deepest point on the Wadden Isle of Ameland over time. The current (2019) prognosis for the eventual subsidence is some 40 cm. During the first decades of production, subsidence prognoses were strongly adjusted after each observation survey (© NAM; source: NAM, 1998; Eysink, 1995).

It was realised that – depending on the underlying (combination of) physical cause(s) – the impact on future subsidence could be very different. In addition, it had become clear that the observed steepness of the edges of the subsidence bowl was much larger than predicted and that the apparent reservoir poromechanical compressibilities required to match the geodetic data substantially differed from laboratory compaction data (NAM, 2017). To assure proper management of the uncertainty, several improvements were requested by the regulator to be developed and implemented by the operator:

• Identification and in-depth study of possible mechanisms that could be the cause of the observed subsidence behaviour;

• Improved data assimilation techniques to test if field data can be used to assess which of the above mechanisms play(s) a role in the field;

• Testing of the newly developed approach against the 30-year-plus production/subsidence history of the Ameland gas field located underneath the Wadden Sea.

The project was coined LTS (Long Term Subsidence). The first phase (LTS-I) focussed on identifying possible and credible physical causes for the observed delayed subsidence. The second phase (LTS-II) focussed on applying this knowledge against an actual Wadden Sea field case (Ameland). A sophisticated data assimilation technique was developed and used to test if – using field and surveyance data – it can be determined which of the remaining post-LTS-I mechanisms factually plays a role in the observed delay in the subsidence above the fields underneath the Wadden Sea. Apart from addressing the three main objectives listed above, the LTS studies generated a number of valuable spin-offs:

• Software to manage large volumes of geodetic data;

• A procedure to combine GPS and benchmark data;

• Improved, more formal and objective methods to identify and manage outliers in the geodetic observations;

• An objective method to take noise in the geodetic data from shallow movements into account;

• The use of spatio-temporal double-differences and the full covariance matrix to confront subsidence predictions against geodetic survey data, eliminating the need for (assumed) stable reference points;

• A solution to prevent ensemble collapse from occurring while testing scenario's against field data using an ensemble-based data assimilation technique;

• An improved workflow to derive in-situ compaction from time-lapse GR-logs of signals from radioactive bullets shot into the formation.

Figure 4 The area around Ameland is measured using a variety of techniques, with different spatial and temporal coverage density (© NAM; source: NAM, 2017).

5 LTS study organisation

The study had to be carried out to the satisfaction of the Inspector General of Mines (the head of the State Supervision of Mines of the Netherlands). Its execution was an official condition for the approval of a number of Wadden Sea Field Development Plan updates. The aim of the study was to improve knowledge of the physical background of the measured time-dependent effects in the subsidence behaviour and its possible influence on expected subsidence in the long term. The study was to lead to a better understanding of the physical processes that explain the subsidence that had already occurred, with the aim of improving the forecasting of future subsidence. The study had to be based on generally accepted rules of physics, objective measurement data and proven scientific methods. The LTS study was carried out on the basis of a proposal by the operator: “Proposed Research Program for Higher Order Subsidence Modelling of The Netherlands Gas fields” (NAM, 2015). The proposal included an overview of the historical development of existing practices, their identified deficiencies and a proposed research program and deliverables for the LTS study. The document resulted from framing and brainstorm sessions within the operator organisation and with experts from SSM and the TNO AGE group (that exclusively advises SSM and the Ministry of Economic Affairs and Climate). Given the importance of the study, it was decided to establish an independent LTS Steering Committee (the SC LTS). The LTS study was carried out by the operator and external parties (public universities and research laboratories) with oversight provided by the SC LTS and SSM. The Wadden Sea Academy – established in 2008 with the task of providing a scientific basis for the management of the natural and social values of the Wadden Sea Region – was asked to create, facilitate and chair the SC LTS. The committee consisted of six internationally renowned experts in the fields of Geodesy, Experimental Rock Mechanics, Theoretical Rock Mechanics and Mining induced Subsidence. The SC LTS members were proposed partly by the operator and partly by the Wadden Sea Society (an independent NGO). They were subsequently also accepted by the Wadden Academy. Chairman and Technical Secretary for the committee were provided by the Wadden Academy. Representatives from SSM and TNO-AGE participated in the SC LTS as observers. The members met on a bi-annual basis between April 2013 and June 2015 to discuss progress and to provide guidance to the first phase of the LTS study. Between meetings, separate discipline meetings took place as well as telephone discussions and emails. Regular meetings were organised by the Wadden Academy to inform stakeholders from nature conservation and public organisations about the project progress and to address their questions and suggestions. For further details see (Wadden Academy, 2015).

At the end of the project a fully independent review was carried out at the request of SSM by an internationally renowned expert not involved in the LTS study (Teatini, 2017).

6 LTS Phase I

During the first phase of the study (LTS-I) the emphasis was on the identification and study of potential mechanisms that could explain the observed non-linear subsidence behaviour. The following six credible explanations were initially identified:

1. The non-linearity is not real but an artefact from noise and uncertainties in the data;

2. Salt flow in thick layers of rock salt overlaying the gas reservoir causes the observed delays;

3. Slow (delayed) depletion in underlaying and adjacent aquifers not captured in the modelling;

4. Non-equilibrium pressure diffusion within the reservoir during the initial production period;

5. Intrinsic non-linear and/or non-elastic in-situ reservoir rock compaction against pressure drop;

6. Collapse of high porosity reservoir rock intervals after reaching a large depletion of the initially strongly over-pressured gas in the reservoir.

The hypotheses were addressed in a number of separate discipline studies. In addition, several other possible influences on the surface subsidence were studied. These included different influence functions relating reservoir compaction to surface subsidence, effects of upscaling and the effects of a strongly heterogeneous and spatially variable overburden. The results of all studies carried out have been published in technical reports and scientific publications (NAM, 2015, 2017) and references therein. The outcomes are briefly summarised below and in Fig. 5.

Figure 5 Four quadrant summary of the LTS elements. Both models and measurements are assigned a functional and a stochastic description. The LTS study aims at confronting these. Images from © NAM; source: NAM (2017).

6.1 Noise and uncertainties

The studies carried out concluded that the observed non-linear subsidence behaviour is real and not an artefact of noise and/or uncertainties in the geodetic and other data. The hypothesis was therefore discarded. The studies did show that subsidence modelling precision can be improved significantly by taking correlation structures in the geodetic data into account. Improved methods were developed for objective outlier identification and handling, data reduction techniques for large geodetic data sets and for the processing and use of GPS data (Samiei-Esfahany and Bähr, 2015), (Williams, 2015). The work was supported by the development of a Bayesian framework to test and validate the quality of subsidence predictions against field measurements (Park and Bierman, 2015).

6.2 Salt flow

Extensive 1D and 3D modelling studies were carried out at the University of Utrecht to assess the potential effect of salt flow in the thick halite layers above the gas reservoirs in the Wadden Sea (Marketos et al., 2015). Results demonstrate that – on its own – such salt flow cannot explain the observed subsidence behaviour and in particular not its temporal behaviour. Also, an observed translation of the subsidence bowl over time is likely the result of changes in the reservoir depletion pattern over time and not a result of salt flow. Nevertheless, the studies demonstrate that salt flow can have an effect on the shape of the subsidence bowl and on the time evolution of subsidence, all very much depending on the values of the in-situ salt material properties. Note that salt flow can influence the shape of the subsidence bowl but not its volume.

6.3 Delayed aquifer depletion

Aquifer pressure depletion and its uncertainty range play a key role in the subsidence in the Wadden Sea. Important uncertainties are how well the aquifers are connected to the gas bearing intervals, the permeability in the aquifers and whether or not residual gas (significantly slowing down depletion rates) is present. Existing reservoir engineering models to calculate the effect of aquifer depletion were improved taking new production and well data into account. Results suggest that aquifer depletion is probably slower than originally assumed. Use of the reservoir engineering models under different parameter assumptions allows the calculation of alternative scenarios that can be tested against field production and subsidence data.

6.4 Anomalous pressure diffusion within the reservoir

The mechanism is based on the assumption that low permeability areas within the reservoir and possibly even at the microscopic pore level could result in a long-term (tens of years) non-equilibrium gas pressure distribution within the reervoir rock that would strengthen its effective macroscopic compressibility (Mossop, 2012, 2015). Although not proven impossible, the mechanism seems unlikely and could not be made more plausible on the basis of the study. The hypothesis was discarded.

6.5 Non-linear and/or non-elastic reservoir rock

Extensive long term (up to 12 weeks) laboratory compaction experiments under simulated in-situ conditions were carried out on Rotliegend reservoir rock samples: the gas-bearing rock in the Wadden Sea gas fields (Hol et al., 2015). The rock was sampled specifically for this study from the Wadden Sea Nes field. Results demonstrate that the compaction of the samples becomes increasingly inelastic (up to 80 %) at higher porosities and that the inelastic component remains considerable (some 50 %) at 20 % porosity (representative for the average in-situ reservoir rock). The laboratory derived inelasticity/creep turned out independent of temperature or the type of pore fluid over the ranges applied during the laboratory experiments. Overall, the compaction numbers measured during the laboratory measurements are comparable to those derived from field data. To capture the observed inelastic creep behaviour a rate-type constitutive model (Pruiksma et al., 2015) was used in the subsequent (LTS2) subsidence modelling to describe the constitutive behaviour of the Rotliegend reservoir rock. Note that the degree of non-linear reservoir rock compaction (with pressure drop) in first order does not influence the shape of the subsidence model while it has a large effect on the volume of the subsidence bowl at a given level of pressure depletion.

6.6 Pore collapse at high pressure depletion

The observed increase of subsidence rate would be due to “pore collapse” of weaker (high porosity) intervals in the reservoir at large amounts of pressure depletion during the later stages of production. Such behaviour has been observed during oil and gas production from chalk fields (Smits et al., 1988), and in the laboratory for high porosity samples from highly over-pressured sandstone gas fields (Schutjens and de Ruig, 1997). At high enough effective stresses the phenomenon is observed for all porous rock but usually only at stress levels far above the range that can occur in de Wadden Sea gas fields. The option could be excluded by long term laboratory experiments on high porosity Rotliegend reservoir samples at effective stresses significantly exceeding those that can occur in-situ. The hypothesis was discarded (Hol et al., 2015).

7 LTS Phase II

During the second phase of the LTS study the emphasis was on testing the relevance of the credible explanations that survived LTS-I in the real world and in particular for the observed subsidence delays in the Wadden Sea. A probabilistic Bayesian framework Ensemble based Subsidence Interpretation and Prediction tool (ESIP), was developed by TNO (Candela et al., 2017) to objectively confront a large set of different scenarios against geodetic observations by calculating the match of the model with the data expressed by a χ2 value. Next to delivering a probabilistic confrontation workflow, the tool also provides an objective statistical description of the outcome, i.e. an expectation case based on the weighted average and a 95 % confidence interval of the posterior model ensemble. TU Delft developed advanced geodetic processing software providing an interface between the geodetic data and the confrontation workflow (van Leijen et al., 2017).

The operator subsequently tested a somewhat modified version of the ESIP tool on the Ameland field underneath the Wadden Sea (NAM, 2017) as recommended by the SC LTS and as requested by SSM. The Ameland case was selected as the first test case given the long (30 year plus) production history and the large set of reservoir, production and geodetic survey data available for the field.

7.1 Scenarios tested

Some 58 history-matched pressure depletion scenarios were tested. Each meets the (gas and water) production data and the pressure data over the full Ameland field history while covering a wide range of different aquifer depletion distributions, with and without the presence of residual gas in the aquifers. For each of these pressure scenarios, parameter values of a generic RTICM subsidence model (Pruiksma et al., 2015) and (van Thienen-Visser et al., 2015) and a (moving) rigid basement influence function (Geertsma, 1966; van Opstal, 1973; Thienen-Visser and Fokker, 2017) were varied in a Monte Carlo simulation. The RTICM model was adapted to be able to handle the large amount of initial gas overpressure in the Ameland field. A time dependent shape factor was added to the influence function with a value dependent on the viscous behaviour of the salt layer above the reservoir. Each set of parameters drawn from the prior distributions during the Monte Carlo simulations results in a subsidence model member with the total set of members defining the ensemble. When confronted against the geodetic data (the full covariant matrix of the spatio-temporal double differences), the resulting fit of each member defines its probability and thereby its weight. The theory used to calculate the test metric accounts for the uncertainties in both the geodetic data and the geo-mechanical model and is based on a further development of (Nepveu et al., 2010). Modifications were successfully made to the applied goodness-of-fit metric to avoid ensemble collapse problems (Snyder et al., 2008) typical for particle filtering methods involving a large number of independent variables (9 independent parameters in the present study).

It is important to not only use the quality of fit of predictions against double difference geodetic data as some simple features of the misfit can then easily be missed. E.g. a relatively small lateral shift of the calculated subsidence bowl relative to its measured position will result in none of the subsidence members achieving a good fit against the double difference geodetic data and ESIP's discriminating capability will be lost. To identify and avoid such issues, visualisation and visual inspection of the predicted and measured subsidence (contours) remains important.

8 Results

Using the ESIP, some 20 000 subsidence members, covering the large range of possible parameter values were tested against the geodetic data for each of the 58 pressure scenarios resulting in a total of more than one million members.

Figure 6 Posterior scenarios within the 95 % confidence band stay below the long-term average subsidence rate limit (© NAM; source: NAM, 2017, Fig. 62).

Figure 7 The Wadden sea tidal flats just after sunset.

The following results were obtained:

1. Likely reservoir and aquifer depletion scenarios and parameter values for the subsidence model, salt flow and the influence function can be identified for Ameland;

2. Posterior parameter uncertainties are much reduced relative to the prior distributions as more geodetic data becomes available over time. In particular early measurement campaigns already significantly narrow the uncertainty for Ameland;

3. Posterior RTICM parameters derived are consistent with those derived from laboratory measurements (de Waal et al., 2015a);

4. The effect of lateral aquifers seems limited;

5. Salt flow alone cannot explain the observed time dependent (delayed) subsidence;

6. Extrapolations show that the likely scenarios within a 95 % confidence band will stay within the defined subsidence capacity including the longer term “acceptable average subsidence rate” limit (Fig. 6);

7. Emergency stop scenarios demonstrate the feasibility of the “Hand on the Tap” approach with its effectiveness obviously decreasing towards the end of field life.