The influence of shallow elastic heterogeneities on the subsidence bowl resulting from deep reservoir depletion

Subsidence due to gas extraction poses significant risks to infrastructure, as even moderate vertical displacements can cause structural damage. Existing models typically predict smooth, broad subsidence bowls, but observational data reveals short-wavelength fluctuations superimposed on this pattern, which are crucial for assessing the potential risk of damage. Current models overlook these variations, leading to an incomplete understanding of subsidence-induced hazards. This research aims to assess the role of shallow elastic heterogeneities in explaining these short-wavelength subsidence fluctuations. We combined InSAR observations with finite element subsidence modelling in the Groningen gas field to analyze the effects of shallow geological structure and elastic moduli variations. Our results show that small-scale elastic heterogeneities cause local displacement changes ranging from 1% to 9%, consistent with InSAR observations, while horizontal strains can vary significantly by 30%–200%. These findings highlight the importance of accounting for short-wavelength fluctuations in subsidence models, as they have direct implications for risk assessment and predicting building damage due to gas extraction.

1 Introduction

It is well known that gas reservoir compaction at depth causes land subsidence (e.g., Van Eijs and van der Wal, 2017; Zoccarato et al., 2018; Gazzola et al., 2023; Candela and Koster, 2022; Fokker et al., 2025). Gas depletion leads to reservoir compaction at depth, and the resulting land subsidence reflects the subsurface elastic response (Fokker et al., 2018; Candela et al., 2022). Most geomechanical models, if not all, indicate a rather smooth subsidence bowl. However, InSAR observations reveal short wavelength variations superimposed on rather smooth subsidence bowls. This, for example, is evident from observations at the Groningen gas field in the Netherlands (Ketelaar et al., 2020; Wouters et al., 2025). Similar subsidence patterns have recently also been observed in other smaller gas extraction locations in the Netherlands (Verberne et al., 2025a), and internationally, for example, in the Ravenna coastal plain (Italy) (Verberne et al., 2025b), and the Hatfield Moors (United Kingdom) (Fibbi et al., 2024).

The short wavelength variations observed by InSAR are of particular importance as they are responsible for localised deformations, such as differential settlements and strains, which are the primary cause of structural damage (Prosperi et al., 2023). It is therefore essential to understand the causes of these local variations in the subsidence bowl and their relative magnitude: they cannot arise from a smooth subsidence bowl. This understanding is the objective of the present contribution.

Multiple processes could lead to short wavelength variations in subsidence bowls. One potential cause for these short wavelength variations is local subsidence caused by phreatic groundwater level changes, indirectly induced by the regional scale subsidence bowl (Kooi et al., 2021). Another candidate is the spatially heterogeneous elastic response of the shallow unconsolidated subsurface to reservoir compaction. Indeed, any elastic heterogeneities in the shallow subsurface can distribute the strain response to deep reservoir compaction in a heterogeneous manner. As a result, the presence of laterally confined small-scale shallow elastic heterogeneities could cause short wavelength variations in the initially smooth subsidence bowl caused by gas extraction. This study explores this hypothesis by integrating numerical models with observational data.

We start with a new post-processing of the InSAR observations to accurately characterise the magnitude and the length scales of the observed subsidence fluctuations. Then, we target a connection to geomechanical modelling. Geomechanical models generally do not account for laterally confined, small-scale shallow elastic heterogeneities. It is challenging to accurately represent these small-scale and finite-size shallow elastic heterogeneities analytically. Numerical codes can effectively model lateral changes in elastic properties, however, most of the time, solely large-scale geological features are accounted for, which again results in a rather smooth subsidence bowl. Therefore, we deploy a numerical solution that incorporates small-scale elastic heterogeneities. A comparison between the observed and modelled fluctuations leads to the proposition that indeed the elastic heterogeneities can be the cause of the observed variability.

2 Groningen gas field: exploitation and typical subsidence modelling

The Groningen gas field is situated in the northern coastal plain of the Netherlands. It was discovered in 1959 and was in production from 1963 until 2024. The field covers an area of approximately 900 km² and is located at depths ranging from 2600 to 3200 m. The exploitation has been causing several issues, the most notable being induced seismicity (Van Thienen-Visser and Breunese, 2015), and land subsidence (Fokker and Van Thienen-Visser, 2016). The subsidence has been observed and modelled, reaching a maximum of 37 cm in the centre of the gas field and approximately 10–20 cm at the edges (Figure 1; Fokker and van Thienen-Visser, 2016). Figure 1 illustrates a typical modelling approach assuming a homogeneous subsurface, which does not account for shallow elastic heterogeneities. As a result, the subsidence map and the subsidence profiles appear to be rather smooth. Figure 1b illustrates a detail of the geology of surficial unconsolidated sediments indicating shallow heterogeneities (TNO-GDN, 2021). The Groningen gas field is delineated in black.

Figure 1

Figure 1. (a) Subsidence calculated with estimated compaction field in 2013 (adapted from Fokker and van Thienen-Visser, 2016), and (b) location of the Groningen gas field. Profiles (a,b) are detailed in Figure 4 for the context of the analysed problem.

3 Methods

3.1 InSAR observations

InSAR observations were used to characterise small-scale subsidence fluctuations superimposed on the large-scale subsidence bowl of the Groningen field. Processed data from raw Sentinel-1 images between 2016 to mid-2020 were used, which were made publicly accessible in the Netherlands by the Department of Waterways and Public Works (Rijkswaterstaat, 2022a; Rijkswaterstaat, 2022b; Verberne et al., 2023). This dataset contains time series data of displacements in the Line Of Sight (LOS) direction for four satellite tracks (two ascending and two descending) in the Groningen area. It includes point-wise time series for both persistent and distributed scatterers, referred to as Level-2 data, which are calibrated to the GNSS network. The uncertainty associated with the estimated InSAR time series is influenced by various factors, including the quality of the scatterers, atmospheric and orbital perturbations, noise in the original radar data, and choices made during the kinematic and geometric parameterization (Brouwer and Hanssen, 2023; Hanssen, 2001). These factors affect the InSAR results, introducing both dispersion and bias, with particular attention given to the estimation of integer phase ambiguities. The technique relies on natural reflections from objects, making it more suitable for urban areas than for rural ones. Specifically, it is designed to optimally monitor coherent scattering objects, such as buildings and infrastructure. Since many of these structures in the Netherlands are founded on deep (Pleistocene) layers, the recorded displacements do not capture shallow soil deformations. However, such deformations can be detected by scatterers with shallow foundations.

The Sentinel-1 satellite radar images have a resolution of approximately 5 m by 20 m. However, this resolution does not directly correspond to the resolution of the InSAR measurement points. The radar images are analyzed on a pixel-by-pixel basis, where the measurement value for each pixel represents the sum of reflections from all objects within a “resolution cell.” A resolution cell is larger than a pixel and contains the raw data from the radar image, which is then combined into a final pixel. Only pixels with stable and relatively strong reflections over time are selected as measurement points. The sampling interval for the data points varies depending on the availability of satellite passes: a 6-day repeat pass up until April 2016, and a 12-day repeat pass from April 2016 onwards (Wegmüller et al., 2015). To ensure consistency, all tracks were resampled onto the same spatial and temporal grid. Temporal resampling was performed using linear interpolation, with the first time epoch of the resampled time series serving as the temporal reference, which was then subtracted from the series. For spatial resampling, two grid sizes, both larger than the Sentinel-1 spatial resolution, were tested: dx = 50 m, dy = 50 m, and dx = 100 m, dy = 100 m. A grid cell averaging algorithm was used to compute the mean value for each non-empty grid cell.

InSAR records the movement of Earth’s surface along the LOS direction. Based on the SAR observation geometry, one can establish the dependence of the LOS displacement on the 3D displacement. To recover the 3D displacement at each grid cell, the four LOS observations with different SAR looking directions are combined. However, due to the near-polar configuration of SAR satellite orbits, the LOS deformation is less sensitive to north-south displacement. By neglecting this component, ascending and descending tracks can be combined to extract vertical and east-west horizontal displacement. For each grid cell, Equation 1 can be written as follows:

$\begin{array}{c}{LOS}^1=e_1{U}_e+u_1{U}_u\\ {LOS}^2=e_2{U}_e+u_2{U}_u\\ \begin{array}{c}{LOS}^3=e_3{U}_e+u_3{U}_u\\ {LOS}^4=e_4{U}_e+u_4{U}_u\end{array}\end{array}(1)$

where subscripts 1, 2, 3, 4 correspond to the four different tracks. $e$ and $u$ represent the east and vertical components of SAR looking vector. $U_e$ and $U_u$ are the east and vertical displacement components. For our analysis, only the vertical displacement time series $,U_u$, has been used. Figure 2 displays the cumulative subsidence from 2016 to mid-2020 on a regular spatial grid of dx = 100 m, dy = 100 m.

Figure 2

Figure 2. InSAR-based observed cumulative subsidence in mm for the Groningen area (that is vertical displacement, $U_u$) from 2016 to mid-2020 and on regular spatial grid with dx = 100 m, dy = 100 m. The large-scale subsidence bowl caused by the Groningen gas field depletion is clearly visible on the North-West part of the map. The elongated blue rectangle corresponds to the selected band for small-scale analysis (see Figure 3).

A rectangular band of the two subsidence gridded maps (dx = 50 m, dy = 50 m and dx = 100 m, dy = 100 m) on top of the Groningen gas field was selected for further analysis of small-scale spatial subsidence fluctuations (see Figures 2, 3). Subsidence profiles were generated by averaging subsidence values along the y-axis (the latitudinal axis) for each longitudinal cell. A Fourier-based low-pass filter was then applied to each profile with a wavelength cut-off of 1000 m. This means that the high-frequency noise, which typically manifests as rapid, short-wavelength fluctuations (potentially caused by atmospheric effects, residual errors from the processing, or other disturbances), will be attenuated, leaving us with a smoother signal that captures the lower-frequency features assumed to be real signal potentially caused by the gas field. Note that the 1000 m wavelength cutoff is deliberately large to exclude potential noise, though the exact value is somewhat arbitrary. The filtered profiles align well for both gridded maps (dx = 50 m, dy = 50 m and dx = 100 m, dy = 100 m), demonstrating that such InSAR post-processing is not affecting the wavelengths longer than 1000 m in the signal. In both filtered profiles, persistent small-scale spatial fluctuations are visible, superimposed over the overall subsidence bowl. The length-scale of these small-scale spatial fluctuations is approximately 3–4 km and their amplitude is around 3–5 mm. After projection onto the vertical, the InSAR resolvability threshold can be as low as 0.5 mm/year (Hanssen, 2001; Wouters et al., 2025); these 3–5 mm are thus considered as real signal. Given that the total subsidence between 2016 and mid-2020 reaches 50 mm, the ratio between the amplitude of small-scale fluctuations and the overall subsidence bowl is approximately 6%–10% (3–5 mm/50 mm). The following section demonstrates that these small-scale fluctuations can be induced by shallow elastic heterogeneities, which locally magnify the depletion-induced subsidence bowl.

Figure 3

Figure 3. Selected band for InSAR small-scale analysis (see Figure 2 for location with respect to the Groningen area). (a) Maps of cumulative subsidence (in mm) from 2016 to mid-2020 and on regular spatial grid of dx = 100 m, dy = 100 m (top) and dx = 50 m, dy = 50 m (bottom). (b) Averaged subsidence profiles for the two gridded maps displayed in (a). (c) Averaged and Fourier-based low-pass filtered subsidence profiles for the two gridded maps displayed in (a).

3.2 Geological framework for subsurface heterogeneities above the Groningen gas field

The objective of this section is to provide a brief and qualitative overview of the geological elastic heterogeneities in the Groningen gas field area. The idea here is that each geological elastic heterogeneity within the subsurface can lead to a deviation from the overall depletion-induced smooth subsidence bowl. Indeed, the assumption of an elastic homogeneous subsurface leading to a smooth subsidence bowl might break down as soon as subsurface elastic heterogeneities are identified.

Directly above the Groningen gas field, deep geological units (below ∼1000 m depth) characterising a typical profile are shown in Figure 4c. This profile is derived from DGMdeep, an open access 3D geological subsurface model of the upper 8500 m of the Dutch subsurface, available both onshore and offshore (TNO-GDN, 2019). The main known deep spatial confined elastic heterogeneities are the presence of salt domes in the Zechstein group. The intermediate depth consists of unconsolidated sediments (roughly between ∼25 and 1000 m depth) which record in the upper 200 m the succession of several glacial cycles (van Dijke and Veldkamp, 1996). Distinctive geological features within this depth domain are glacial tunnel valley infills that intersect the subsurface, locally with base depths ranging from 30 m to 800 m. The width of the mapped valleys is approximatively 1.5–2 km, although near surface valley infills are known to be broader. The depth of their most contrasting lithological infill, i.e., clay, within the area of the subsidence bowl of the Groningen gas field typically lies between 30–190 m (Figure 4b; TNO-GDN, 2025a). The lithological infill of the tunnel valleys consists from bottom to top of coarse sand, medium coarse sand, clay, fine sand. The clay within the tunnel valleys has peculiar geomechanical properties owing to their extreme high content of clay particles (deep lacustrine deposit), and subsequent loading and shearing by land ice during later glaciations (Schokking, 1998).

Figure 4

Figure 4. Geological heterogeneities on top of the Groningen gas field. (a) Soil-type heterogeneities at shallow depth in north-south direction (derived from TNO-GDN’s GeoTOP model; TNO-GDN, 2025b), with the blue-line indicating the base of the Holocene coastal deposits, and the red-line the extend of the upper part of the glacial valleys (cf. Kruiver et al. 2017a). This profile is observed in Figure 1b represented with a blue line (profile a). (b) Depth of the base of the clayey infill of the glacial tunnel valleys (TNO-GDN, 2025a), with the contour of the Groningen gas field provided in black for spatial reference, and (c) salt domes in the Zechstein group (TNO-GDN, 2019), which profile is observed in Figure 1b represented with a red line (profile b). The arrow on the left side of the figure indicates the transition from deep-scale strata to shallow-scale strata.

At the shallowest depths (ranging between 0 and ∼25 m), the GeoTOP subsurface model indicates small-scale (hectometer-to kilometer-scale) soil-type heterogeneities (variations of peat, clay, and sand) (Figure 4a). GeoTOP is a 3D geological subsurface model of the upper 50 m of the Netherlands. It is open access and provides information on lithoclasses (soil types) and lithostratigraphy (geological units subdivided based on lithoclasses and vertical position) (Stafleu et al., 2021; TNO-GDN, 2025b). The deepest occurrences of heterogeneities primarily consist of clay and peat infills formed during the last interglacial, which are confined in depressions carved into underlying stiff sand (Peeters et al., 2015). In the area of the Groningen gas field, the clayey and peaty valley infills reaches thicknesses up to 10 m with a width of 5 km. The most surficial soft soil deposits in the area are alternating tidal clay and peat beds formed during the Holocene epoch (Figure 1b). Here, they formed on a regional scale with a thickness ranging between >1 and ∼15 m in tidal basin systems grading landwards into drowned brook valleys (Vos and Knol, 2015). The soft soil beds are locally alternating with more sandy tidal channel deposits.

Our focus is on explaining the small-scale (kilometer-scale) fluctuations observed with InSAR (see Figure 3). The deeper the elastic heterogeneities, the smoother their potential impact on the subsidence bowl. Indeed, a preliminary geomechanical analysis confirms that deep salt domes can only lead to very diffuse and large-scale (multi-kilometer-scale) spatial fluctuations of the subsidence bowl. Kilometer-scale glacial valleys might provide a suitable explanation for InSAR observations, but because of their depth, preliminary geomechanical analysis confirms that their influence on the subsidence bowl is smoothed down and distributed over larger scales. Only the very shallow elastic heterogeneities identified in the GeoTOP model (Figure 4a), directly below the ground surface, remain as a viable candidate to produce hectometer-to kilometer-scale subsidence fluctuations of similar magnitude to those observed in InSAR. The subsequent analysis will focus on these shallow elastic heterogeneities.

3.3 Subsidence modelling considering elastic heterogeneities

Most analytical and semi-analytical methods used to model subsidence due to reservoir depletion rely on simplifying assumptions. These approaches typically produce smooth, bowl-shaped subsidence patterns, as localised elastic heterogeneities cannot be incorporated into the boundary conditions of a semi-analytical solution (Fokker and Orlic, 2006; Paullo Muñoz and Roehl, 2017; Weijermars, 2023). Additionally, while it is common to find studies on land subsidence using more advanced numerical frameworks, the importance of small-scale shallow subsurface elastic heterogeneities is often neglected (Gambolati et al., 2001).

To account for the effect of small-scale shallow subsurface elastic heterogeneities, the following numerical approach was implemented. A finite-element strategy was selected to analyse the current problem, conducted using Diana Fea (2017), which has been validated in the simulation of reservoir depletion analyses (e.g., Buijze et al., 2019), and heterogeneous small-scale problems (e.g., Buijze et al., 2020). The aim was not to replicate the entire complexity of the Groningen subsurface (Figure 4) but to focus on the effect of a single shallow elastic heterogeneity. To achieve this, a generic synthetic model, based on the Groningen gas field setting, was developed.

The applied numerical model consists of two elastic layers (Figure 5). The deeper unit is located between 1000 m and 2300 m; while the shallower unit is located between 0 m and 1000 m. A localised surface elastic heterogeneity is situated at a horizontal distance (from the left boundary to the centre of the elastic heterogeneity) of 3,000 m. Its depth is d = 50 m, and its width w is varied to examine its effects on the ground surface. The values considered for w are 25 m, 50 m, 100 m, 200 m, 300 m, 400 m, 500 m, 700 m, 1000 m, 1500 m, 2000 m, and 2500 m.

Figure 5

Figure 5. The applied numerical model considering a shallow elastic heterogeneity. Reservoir subject to depletion in green. Surface elastic heterogeneity H in black, its dimensions w is varied in the simulations.

The elastic properties of the deeper unit and the reservoir are based on literature and were fixed in our analysis, with a Young’s modulus of E = 12.85 GPa and a Poisson’s ratio of ν = 0.28. The elastic properties of the shallower layer unit were derived from (i) shear wave velocities and densities of 29 different lithostratigraphic soil units, and (ii) from literature, following the work of Orlic (2016), and Kruiver et al. (2017a), Kruiver et al. (2017b). Note that the complete delineation of the soil units can be found in Kruiver et al. (2017a), while the input for the derivation of the elastic properties can be found in the respective GeoTOP model dataset (Ntinalexis et al., 2022). Figure 6a displays the Poisson’s ratio distribution of two wells situated above the Groningen gas field, while Figure 6b illustrates the best estimate elastic modulus for various soil types. Based on these observations, we selected an average elastic modulus E = 0.2 GPa and an average Poisson’s ratio of ν = 0.45 for the shallow layer unit. Lastly, the elastic modulus of the elastic heterogeneity is based on the estimated upper and lower limit values in Figure 6b, specifically, E_soft = 0.025 GPa and E_stiff = 1.0 GPa.

Figure 6

Figure 6. (a) Poisson’s ratio for two deep wells in the Groningen field, and (b) best estimated elastic modulus of different soil types.

The model is fixed in the normal direction at the bottom and both vertical boundaries to allow for realistic deformation of the domain. To induce depletion, a pressure drop of 5 MPa was applied to the reservoir. This pressure drop magnitude mimics the average pressure drop at Groningen from 2016 to mid-2020 in order to align with the InSAR observation time window.

Due to the considerable size of the spatial model extent and the necessity to use relatively small elements to accurately capture the elastic response of the subsurface around the shallow elastic heterogeneity, a two-phase analysis approach was implemented. The first phase involved a field-scale analysis. In this phase, the entire domain was simulated using a coarse mesh (with elements of approximately 100 m × 100 m) to capture large-scale deformations (Figure 7a), resulting in a set of deformation contours (Figure 7b). The second phase focused on a reduced domain, which includes only the shallow subsurface and the elastic heterogeneity. This reduced domain was re-meshed with a fine mesh to capture short-wavelength variations. Note that the reduced domain in Figure 7c corresponds to the red dotted line in Figure 7b. The deformations from the first phase were then applied as prescribed displacements at the boundaries of the reduced domain, replicating the same deformation patterns. The final small-scale mesh used to analyse ground displacements and deformations consists of square elements, each measuring 2.5 m per side. Finally, Figure 7d shows a comparison between the surface deformation using the large- and small-scale models with generic material properties and a small-scale mesh of elements approximately 10 times smaller than those from the large-scale model. As observed, the two-phase analyses approach can interpolate the large-scale results into the small-scale model to reproduce ground deformations accurately, with differences being virtually zero. This comparison also serves as a mesh refinement validation. The surface displacement curves obtained with the coarse and refined meshes are indistinguishable, confirming that the solution is numerically stable and mesh-converged. Note that in the large-scale model, the boundaries are fixed using rollers. This prevents normal-to-the-edge displacements and only tangential movements are allowed, thereby avoiding the development of localised unfeasible stresses. In contrast, all boundaries in the small-scale model are prescribed, meaning that they are rigid boundaries with positions taken directly from the large-scale model. It is important to ensure that these prescribed conditions coincide with locations where displacements have been evaluated, such as edges with explicit nodes, so that interpolation remains accurate. Although this rigid prescription could, in principle, restrict boundary relaxation, the imposed displacements are taken directly from the large-scale solution and are therefore expected to introduce virtually no differences.

Figure 7

Figure 7. Two-phase analysis. (a) Large-scale model with a coarse mesh, (b) Illustrative results from the large-scale model, (c) segment of the large-scale model with a fine mesh and prescribed displacements from the large-scale model, and (d) validation between the large- and the small-scale models results. Sketches not to scale.

Finally, it is important to emphasise that the present research serves as an introduction to the effects of elastic heterogeneities on the subsidence bowl, focusing only on elastic heterogeneities. Other important sources of variability, such as poroelastic effects, groundwater changes and plasticity, are certainly worth exploring but remain beyond the scope of this study, for clarity and conciseness.

4 Results

Figure 8 presents an example of a subsidence simulation due to reservoir depletion, considering an elastic heterogeneity of width w = 500 m. Figure 8a illustrates the large-scale simulation considering a stiff elastic heterogeneity (E_stiff = 1.0 GPa). Although the elastic heterogeneity is not the smallest that can be analysed, its size is insufficient to capture its effects at the ground surface. Furthermore, to illustrate the element scaling issue, one element of the mesh with dimensions 100 m × 100 m has been highlighted in red (next to the elastic heterogeneity). Clearly, this element, despite appearing being small, is too large in relation to the elastic heterogeneity, therefore insufficient to yield accurate deformation results. However, after the remeshing step in the area highlighted in red around the elastic heterogeneity in Figure 8a (second phase), the domain is refined in terms of element size. In this case, an arrangement of 3 × 3 elements of size 6 m × 6 m is highlighted within the elastic heterogeneity. As observed, the elements are now small enough, and suitable to capture the impact of the elastic heterogeneity at the ground surface. Note that the horizontal and vertical scales of the figures are not the same, causing a change in the dimensions of the elastic heterogeneity.

Figure 8

Figure 8. (a) Field-scale simulation using a coarse mesh and one element of the mesh of dimensions 100 m × 100 m highlighted in red, and (b) small-scale simulation using a fine mesh and an arrangement of 3 × 3 elements of size 6 m × 6 m highlighted in red. The black edge shows the position of the elastic heterogeneity.

Figure 9 shows the ground surface vertical displacement considering both soft and a stiff elastic heterogeneities. Figures 9a,b show the displacements when a soft elastic heterogeneity is present. When w is small, a sharp upward vertical displacement occurs. As w increases, the displacement at the centre of the elastic heterogeneity diminishes and redistributes towards the edges. In contrast, when the elastic heterogeneity is stiffer than the surrounding material, a downward vertical deformation occurs (Figures 9c,d). As w decreases, the vertical deformation also decreases. However, this reduction is not substantial, and the vertical deformation persists despite changes in w.

Figure 9

Figure 9. (a) Vertical displacements considering a soft elastic heterogeneity of varying width w, (b) close view of results, (c) vertical displacements considering a stiff elastic heterogeneity of varying width w, and (d) close view of results. Red curve is for the case without an elastic heterogeneity. The red curve illustrates the case without heterogeneity.

In contrast, Figure 10 shows the effects of the elastic heterogeneity in the horizontal displacement. Figures 10a,b illustrate the effects of a soft elastic heterogeneity, while Figures 10c,d illustrates the effects of a stiff elastic heterogeneity.

Figure 10

Figure 10. (a) Horizontal displacements considering a soft elastic heterogeneity of varying width w, (b) close view of results, (c) horizontal displacements considering a stiff elastic heterogeneity of varying width w, and (d) close view of results.

Note that the opposite behaviour of soft and stiff heterogeneities (i.e., uplift for the soft case and downward deflection for the stiff case) can be explained by stress redistribution within the shallow subsurface. When the heterogeneity is softer than its surroundings, it is compressed laterally by the stiffer material, which results in a relative upward deflection at the surface above the heterogeneity. In contrast, when the heterogeneity is stiffer than its surroundings, it resists compression and instead deforms laterally, leading to a relative downward deflection at the surface.

Regarding the relative change of the vertical displacement due to the elastic heterogeneities, this is defined as the ratio between the magnitude of the local deflection of the vertical displacement caused by the elastic heterogeneity and the magnitude of the overall subsidence bowl at the location of the elastic heterogeneity without the presence of the elastic heterogeneity (Figure 11a). Note that this corresponds to the ratio between the amplitude of the small-scale fluctuations and the overall subsidence bowl derived from InSAR observations.

Figure 11

Figure 11. (a) relative change of vertical movement as a function of w for both soft and stiff elastic heterogeneity, and (b) relative change of the horizontal displacement as a function of w, considering both soft and stiff elastic heterogeneity.

The dependence of this relative change on the width of the elastic heterogeneity, for both a soft (E_soft = 0.025 GPa) and a stiff elastic heterogeneity (E_stiff = 1.0 GPa), is displayed in Figure 11a. Overall, the relative change is between 2% and 9%; it is thus of the same order of magnitude as observed in InSAR (see Figure 3). Also, the relative change of the horizontal displacement as a function of the width of the elastic heterogeneity and for both a soft and stiff elastic heterogeneity is displayed (Figure 11b). In both cases, stiff and soft, the relative change increases as w increases. Overall, the change is between 1% and 9%. Additionally, the results indicate that surface displacements and deformations consistently occur in the presence of an elastic heterogeneity, but the amplitude of these changes does not vary significantly with heterogeneity width (w). Instead, the primary role of w is to control how and where the deformation is expressed. For small heterogeneities, the response is concentrated directly above the inclusion, whereas for large heterogeneities, the deformation shifts toward the edges. An additional observation is that vertical and horizontal relative changes behave in a contrasting manner, when the vertical response is strong, the horizontal response is weak, and vice versa.

To conclude the analysis of the surface elastic heterogeneity, we extend the 2D model to a simplified 3D setting in order to examine the effect of elastic heterogeneity thickness. As illustrated in Figure 12, the elastic heterogeneity in the 3D model is located at the centre of the domain. In this case, the length w of the elastic heterogeneity is kept constant at 700 m, while the thickness η is varied. The values considered for η are 50, 100, and 300 m. The reservoir characteristics in the extra dimension are assumed to be constant across its extent.

Figure 12

Figure 12. Three-dimensional (3D) sketch of the model (not to scale).

Figure 13 presents an example of the 3D analysis. Figure 13a shows the mesh of the model, with the small elastic heterogeneity highlighted in red. Figure 13b displays the vertical deformation. Due to the large scale of the model, the effects of the elastic heterogeneity are barely visible. However, Figure 14 illustrates the effects of the elastic heterogeneity at a smaller scale. Similar to the 2D framework, due to the presence of the small elastic heterogeneities, the subsidence profile is deflected. The difference between 3D and 2D cases is not substantial, and as expected when the thickness η increases, the 3D case gradually converges to the 2D case. It is observed that, in the case of the soft elastic heterogeneity, a thickness η matches the 2D solution quickly (that is, after η reaches 300 m). However, a stiff elastic heterogeneity need a larger thickness to match the 2D scenario (that is, after η reaches 1000 m).

Figure 13

Figure 13. Example of the 3D simulation. (a) full mesh with the elastic heterogeneity colored in red and the reservoir in blue, and (b) vertical displacements. Note that the contrast between the different grey zones is due to the meshing and not to use of different materials.

Figure 14

Figure 14. Vertical displacement above a 3D elastic heterogeneity along the w-direction. Effect of the thickness η considering a (a) soft elastic heterogeneity, and a (b) stiff elastic heterogeneity. Red curve is the reference 2D case for w = 700 m.

Although it is well known that 3D simulations converge to the 2D case when the geometry remains constant in the out-of-plane direction, it is useful to highlight this behaviour in the present context. The convergence occurs because the 2D model inherently assumes that the geometry in the third direction is infinite. As the thickness of the elastic heterogeneity in the 3D model increases, the results progressively approximate those of the 2D case, until the heterogeneity is sufficiently thick that the differences become negligible. This behaviour supports the use of 2D models for the analysis of real-case problems where the geometry can be reasonably considered extensive in one direction, thereby reducing computational cost and time.

5 Discussion and conclusion

Our use-case regards the Groningen gas field, but the generic set-up of the modelling allows its application to any gas field in the Netherlands or elsewhere. First, we demonstrated that, using InSAR observations, small-scale fluctuations of the subsidence bowl with amplitude higher than the expected noise level can be clearly observed. The magnitude of these small-scale (kilometer-scale) fluctuations is between 6% and 10% of the background signal. Second, we assessed the expected depth and size of the subsurface elastic heterogeneities above the Groningen gas field. Considering that the effect of deeper elastic heterogeneities on the subsidence bowl is smoothed out by the subsurface elastic response, we focused our modelling efforts on the shallower elastic heterogeneities. Indeed, shallow lithological elastic heterogeneities are clearly present in our subsurface model (see Figure 4). Results of our 2D numerical models indicated that the relative changes in vertical and horizontal movements are between 1% and 9%. Results of the 3D numerical models were very close to those obtained with our 2D numerical models; only slightly larger for horizontally constrained elastic heterogeneities. The modelled relative changes are thus in the range of those observed with InSAR. Physically, a 1%–9% relative change corresponds to local steepening or flattening of the ground-surface displacement profile over short (kilometre-scale) distances. Because horizontal strain is approximately the spatial gradient of horizontal displacement, such localised perturbations amplify differential settlement and strain demand in structures, thereby elevating damage risk. Of course, this correlation does not necessarily imply that the observed fluctuations are caused by the shallow elastic heterogeneities. However, our numerical models have been deployed with all known realistic ingredients: the pressure depletion, reservoir thickness, expected ranges for elastic moduli. It should also be noted that our modelling does not include other relevant processes, such as plastic deformation, hydromechanical coupling, time-dependent effects, and stochastic heterogeneities, are also relevant and warrant comprehensive investigation. These processes were not included here because their treatment requires a substantially broader modelling framework, and their inclusion would dilute the specific focus on elastic heterogeneities and complicate the delineation of their effects. A full investigation of such processes would require an extensive parametric study, which lies beyond the scope of the present work but represents an important direction for future research. We conclude that our approach has demonstrated that shallow elastic heterogeneities, effectively present in the shallow subsurface, are a candidate to explain the small-scale variations of the subsidence bowl observed in InSAR.

A pressing question concerns the effect of these fluctuations in the displacement profiles on the potential induced damages to buildings and/or infrastructures. It is well known that the spatially localised variations of the subsidence profile can potentially lead to damages to housings and/or infrastructures (Prosperi, 2025). For a deep source of subsidence, as it is the case for a gas depleting reservoir, the induced horizontal strains at the ground surface, proportional the local spatial gradient of the horizontal displacement, are expected to be the main drivers of damage to housings and/or infrastructures (Geurts et al., 2021). Figure 15 presents the induced horizontal strains derived from the profiles of horizontal displacements of Figures 10a,b.

Figure 15

Figure 15. Horizontal strains considering a soft elastic heterogeneity.

The induced horizontal strains at the ground surface increase with the decrease of the width of the elastic heterogeneity. The strain peaks at the edge of the elastic heterogeneity are caused by the sharp transition in Young’s modulus between the elastic heterogeneity and its surroundings. In reality, one expects plasticity to cap these sharp peaks in elastic strains; therefore, we better focus our analysis and interpretation on the centre of the elastic heterogeneity, which is unaffected by this edge effect. The relative change of the horizontal strain at the centre of the elastic heterogeneity is between 30% and 200%. It means that moderate changes in horizontal displacements (1%–9%) distributed over a short distance can lead to substantial changes of the horizontal strain (30%–200%). Determining how changes in horizontal strain translate to the probability of damage to housing and infrastructure is beyond the scope of this study. However, we note here that the building response needs to be considered, presumably with the use of a fragility curve that translates induced horizontal strains to a probability of exceedance of a certain damage state.

The three main takeaways from our study are the following:

1. The presented numerical model, inspired by the Groningen gas field and populated with realistic model parameters, indicates that shallow and small-scale elastic heterogeneities can result in local relative changes of approximately 1% and 9% in the vertical and horizontal displacements at the ground surface.

2. The modelled fluctuations in the vertical displacements fall within the same range as those observed with InSAR.

3. The moderate changes of approximately 1% and 9% in the horizontal displacements can lead to substantial changes of around 30% and 200% in the horizontal strain.

Data availability statement

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

Author contributions

JG: Validation, Writing – review & editing, Investigation, Formal Analysis, Methodology, Writing – original draft, Visualization. TC: Investigation, Conceptualization, Methodology, Supervision, Project administration, Writing – review and editing, Writing – original draft. MP: Writing – review and editing, Supervision, Project administration. GZ: Validation, Writing – review and editing, Visualization, Formal Analysis. KK: Writing – review and editing, Supervision. PF: Investigation, Writing – review and editing, Formal Analysis, Methodology. FA: Investigation, Writing – review and editing.

Funding

The author(s) declared that no financial support was received for this work and/or its publication.

Acknowledgements

This study was part of the KEM-program of the Dutch Ministry of Climate Policy and Green Growth. We thank the reviewers for their useful comments. Also, we are indebted to Marcel Bakker for kindly providing us maps and metrics of the subglacial tunnel valley deposits.

Conflict of interest

The author(s) declared that this work was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

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

References

Buijze, L., van den Bogert, P. A., Wassing, B. B. T., and Orlic, B. (2019). Nucleation and arrest of dynamic rupture induced by reservoir depletion. J. Geophys. Res. Solid Earth 124 (4), 3620–3645. doi:10.1029/2018jb016941

CrossRef Full Text | Google Scholar

Brouwer, W.S., and Hanssen, R.F. (2023). A treatise on InSAR geometry and 3-D displacement estimation. IEEE Trans. On Geosc. and Rem. Sens.61, 1–11. doi:10.1109/TGRS.2023.3322595

CrossRef Full Text | Google Scholar

Buijze, L., Guo, Y., Niemeijer, A. R., Ma, S., and Spiers, C. J. (2020). Nucleation of stick-slip instability within a large-scale experimental fault: effects of stress heterogeneities due to loading and gouge layer compaction. J. Geophys. Res. Solid Earth 125 (8), e2019JB018429. doi:10.1029/2019JB018429

PubMed Abstract | CrossRef Full Text | Google Scholar

Candela, T., and Koster, K. (2022). The many faces of anthropogenic subsidence. Science 376 (6600), 1381–1382. doi:10.1126/science.abn3676

PubMed Abstract | CrossRef Full Text | Google Scholar

Candela, T., Chitu, A. G., Peters, E., Pluymaekers, M., Hegen, D., Koster, K., et al. (2022). Subsidence induced by gas extraction: a data assimilation framework to constrain the driving rock compaction process at depth. Front. Earth Sci. 10, 713273. doi:10.3389/feart.2022.713273

CrossRef Full Text | Google Scholar

Diana Fea, B. V. (2017). User's manual. Netherlands: Delft. Available online at: https://www.dianafea.com.

Google Scholar

Fibbi, G., Novellino, A., Bateson, L., Fanti, R., and Dels Soldato, M. (2024). Multidisciplinary assessment of seasonal ground displacements at the hatfield moors gas storage site in a peat bog landscape. Sci. Rep. 14, 22521. doi:10.1038/s41598-024-73548-9

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Fokker, P. A., and Van Thienen-Visser, K. (2016). Inversion of double-difference measurements from optical leveling for the Groningen gas field. Int. J. Appl. Earth Observation Geoinformation 49, 1–9. doi:10.1016/j.jag.2016.01.004

CrossRef Full Text | Google Scholar

Fokker, P. A., Van Leijen, F. J., Orlic, B., Van Der Marel, H., and Hanssen, R. F. (2018). Subsidence in the Dutch Wadden Sea. Neth. J. Geosciences 97, 129–181. doi:10.1017/njg.2018.9

CrossRef Full Text | Google Scholar

Fokker, P. A., Candela, T. G. G., Erkens, G., Hanssen, R. F., Kooi, H., Koster, K., et al. (2025). “Subsidence,” in Geology of the Netherlands. Editors J. H. Ten Veen, G.-J. Vis, J. De Jager, and Th.E. Wong second edition (Amsterdam: Amsterdam University Press), 825–847.

Google Scholar

Gambolati, G., Ferronato, M., Teatini, P., Deidda, R., and Lecca, G. (2001). Finite element analysis of land subsidence above depleted reservoirs with pore pressure gradient and total stress formulations. Int. J. Numer. Anal. Methods Geomechanics 25 (4), 307–327. doi:10.1002/nag.131

CrossRef Full Text | Google Scholar

Gazzola, L., Ferronato, M., Teatini, P., Zoccarato, C., Corradi, A., Dacome, M. C., et al. (2023). Reducing uncertainty on land subsidence modeling prediction by a sequential data-integration approach. Application to the arlua off-shore reservoir in Italy. Geomechanics Energy Environ. 33, 100434. doi:10.1016/j.gete.2023.100434

CrossRef Full Text | Google Scholar

Geurts, C. P. W., Pluymaekers, M. P. D., and Rots, J. G. (2021). Report No.: TNO 2021 R10325B. Schade aan gebouwen door diepe bodemdaling en -stijging. Delft: TNO report.

Google Scholar

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

Google Scholar

Ketelaar, G., Bähr, H., Liu, S., Piening, H., van der Veen, W., Hanssen, R., et al. (2020). Integrated monitoring of subsidence due to hydrocarbon production: consolidating the foundation. PIAHS 382, 117–123.

CrossRef Full Text | Google Scholar

Kooi, H., Landwehr, J. C., Stuurman, R., van Meerten, J. J., Levelt, O., and Korff, M. (2021). Indirecte schade-effecten van diepe bodemdaling en -stijging bij het Groningen gasveld en gasopslag Norg. Deltares rapport 11207096-002-BGS-0001, 67.

Google Scholar

Kruiver, P. P., Wiersma, A., Kloosterman, F. H., de Lange, G., Korff, M., Stafleu, J., et al. (2017a). Characterisation of the Groningen subsurface for seismic hazard and risk modelling. Neth. J. Geosciences 96, s215–s233. doi:10.1017/njg.2017.11

CrossRef Full Text | Google Scholar

Kruiver, P. P., van Dedem, E., Romijn, R., de Lange, G., Korff, M., Stafleu, J., et al. (2017b). An integrated shear-wave velocity model for the Groningen gas field, the Netherlands. Bull. Earthq. Eng. 15, 3555–3580. doi:10.1007/s10518-017-0105-y

CrossRef Full Text | Google Scholar

Ntinalexis, M., Kruiver, P. P., Bommer, J. J., Ruigrok, E., Rodriguez-Marek, A., Edwards, B., et al. (2022). A database of ground motion recordings, site profiles, and amplification factors from the Groningen gas field in the Netherlands. Earthq. Spectra 39, 687–701. doi:10.1177/87552930221140926

CrossRef Full Text | Google Scholar

Orlic, B. (2016). Geomechanical effects of CO2 storage in depleted gas reservoirs in the Netherlands: inferences from feasibility studies and comparison with aquifer storage. J. Rock Mech. Geotechnical Eng. 8 (6), 846–859. doi:10.1016/j.jrmge.2016.07.003

CrossRef Full Text | Google Scholar

Paullo Muñoz, L. F., and Roehl, D. (2017). An analytical solution for displacements due to reservoir compaction under arbitrary pressure changes. Appl. Math. Model. 52, 145–159. doi:10.1016/j.apm.2017.06.023

CrossRef Full Text | Google Scholar

Peeters, J., Busschers, F. S., and Stouthamer, E. (2015). Fluvial evolution of the rhine during the last interglacial-glacial cycle in the southern north Sea basin: a review and look forward. Quat. Int. 357, 176–188. doi:10.1016/j.quaint.2014.03.024

CrossRef Full Text | Google Scholar

Prosperi, A. (2025). Modelling of subsidence induced damage to masonry buildings. Influence of soil heterogeneity and settlement and development of fragility curves. (doctoral thesis). Delft University of Technology, Delft, Netherlands.

Google Scholar

Prosperi, A., Longo, M., Korswagen, P. A., Korff, M., and Rots, J. G. (2023). Sensitivity modelling with objective damage assessment of unreinforced masonry facades undergoing different subsidence settlement patterns. Eng. Struct. 286, 116113. doi:10.1016/j.engstruct.2023.116113

CrossRef Full Text | Google Scholar

Rijkswaterstaat (2022a). Dataregister rijkswaterstaat. Available online at: https://maps.rijkswaterstaat.nl/dataregister/srv/eng/catalog.search#/metadata/d9c04bc0-d8b2-4791-a072-8b46c25166ac?tab=relations.

Google Scholar

Rijkswaterstaat (2022b). FAQ InSAR based deformation service for the Dutch built environment - rijkswaterstaat publicatie platform. Available online at: https://open.rijkswaterstaat.nl/overige-publicaties/2020/faq-insar-based-deformation-service-for.

Google Scholar

Schokking, F. (1998). Anisotropic strength behaviour of a fissured overconsolidated clay in relation to Saalian glacial directions. Eng. Geol. 49, 31–51. doi:10.1016/s0013-7952(97)00035-5

CrossRef Full Text | Google Scholar

Stafleu, J., Maljers, D., Busschers, F. S., Schokker, J., Gunnink, J.L., and Dambrink, R.M. (2021). “Models created as 3-D cellular voxel arrays,” in Applied Multidimensional Geological Modeling: Informing Sustainable Human Interactions with the Shallow Subsurface. Editors A. K. Turner, H. Kessler, and M. J. van der Meulen (John Wiley and Sons, Ltd), 247–271. doi:10.1002/9781119163091.ch11

CrossRef Full Text | Google Scholar

TNO – GDN (2019). DGM-diep v5.0. TNO - Geol. Dienst Nederl. Available online at: https://www.dinoloket.nl/en/subsurface-models/map.

Google Scholar

TNO – GDN (2021). Geological map 2021. Available online at: https://www.dinoloket.nl/en/subsurface-models/map.

Google Scholar

TNO – GDN (2025a). BRO REGIS II v2.2.3 TNO - geologische dienst nederland. Available online at: https://www.dinoloket.nl/en/subsurface-models/map.

Google Scholar

TNO – GDN (2025b). BRO GeoTOP v1.6.1 TNO - Geologische dienst nederland. Available online at: https://www.dinoloket.nl/en/subsurface-models/map.

Google Scholar

Van Dijke, J. J., and Veldkamp, A. (1996). Climate-controlled glacial erosion in the unconsolidated sediments of northwestern Europe, based on a genetic model for tunnel valley formation. Earth Surf. Process. Landforms 21 (4), 327–340. doi:10.1002/(sici)1096-9837(199604)21:4<327::aid-esp540>3.0.co;2-p

CrossRef Full Text | Google Scholar

Van Eijs, R., and van der Wal, O. (2017). Field-wide reservoir compressibility estimation through inversion of subsidence data above the Groningen gas field. Neth. J. Geosciences 96, 117–129. doi:10.1017/njg.2017.30

CrossRef Full Text | Google Scholar

Van Thienen-Visser, K., and Breunese, J. N. (2015). Induced seismicity of the Groningen gas field: history and recent developments. Lead. Edge 34, 664–671. doi:10.1190/tle34060664.1

CrossRef Full Text | Google Scholar

Verberne, M., Koster, K., Lourens, A., Gunnink, J., Candela, T., and Fokker, P. A. (2023). Disentangling shallow subsidence sources by data assimilation in a reclaimed urbanized coastal plain, south Flevoland polder, the Netherlands. J. Geophys. Res. Earth Surf. 128, e2022JF007031. doi:10.1029/2022jf007031

CrossRef Full Text | Google Scholar

Verberne, M. A. M., Koster, K., De Bresser, H., and Fokker, P. A. (2025a). Disentangling subsidence from shallow soil processes and gas extraction in a Dutch UNESCO world heritage polder with InSAR and data assimilation. Earth Syst. Environ. doi:10.1007/s41748-025-00885-8

CrossRef Full Text | Google Scholar

Verberne, M. A. M., Teatini, P., Koster, K., Fokker, P. A., and Zoccarato, C. (2025b). An integral approach using InSAR and data assimilation to disentangle and quantify multi-depth driven subsidence causes in the Ravenna coastland, northern Italy. Geomechanics Energy Environ. 43, 100710. doi:10.1016/j.gete.2025.100710

CrossRef Full Text | Google Scholar

Vos, P. C., and Knol, E. (2015). Holocene landscape reconstruction of the Wadden Sea area between marsdiep and weser. Neth. J. Geosciences 94, 157–183. doi:10.1017/njg.2015.4

CrossRef Full Text | Google Scholar

Wegmüller, U., Santoro, M., Werner, C., and Cartus, O. (2015). “On the estimation and interpretation of Sentinel-1 tops InSAR coherence,” in Proc. ‘Fringe 2015 workshop’.

Google Scholar

Weijermars, R. (2023). Surface subsidence and uplift resulting from well interventions modeled with coupled analytical solutions: application to Groningen gas extraction (netherlands) and CO2-EOR in the kelly-snyder oil field (west Texas). Geoenergy Sci. Eng. 228, 211959. doi:10.1016/j.geoen.2023.211959

CrossRef Full Text | Google Scholar

Wouters, M. C., Govers, R., and Hanssen, R. (2025). Development of an efficient model to calculate subsidence above the Groningen gas field. Neth. J. Geosciences 104, e12151. doi:10.70712/njg.v104.12151

CrossRef Full Text | Google Scholar

Zoccarato, C., Ferronato, M., and Teatini, P. (2018). Formation compaction vs land subsidence to constrain rock compressibility of hydrocarbon reservoirs. Geomechanics Energy Environ. 13, 14–24. doi:10.1016/j.gete.2017.12.002

CrossRef Full Text | Google Scholar