History match of the Groningen field dynamic reservoir model to subsidence data and conventional subsurface data

Abstract

This paper presents the method applied to history-match the Groningen field dynamic model to conventional data (pressure data and water influx data) and to subsidence data, which is a novelty. Modelled subsidence is matched to subsidence data based on a simplified geomechanical model, which was built into the dynamic simulator.

A two-tier method was used to first achieve a match on a field-wide scale using field-average history-match quantifiers, which was then further improved at a regional/well level using regional history-match quantifiers. The history match was assisted by a space-filling experimental design. The resulting model has a field-average match to pressure data of ±2.17bar with a measurement uncertainty of ±0.4bar, to water influx data of ±2m with a measurement uncertainty of ±0.5m, and to subsidence data of ±4cm with a measurement uncertainty of ±1cm.

The output from this model is used as input for compaction, subsidence and production forecasts feeding into the hazard and risk assessment completed by NAM for the Groningen Winningsplan 2016.

Introduction

An update of the history matching of the Groningen field dynamic reservoir model was done as part of the 2015 Groningen Field Review (GFR2015) which was input to the Groningen Winningsplan update in 2016. Although the previous GFR2012 model vintage was still well suited for field development purposes, this update was instigated in light of the production-induced seismicity.

Results from the dynamic reservoir model are among the main inputs used to predict future subsidence. Subsidence is the surface imprint of reservoir compaction, which in turn is caused by pressure depletion due to gas production. Compaction is also believed to be the driving energy source for the seismicity observed in Groningen and therefore of special interest for this model update. Since direct measurements of compaction are sparse, being only available from five wells in the reservoir, the compaction dataset is insufficient to constrain the full reservoir model. Subsidence data are, however, readily available from levelling and satellite surveys and can thus be used to constrain pressure depletion. The previous model, GFR2012, had not been constrained to subsidence data. Calculated subsidence based on the results of that model resulted in notable local mismatches to observed subsidence. Aquifers adjacent to the field (e.g. underneath the city of Groningen) had not been incorporated in the GFR2012 model area. Pressure depletion is expected in these aquifers and therefore the model area has been extended to capture subsidence effects at the edges of the Groningen closure (Fig. 1). This required the inclusion of nine peripheral fields and their respective histories of pressure depletion.

Fig. 1. GFR2012 (left) and GFR2015 (right) dynamic simulation grid. Gas saturation in blue, and water in red.

The objective of the paper is to show how subsidence data can be used to constrain dynamic reservoir models additional to pressure and water ingress data. We explain the history-matching workflow applied to match multiple model realisations, composed of the Groningen field and peripheral fields, to the available subsurface data and the subsidence data. The final product of the workflow is a single best-matched model realisation. For forecasting, a workflow was used to generate a set of 1000 model realisations based on the uncertainty ranges as established in the history-matching methodology. From this set three models are selected representing the P10/P50/P90 realisations with respect to Ultimate Recovery. The applied forecast uncertainty workflow is beyond the scope of this paper. These models are used as the basis for forecasts of gas recovery, pressure distribution and subsidence.

Dynamic data

In total, 347 wells have been drilled in the Groningen field throughout its field life, and these have been subject to surveillance activities (measurements) over time. These data are used to calibrate the Groningen field dynamic model by minimising the mismatch of model output to field measurements; in previous model vintages reservoir pressure and water rise measurements at well locations have been used. The resulting dataset allows for good control on the dynamic behaviour of the reservoir near these wells. In areas where wells are sparse, temporal information to constrain the model is thus also scarce – this includes the north of the field and the aquifers surrounding the Groningen field.

Reservoir pressure data (SPG and RFT)

Reservoir pressure is routinely measured at the various production clusters (generally located towards the centre and the south of the field) and in the observation wells (typically located towards the periphery of the field). In total, some 1800 downhole static pressure gradient (SPG) measurements have been recorded to date in the Groningen field (see Fig. 2). Due to the typical accuracy of the gauges, multiple gauge types have been used; the uncertainty in the pressure measurements is ±0.4bar.

Fig. 2. Down-hole pressure measurements in the Groningen field over time converted to datum depth (2875m TVD NAP) in the gas phase, from SPG surveys.

For the Groningen and adjacent fields, 41 repeat formation tests (RFT) are available. RFT logging is done shortly after drilling the well and provides pressure measurements as a function of depth. Model output is generated to compare to both RFT and SPG measurements.

Aquifer ingress data

Aquifer ingress resulting from reservoir pressure depletion has been monitored in 30 wells by time-lapse Pulsed Neutron Log (PNL) surveys, with a total of 217 measurements to date. The maximally interpreted heights of change due to water saturation increase are collected in a dataset, representing water rise at well locations as a function of time. (PNL measures the capture cross section, also called sigma, which represents the ability of a given formation to capture thermal neutrons. Sigma will be affected by water movement, because chlorine in the salty water is the most effective neutron absorber in the reservoir. A gamma-ray tool connected to the PNL tool allows for depth correlation by log comparison. The water ingress height is determined by comparing the sigma log measurement of different surveys over time.) The interpreted depth of the gas–water contact is compared to model output. Due to the resolution of the gauge and resolution in the depth correlation the uncertainty in the interpreted height is ± 0.5m.

Subsidence data

Subsidence data are readily available from satellite and levelling surveys and provide the imprint of pressure depletion over time for the entire onshore part of the Groningen field and its surroundings, also in places where there is no well control. Subsidence has been measured throughout the production life of Groningen. It therefore provides additional information that has been used to constrain this update of the dynamic model.

The total dataset of subsidence measurements is filtered by the NAM (Nederlandse Aardolie Maatschappij B.V.) Geomatics department for stable benchmarks. Stable means the benchmark is both physically stable and the observations are free of human measurement and recording error. A benchmark is classified as stable when the estimated deviation from the trend obtained from proximal benchmarks is no greater than 20mm. Benchmarks may be excluded from the set, because abnormally fast settlement or uplift occurred (greater than 20mm with respect to a previous survey) caused by physical displacement of the benchmark (such as by heavy vehicles), or because the benchmark was accidentally reinstated at a wrong location after construction work. Furthermore, abnormally large deviations from the proximal observed trend (greater than 2mma⁻¹) typically related to shallow effects (e.g. soil compaction, solution salt mining, building settlement, groundwater level changes) may render benchmarks unstable. Of all benchmarks, 85% are classified as stable.

The first subsidence (levelling) survey that covered the full field was taken in 1972. Earlier surveys only partially covered the Groningen field, as it took some years after the initial discovery to establish the full extent of the field. The most recent full field (levelling) survey was acquired in 2013 (and the next repeat survey is currently planned for 2018). The measurement error over the 1973–2013 epoch is in the order of 1cm (NAM, 2014).

History-matching the dynamic model

Upscaling

Shell's in-house 3D numerical, fully implicit, finite volume reservoir modelling simulator is used to simulate the dynamic behaviour of the Groningen field and peripheral fields. The static geological subsurface pillar grid model (Visser & Solano Viota, 2017) was upscaled from 50,714,125 voxels (479×605×175) to 543,600 grid blocks (120×151×30) of which 317,896 are active. Each grid block measures approximately 400m by 400m, with an average thickness of 6m.

The static properties for the entire numerical grid (not just within the Groningen closure) after upscaling are as follows. Net-to-gross ranges from 0.006m³m⁻³ to 1m³m⁻³, with an average of 0.92m³m⁻³. Porosity is modelled compressible and compacts under pressure depletion; it initially ranges from 0.04m³m⁻³ to 0.27m³m⁻³, with an average of 0.14m³m⁻³. Horizontal permeability is modelled isotropic and ranges from 0.25mD to 1845mD, with an average of 122mD; vertical permeability is lower and ranges from 0.0014mD to 1043mD with an average of 33mD.

Initialisation

The model is divided over 45 initialisation regions of which 17 are within the Groningen field. Regions are based on different observed initial free water levels or observed pressure lags in regions separated by partly sealing faults. Within the Groningen closure, the temperature at datum depth varies and ranges from 87°C in the southeast to 116°C in the northwest with an average of 101°C. The model simulates a gas and a water phase, and the pressure–volume–temperature behaviour of the gas is based on the analysis of multiple fluid samples (Burkitov et al., 2016). Capillary pressure is modelled using Brooks-Corey and is based on saturation logs which compare to core flood experiment results. Relative permeability is modelled based on 19 special core experiments from core plugs collected across the Groningen and adjacent fields; no upscaling of these properties has been applied because the distribution of properties in Groningen is relatively homogeneous (Burkitov et al., 2016). Nine finite linear analytical aquifers are attached at the lateral boundaries of the numerical model.

Over the entire history, 347 wells have been active, including producers in the Groningen field and in peripheral fields, observation wells and injection wells. Alvestad inflow equations are used to model the wells.

Subsidence modelling

Reservoir compaction in each grid cell is calculated based on reservoir pressure and matrix compressibility. Matrix compressibility is modelled as a function of porosity based on a trend observed in core experiments, and ranges from 1.7×10⁻⁶bar⁻¹ to 1.9×10⁻⁵bar⁻¹, with an average of 4.3×10⁻⁶bar⁻¹. Using the Geertsma and van Opstal equation (Geertsma & van Opstal, 1972), the calculated reservoir compaction is converted to an areal grid of surface subsidence. (This method is set up as a proxy to constrain dynamic modeling, and significantly simplifies the overburden (Burkitov, 2016). For NAM's official subsidence predictions a more detailed analysis is done by the Geomechanics Department (van der Wal & van Eijs, 2016).) Reservoir compaction at any given location yields a dedicated bowl of surface subsidence, with its radius related to the reservoir depth. Hence, the ability to laterally resolve reservoir compaction based on surface subsidence is limited to the reservoir depth. In the discrete models used, the total subsidence at surface is calculated as the sum of the subsidence bowls associated with the compaction of each individual grid cell. The approach assumes a homogeneous half-space with a Poisson's ratio of 0.25. Both the modelled subsidence and the measured subsidence at stable benchmarks are projected on a large grid of approximately 4km by 4km which is used for visual comparison, and over which a root-mean-squared error (RMSE) is calculated. The assumption in the workflow is that all subsidence is caused by compaction in the Groningen field or adjacent fields only. Even though other sources impact subsidence in a small fraction of the surface area (e.g. solution salt mining), it is beyond the scope of this work to model these impacts. Despite the uncertainty in the modelled subsidence, the method does reveal spatial differences in subsidence related to differential pressure depletion in areas outside of well control.

History-matching methodology

The upscaled dynamic model was initialised in the dynamic simulator and run in history-match mode: the 60 years of historic production is withdrawn in monthly time steps from the individual wells. To improve the initial match to the four available historical data types, a method of history matching assisted by a space-filling experimental design was applied (Burkitov et al., 2016). This is a design of experiments (DOE) technique, which has the advantage of being very flexible and quick to implement (Wilson, 2015), and it uses full field models, instead of the typically used proxy models. The following six steps explain the process.

1. 1. Set-up of the variable model parameters. The multidisciplinary subsurface team determines the relevant variable model parameters that are uncertain and likely to impact the mismatch to data. This modelling exercise is by no means the first attempt to history-match Groningen, and in consequence the process was significantly steered by the history-matching parameters used in GFR2012. The applied method relies on continuous variables and is not suitable for categorical parameters, such as different static models, without repeating the entire workflow for every realisation of such a categorical parameter.

   Certain variable parameters were introduced to improve the match to data that had not been incorporated before, including parameters to match subsidence data. In total, 96 parameters are used in the history-matching workflow. While this is a large number, it allows for achieving a match in all regions of the field. Geologists, geomechanicists and petrophysicists have been continuously involved in constraining the parameter ranges within realistic bounds, to filter out unexplainable, unrealistic or inconsistent combinations of parameters. The set of variable parameters is selected to impact the history match for both the full field and the output of regions or wells and comprises the following types: gross bulk volume, permeability multiplier, fault seal factors, aquifer length, free water level, saturation functions and compressibility uncertainty. The most sensitive parameter varies per observable and has been determined using a tornado design. For subsidence, a fault separating the Groningen field from an adjacent aquifer is most impactful, for SPG the permeability multiplier in the Rotliegendes is most impactful, and for PNL the residual gas uncertainty is most impactful.

2. 2. Experimental design. The next step is to establish how (combinations of) settings for the variable model parameters can improve the match. Given the large set of variable parameters, certain experimental designs would require an unrealistically large number of simulations. A space-filling experimental design can be used, as the number of simulations required is a choice of the user. A modified Van Der Corput space-filling approach (Van Der Corput, 1935) was used to generate a set of 1000 different models. The Van Der Corput space-filling approach was modified to be used in multiple dimensions, instead of one, by randomly ordering the sampled parameter value of the first dimension for each higher dimension. Essentially, for each model realisation, every variable model parameter is sampled within its allocated uniform range, which is set symmetrically around its base case. The workflow results in models with parameters symmetrically around each base case parameter, but does not penalise deviations from the prior model base case parameter setting.

3. 3. Quantify mismatch. To quantify the mismatch between model output and the dynamic data (SPG, RFT, PNL and subsidence data), a RMSE is calculated for each data type. PNL and subsidence mismatches are not weighted differently. Pressure mismatches are weighted based on reliability, which is quantified 1 (poor) to 10 (good) depending on the tools used and the measurement protocol applied. Over the field life of Groningen, pressure gauges have improved in accuracy and the measurement protocols were improved, allowing longer stabilisation periods prior to taking pressure measurements. A field-average RMSE and local RMSEs per well, cluster and region are calculated.

4. 4. Model selection and improvement. Out of the set of models resulting from the space filling, a selection is made based on a low field average mismatch to subsidence, pressure and water influx. Because of practical reasons, the unfamiliarity with the subsidence model response and a time constraint to deliver a match, a graphical method was used to find an optimal solution, instead of an iterative method which converges to a minimum using a numerical solver. The three separate mismatches are plotted in 3D with RMSE values increasing away from the origin (see Fig. 3, left). The models closest to the origin are in theory the ones with the best match to the three constraining datasets.

   Although the set of optimal models towards the origin have a good overall history match, they will not necessarily match all individual wells. The models selected in Figure 3 (left) are further examined for their quality of match at a well level.

   To improve the match locally, local RMSE functions are compared to all variable model parameters to find possible correlations within the full set of model realisations. When there is an apparent correlation between a parameter and a local mismatch, it can be used to manually set the value of that parameter at its optimum value within its defined range, to obtain a model with a lower mismatch. An example is given in Figure 4, which shows a pressure match to the Harkstede-2A observation well. This well was drilled in a compartment separated from a nearby production cluster by a fault group, and the pressure mismatch was found to have a clear correlation with its associated fault seal factor. Consequently, the fault seal factor was fixed at 10^−1.2.

   As mentioned, subsidence was partly included in the workflow to improve pressure match in the periphery. An example that shows that an improved match to subsidence results in an improved match to pressure is given in Figure 5 for the Warffum-2b well. Two models are shown, where the only difference is the sealing behaviour of a fault separating the Warffum-2b well from the aquifer. Improving the pressure match improves the subsidence match at this location. This example does give some confidence that subsidence data can be used to constrain the model outside of well control, for instance in the aquifers to the west of the field where no well information is available.

5. 5. Improve definition of variable parameters. The set of variable model parameters defined at step 1 did not always deliver a match to the dynamic data locally. The occurrence of certain combinations of matches could be mutually exclusive in the ensemble of generated models and a match may not be possible without altering the variable parameter set-up. Inconsistent variable parameter settings are the result of effective variable parameter choices, based on engineering/ geoscience judgement, mostly related to uncertainty in fault transmissibility. Given a large number of over 1600 interpreted faults within the Groningen closure and a preference to limit the number of variable parameters, only a limited selection of faults is modelled to be sealing. When the selected set causes inconsistent results, revision of the effective parameters is done by the subsurface team, with the relevant disciplines involved. An example is illustrated in Figure 6 for the Ten Boer-4 observation well, which is separated by multiple faults from the nearby Eemskanaal cluster of production wells. Ten Boer-4 was drilled in a compartment which is lagging in pressure and where no water rise was observed historically. Certain faults were selected to vary in transmissibility near this well. For an ensemble of 1000 models, the pressure and PNL matches of the Ten Boer-4 observation well were found to be mutually exclusive (Fig. 6, left). A relatively high transmissible fault group improves the pressure match for this well because it improves pressure communication to Eemskanaal, although, conversely, higher transmissibility in this specific group of faults results in increased water influx which had not been measured by PNL surveys. With the subsurface team an alternative set of faults was set to vary in transmissibility, and in a repeat of the workflow, no inconsistent results occurred for this example (Figure 6, right).

6. 6. Repeat steps 2–5 until the match is consistent. When the set of variable model parameters is adjusted as per step 5, the workflow steps 2–5 are repeated until a consistent match has been achieved. A match has minimal local and global RMSE for different data types in the ensemble of models.

Fig. 3. Assisted History Matching workflow from the space-filling exercise towards the manual tuning.

Fig. 4. Local history-match improvements using the ensemble of 1000 models. (A) A prior model selected in Figure 3. (B) The correlation between a certain variable fault-group transmissibility and the pressure match of a local well. Highlighted are all models selected in Figure 3, (C) applying the best-match sealing factor.

Fig. 5. Pressure match and subsidence match in the northwestern part of the model.

Fig. 6. The local pressure and PNL match of the Ten Boer 4 observation well. (A) Initial variable fault group: no simultaneous match for both pressure and PNL. (B) Alternative variable fault group: PNL and pressure match converge to the origin when transmissibility is decreased.

Results

By applying the workflow as described in this paper, a best-matched model was established. Using the RMSE as a mismatch function, the following results were obtained.

• • The average mismatch to SPG pressures is ±2.17bar, with a data uncertainty of ±0.4bar and no significant trend of pressure mismatch over time (SGS Horizon, 2016), as shown in Figure 7 Early in the field life, the mismatch to pressure is higher. This is partly due to operational procedures, since the well shut-in duration prior to measurement, which directly impacts the measured pressure, was not recorded. Later in field life, the mismatch stabilises.

• • The average mismatch of the selected model to subsidence is ±4cm, with a data uncertainty of ±1cm. The results are shown in Figure 8. In the northeast of the field, the model predicts too much subsidence, while in the south it predicts too little. The best match is found in the west and in the central part.

  Residuals from the subsidence match have a clear areal trend (Fig. 8, right). This cannot be attributed to a mismatch in reservoir pressure, which has been matched accurately. Physical elements not included in the model could help explain the residuals in the subsidence. For instance, the southern mismatch could be related to pressure depletion of the Carboniferous formation underlying the Groningen reservoir. The resulting compaction could improve the match. The mismatch in the north of the field could be related to gas saturation trapped in the aquifer below the gas–water contact. The trapped gas would reduce the pressure drop in the aquifer and this would, in turn, reduce compaction, thus improving the over-prediction of subsidence in the north.

• • The mismatch to water rise is ±2m, with a data uncertainty of ±0.5m, which is considered to be a good match given the coarse vertical resolution of the dynamic model (SGS Horizon, 2016).

Fig. 7. Pressure mismatch as a function of time for all SPG data points; the colour indicates the shut-in time in days prior to the survey.

Fig. 8. Subsidence proxy output (A), measurement output (B) and the difference between model and measurement (C), all in cm. No subsidence measurements are available for the northern and eastern part of the grid, which lies offshore.

Simulation of a single model takes approximately 6 hours including the subsidence calculation.

Conclusions and recommendations

In light of the production-induced seismicity, the modelling of the Groningen reservoir is extended far beyond what is required in a conventional setting for making business decisions about the development of the field. A novel approach to history matching was introduced which assimilates subsidence data into a dynamic model by using a simplified geomechanical model. At well locations, a good match to subsidence data corresponds to a good match to pressure data; subsidence data can help constrain the model outside of well control. Residuals in the subsidence match could be caused by a depleting Carboniferous and residual gas in the aquifer.

The history-matching method applied has proved to be a relatively quick method to achieve a good match to multiple types of data. This allowed for assimilation of large amounts of available data within a relatively tight schedule of 3 months. It also proved useful in indicating inconsistencies in the definition of variable model parameters.

Future studies to further increase the dynamic understanding of the Groningen field and to improve the predictive modelling capability are outlined in the Study and Data Acquisition Plan (NAM, 2016). The dataset constraining the dynamic model will be enhanced by including static reservoir pressure measurements based on shut-in tubing head pressure data. Reservoir depletion and aquifer ingress can be constrained by including the historic gravity surveys, including the recent 2015 survey. Potentially, depletion of the underlying gas-bearing Carboniferous formation may explain the observations from the gravity data. A petrophysical study is being undertaken to investigate the occurrence of (residual) gas below the gas–water contact. Rock compressibility was populated in the dynamic model based on a trend observed in core experiments, by means of a single porosity transform throughout the entire grid. Areal trends that are not honoured by this transform will be captured by updating the compressibility grid using a model-based inversion of subsidence data and calculated reservoir pressure. This inversion will be done with the dedicated geomechanical model and will honour the full subsidence dataset instead of only using the 1972 and 2013 sets.