Original Article

Development of an efficient model to calculate subsidence above the Groningen gas field

Marius C. Wouters¹, Rob Govers¹ and Ramon F. Hanssen²

¹Tectonophysics Group, Department of Earth Sciences, Faculty of Geosciences, Utrecht University, Utrecht, The Netherlands; ²Department of Geoscience and Remote Sensing, Delft University of Technology, Delft, The Netherlands

Abstract

In order to constrain different drivers of subsidence in the Groningen gas field region, the integration of geomechanical simulations into a data assimilation procedure is crucial. Existing geomechanical models vary in complexity depending on their implementation of the available input data of the subsurface geometry and properties and reservoir pressure. High-complexity models are associated with many parameters to be estimated and tend to be computationally expensive, hindering their practical use in data assimilation. We develop a mechanical model that is optimised in terms of model complexity for the context of simulating surface deformation above the Groningen gas field. The reservoir discretisation and vertical elastic layering are simplified such that model details that are unlikely to be generating surface signals resolvable in geodetic data are eliminated. We demonstrate that the optimised model is ~100 times more numerically efficient than complete models. We also determine the sensitivity of subsidence to the lateral compaction resolution and the elastic layering of our efficient model, to constrain the model resolution in future data assimilation applications for Groningen.

Keywords: Subsidence; reservoir compaction, elastic layering; Groningen; The Netherlands

 

 

Introduction

Gas production from the large Groningen field (located in the north-east of the Netherlands) between 1963 and 2023 has induced land subsidence and seismicity, which are both continuing after production termination. In contrast to the extensive damages caused by the seismicity, production-induced subsidence has most likely not led to damage to houses (Geurts et al., 2023). Still, the subsidence presents sea defence and water management challenges, as the area is near (and partly even below) sea-level (De Waal et al., 2015). In addition, subsidence occurring below the UNESCO (United Nations Educational, Scientific and Cultural Organization) world heritage listed Wadden Sea could threaten the sensitive coastal wetland habitat (De Waal & Schouten, 2020; Fokker et al., 2018). Subsidence rates have been highest (~8 mm/year) above the central part of the field and decreasing towards the edges of the field (see Figure 1), which has resulted in a maximum cumulative subsidence that is currently over 40 cm (NAM, 2021).

Figure 1. Average InSAR velocities projected onto the vertical (Brouwer & Hanssen, 2023) in the Groningen gas field region for the period January 2015 – June 2020 (modified from Bodemdalingskaart 2.0 (2020)). The outline of the reservoir is shown by the thick black line, with up to ~8 mm/year above the centre of the field. Subsidence due to salt mining in Veendam and Winschoten (not part of this study) is also clearly visible.

Previous studies have focussed on the pressure depletion and resulting reservoir thinning (compaction) as the cause of subsidence above the field (e.g. Bierman et al., 2015; Fokker & Van Thienen-Visser, 2016; NAM, 2020b; Smith et al., 2019; Van Eijs & Van der Wal, 2017). However, significant subsidence is also observed in other parts of the Netherlands, away from hydrocarbon production or salt mining, driven by deformation of near-surface soil layers (clay and peat). In Groningen, subsidence resulting from these shallow processes is heterogeneously distributed over the region, with average rates of several mm/year locally (Fokker et al., 2018), potentially up to ~15 mm/year (Erkens et al., 2017). This shallow subsidence is only partially represented in conventional geodetic datasets, including the Interferometric Synthetic Aperture Radar (InSAR) data shown in Figure 1, because the detected signals originate mostly from deeply founded structures, which are insensitive to shallow deformation. Following periods of drought and precipitation, shallow subsidence signals can be up to ~10 cm in amplitude in peatlands in the Netherlands, as observed by in-situ extensometers and with the advanced distributed scatterer InSAR processing of Conroy et al. (2022, 2023).

Our overall project aims to use data assimilation techniques to disentangle these various processes from the geodetic (mostly InSAR) observations and constrain subsurface properties with associated uncertainties (see Hanssen et al., 2020; Kim & Vossepoel, 2023). Candela et al. (2022) find that disentanglement is theoretically possible. Our data assimilation approach will require physical models for all relevant processes. In the present study, we focus on the physical model representation of the response Groningen gas field region to reservoir compaction.

Detailed 3D descriptions of the stratigraphy, fault structure and material properties (e.g. De Jager & Visser, 2017; NAM, 2020c; Visser & Solano Viota, 2017) and high-resolution models of reservoir pressure evolution (Landman & Visser, 2023; Van Oeveren et al., 2017) in the Groningen subsurface are available. Existing geomechanical models have incorporated the details of this input data to a varying degree. The simplest models employ an analytical solution of a uniform elastic half-space to represent the subsurface (e.g. Bierman et al., 2015). Others incorporate information on reservoir pressures and material properties using either a half-space with a rigid basement (e.g. NAM, 2020a, 2020b; Van Eijs & Van der Wal, 2017) or a layered half-space (Fokker & Van Thienen-Visser, 2016). The most complex models use a fully 3D finite element approach (Guises et al., 2015; Lele et al., 2016). Whilst the simplest models potentially ignore subsurface aspects with significant surface expressions (e.g. the elastic profile, viscous behaviour of the caprock and strong lateral compaction variations not captured by low model resolution), the more detailed models could include many (small-scale) features without significant surface expressions. The high model resolution associated with the more detailed models leads to high computational costs, especially given the (many) forward calculations that are needed for data assimilation. Additionally, a large number of parameters to be estimated can lead to data assimilation problems like weight collapse or highly non-unique unrealistic solutions when most of the surface information is projected onto only a small subset of parameters.

We, therefore, develop a new model, with a complexity that is optimised for use with InSAR data in a data assimilation context. Here, we have two goals. The first goal is for the model to be able to reproduce the surface deformation resulting from reservoir compaction with an accuracy that suffices given the uncertainties of modern InSAR data. This results in a numerically efficient forward model. The second goal is to constrain the resolvability from InSAR data of lateral variations in reservoir compaction and vertical variations in the elastic response of the subsurface to compaction. The results will be applied in our follow-up study to restrict the number of parameters that need to be estimated by data assimilation.

First, we review information about the Groningen gas field region that is relevant for subsidence there. Then, we present the overall setup of our numerical model. Next, we show results for the optimal discretisation of the reservoir, sensitivity to lateral variations in reservoir compaction and sensitivity to layering of elastic properties. Finally, we describe the implications of our results for the planned compaction data assimilation study and discuss the impact of our model simplifications.

Subsidence in the Groningen gas field region

The Groningen gas field was discovered in 1959 and was in operation by the Nederlandse Aardolie Maatschappij (NAM) since 1963. At the time, it stood as the largest known gas field in the world (De Jager & Visser, 2017). The reservoir is roughly 30 × 30 km in size, with a geometry that has been mapped in great detail (NAM, 2020c; Figure 2). The top of the Rotliegend sandstone reservoir layer (Permian age) is located between 2.6 and 3.2 km. The gas-bearing reservoir layer, which varies in thickness between 0.15 and 0.25 km, is topped by the Zechstein evaporite caprock. Since the start of production, the reservoir pressure decreased from approximately 35 MPa to an average of 7 MPa by 2020 (NAM, 2020a). This pressure reduction has led to reservoir compaction, causing subsidence at the surface. Since 1964, periodic levelling surveys were conducted to monitor subsidence above the field (NAM, 2020b).

Figure 2. Example tessellation of the reservoir for TriaMax = 20 km². The original Groningen reservoir outline is shown in red, and the simplified outline and mesh are shown in black (ClipMax = 0.94 km²). Centres of dilatation are located in the centroids of mesh triangles. The mesh is finer along the reservoir edge to fill the shape. Plotted in the local Dutch RD coordinate system.

Subsidence monitoring efforts were expanded with InSAR measurements, starting in 1992 (Ketelaar, 2009; NAM, 2015). Additionally, since 2013, a network of permanent Global Navigation Satellite System (GNSS) stations has been installed above the field (Van der Marel, 2020). InSAR has become a leading technique for semi-continuous monitoring of surface deformation, particularly the PS-InSAR technique (Ferretti et al., 2001) based on coherent point scattering objects, like buildings and infrastructure. As many of these objects in the Netherlands are founded on relatively deep (Pleistocene) and stable layers (see Dheenathayalan, 2019; Ketelaar, 2009; Van Leijen et al., 2020), their velocities are predominantly driven by reservoir compaction. Scatterers with a shallow foundation also exhibit subsidence signals caused by soil deformation.

Little subsidence was observed during the first years of production (e.g. Schoonbeek, 1976; Van Eijs & Van der Wal, 2017). As the production increased in the early 1970s, the subsidence rate increased as well, albeit with a delay of ~3 years (Hettema et al., 2002). Multiple processes have been proposed to explain this delay in the onset of subsidence, including a rate-dependent inelastic reservoir response (De Waal, 1986; Hol et al., 2015; Pijnenburg et al., 2019; Pruiksma et al., 2015), a heterogeneous permeability distribution (Mossop, 2012) and flow of the viscoelastic evaporite caprock (Marketos et al., 2015). Since around 1980, the average subsidence rate above the centre of the field was ~8 mm/year (e.g. Van Eijs & Van der Wal, 2017).

Gas production varied greatly during the Groningen field life, even within seasons, as production peaked during winter (see e.g. Hettema et al., 2017). Pressure equilibration after these seasonal production swings occurred within half a year, whilst equilibration of regions along the reservoir edge may have taken a few years (Van Thienen-Visser et al., 2014; Van Wees et al., 2017). Following decades of impactful production-induced earthquakes and a systematic inadequate response to the earthquakes (Parliamentary Committee of Inquiry into Natural Gas Extraction in Groningen, 2023), production was sharply reduced in 2014 and terminated in 2023. Pijpers and Van Der Laan (2016) see a deceleration of subsidence rate in GNSS data approximately 9 weeks after the production reduction, although Qin et al. (2019) find no significant deceleration in InSAR data.

According to NAM (2021), the maximum total amount subsidence since the start of production is over 40 cm, just northwest of the centre of the field. NAM (2020a, 2020b) and Wildenborg et al. (2022) forecast subsidence to continue well into the 21^st century, accumulating an additional 10 cm. As the region above the field is near or below the sea level, surface subsidence can potentially lead to water management problems (e.g. groundwater rise and salination). To mitigate such problems, reliable subsidence forecasts are critical. However, past forecasts for the total subsidence have varied greatly (30–100 cm), due to large uncertainties in the field observations and variations in modelling choices (De Waal et al., 2015). Incorporating these uncertainties in a subsidence inversion and prognosis approach (e.g. by data assimilation) enables reporting model results of (future) subsidence with associated uncertainty ranges whilst also leading to a more comprehensive picture of the reservoir characteristics in general. Van Thienen-Visser and Fokker (2017) propose this approach as a future research topic, but currently, this has not been published. We plan to perform such a data assimilation study. A data assimilation approach involves many model realisations and thus requires an efficient and flexible model, which is the subject of this study. Whereas Fokker and Van Thienen-Visser (2016) and Van Thienen-Visser and Fokker (2017) used levelling data, we focus mostly on InSAR data, which could be sensitive to some near-surface subsidence signals too.

In 1973, Geertsma published an influential paper presenting a semi-analytical solution of the response of an elastic half-space to a depleting disc-shaped reservoir. Geertsma and Van Opstal (1973) adapted the method to a reservoir with an arbitrary shape and applied it to the Groningen field. Van Opstal (1974) added a rigid basement. These semi-analytical solutions have been heavily used in Groningen to estimate reservoir compaction parameters and the compaction distribution from subsidence data (e.g. Bierman et al., 2015; Smith et al., 2019; Van Eijs & Van der Wal, 2017). However, the Groningen subsurface is strongly layered with significant variations in elastic parameters and a prominent viscoelastic evaporite caprock layer. Fokker and Van Thienen-Visser (2016) used a layered version of Geertsma’s solution by Fokker and Orlic (2006) to find that the subsidence response of a model with an elastic layering based on the Groningen subsurface is significantly different from that of a uniform half-space. However, they did not examine the sensitivity of the surface deformation to the elastic contrasts, layer depths and layer thicknesses. In this study, we investigate this sensitivity to the elastic profile further, so we also develop layered models.

Numerical model and research approach

We model the volumetric strain associated with the depleting reservoir using the centre of dilatation approach, which has been developed in the context of surface deformation above expanding magma chambers (McTigue, 1987; Mogi, 1958), also referred to as a point source, Mogi source or nucleus of strain. The Mogi’s model consists of a pressurised spherical cavity in an elastic half-space. We model the centre of dilatation using the equivalent representation of a superposition of three orthogonal tensile dislocations. Given the appropriate amount of opening (or closing, in our case of a depleting reservoir) applied to the dislocations, this representation is identical to Mogi’s solution (Kumagai et al., 2014).

A single centre of dilatation accurately represents a volume undergoing volumetric strain (e.g. magma chamber or reservoir) when its radius is small compared to its depth. Consequently, we use a distribution of centres of dilatation to properly model the laterally extensive Groningen reservoir. The surface deformation is computed as a superposition (numerical integration) of the deformation resulting from individual centres of dilatation. Similar setups have been used in the context of magmatic intrusions (Vasco et al., 1988), depleting aquifers (Carlson et al., 2020) and hydrocarbon reservoirs (Du et al., 2005; Fokker & Van Thienen-Visser, 2016; Fokker et al., 2016; Geertsma & Van Opstal, 1973). To keep our analysis focused on the subsurface response to compaction, the model is driven kinematically, directly by volumetric strain representing the reservoir compaction. Our approach differs significantly from a dynamic approach, where reservoir compaction is driven by pressure depletion (e.g. Fokker & Van Thienen-Visser, 2016; Lele et al., 2016; NAM, 2020a, 2020b; Van Eijs & Van der Wal, 2017).

The outline of the Groningen reservoir is highly irregular in map view (Figure 2). We opt to discretise the reservoir with triangles, as they can be fitted easily to the reservoir shape. We use Triangle (Shewchuk, 1996) to generate high-quality Delaunay triangular tessellations (or meshes) of the reservoir. To complete the 3D geometry of the reservoir, the triangles are taken to represent the top of the reservoir by attributing the local reservoir thickness to them, that is, they are non-uniform triangular prisms. The centres of dilatation are positioned in the centroids of the prisms (Figure 2), with the dilatation magnitude representing the total volumetric strain of the prismatic volume. Selection of the number of centres of dilatation is discussed later.

We compute the response of a horizontally layered half-space to reservoir compaction using the linear viscoelastic semi-analytical code of Wang et al. (2006). This code is available online and has been widely cited and used by others, also in the context of hydrocarbon reservoirs (e.g. Fouladi Moghaddam et al., 2021; Vasco et al., 2017). Even though we focus on the elastic response to compaction in this study, follow-up studies could use other features of the code to expand on this work, like using viscoelastic behaviour in the reservoir and overburden or incorporating slip on major reservoir faults. Geoid and gravity solutions may help discriminate between shallow and deep subsidence. The code operates in two steps. The initial program, PSGRN, is designed to compute the Green’s functions for fundamental dislocations. In our case, we use the closing tensile dislocations to simulate centres of dilatation. The free surface is traction free, and displacements are zero at an infinite distance from the dislocations. The Green’s functions are used by the subsequent program, PSCMP, which calculates the surface deformation resulting from the distribution of centres of dilatation.

In line with the horizontal layering that is assumed in the model, we assume that the reservoir is situated between 2.8 and 3.0 km depth, with the centres of dilatation situated in the middle at a depth of 2.9 km. In reality, the reservoir is not flat, with depth variations on the order of 0.1 km throughout the field in reality (NAM, 2020c). These 0.1 km deviations from our modelled reservoir depth result in surface velocities that differ less than 0.12 mm/year, which we see as insignificant in this study. Initially, we apply a volume change to each centre of dilatation that is proportional to the size of the corresponding prismatic volume, such that the volumetric strain is uniform throughout the reservoir. The present study focusses on subsidence rates averaged over multiple years, so we interpret this volumetric strain as a strain rate and choose its value such that the resulting subsidence rate is representative for the approximate average subsidence rate above the Groningen field before the 2014 production reduction. Given the large lateral extent of the Groningen reservoir (~30 km) relative to its thickness (~0.2 km), the reservoir deforms mostly by uniaxial vertical compaction except near reservoir edges (e.g. Paullo Muñoz & Roehl, 2017; Tempone et al., 2010). We apply a uniform volumetric strain rate of 4.2 ∙ 10^-5 year^-1, which is equivalent to a compaction rate of 8.4 mm/year, which is on the high end of the range of observed compaction rates above the Groningen field (e.g. Van Eijs & Van der Wal, 2017). We also explore laterally varying compaction rates.

We assess the consequences of a variety of model simplifications by comparing the surface velocity results of different models. The velocity difference in direction i (i = eastward, northward, or up) at the j^th surface location r_j is computed from the surface velocities v¹ and v² of two models (1 and 2):

To determine whether the model simplifications have an imprint on velocity components above or below the observational uncertainty (see next section), we compute the maximum velocity difference in any direction and at any surface location:

We sample the modelled surface deformation every 1.2 km. We select this spacing such that the shape of the surface deformation is sufficiently well represented, which is essential for our comparison with unevenly spaced InSAR observations in our follow-up study. We verify that decreasing the spacing further does not affect our resolvability analysis results, as we show below using denser sampling. When actual InSAR observations are used in our follow-up study, modelled surface deformation can be computed directly on the unevenly distributed scatterer locations or by interpolation We verified our implementation of the PSGRN/PSCMP code against the response of a uniform elastic half-space to a pressure drop in a disc-shaped reservoir (Geertsma, 1973). We confirmed that differences with respect to Geertsma’s solution converged to the machine precision with increasing mesh density, and we concluded that our approach works properly.

Results

Observational uncertainty

Observational uncertainty plays a central role in our study. The uncertainty associated with estimated PS-InSAR velocities is defined by factors such as the quality of the scatterer, atmospheric and orbital perturbances, noise in the original radar observations, and choices in the kinematic and geometric parameterisation (Brouwer & Hanssen, 2023; Hanssen, 2001). These factors and choices influence the InSAR results in terms of dispersion and bias, with a particular emphasis on the estimation of the integer phase ambiguities. Ideally, the precision of InSAR is characterised by the stochastic model of the estimates (Van Leijen et al., 2020). Whilst the stochastic model is specific to each dataset (satellite mission, spatial and temporal domains, and processing choices), here, we aim to design a mechanical model that is optimised for application with a typical representative InSAR product over the Groningen study area in general. Thus, we opt for a pragmatic solution of a resolvability threshold of 0.5 mm/year, above which we see velocity signals to be potentially resolvable. We arrive at this value by assuming that the standard deviation of the estimate for the average velocities of a single InSAR scatterer, after projection onto the vertical (PoV) or onto the horizontals (PoH), is 0.18 mm/year. Furthermore, we assume the estimates to be normally distributed, and in this study, we regard the estimated velocities of different scatterers to be spatially uncorrelated. Consequently, with a confidence level of 95%, we assume that for velocity signals larger than 0.36 mm/year, the null-hypothesis of stability is rejected. Thus, with a discriminatory power of 80%, velocity signals greater than 0.5 mm/year are potentially resolvable by InSAR, which we refer to as the resolvability threshold.

Simplifying the reservoir discretisation

By using a dense mesh with closely spaced centres of dilatation, the surface deformation solution converges to that of a continuous reservoir, albeit at a significant computational cost. A coarse mesh introduces artefacts that become more prominent as the horizontal distance between centres of dilation increases, causing responses of individual centres of dilatation to become visible at the surface (Geertsma & Van Opstal, 1973). Artefacts are smaller for a deeper reservoir than for a shallower reservoir because the overburden response smooths out the reservoir signals. Some studies mention a rule of thumb, based on reservoir depth, for either the minimum lateral wavelength of the surface deformation (e.g. Smith et al., 2019) or the maximum horizontal distance between centres of dilation (e.g. Dusseault & Rothenburg, 2002). However, besides being approximate, these rules of thumb depend on the assumed or desired accuracy of the surface deformation. We therefore determine the required spacing of the centres of dilation ourselves based on the resolvability of the InSAR data. We first investigate what minimal interior mesh density is needed to avoid the artefacts exceeding the resolvability threshold of 0.5 mm/year. An example mesh is displayed in Figure 2.

This example also illustrates the need for the mesh density to increase towards the reservoir edges, as a consequence of fitting triangles to a complex reservoir shape. A simplified reservoir shape requires less mesh refinement towards the edge and is thus computationally more efficient. However, the surface deformation is altered as the reservoir shape is simplified. So, in the second step, we determine the maximum level of reservoir shape simplification, for which artefacts introduced by the simplification remain below the resolvability threshold.

In the third step, we combine our findings on the interior mesh density and shape simplification and present an efficient reservoir discretisation for computing the total surface deformation. In all three steps, we assume uniform elastic half-space properties. Deformation induced by any arbitrary dislocation source in an elastic half-space is only dependent on the Poisson’s ratio (c.f. equation 3.107 in Segall (2010)). Hence, the response to our applied volumetric strain rate also only depends on Poisson’s ratio. We use a Poisson’s ratio of 0.32 for our elastic half-space, which is the average of the Zeerijp-2 (ZRP-2) sonic well log (Romijn, 2017). Whilst a different choice for the Poisson’s ratio of the half-space, like the 0.25 of Bierman et al. (2015), would alter surface velocities slightly, it does not significantly affect our results for the resolvability of reservoir features. Later, we specifically investigate the influence of the elastic structure of the half-space.

Interior mesh density

Because we focus here on the influence of the interior mesh density, not the increased mesh density around a complex reservoir outline, we use a synthetic 30 × 30 km square reservoir here because no refinement around the edge is needed. As its dimensions agree roughly with those of the Groningen reservoir, we may expect semi-realistic surface velocities. The mesh density is controlled by the triangles sizes, which are limited by the maximum triangle surface area (TriaMax) set in the tessellation. This TriaMax is a useful metric in our context because the centres of dilatation represented by the largest triangles contribute most to the surface deformation and, therefore, also introduce the strongest surface artefacts. We compute the surface deformation for models with successively finer meshes, with TriaMax values ranging from 225 to 0.3 km². All model results are compared to the results of an ultrafine reference model. This reference model has a TriaMax of 0.025 km², and its results are insensitive to further refinement. We interpret the differences with the reference model as artefacts introduced by coarser meshes.

The middle row of Figure 3 shows the computed surface deformation for the ultrafine reference model. As the half-space deforms in response to the compacting reservoir, subsidence rates range from 9.4 mm/year above the centre of the reservoir, to 4.9 mm/year above the reservoir edges, to less than 0.5 mm/year beyond 10 km distance from the reservoir. This pattern is similar to the solution of Geertsma (1973) for a disc-shaped reservoir. He also saw a maximum subsidence signal being larger than the applied reservoir compaction (in our case 8.4 mm/year), as the reservoir is deflected downwards slightly. Horizontal motion is towards the centre of the reservoir and is most prominent above the reservoir edges, again similar to Geertsma’s solution. The results for a coarse mesh (TriaMax = 25 km²) in Figure 3 (top row) are similar to the reference model on the scale of the entire reservoir, but the imprint of individual centres of dilatation causes local deviations. Differences between the results for the coarse grid and the reference model (bottom row) are most prominent in the verticals and show the imprint of the individual centres of dilatation. The maxΔv of this coarse mesh density is 1.6 mm/year, which is significantly larger than the assumed resolvability threshold of 0.5 mm/year.

Figure 3. Results of models aimed at constraining the mesh density within the interior of the reservoir. The outline of the square reservoir is shown in black. The top row shows surface velocities for a course mesh (TriaMax = 20 km²). The middle row shows surface velocities for a model with an ultrafine mesh (TriaMax = 0.025 km²) that we use as a reference. The bottom row shows the Δv_i differences of the results of the coarse model relative to the ultrafine model.

Figure 4 shows the velocity differences as a function of TriaMax. The maxΔv values decrease with decreasing TriaMax and approximately follow a linear fit in our log-log plot, suggesting that maxΔv ≈ b TriaMax^a, with constants а ~ 1.4 and b ~ 0.026 (with maxΔv in mm/year and TriaMax in km²). Since the number of centres of dilatation (N) is inversely proportional to TriaMax, maxΔv decreases rapidly with increasing N (N^-1.4). The fit intersects the resolvability threshold at TriaMax = 8.7 ± 0.4 km². So, for an efficient computation of the surface response to compaction, a TriaMax < 8.7 km² in the interior of the reservoir suffices.

Figure 4. Largest surface velocity differences between models with a coarse mesh and the ultrafine reference model (maxΔv) as a function of TriaMax of the coarse meshes. The red-dashed line indicates a linear fit. The assumed resolvability threshold of InSAR observations is shown as the horizontal black-dashed line.

Avoiding refinement around complex reservoir outline

Applying the constraint of TriaMax < 8.7 km² to the actual Groningen reservoir shape (NAM, 2020c) leads to a mesh with over 24,500 triangles. This large number is mostly due to the intricate outline of the Groningen reservoir requiring many small triangles. However, compacting small-scale geometric features of the reservoir shape do not necessarily lead to resolvable signals at the surface. Moreover, the exact shape of these small-scale features is likely uncertain, particularly where the reservoir is thin. Therefore, we simplify the reservoir shape using the algorithm of Visvalingam and Whyatt (1993), which minimises the number of points on a curve (in our case the reservoir outline) whilst trying to retain shape. Points are eliminated based on the triangular area that is added or removed by their elimination. The algorithm uses a tolerance value, specified by the user, to set the simplification level. The red line in Figure 2 shows an example. Two types of regions emerge: regions that are added to the original reservoir by the simplification and regions that are removed from the original reservoir. These regions are found by clipping the original with the simplified outline polygon using the Pyclipper package by Johnson et al. (2022), which is based on Vatti (1992).

As the level of simplification increases, the size of both kinds of clipped regions increases too. Since the compaction distribution here is assumed to be uniform, the magnitude of volumetric strain in a reservoir region is directly proportional to the size of the region, and the largest clipped regions lead to the largest magnitude artefacts at the surface. We, therefore, use the size of the largest clipped region (ClipMax) as a metric for simplification. Note that this value is generally different from the tolerance value (see Figure 5). We compare the results for models with simplified reservoir shapes with a reference model based on the original reservoir shape from NAM (2020b) to determine maxΔv. To isolate the imprint of simplifications of the outline, we run these experiments with a fine interior mesh (TriaMax = 0.2 km²) in all models, including in the reference model.

Figure 5. Tolerance value applied in the Visvalingam & Whyatt algorithm for reservoir shape simplification and the resulting largest clipped region (ClipMax). Dashed line shows direct proportionality.

Figure 6 shows that surface velocities for the model with the simplified reservoir shape are similar to velocities for the reference model with the original reservoir shape. Unsurprisingly, the surface velocity differences between models with the simplified and original reservoir are largest, where the two shapes differ most. Along the north-western part of the reservoir, the vertical velocity differences are largest at 1.5 mm/year, which is above our assume resolvability threshold.

Figure 6. Results of models that focus on the complexity of the reservoir outline. The top row shows surface velocities for a simplified reservoir shape (ClipMax = 10.2 km²). The middle row shows velocities for the original reservoir shape. The bottom row shows the Δv_i differences between results of the simplified model relative to the original model. The original outline has a red colour, and the simplified outline is black. Plotted in the Dutch RD coordinate system.

Figure 7 summarises the results of model experiments where we varied ClipMax. Velocity differences between simplified and reference models increase with ClipMax. In this test with a relatively fine mesh, most clipped regions consist of multiple triangles and are thus represented by multiple centres of dilatation. Still, the surface impact of clipped regions that are modest in size (ClipMax < 10 km²) is similar to that of single centre of dilatation, where the total volumetric strain corresponding to the entire clipped region is applied to the single centre (see red dashed line in Figure 7). When the level of shape simplification is small (ClipMax < 0.02 km²), artefacts due to variations in the meshing of the interior dominate over the artefacts introduced by the shape simplification, and, thus, the maxΔv results deviate from the dashed red line. The intersection with the assumed resolvability threshold is around a ClipMax of 2.3 km². So, for an efficient computation of the surface response to compaction, a reservoir shape simplification level with a ClipMax < 2.3 km² suffices.

Figure 7. Largest surface velocity differences between models with a simplified reservoir shape and the reference model (maxΔv) as a function of ClipMax of the simplified reservoirs. For reference, the red-dashed line shows the largest surface velocity above a single centre of dilatation representing a reservoir region of size ClipMax. Results are shown for both the original sampling of the surface velocities with a 1.2 km spacing and for a denser sampling with a 0.6 km spacing.

In Figure 7, we also show the results for a denser surface deformation sampling, a spacing of 0.6 km instead of the 1.2 km used everywhere else, as this test turned out to be most affected by the denser sampling. Still, the results are practically identical, especially around the resolvability threshold, and our conclusions are unaffected by a denser sampling. So, we consider the 1.2 km spacing sufficient for the purpose of this study.

Optimised reservoir discretisation

Combining both the constraint on the mesh density (TriaMax ≤ 8.7 km²) and the constraint on the reservoir shape complexity (ClipMax ≤ 2.3 km²) leads to a mesh of just 196 triangles. However, comparing the results of this model to a model using the original reservoir shape and an ultrafine interior mesh density (TriaMax = 0.025 km²) shows differences up to 0.6 mm/year, which is above the resolvability threshold of 0.5 mm/year. The reason is that the two simplifications lead to a superposition of artefacts near the reservoir edge so that the net artefacts become too large. To avoid this, we aim to keep the artefacts introduced by the individual simplifications below 0.25 mm/year so that their combined imprint never exceeds 0.5 mm/year. Figure 4 shows that now the internal mesh density needs TriaMax ≤ 5.2 km², and Figure 7 shows that the shape simplification requires ClipMax ≤ 1.15 km². Note that the tolerance value of the Visvalingam–Whyatt algorithm does not directly correspond to ClipMax, and minor manual tuning is, therefore, required (see Figure 5).

A simplified model for Groningen constructed with a mesh using these constraints consists of 298 centres of dilatation. A model without these simplifications, with the original reservoir shape and a TriaMax of 0.25 km² (which leads to a spacing of centres of dilatation that is similar to the reservoir resolution of the operator’s geomechanical model), consists of 29,493 centres of dilatations. Given the linearity of our model, this two orders of magnitude reduction in model size leads to a two orders of magnitude gain in numerical efficiency (for details, see Table 1).

Table 1. Computation times for the Green’s functions (PSGRN) and for the surface deformation (PSCMP)

	Number of elastic layers	PSGRN runtime (s)	TriaMax (km2)	ClipMax (km2)	Number of centres	PSCMP runtime (s)
Uniform half-space, complex mesh	1	0.5	0.25	0	29,493	2,491
Complex layering, optimised mesh	31	13	5.2	1.15	298	23
Preferred layering, optimised mesh	6	3.5	5.2	1.15	298	23
Results are shown for a model with a uniform half-space, a model with the 31 layers based on ZRP-2, and a model with our preferred six layers. The first model uses a very complex mesh, with no reservoir shape simplification and a spacing of centres of dilatation that is similar to the reservoir resolution of the operator’s geomechanical model. The second and third models use the optimised reservoir discretisation. All computations are performed on an AMD EPYC 7451 processor.

Sensitivity of surface deformation to lateral variations in reservoir compaction

Reservoir compaction is not uniform throughout the reservoir. There are several reasons for lateral variations. First of all, reservoir compressibility (a material property) may vary up to a factor two throughout the reservoir (NAM, 2020b; Van Eijs & Van der Wal, 2017). In addition, the thickness of the gas-bearing layer shows strong lateral variations, especially across faults and towards the reservoir edges (De Jager & Visser, 2017). Furthermore, whilst the production strategy of equalising the reservoir pressure has led to a small spread in pressures between production clusters (Geurtsen & De Zeeuw, 2017; Van Oeveren et al., 2017), some areas near the edges of the reservoir, distant from the wells, have experienced only half the pressure drop of the overall field average (NAM, 2020b, 2021). And finally, the Loppersum clusters have seen a small pressure lag relative to the mean field pressure following an imposed production reduction in 2014 and eventual termination in 2018 (NAM, 2022). Previous studies indeed found strong variability in compaction (Bierman et al., 2015; De Zeeuw & Geurtsen, 2018; Smith et al., 2019; Van Thienen-Visser & Fokker, 2017). The result of Van Thienen-Visser and Fokker (2017) shows the strongest lateral variability, with regions of higher-than-average compaction directly bordering almost undeformed regions along some reservoir edges. Thus, we see a lateral compaction rate contrast of twice the reservoir average to be at the high end of what may realistically be expected.

In an inversion context (like data assimilation), the question arises: what lateral variations in the reservoir compaction rate are resolvable in surface deformation observations? We are particularly interested in the minimum length-scale in the reservoir, for which the compaction rate contrasts can be resolvable at the surface. To investigate this, we use models with reservoir compaction rates applied as a checkerboard pattern of alternating uniformly strongly and weakly compacting reservoir segments. We use a square reservoir and a fine mesh (TriaMax = 0.1 km²) so that we can also test detailed checkerboards, with small tiles. The half-space is again uniform elastic, with a Poisson’s ratio of 0.32.

Figure 8 shows the surface deformation for a checkerboard with 5 × 5 km tiles and alternating 16.8 mm/year (twice the reservoir average) and zero compaction rates. The differences with the results of the model with a uniform compaction rate represent the impact of the lateral compaction variations. Whilst the imprint of the checkerboard pattern is clearly visible at the surface, the pattern is smoother because of the elastic response (bending) of the over- and underburden. Here, maxΔv is 1.4 mm/year, which would be resolvable from the observations as it exceeds the assumed resolvability threshold.

Figure 8. Results of models that focus on compaction variations. The top row shows surface velocities for a model with a checkered compaction rate consisting of 25 km² tiles. The middle row shows velocities for a uniform compaction rate. The bottom row shows the Δv_i differences of the results for the checkered relative to the uniform compaction rate. Applied compaction rates are plotted on the right.

Figure 9 shows that maxΔv gets smaller with decreasing tile size as the smoothing effect becomes more influential. For the compaction rate contrast of 16.8 mm/year (blue line), the resolvability threshold is exceeded when the tiles are larger than ~12 km², that is, 3.4 × 3.4 km. In other words, the smallest lateral variation with a contrast of 16.8 mm/year that can be resolved is 3.4 km.

Figure 9. Largest surface velocity differences between models with checkered compaction and the uniform compaction model (maxΔv) as a function of the size of the checkerboard tiles. Both the results for the compaction rate contrast of 16.8 and 4.2 mm/year are plotted.

The grey line in Figure 9 shows maxΔv as a function of tile size for a smaller compaction rate contrast of 4.2 mm/year (half the reservoir average). The surface deformation exceeds the observational resolvability threshold for a tile size of ~36 km² (lateral resolution of 6.0 km). Figure 10 shows the smallest resolvable compaction rate contrasts as a function of the checkerboard tile size. As the tile size decreases, the magnitude of the contrasts needed for resolvable surface velocities increases. The rightmost point in the plot comes from an experiment, where tiles have the dimension of a quarter of the field. The result suggests that fieldwide compaction rate contrasts smaller than 1.3 mm/year are not resolvable from the observations, with smaller contrasts being indistinguishable from a situation with no lateral variations on this scale.

Figure 10. Smallest compaction rate contrasts between tiles for which the maxΔv with a uniform model exceeds the resolvability threshold as a function of tile size. The blue and grey squares correspond to the 16.8 and 4.2 mm/year contrasts of Figure 9, respectively.

Sensitivity of surface deformation to elastic layering

Elastic properties in the Groningen gas field region vary significantly with depth (Figure 11). Since the elastic response of a layered half-space can differ from that of a uniform half-space in both amplitude and pattern (e.g. Fokker & Orlic, 2006; Mehrabian & Abousleiman, 2015; Segall, 2010; Van Opstal, 1974), the surface velocities in Groningen are likely affected by the elastic structure too. Fokker and Van Thienen-Visser (2016) computed the surface deformation above a single centre of dilatation for both a uniform half-space and an elastic layering based on the Groningen subsurface and, indeed, found surface deformation to differ significantly between the two (up to ~100%). In preliminary tests using the Groningen reservoir shape and an elastic layering based on the ZRP-2 profile (Figure 11), we find surface velocities that differ more than 1 mm/year from a model using a uniform half-space. We also tested a model using a half-space with a sub-reservoir basement, with the basement’s Young’s modulus being 3 orders of magnitude higher than that of the top layer, analogous to the rigid basement of Van Opstal (1974). Regardless of the depth of the basement (3, 4 or 7 km), we found that surface velocities differ more than 1.7 mm/year from the layered model. So, for both the uniform and basement layer models, the difference with a layered model exceeds the uncertainty threshold of 0.5 mm/year, meaning that some information about elastic layering may be resolvable from the surface observations in a data assimilation procedure. Here, we aim to further explore the sensitivity of surface velocities to elastic layering in the Groningen gas field area. We seek to identify the depth and thickness of layers for which their elastic properties are potentially individually resolvable and build a preferred layering consisting of such layers.

Figure 11. Profiles of Young’s modulus (a) and Poisson’s ratio (b) from the ZRP-2 sonic well log data (Romijn, 2017). An example layering with six layers is plotted in red, for both the reference layering (solid) and a perturbed profile with the Young’s modulus of middle layer doubled (dashed). The blue line depicts our preferred layering. The arrows indicate that the base layers continue to infinity. Main stratigraphic units in the well are illustrated as coloured bands in the background. A slight vertical scaling is applied to the ZRP-2 data to align the middle of the Slochteren reservoir formation of the well log with the 2.9 km deep middle of the reservoir of our numerical model setup.

As thicknesses, depths and elasticities of the strata vary throughout the field, the elastic profile also varies laterally (Romijn, 2017). In contrast, elasticity in our approach of a layered half-space is uniform within a layer, also laterally. We thus seek an elastic layering that best represents the overall elastic structure of the entire field. Whilst the ZRP-2 profile is not necessarily exactly representative for entire field, we assume that an average elastic profile of the entire field is likely similar in shape. It is, therefore, reasonable to use the profile as a basis and search for elastic properties best representing the layer in a range around the profile. This range would be the same as that of the prior probability function in a data assimilation approach.

Here, we explore the effect of changing the elastic properties within this range, using reference profiles directly based on the ZRP-2 profile and perturbing the elastic properties of each individual layer. We perform a series of tests, each test with a different number of equally thick layers within the depth interval of the ZRP-2 profile. For each test, the reference profile is constructed by averaging the ZRP-2 data of Young’s modulus and Poisson’s ratio over the layer intervals. For lack of observations of elasticity variations deeper below the reservoir, we use a single base layer, with elastic properties corresponding to the bottom of the ZRP-2 profile. Because of our layered half-space approach, this base layer is infinitely thick. For example, the solid red lines in Figure 11 represent the elastic profile for five equally thick layers above the base layer, that is, a 6-layer model.

The dotted line in Figure 11 shows one of the perturbed profiles based on the 6-layer model, with an increased Young’s modulus of the middle layer. The applied perturbations are based on the range of elastic property values throughout single stratigraphic units in the field-wide velocity model of Romijn (2017). We perturb each layer’s Young’s modulus by doubling and halving the reference value, and we perturb the Poisson’s ratio by adding and subtracting 0.1 from the reference value, so running four perturbed models per layer. For the 6-layer profile, this amounts to 24 tested models. All models use a high mesh density (TriaMax = 0.2 km²) and little reservoir shape simplification (ClipMax = 0.015 km²) to avoid numerical artefacts > 0.01 mm/year (see Figure 7).

Surface velocities resulting from the reference and perturbed 6-layer models of Figure 11 are shown in Figure 12. This specific perturbation of the Young’s modulus leads to a decrease in horizontal velocity magnitudes and a slight widening of the subsidence bowl. Maximum differences between these 6-layer models are 0.77 mm/year, which is above the 0.5 mm/year threshold, so we deem this specific perturbation resolvable. Whilst the ratio between horizontal and vertical velocities is relatively constant in all previous models (see Figures 3, 6 & 8), we see here (Figure 12), and in general for most elasticity perturbations, that the perturbations lead to a significant change of the ratio between horizontal and vertical velocities. This demonstrates that this ratio contains information regarding the elasticity structure.

Figure 12. Results of models that focus on elasticity layering. The top row shows surface velocity components for the perturbed 6-layer elasticity profile shown by the dashed line in Figure 11. The middle row shows surface velocities for the 6-layer reference model with elasticities of the solid red line in Figure 11. The bottom row shows the Δv_i differences of the perturbed relative to the reference model.

Figure 13a summarises the results for the perturbations of the 6-layer model, showing the largest maxΔv for each layer. Young’s modulus perturbations are resolvable in all layers, but the top two with the bottom two layers (base layer and layer containing the reservoir) lead to the strongest signals. Perturbing the Poisson’s ratio instead of the Young’s modulus generally results in smaller maxΔv values, with only the bottom two layers being resolvable here.

We perform this procedure for profiles with two through 36 sublayers. Figures 13b and 13c summarise the results of all these tests, for perturbations of Young’s modulus and Poisson’s ratio, focussing on the resolvability per layer. On the right of the graphs, the results for models with two layers above the base layer are shown. All dots are filled with black, indicating that elasticity perturbations are resolvable. Red dots further to the left highlight the results for the models of Figure 13a. The models with more (thinner) layers are plotted to the left, showing generally even lower resolvability. We find that the base layer is resolvable, regardless of the number of overlying elastic layers.

Figure 13. (a) The maxΔv values for the six layers of the example profile shown in Figure 11. (b, c) Resolvability as a function of depth and thickness of the perturbed layer, for perturbations of the Young’s modulus (b) and Poisson’s ratio (c). The models corresponding to (a) are plotted in red. Approximate boundaries between resolvable and not resolvable models are drawn as black dashes. Results for the infinitely thick base layer are not shown. In all tested profiles, the base layer perturbations are resolvable. The layer depths and thicknesses of our preferred layerings are plotted as blue crosses.

The dashed black lines in Figures 13b and 13c approximately separate non-resolvable thin layers from resolvable thicker layers. They show that elastic properties of the top layer can only be resolved from observations when it is 1.5 km thick, whilst at depth, thinner layers can also be resolved. Especially around the reservoir depth, we see thin (<0.1 km) layers that are resolvable. This sensitivity to the reservoir elasticity originates from the response to the applied volumetric strain there. When the elastic properties of the reservoir layer are perturbed, the applied strain is kept constant. So, the stresses in the reservoir and connected over- and underburden, induced by the applied strain, change proportionally to the elasticity perturbation. With this significant change of stress state in the half-space, the deformation pattern also differs significantly.

For some layers around the middle of the profile, Poisson’s ratio perturbations are not resolvable, whilst Young’s modulus perturbations are. In addition, some isolated thin layers are resolvable, whilst similar layers just above or below are not. These thin layers lie around depths, where the elastic properties in the ZRP-2 data (Figure 11) exhibit sharp contrasts: for example, around 1.6 and 2.1 km for the Young’s modulus and 0.8 km for the Poisson’s ratio. As this seems to be related to the detailed shape of the ZRP-2 profile, it is likely not representative of the resolvability of elastic layering in the field in general.

We use the results of Figures 13b and 13c to construct our preferred layering (shown in Figure 11), with layers that are all individually resolvable (blue crosses in Figures 13b and 13c lie to the right of the dashed lines). This layering consists of six layers for the Young’s modulus, and only 5 distinct Poisson’s ratio values because of its lower resolvability. The Poisson’s ratio value for the second and third layer from the top is set as identical. The top layer is 1.5 km thick, with the layers below decreasing in thickness towards the reservoir. The vertical resolution of this preferred profile is similar to that of NAM (2013) and Fokker and Van Thienen-Visser (2016), with both studies using seven elastic layers.

Computational efficiency of our proposed model

Our workflow consists of two steps: Green’s functions of the layered half-space response to a centre of dislocation are computed by the PSGRN program, and the deformation of the 2,601 points in our employed surface grid resulting from the collection of centres of dislocation representing the reservoir is computed by the PSCMP program. So, whilst the PSCMP program is run for each new distribution of centres of dilatation, the PSGRN does not need to be rerun as long as the elastic layering and elastic properties remain the same.

The Green’s function computation time increases with the number of elastic layers (Table 1). For our preferred six layers (Figure 11), PSGRN computation is almost four times as fast as for a model using all 31 stratigraphic (sub)layers identified in the ZRP-2 well. The PSCMP computation time of the surface velocities is linearly related to the number of centres of dilatation in the reservoir. Our reservoir discretisation simplifications lead to a gain of two orders of magnitude in computational efficiency. Given the linearity of the problem, the workload can be divided between different nodes, to speed up computation even more. We conclude that our preferred model is efficient. It has 11 elastic parameters and 294 centres of dilatation (i.e. 305 parameters in total) that need to be estimated by the data assimilation procedure.

Discussion

Implications for a compaction inversion

In a follow-up study, we plan to use InSAR data to invert for the reservoir compaction evolution in Groningen using data assimilation. We propose to use our optimised reservoir discretisation for the first version of the numerical model used in this inversion. If the results of the first inversion suggest that compaction behaviour within some sections of the reservoir is relatively uniform, our findings of Figure 10 can be used to divide the reservoir into groups of centres of dilatation (compartments), where the applied compaction is uniform within each compartment (like the checkerboard compartments). This reduces the number of compaction parameters that need to be estimated in the data assimilation.

Previous inversions for reservoir compaction and compressibility have produced results with a range of lateral reservoir resolutions. Where the inversions of Van Eijs and Van der Wal (2017) and Smith et al. (2019) use a 1 km lateral resolution, Fokker and Van Thienen-Visser (2016) involves a 3.2 km resolution and Bierman et al. (2015) report results both at a 2.5 km and a 5 km resolution. For reference, a 5 km lateral resolution corresponds to ~50 reservoir compartments, which should be individually resolvable when the compaction rate contrasts between the compartments are on the order of 6 mm/year (see Figure 10).

In addition to inverting for compaction, we also plan to invert for the elastic profile. We propose to use our preferred layering (blue lines in Figure 11) as the base (prior) and let the inversion modify the elastic properties per layer to estimate the profile best representing the Groningen subsurface (posterior). Our preferred layering consists of six layers (5 distinct Poisson’s ratio values). If actual inversion results would suggest a low sensitivity to certain layers, or if elastic properties of different layers would be estimated to be similar, the number of layers could be decreased further, reducing the Green’s function computation time. More layers are likely not required.

Impact of model simplifications

Our model is driven kinematically (directly by reservoir compaction), as opposed to a model that is driven dynamically (by pressure depletion). A dynamic approach requires a choice for a reservoir rheology, for the relationship between pressure depletion and compaction. In addition, a dynamic approach requires model inputs such as distribution of rheological properties (e.g. uniaxial compaction coefficient; NAM, 2020b) and reservoir pressures (Landman & Visser, 2023; Van Oeveren et al., 2017) to be selected. A kinematic approach is unaffected by possible biases and uncertainties in these model inputs. However, a downside of a kinematic approach is the insensitivity to the Young’s modulus magnitudes, both in a uniform (c.f. equation 3.107 in Segall (2010)) and a layered elastic half-space, with the latter only being dependent on the shape of the profile, not the actual Young’s modulus magnitudes. So, even though we have used the ZRP-2 well data as a basis for our layering, scaling the Young’s modulus values of each layer by the same constant results in an elastic response that is identical to the original layering. This means that in order to convert our current model to a dynamic model, an additional inversion step to obtain information on the reservoir Young’s modulus magnitude is required.

In this study, we have used the outline of the reservoir to delimit the extent of compaction. However, other studies have attributed a portion of the observed surface deformation to compaction outside of the reservoir, from pressure depletion in laterally connected aquifers (e.g. NAM, 2020b; Van Thienen-Visser & Fokker, 2017). If we would want to study the potential activity of connected aquifers, our distribution of centres of dilatation could be extended beyond the reservoir outline. As the compaction rates in the aquifers are most likely low, lateral contrasts are likely small too, and, thus, a relatively low lateral compaction resolution in the aquifers would probably suffice.

Another model complexity to potentially analyse in a future study is the impact of viscous behaviour, specifically salt creep (e.g. Spiers et al., 1990) in the thick Zechstein evaporite caprock. As the PSGRN/PSCMP code is capable of simulating layers with a Burgers’ rheology, both transient creep (Kelvin-Voigt) and steady-state creep (Maxwell) can be considered for the caprock, similar to Marketos et al. (2016).

Choice for resolvability threshold

We evaluate the potential resolvability of signals in InSAR data with the simplified resolvability threshold approach. We base the threshold value on an estimate for the precision of velocities of single scatterers (from measurement and processing noise). Here, we assume that velocities of different scatterers are spatially uncorrelated. However, in reality, transmission of compaction signals trough the overburden causes spatial smoothing of the surface signals (regardless of the exact elastic structure, reservoir geometry and compaction). Integrating information of these spatial correlations could improve precision of the velocity estimates, potentially improving resolvability. Additionally, whilst we simulate a single velocity result per location, estimating single linear trends per scatter from the InSAR data may lead to relevant temporal variations being missed, especially for a long dataset. Hence, we will estimate variable trends or piece-wise velocities, for which the precision could be lower than for single linear trends. So, for our follow-up inversion study, we propose for the resolvability and model complexity to be re-evaluated considering the specific (stochastic) characteristics of the selected InSAR dataset.

Conclusions

Here, we present a geomechanical model optimised for the context of simulating surface velocities due to reservoir compaction of the Groningen gas field. We apply simplifications of both the reservoir shape and the internal meshing of the reservoir without significantly affecting the calculated surface velocities, keeping the numerical artefacts introduced by the simplifications below the assumed resolvability threshold. Compared to a model without simplification of the reservoir shape and with an internal meshing similar to the operator’s geomechanical model, the applied simplifications enable a two orders of magnitude reduction in model complexity and, consequently, a two orders of magnitude gain in computational efficiency.

Resolvability of reservoir compaction in geodetic data depends both on the lateral reservoir resolution of the model and on the degree of lateral compaction variability. Checkerboard tests show that lateral compaction contrasts of 16.8 mm/year in magnitude can be resolved when the lateral length scale of the variations is at least 3.4 km. Resolving features within a smoother compaction rate distribution is more challenging. For contrasts of 4.2 mm/year in magnitude, variations need to be at least 6.0 km in size. In an inversion study, the expected level of compaction variability may not be known a priori. Therefore, we propose to carry out an initial inversion where compaction is estimated at a high lateral resolution. If this initial inversion indicates that the compaction behaviour within some reservoir sections is relatively uniform, the resolution of these sections can be decreased, guided by our findings on the resolvability of lateral compaction variations.

The elastic profile of the Groningen subsurface impacts both the pattern and magnitude of the surface velocity response to reservoir compaction. The results of our layered models show that the surface velocities are more sensitive to the Young’s modulus values of the layers than to their Poisson’s ratios. Sensitivity is highest for layers close to the reservoir depth, where layers <0.1 km thick are resolvable, and lowest for the top layer, which needs to be 1.5 km thick to be resolvable. We introduce a preferred elastic layering consisting of six layers (5 distinct Poisson’s ratio values), with thicknesses such that all layers are potentially individually resolvable at the surface. More layers are likely not required.

Acknowledgements

This work was funded by the Dutch Research Council (NWO) grant NWO DEEP.NL.2018.052. All map figures, except for Figure 1, were generated using Generic Mapping Tools version 6 (Wessel et al., 2019). Input, output and scripting files that were used for this paper are digitally stored in the Yoda repository of the Utrecht University (https://doi.org/10.24416/UU01-W3PCRY). We thank Lukas van der Wiel for his continuous numerical support and Rob van Eijs and an anonymous reviewer for their incisive reviews. Author contributions are following the CReDiT taxonomy: Marius Wouters: Conceptualisation (equal), Methodology (lead), Software (lead), Investigation (lead), Data curation (lead), Writing – Original draft (lead), Writing – Review & Editing (equal), Visualisation (lead); Rob Govers: Conceptualisation (equal), Methodology (supporting), Software (supporting), Resources (lead), Writing – Review & Editing (equal), Supervision (lead), Project administration (lead), Funding acquisition (supporting); Ramon Hanssen: Conceptualisation (supporting), Methodology (supporting), Writing – Review & Editing (supporting), Supervision (supporting), Project administration (supporting), Funding acquisition (lead).

References

Bierman, S.M., Kraaijeveld, F. & Bourne, S.J., 2015. Regularised direct inversion to compaction in the Groningen reservoir using measurements from optical leveling campaigns. Shell Global Solutions International B.V. https://nam-feitenencijfers.data-app.nl/download/rapport/cc5ea278-c093-457b-b930-1869a3c26c21?open=true.

Bodemdalingskaart 2.0, 2020. Bodemdalingskaart.nl.

Brouwer, W.S. & Hanssen, R.F., 2023. A treatise on InSAR geometry and 3-D displacement estimation. IEEE Transactions on Geoscience and Remote Sensing 61: 1–11. DOI: 10.1109/TGRS.2023.3322595.

Candela, T., Chitu, A.G., Peters, E., Pluymaekers, M., Hegen, D., Koster, K. & Fokker, P.A., 2022. Subsidence induced by gas extraction: a data assimilation framework to constrain the driving rock compaction process at depth. Frontiers in Earth Science 10: 1–18. DOI: 10.3389/feart.2022.713273.

Carlson, G., Shirzaei, M., Ojha, C. & Werth, S., 2020. Subsidence-derived volumetric strain models for mapping extensional fissures and constraining rock mechanical properties in the San Joaquin Valley, California. Journal of Geophysical Research: Solid Earth 125(9): 1–16. DOI: 10.1029/2020JB019980.

Conroy, P., Van Diepen, S.A.N., Van Asselen, S., Erkens, G., Van Leijen, F.J. & Hanssen, R.F., 2022. Probabilistic Estimation of InSAR Displacement Phase Guided by Contextual Information and Artificial Intelligence. IEEE Transactions on Geoscience and Remote Sensing 60: 1–11. DOI: 10.1109/TGRS.2022.3203872.

Conroy, P., Van Diepen, S.A.N., Van Leijen, F.J. & Hanssen, R.F., 2023. Bridging Loss-of-Lock in InSAR Time Series of Distributed Scatterers. IEEE Transactions on Geoscience and Remote Sensing 61: 1–11. DOI: 10.1109/TGRS.2023.3329967.

De Jager, J. & Visser, C., 2017. Geology of the Groningen field – an overview. Geologie En Mijnbouw/Netherlands Journal of Geosciences 96(5): s3–s15. DOI: 10.1017/njg.2017.22.

De Waal, J.A., 1986. On the rate type compaction behaviour of sandstone reservoir rock (Delft University of Technology). Delft University of Technology. https://repository.tudelft.nl/islandora/object/uuid%3Ab805782b-2eb4-4f72-98f4-f727c4ea9df0.

De Waal, J.A., Muntendam-Bos, A.G. & Roest, J.P.A., 2015. Production induced subsidence and seismicity in the Groningen gas field-can it be managed? Proceedings of the International Association of Hydrological Sciences 372: 129–139. DOI: 10.5194/piahs-372-129-2015.

De Waal, J.A. & Schouten, M.W., 2020. Regulating subsidence and its uncertainty in the Dutch Wadden Sea. Proceedings of the International Association of Hydrological Sciences 382: 63–70. DOI: 10.5194/piahs-382-63-2020.

De Zeeuw, Q. & Geurtsen, L., 2018. Groningen dynamic model update 2018 – V6. NAM. https://nam-onderzoeksrapporten.data-app.nl/reports/download/groningen/en/b8a36d79-118e-4c2b-8380-ea338c9aafc8.

Dheenathayalan, P., 2019. Optimizing the exploitation of persistent scatterers in satellite radar interferometry. Delft University of Technology. DOI: 10.4233/uuid:aa1ef96f-4da9-41ff-bff8-30186ef2a541.

Du, J., Brissenden, S.J., McGillivray, P., Bourne, S., Hofstra, P., Davis, E.J., Roadarmel, W.H., Wolhart, S.L., Marsic, S., Gusek, R.W. & Wright, C.A., 2005. Mapping reservoir volume changes during cyclic steam stimulation using tiltmeter-based surface-deformation measurements. SPE International Thermal Operations and Heavy Oil Symposium. (Calgary, Alberta): 1–12. DOI: 10.2118/97848-PA.

Dusseault, M.B. & Rothenburg, L., 2002. Analysis of deformation measurements for reservoir management. Oil & Gas Science and Technology 57(5): 539–554. DOI: 10.2516/ogst:2002036.

Erkens, G., Stafleu, J. & Van den Akker, J., 2017. Bodemdalingvoorspellingskaarten van Nederland, versie 2017. Deltares rapport klimaateffectatlas (Utrecht).

Ferretti, A., Prati, C. & Rocca, F., 2001. Permanent scatterers in SAR interferometry. IEEE Transactions on Geoscience and Remote Sensing 39(1): 8–20. DOI: 10.1109/36.898661.

Fokker, P.A. & Orlic, B., 2006. Semi-analytic modelling of subsidence. Mathematical Geology 38(5): 565–589. DOI: 10.1007/s11004-006-9034-z.

Fokker, P.A., Van Leijen, F.J., Orlic, B., Van der Marel, H. & Hanssen, R.F., 2018. Subsidence in the Dutch Wadden Sea. Geologie En Mijnbouw/Netherlands Journal of Geosciences 97(3): 129–181. DOI: 10.1017/njg.2018.9.

Fokker, P.A. & Van Thienen-Visser, K., 2016. Inversion of double-difference measurements from optical levelling for the Groningen gas field. International Journal of Applied Earth Observation and Geoinformation 49: 1–9. DOI: 10.1016/j.jag.2016.01.004.

Fokker, P.A., Wassing, B.B.T., Van Leijen, F.J., Hanssen, R.F. & Nieuwland, D.A., 2016. Application of an ensemble smoother with multiple data assimilation to the Bergermeer gas field, using PS-InSAR. Geomechanics for Energy and the Environment 5: 16–28. DOI: 10.1016/j.gete.2015.11.003.

Fouladi Moghaddam, N., Nourollah, H., Vasco, D.W., Samsonov, S.V. & Rüdiger, C., 2021. Interferometric SAR modelling of near surface data to improve geological model in the Surat Basin, Australia. Journal of Applied Geophysics 194: 1–17. DOI: 10.1016/j.jappgeo.2021.104444.

Geertsma, J., 1973. Land subsidence above compacting oil and gas reservoirs. Journal of Petroleum Technology 25: 734–744. DOI: 10.2118/3730-PA.

Geertsma, J. & Van Opstal, G., 1973. A numerical technique for predicting subsidence above compacting reservoirs, based on the nucleus of strain concept. Verhandelingen van Het Koninklijk Nederlands Geologisch Mijnbouwkundig Genootschap 28: 63–78.

Geurts, C.P.W., Pluymaekers, M.P.D., Rots, J.G., Korswagen, P.A. & TU Delft, 2023. Response to the review comments: summarizing report on the additional studies into the direct effects of deep subsidence. TNO. https://www.schadedoormijnbouw.nl/media/5vcfs4g3/20231113-tno-onderzoek-diepe-bodemdaling-r12185.pdf.

Geurtsen, L. & De Zeeuw, Q., 2017. Monitoring reservoir pressure in the Groningen gasfield. NAM. https://nam-feitenencijfers.data-app.nl/download/rapport/7e818de8-387a-4ab9-a3f4-21c8211e3d5b?open=true.

Guises, R., Embry, J. & Barton, C., 2015. Dynamic geomechanical modelling to assess and minimize the risk for fault slip during reservoir depletion of the Groningen field. Baker RDS. https://nam-onderzoeksrapporten.data-app.nl/reports/download/groningen/en/d32ec1fd-1d59-4f6c-a462-93e2515e9fd9.

Hanssen, R.F., 2001. Radar interferometry: data interpretation and error analysis. Springer.

Hanssen, R.F., Wouters, M.C., Amootaghi, A., Asopa, U., Bruna, M., Janssen, C., Kim, S.S.R., Vossepoel, F.C., Stouthamer, E. & Govers, R., 2020. Monitoring and modeling land subsidence due to hydrocarbon production integrating geodesy and geophysics. AGU Fall Meeting Abstracts. AGU (Washington, DC). https://agu.confex.com/agu/fm20/meetingapp.cgi/Paper/754443.

Hettema, M.H.H., Jaarsma, B., Schroot, B.M. & Van Yperen, G.C.N., 2017. An empirical relationship for the seismic activity rate of the Groningen gas field. Geologie En Mijnbouw/Netherlands Journal of Geosciences 96(5): s149–s161. DOI: 10.1017/njg.2017.18.

Hettema, M.H.H., Papamichos, E. & Schutjens, P., 2002. Subsidence delay: field observations and analysis. Oil and Gas Science and Technology 57(5): 443–458. DOI: 10.2516/ogst:2002029.

Hol, S., Mossop, A.P., van der Linden, A.J., Zuiderwijk, P.M.M. & Makurat, A.H., 2015. Long-term compaction behavior of Permian sandstones – an investigation into the mechanisms of subsidence in the Dutch Wadden Sea. 49th US Rock Mechanics/Geomechanics Symposium. American Rock Mechanics Association (San Fransisco, California): 1–8.

Johnson, A., Charlton, M., Treyer, L. & Ratajc, G., 2022. Pyclipper – Cython wrapper for the C++ translation of the Angus Johnson’s Clipper library (C. Lupo, ed.). ver. 1.3.0. https://pypi.org/project/pyclipper.

Ketelaar, V.B.H., 2009. Satellite Radar interferometry – subsidence monitoring techniques (F. D. Van der Meer, ed.). Springer (Dordrecht): Vol. 14. DOI: 10.1007/978-1-4020-9428-6.

Kim, S.S.R. & Vossepoel, F.C., 2023. On spatially correlated observations in importance sampling methods for subsidence estimation. Computational Geosciences 28: 91–106. DOI: 10.1007/s10596-023-10264-9.

Kumagai, H., Maeda, Y., Ichihara, M., Kame, N. & Kusakabe, T., 2014. Seismic moment and volume change of a spherical source. Earth, Planets and Space 66(1): 1–10. DOI: 10.1186/1880-5981-66-7.

Landman, A.J. & Visser, C., 2023. Groningen Dynamic Model Update 2023 – V7. NAM. https://nam-onderzoeksrapporten.data-app.nl/reports/download/groningen/en/b77a608e-a83c-4857-a75b-1da910bdc4b.

Lele, S.P., Hsu, S., Garzon, J.L., DeDontney, N., Searles, K.H., Gist, G.A., Sanz, P.F., Biediger, E.A.O. & Dale, B.A., 2016. Geomechanical modeling to evaluate production-induced seismicity at Groningen field. Abu Dhabi International Petroleum Exhibition and Conference (Abu Dhabi): 1–18. DOI: 10.2118/183554-ms.

Marketos, G., Govers, R. & Spiers, C.J., 2015. Ground motions induced by a producing hydrocarbon reservoir that is overlain by a viscoelastic rocksalt layer: a numerical model. Geophysical Journal International 203: 228–242. DOI: 10.1093/gji/ggv294.

Marketos, G., Spiers, C.J. & Govers, R., 2016. Impact of rock salt creep law choice on subsidence calculations for hydrocarbon reservoirs overlain by evaporite caprocks. Journal of Geophysical Research: Solid Earth 121: 4249–4267. DOI: 10.1002/2016JB012892.

McTigue, D.F., 1987. Elastic stress and deformation near a finite spherical magma body: resolution of the point source paradox. Journal of Geophysical Research 92(B12): 12931–12940. DOI: 10.1029/jb092ib12p12931.

Mehrabian, A. & Abousleiman, Y.N., 2015. Geertsma’s subsidence solution extended to layered stratigraphy. Journal of Petroleum Science and Engineering 130: 68–76. DOI: 10.1016/j.petrol.2015.03.007.

Mogi, K., 1958. Relations between the eruptions of various volcanoes and the deformations of the ground surfaces around them. Bulletin of the Earthquake Research Institute 36: 99–134.

Mossop, A., 2012. An explanation for anomalous time dependent subsidence. 46th U.S. Rock Mechanics/Geomechanics Symposium. https://onepetro.org/ARMAUSRMS/proceedings-abstract/ARMA12/All-ARMA12/ARMA-2012-518/120766.

NAM, 2013. Technical addendum to the Winningsplan Groningen 2013 subsidence, induced earthquakes and seismic hazard analysis in the Groningen Field. https://www.sodm.nl/documenten/publicaties/2015/06/23/18.-technical-addendum-to-the-winningsplan-groningen-2013-subsidence-induced-earthquakes-and-seismic-hazard-analysis-in-the-groningen-field-nam-november-2013.

NAM, 2015. Meetregister Noord Nederland 2014. https://www.nlog.nl/sites/default/files/meetregisternoordnederland2014.pdf.

NAM, 2020a. Bodemdaling door aardgaswinning in Groningen, Friesland en het noorden van Drenthe. https://nam-feitenencijfers.data-app.nl/download/rapport/aca1b62a-12c5-4918-9ec9-0f0b17aebe56?open=true.

NAM, 2020b. Groningen long term subsidence forecast. https://nam-onderzoeksrapporten.data-app.nl/reports/download/groningen/en/d8970d78-f51a-4a3b-85d4-f80f42d055af.

NAM, 2020c. Petrel geological model of the Groningen gas field, the Netherlands. Open access through EPOS-NL. Yoda data publication platform. Utrecht University. DOI: 10.24416/UU01-1QH0MW.

NAM, 2021. Reservoir pressure and subsidence groningen field update for production profile GTS – raming 2021. https://nam-onderzoeksrapporten.data-app.nl/reports/download/groningen/en/354157f5-5f0b-4fe1-a95a-7da3d2858686.

NAM, 2022. Rapportage Seismiciteit Groningen – November 2022. https://www.rijksoverheid.nl/documenten/rapporten/2022/11/30/nam-rapportage-seismiciteit-groningen-november-2022.

Parliamentary Committee of Inquiry into Natural Gas Extraction in Groningen, 2023. Conclusions and recommendations. In Groningers before Gas. Tweede Kamer der Staten-Generaal. https://www.tweedekamer.nl/Groningen/rapport/boek1_Engelstalig.

Paullo Muñoz, L.F. & Roehl, D., 2017. An analytical solution for displacements due to reservoir compaction under arbitrary pressure changes. Applied Mathematical Modelling 52: 145–159. DOI: 10.1016/j.apm.2017.06.023.

Pijnenburg, R.P.J., Verberne, B.A., Hangx, S.J.T. & Spiers, C.J., 2019. Inelastic deformation of the slochteren sandstone: stress-strain relations and implications for induced seismicity in the Groningen Gas Field. Journal of Geophysical Research: Solid Earth 124: 5254–5282. DOI: 10.1029/2019JB017366.

Pijpers, F. & Van der Laan, D.J., 2016. Trend changes in ground subsidence in Groningen. Statistics Netherlands. https://www.cbs.nl/-/media/_pdf/2016/24/subsidence_may2016.pdf.

Pruiksma, J.P., Breunese, J.N., Van Thienen-Visser, K. & De Waal, J.A., 2015. Isotach formulation of the rate type compaction model for sandstone. International Journal of Rock Mechanics and Mining Sciences 78: 127–132. DOI: 10.1016/j.ijrmms.2015.06.002.

Qin, Y., Salzer, J., Maljaars, H. & Leezenberg, P.B., 2019. High resolution InSAR in the Groningen area. SkyGeo. https://nam-onderzoeksrapporten.data-app.nl/reports/download/groningen/en/26e6a722-b7c2-4fc1-a877-38066df51f14.

Romijn, R., 2017. Groningen Velocity Model 2017 – Groningen full elastic velocity model. NAM. https://nam-feitenencijfers.data-app.nl/download/rapport/9a5751d9-2ff5-4b6a-9c25-e37e76976bc1?open=true.

Schoonbeek, J.B., 1976. Land subsidence as a result of gas extraction in Groningen, the Netherlands. Proceedings of the Second International Symposium on Land Subsidence. International Association of Hydrological Sciences: 267–284. https://www.landsubsidence-unesco.org/wp-content/uploads/2019/03/Proceedings_Anaheim.pdf.

Segall, P., 2010. Earthquake and volcano deformation. Princeton University Press (Princeton, New Jersey) DOI: 10.5860/choice.48-0287.

Shewchuk, J.R., 1996. Triangle: engineering a 2D quality mesh generator and delaunay triangulator. In Applied computational geometry: towards geometric engineering. Series: Lecture Notes in Computer Science. 203–222 pp. https://people.eecs.berkeley.edu/~jrs/papers/triangle.pdf.

Smith, J.D., Avouac, J.P., White, R.S., Copley, A., Gualandi, A. & Bourne, S., 2019. Reconciling the long-term relationship between reservoir pore pressure depletion and compaction in the groningen region. Journal of Geophysical Research: Solid Earth 124(6): 6165–6178. DOI: 10.1029/2018JB016801.

Spiers, C.J., Schutjens, P.M.T.M., Brzesowsky, R.H., Peach, C.J., Liezenberg, J.L. & Zwart, H.J., 1990. Experimental determination of constitutive parameters governing creep of rocksalt by pressure solution. In: R.J. Knipe & E.H. Rutter (eds.): Deformation mechanisms, rheology and tectonics. Geological Society Special Publication (London). Vol. 54, pp. 215–227. DOI: 10.1144/GSL.SP.1990.054.01.21.

Tempone, P., Fjær, E. & Landrø, M., 2010. Improved solution of displacements due to a compacting reservoir over a rigid basement. Applied Mathematical Modelling 34(11): 3352–3362. DOI: 10.1016/j.apm.2010.02.025.

Van der Marel, H., 2020. Comparison of GNSS processing methodologies for subsidence monitoring: NAM GNSS alternative processing method project. Delft University of Technology. http://resolver.tudelft.nl/uuid:d7a427a2-9db4-4468-9a83-5bc18e64bf47.

Van Eijs, R.M.H.E. & Van der Wal, O., 2017. Field-wide reservoir compressibility estimation through inversion of subsidence data above the Groningen gas field. Geologie En Mijnbouw/Netherlands Journal of Geosciences 96(5): s117–s129. DOI: 10.1017/njg.2017.30.

Van Leijen, F.J., Samiei-Esfahany, S., Van der Marel, H. & Hanssen, R.F., 2020. Improving the functional and stochastic model of InSAR. Delft University of Technology. https://nam-onderzoeksrapporten.data-app.nl/reports/download/groningen/en/a040d54f-da6a-4a9d-aef6-1f1fedde4627.

Van Oeveren, H., Valvatne, P., Geurtsen, L. & Van Elk, J., 2017. History match of the Groningen field dynamic reservoir model to subsidence data and conventional subsurface data. Geologie En Mijnbouw/Netherlands Journal of Geosciences 96(5): s47–s56. DOI: 10.1017/njg.2017.26.

Van Opstal, G.H.C., 1974. The effect of base-rock rigidity on subsidence due to reservoir compaction. 3rd Congress of the International Society of Rock Mechanics (Denver, Colorado): 1102–1111.

Van Thienen-Visser, K. & Fokker, P.A., 2017. The future of subsidence modelling: compaction and subsidence due to gas depletion of the Groningen gas field in the Netherlands. Geologie En Mijnbouw/Netherlands Journal of Geosciences 96(5): s105–s116. DOI: 10.1017/njg.2017.10.

Van Thienen-Visser, K., Nepveu, M., Van Kempen, B., Kortekaas, M., Hettelaar, J., Peters, L., Van Gessel, S. & Breunese, J., 2014. Recent developments of the Groningen field in 2014 and, specifically, the southwest periphery of the field. https://www.nlog.nl/sites/default/files/finaltnoreportekl.pdf.

Van Wees, J.D., Fokker, P.A., Van Thienen-Visser, K., Wassing, B.B.T., Osinga, S., Orlic, B., Ghouri, S.A., Buijze, L. & Pluymaekers, M., 2017. Geomechanical models for induced seismicity in the Netherlands: inferences from simplified analytical, finite element and rupture model approaches. Geologie En Mijnbouw/Netherlands Journal of Geosciences 96(5): s183–s202. DOI: 10.1017/njg.2017.38.

Vasco, D.W., Harness, P., Pride, S. & Hoversten, M., 2017. Estimating fluid-induced stress change from observed deformation. Geophysical Journal International 208: 1623–1642. DOI: 10.1093/gji/ggw472.

Vasco, D.W., Johnson, L.R. & Goldstein, N.E., 1988. Using surface displacement and strain observations to determine deformation at depth, with an application to Long Valley Caldera, California. Journal of Geophysical Research 93(B4): 3232–3242. DOI: 10.1029/JB093iB04p03232.

Vatti, B.R., 1992. A generic solution to polygon clipping. Communications of the ACM 35(7): 56–63. DOI: 10.1145/129902.129906.

Visser, C.A. & Solano Viota, J.L., 2017. Introduction to the Groningen static reservoir model. Geologie en Mijnbouw/Netherlands Journal of Geosciences 96(5): s39–s46. DOI: 10.1017/njg.2017.25.

Visvalingam, M. & Whyatt, J.D., 1993. Line generalisation by repeated elimination of points. The Cartographic Journal 30(1): 46–51. DOI: 10.1179/000870493786962263.

Wang, R., Lorenzo-Martín, F. & Roth, F., 2006. PSGRN/PSCMP – a new code for calculating co- and post-seismic deformation, geoid and gravity changes based on the viscoelastic-gravitational dislocation theory. Computers and Geosciences 32: 527–541. DOI: 10.1016/j.cageo.2005.08.006.

Wessel, P., Luis, J.F., Uieda, L., Scharroo, R., Wobbe, F., Smith, W.H.F. & Tian, D., 2019. The Generic Mapping Tools Version 6. Geochemistry, Geophysics, Geosystems 20: 5556–5564. DOI: 10.1029/2019GC008515.

Wildenborg, T., Peters, L., Moghadam, A., Fokker, P., Geel, K., Nelskamp, S., Bottero, S., Wiersma, A. & Marsman, A., 2022. KEM-19 – evaluation of post-abandonment fluid migration and ground motion risks in subsurface exploitation operations in the Netherlands. TNO and Deltares. https://kemprogramma.nl/file/download/5d5758ed-189e-4cc3-b7dc-8185627c99ac/kem-19-d1-literature-review_deliverable-d1_v2022.pdf.