Comparison of hydrocarbon and geothermal energy production in the Netherlands: reservoir characteristics, pressure and temperature changes, and implications for fault reactivation

Abstract

The Netherlands is in the midst of an energy transition with hydrocarbon production gradually declining, whereas the role of sustainable energy technologies is on the rise. One of these technologies is geothermal energy production from porous reservoirs at 1.5–3 km depth. As the number of geothermal projects increases, there is a growing concern that felt and/or damaging induced seismic events could occur as a result of geothermal operations. Over the last two decades, such events have occurred in the Netherlands due to gas production, notably in the Groningen gas field. However, the occurrence of felt events is limited to hydrocarbon fields in certain regions or reservoirs. Understanding where and for which plays these events are observed helps to estimate seismogenic potential for geothermal operations and other sustainable subsurface activities. Here, we summarise and review the main similarities and differences in terms of geological and geomechanical characteristics between the hydrocarbon and geothermal plays in the Netherlands, and we consider the differences in pressure and temperature changes. By doing so, we provide better insights into the factors that could play a role for fault reactivation and induced seismicity, and how these differ for hydrocarbon production and geothermal operations in the Netherlands. The review shows that geological characteristics for most geothermal target reservoirs are similar to those of hydrocarbon, albeit geothermal projects so far target higher porosity rocks than hydrocarbon reservoirs. On the other hand, pressure and temperature changes are very different, with significant depletion for hydrocarbon fields vs significant cooling around geothermal injection wells. The different operations result not only in different expected stress change magnitudes but also in a distinct spatio-temporal stress build-up on faults, which has implications for seismogenic potential and monitoring of these different operations.

Introduction

For decades, natural gas has been one of the main energy sources in the Netherlands. Since the discovery of the Groningen gas field in 1959, gas was produced from over 350 different gas fields, with an annual production of 19.1 BCM from 217 producing fields in 2021 (Ministry of Economic Affairs & Climate, 2021). Presently, many smaller gas fields are nearing the end of their productive lifespan and gas production is declining. Also, the gas production from the giant Groningen gas field is currently being phased out due to the occurrence of numerous induced seismic events that caused damage to infrastructure and housing (Muntendam-Bos et al., 2022). At the same time, there is a strong demand for more sustainable forms of energy. One of the sustainable energy technologies in the energy transition is geothermal energy production. In the Netherlands, medium temperature (50–100°C) fluids can be produced by circulating fluid through porous and/or permeable sedimentary reservoirs at depth (1.5–3.5 km). With the help of a heat pump, even shallower, lower temperature reservoirs could also be exploited. Around 20 geothermal projects are currently operational within sedimentary formations with primary porosity and permeability, having produced 6.2 PJ in the year 2020 (Mijnlieff, 2020). This is expected to increase to ∼15 PJ/year in 2030 (PBL, 2020). Deeper, tighter and/or fractured sedimentary formations requiring stimulation are not yet being exploited, but may be in the future (ter Heege et al., 2020).

One of the concerns related to geothermal energy production is the occurrence of seismic events large enough to be felt at the surface (Buijze, et al., 2019a; Evans et al., 2012; Foulger et al., 2018). To date however, no event with a magnitude large enough to be felt (∼M > 1.5–2) has been recorded near any of the operational geothermal sites in the Netherlands that target sandstone reservoirs (Buijze et al., 2019a; Muntendam-Bos et al., 2022). However, microseismicity with local magnitudes up to M_L 1.7 has been recorded near two doublets targeting fractured carbonates in the south-east, resulting in the suspension of one of the doublets (Vörös & Baisch, 2022). Some concerns about the occurrence of induced seismicity near the other doublets remain; the operation time of most doublets is still relatively short, and the number of geothermal projects is limited. Moreover, concerns remain because of the country’s history with induced seismicity related to gas production, notably the Groningen field but also ∼ 20 other fields with maximum magnitudes between 2.0 and 3.5 (Dost et al., 2012; Muntendam-Bos et al., 2022; Van Eck et al., 2006; Van Thienen-Visser & Breunese, 2015). Geothermal projects operate in the same sandstone reservoirs that host several of the main hydrocarbon plays (Mijnlieff, 2020). The targeted reservoir characteristics might be different due to different economic boundary conditions required for geothermal operations. Also, geothermal operations result in different pressure and temperature changes compared to hydrocarbon production, which leads to different stress changes. In this article, we consider the similarities and differences between hydrocarbon and geothermal energy production and list the implications for induced seismicity. In the first part, we review the main reservoir characteristics, as well as pressure and temperature changes. In the second part, we summarise the observed induced seismic events and discuss differences in the stressing mechanisms expected for gas production and geothermal energy production. We limit the comparison to the current producing or abandoned onshore hydrocarbon fields and current geothermal targets – that is, geothermal operations in permeable sedimentary rocks not requiring stimulation.

Review of geological and operational characteristics of hydrocarbon and geothermal plays

Here, we summarise the main geological characteristics of the hydrocarbon plays and the geothermal plays in the Netherlands. Because geothermal production is restricted to the onshore, emphasis is placed on the onshore and near-shore fields and characteristics.

Reservoir rocks and caprocks

In the Netherlands, both gas and oil plays are present, with gas plays being the most prolific. The main source rocks of the gas plays in the Netherlands are predominantly coals of the Upper Carboniferous Limburg Group, and the main source rocks for the oil plays are the younger Jurassic shales, notably the Posidonia shale (Fig 1). The extent of the hydrocarbon plays is controlled by the maturity of the source rock but also by the presence of good reservoir and caprocks. Most geothermal plays target the same formations that also host the hydrocarbon fields. The geothermal play types in the Netherlands are so-called conductive intra-cratonic basin play types (Moeck, 2014). The Netherlands is a tectonically quiet area with an average thermal gradient of 31°C/km at depth (Békési et al., 2018; Bonté et al., 2012) and a lower gradient of 20°C/km in the shallowest 400 m (Gies et al., 2021). Current geothermal operations target sedimentary aquifers which have temperatures suitable for horticulture and district heating (∼50–100°C). Geothermal heat is produced via balanced circulation through these formations. Productivity of the geothermal systems is controlled by sufficient transmissivity or permeability-thickness (∼permeability × thickness) and a sufficiently high reservoir temperature (Van Wees et al., 2012). Many operational geothermal projects are situated in regions which also host hydrocarbon fields, as ample seismic data are available (e.g. Mijnlieff, 2020), although projects must remain outside of the direct influence area of hydrocarbon fields. However, due to the different economic requirements for geothermal production, also other areas without hydrocarbon fields are also of interest for geothermal (e.g. van Wees et al., 2020)

Fig. 1. Overview of main lithostratigraphic groups, tectonic events, lithologies and gas, oil and geothermal plays in the Netherlands. (based on: de Jager & Geluk, 2007; Duin et al., 2006; Geluk, 2005; Van Adrichem Boogaert & Kouwe, 1997, www.nlog.nl, www.dinoloket.nl/nomenclator).

Fig 2 shows the location of hydrocarbon fields and geothermal doublets, as well as the main structural elements present during Mesozoic rifting. The structural elements are divided in three types; highs, where Permian and/or Carboniferous rocks have not been deposited or have been eroded, platforms, which are characterised by the absence of the Jurassic formations due to erosion, and (inverted) basins, where the Jurassic has been preserved (Kombrink et al., 2012). In the following, we briefly summarise the main hydrocarbon and geothermal plays in the Netherlands; for an extensive overview see for example, Wong et al. (2007), Mijnlieff (2020), van Wees et al., (2020) and Willems et al. (2020), as well as www.nlog.nl for production licences and geological information of hydrocarbon fields.

Fig. 2. Location of hydrocarbon fields and geothermal projects in the Netherlands. Background shading indicates structural elements present at Upper Jurassic – Lower Cretaceous time (after: Békési et al., 2018; Duin et al., 2006; Kombrink et al., 2012). Dotted pink line: Zechstein salt 100 m isopach indicating the southern boundary of the Zechstein salt deposits. Dotted brown line; southern extent of the Ten Boer Claystone Member (DGM-diep v5 https://www.nlog.nl/dgm-diep-v5-en-offshore). Colored shapes: hydrocarbon fields, circles: geothermal projects. Hydrocarbon fields and geothermal projects are colored according to the stratigraphic group of the main reservoir formation. Red outlines indicate oil production, black outlines gas production (field info from www.nlog.nl). N: North Sea Supergroup, NU: Upper North Sea Group, NL: Lower North Sea Group, CK: Chalk, KN: Rijnland Group, SG: Scruff Group, SL: Schieland Group, AT: Altena group, RN: Upper Germanic Trias Group, RB: Lower Germanic Trias Group, ZE: Zechstein Group, RO: Upper Rotliegend, DC: Limburg Group, CL: Carboniferous Limestone.

Kolenkalk Group, Lower Carboniferous age

Lower Carboniferous carbonates do not host hydrocarbon fields but are a geothermal target. Two doublets targeted the Zeeland Formation of the Kolenkalk Group which, in the few locations where they lie at shallower depth and have been drilled, are comprised of shales and (locally) of karstified platform carbonates (Kombrink, 2008, Mozafari et al., 2019). These platforms are located at accessible depths along the northern and southern fringe of the Ruhr Valley Graben and around the London Brabant Massif (Bouroullec et al., 2019; ter Heege et al., 2020). Instead of circulation through a porous matrix, circulation occurs through karsts, fractures and faults. Unfortunately, operations of the Dutch doublets have been terminated, in one case due to the occurrence of seismic events (Muntendam-Bos et al., 2022; Vörös & Baisch, 2022). Just south of the Dutch border, the same carbonates are targeted at Balmatt, near Mol in Belgium (Broothaers et al., 2021). Here also, seismic events up to M_L 2.2 occurred, which could be felt in the villages nearby (Kinscher et al., 2023). The events led to suspension of operations.

Limburg Group, Upper Carboniferous age

A minor gas play is found in the Limburg Group, notably in the fluvial sandstones of the Tubbergen Formation. These reservoirs occur mostly where the Upper Rotliegend is absent, like in the Lower Saxony Basin in the east of the Netherlands, where various formations from the Limburg Group are directly overlain by the Zechstein salt (NAM, 2013), for example, the Coevorden Field (Kombrink et al., 2007), as well as on the Cleaver Bank High in the Dutch offshore. The sands of the Limburg Group are of fluvial origin and are intercalated with shale, silt and coal layers, with net-to-gross 30–70%. No geothermal projects currently target these formations.

Upper Rotliegend Group, Permian age

The Upper Rotliegend is the largest gas play in the Netherlands, hosting the giant Groningen gas field in the north of the country, as well as numerous other fields (de Jager & Visser, 2017). The producing formation is the Slochteren Formation, which consists of both aeolian and fluvial sandstones with high net-to-gross ratios. It is present across most of the onshore area of the Netherlands, apart from the Lower Saxony Basin and the various highs (Figs 2 and 3), where Cretaceous uplift resulted in the erosion of the Permian of where no Rotliegend was deposited in the first place (Duin et al., 2006; Geluk, 1999). The southern extent of the Upper Rotliegend gas play is controlled by the presence of the Zechstein Group. This thick sequence of salt overlying the Upper Rotliegend forms an excellent caprock to the Rotliegend in the north of the Netherlands. However, towards the south, the salt pinches out and the Zechstein no longer consists of salt but of limestone and shale, and therefore does not constitute a proper seal. In the southernmost parts of the country, it is absent due to erosion in the Kimmeridgian, limiting the southern extent of this play (Fig 2). Towards the north, the clay content of the Slochteren Formation progressively increases, transitioning into the Silverpit Formation, with the Ten Boer Claystone Member overlying the Slochteren sandstone (De Jager & Visser, 2017), and the Ameland Claystone Member separating the Upper and Lower Slochteren sandstones. The claystone members may contain gas, but their contribution is minor compared to the production from the Slochteren Formation. Further south, these members are absent and the Slochteren Formation directly underlies the Basal Zechstein. North-northwest of the Dutch onshore also many gas fields are found, but north of onshore Netherlands the Slochteren Formation is largely absent, thereby limiting the widespread occurrence of Upper Rotliegend gas fields, apart from small area in the northermost part of the Dutch offshore. The Slochteren Formation is unconformably underlain by silts, shales, sandstones, interbedded with coal seams, of the Limburg Group.

Fig. 3. Depths of main litho-stratigraphic groups targeted by hydrocarbon or geothermal. Location of hydrocarbon fields (blue polygons) and geothermal doublets (red stars) are indicated. Depths are retrieved from DGM 5.0 (www.nlog.nl). KN: Rijnland, S: Schieland + Scruff, RN: Upper Germanic Trias, RB: Lower Germanic Trias, ZE: Zechstein, RO: Upper Rotliegend (see also Fig 1).

The Slochteren Formation is also a prime geothermal target. In 2021, nine doublets circulated through the Slochteren (Ministry of Economic Affairs & Climate, 2021) with several more planned to start soon. Whereas the extent of the hydrocarbon play is limited by the presence of the Zechstein salt caprock, the geothermal play extends further south. Favorable conditions are found in the centre of the Netherlands, in the Central Netherlands Basin, North Holland Platform and Friesland Platform around the Texel-IJsselmeer High (Fig. 2) (Vrijlandt et al., 2019). At these locations, the Ten Boer and Ameland Claystone Members are mostly absent, and the Slochteren Formation is present as a single interval, directly overlain by the carbonates and anhydrites of the basal Zechstein, with only a thin sequence of Zechstein salt present. Around the Texel-IJsselmeer High, the Slochteren Formation is primarily of aeolian origin (Mijnlieff, 2020). Occasionally anhydrite beds are present. The Slochteren Formation is also present further south, but the combination of deeper burial depth and southward thinning has led to poor transmissivity in these areas.

Zechstein Group, Permian age

A different hydrocarbon play type is found in the Zechstein, where the main reservoir rocks are within marginal platform carbonates embedded in salt. The main producing members are the Z2 and Z3 carbonates, which host several fields located in the east, centre and west Netherlands and the central offshore, along the southern boundary of the Zechstein salt. Locally, the gas is enriched in condensates derived from source rocks present within the Zechstein. A number of fields form stacked fields with producing horizons both in the Zechstein and the underlying Slochteren Formation and/or Upper Carboniferous rocks.

Main Buntsandstein Group, Triassic age

The Triassic play is the second-largest gas play and is mainly present in the central and southern Dutch onshore, the Lower Saxony Basin and the offshore; in the northern onshore, it has largely been eroded during the Jurassic (Figs 2 and 3). In the northern half of the Netherlands, only few gas fields are found in the Triassic as most gas from the Carboniferous source rocks is trapped below the Zechstein salt; only locally leakage through the Zechstein caprock via, for example, faults allowed for the charge of Triassic gas fields. Several gas shows in the Triassic above the Zechstein salt are found in the Lower Saxony Basin (e.g. the Roswinkel field, Roest et al., 1998) or in the west of the Central Netherlands Basin, where several gas fields produce(d) from both the Rotliegend and the Triassic (e.g. Bergen, Bergermeer). In the West Netherlands and Broad Fourteens Basins, the absence of the Zechstein salt allowed upwards migration of gas from the Westphalian source rock, trapping it within sandstones of the Lower Germanic Triassic. Numerous gas fields are found along the fault-bounded margins of both basins. Occasionally oil is also found, likely derived from the overlying Jurassic shales. The main producing horizons are found in the Main Buntsandstein Group; the Volpriehausen, Detfurth and Hardegsen Formations of fluvial and aeolian origin. These are separated by claystones but towards the south they form one massive stack of sandstone. In addition, several productive units are found in the Upper Germanic Trias Solling and Röt Formations. The formations are capped by Upper Germanic Triassic shales, evaporites and limestones (Muschelkalk, Keuper, Röt and/or Solling Formations) or Lower Cretaceous shales (Vlieland Claystone) and are underlain by, for example, the Rogenstein Formation of the Lower Buntsandstein Subgroup, an oolitic siltstone.

For geothermal operations on the other hand, the Triassic is only a minor target. The sandstones of the Main Buntsandstein Subgroup are the main target for geothermal circulation. Only one doublet operates in the Triassic (Vierpolders), and a second is being realised, both located at the southern margins of the West Netherlands Basin. The porosity and permeability of the Triassic is variable and can be low in several locations because of its relatively deep palaeo-burial depth and cementation, as also demonstrated by the Naaldwijk geothermal well NLW-GT-01 (Boersma et al., 2021). Along the southern edge of the West Netherlands Basin favourable transmissivity-temperature conditions are proven to exist (Mijnlieff, 2020). In the southeastern part of the country, the poorly known Lower Buntsandstein Nederweert Sandstone Formation is considered a potential target for geothermal exploration, although it has not been drilled for this purpose yet.

Schieland & Rijnland Group, Upper Jurassic and Lower Cretaceous age

Most Dutch oil fields are found in Upper Jurassic and Lower Cretaceous formations, which overlie the Lower Jurassic source rock, the Posidonia Shale. These fields are located mostly in the Mesozoic basins (Duin et al., 2006) and the Dutch Central Graben (Fig 2). Many oil fields also have a gas cap, with gas derived mostly from the Carboniferous source rocks. Most of the oil fields in the West Netherlands Basin have ceased production. The main producing reservoir units were the Rijswijk, Berkel, IJsselmonde and De Lier Members of the Lower Cretaceous Vlieland Sandstone Formation, host to for example, the Berkel, De Lier and Rotterdam fields. The reservoir rocks are mostly shallow marine sandstones with significant variability in thickness due to syn-tectonic sedimentation, and variability in facies and hence transmissivity (Willems et al., 2020). The main caprock of these reservoirs is the Vlieland Claystone Formation. Some reservoirs are also found in the Upper Cretaceous Holland Formation, and a few oil fields are found in Upper Jurassic Schieland Group (e.g. Wassenaar), with the fluvial Delft Sandstone Member forming the main reservoir. Another Lower Cretaceous oil play is found in the Lower Saxony Basin, with the Schoonebeek field being the biggest onshore oil field in Europe. The main reservoir here are the Bentheim sandstone of the Vlieland Formation, overlain by claystones and marls. The Jurassic/Cretaceous formations are also host to several gas fields, most of which are found in the Friesland platform. The main reservoir formation here is the Friesland Member of the Vlieland Sandstone Formation, capped by the Vlieland Claystone Formation and Holland Marl Members.

The Jurassic/Cretaceous play is currently the largest geothermal play, with most operational doublets targeting the Delft Sandstone Member of the Schieland Group (Willems et al., 2020, Ministry of Economic Affairs & Climate, 2021). The Delft Sandstone is of fluvial origin and has a high net-to-gross ratio. Favorable conditions are predominantly found in the West Netherlands Basin. Because the Schieland Group is composed of syn-rift deposits, geothermal operations target the higher thickness formations in the graben structures, as opposed to hydrocarbon fields which are mainly found on the highs having a relatively thin reservoir section (e.g. Donselaar et al., 2015). The underlying Alblasserdam Member, which has a lower net-to-gross ratio, is also targeted by various doublet systems. A few doublets also (co-)produce from the marine sands of the Vlieland Formation, which is considered less productive due to its aforementioned variability in facies (Fig 1). In the north, the Vlieland Formation is also present, but its conditions are currently considered to be unfavourable for geothermal operations (Mijnlieff, 2020).

Chalk Group

Contrary to Denmark and Norway where multiple hydrocarbon fields are found in the Chalk, only a few hydrocarbon fields are found in the Chalk in the Netherlands. These include an oil field in the northern offshore (Hanze F2A Field) (Hofmann et al., 2002) and the Harlingen Field (van den Bosch, 1983) and De Wijk Field onshore.

North Sea Supergroup, Cenozoic age

At two locations, gas is found in the Basal Dongen Formation of the North Sea Group, in the De Wijk and Wanneperveen gas fields (Bruijn, 1996). Both produce from multiple intervals, including also the Triassic. The Cenozoic deposits may constitute a large geothermal play, as they are present throughout most of the Netherlands. Temperatures are relatively low (10–50°C), but with the use of a heat pump, temperatures can be raised to temperatures that can be used for greenhouse heating. In fact, a first geothermal exploration well was drilled in 1987 to explore geothermal potential of the Cenozoic, but the reservoir quality encountered was too low. Currently, a single doublet is operational, targeting the Brussel Sandstone Member, but in the future more doublets may be developed in this play (Geel et al., 2022; Veldkamp et al., 2022). Apart from the Brussel Sandstone Member, various other Cenozoic Formations have been identified as potential geothermal target, but they have not been drilled yet.

Depth, porosity, permeability and thickness

The different requirements for economically viable hydrocarbon vs geothermal operations are not only reflected in different project locations (Fig 2) but also in the targeted depth ranges, the porosity and the permeability which relate to the project locations. As discussed earlier, hydrocarbon reservoirs require the presence of a mature source rock, a structural or stratigraphic trap, and a caprock. In addition, the reservoir permeability needs to be such that commercial flow rates can be attained at the wells. Gas reservoirs with permeabilities lower than 0.1–1 mD and porosities lower than 10% are considered tight and will require extra stimulation to achieve economic flow rates (e.g. Herber & de Jager, 2010). A number of gas fields in the (offshore) Netherlands, predominantly in the Slochteren Formation, are not developed because of their tight nature or because so-called cut-off volumes are too small; these are considered ’stranded assets' (Herber & De Jager, 2010). Most of the producing and abandoned fields have permeabilities between 1 and 1000 mD and average porosities between 5 and 25%, with depths ranging from 1000 to 4000 m depth. Porosities and permeabilities generally decrease with depth and hence typically also with formation age, though significant scatter is present in the data, partly because of inversion following deep burial. Fair permeability and porosity may still be found at depths down to 4000 m, depending on the amount of inversion, the occurrence of diagenetic processes during the geological history in combination with the timing of the gas charge (e.g. Van Kempen et al., 2018).

The produced geothermal power is a function of the temperature difference between the formation temperature and the injected cold water and the heat capacity, as well as the flow rate, which depends strongly on the reservoir transmissivity, as well as fluid viscosity (van Wees et al., 2012). Transmissivity values exceeding ∼ 5 Dm are considered a requirement for commercially attractive operations. The temperature increases with depth, but permeability (and transmissivity) decreases, resulting in an optimal depth window for geothermal energy production. This window is found at depths between 1500 and 3000 m, a narrower range compared to hydrocarbon field depths. The producing geothermal doublets in the Netherlands target relatively permeable (>50 mD) reservoirs with thicknesses of 50 m to > 200 m, translating into transmissivities from 5 to > 100 Dm. Average porosities often exceed 15% and are on average higher than those of the hydrocarbon fields. Note that many of the doublets currently target ‘sweet spots’ in term of transmissivity. In the future, more doublets may target lower transmissivity ranges, depending on economic conditions such as the gas price but also the costs of drilling and heat pumps and the ability to stimulate the reservoir. Operations in deeper, presumably tighter (sandstone) formations do require higher injection pressure or may require stimulation, as in for example, Enhanced Geothermal Systems, which are currently not considered in the Netherlands. Also further exploitation of fractured carbonates and tight siliciclastic rocks may be considered (ter Heege et al., 2020). Operations in shallower formations may become economically attractive with the help of a heat pump.

Elastic properties of the reservoir rocks and caprock lithologies

Elastic properties are important for the magnitude of stress changes that can occur due to pressure and temperature changes. With increasing depth and decreasing porosity, the elastic properties change due to mechanical and chemical compaction. In Fig 5, elastic properties obtained from well log measurements (Hunfeld et al., 2021) are shown for different litho-stratigraphic groups, with static Young’s modulus E _stat recalculated from the dynamic Young’s modulus (Eissa & Kazi, 1988). Logs were taken within hydrocarbon reservoirs (squares) or along exploration wells targeting hydrocarbons. No log data for any of the producing geothermal wells was included in the well log dataset, but several of the wells are located close to and in similar formation as the geothermal doublets (Hunfeld et al., 2021).

Fig. 4. Porosity, permeability and reservoir thicknesses for hydrocarbon reservoirs (squares) and geothermal doublets (circles). Hydrocarbon reservoirs include producing and abandoned fields, but not undeveloped fields. Porosity and permeability are formation-averaged values, from petrophysical and core analyses of wells within hydrocarbon fields (www.nlog.nl), or well tests within geothermal wells (Supplementary Materials 1). Reservoir thickness is the net gas-bearing reservoir section for hydrocarbon reservoirs (Van Thienen-Visser et al., 2012), and the net productive reservoir thickness for geothermal. In c dashed lines indicate temperatures for an average geothermal gradient of 31 °C/km, and in d dashed lines depict permeability x thickness.

Fig. 5. Elastic properties of main litho-stratigraphic groups obtained from formation-averaged well log data (Hunfeld et al., 2021) against depth. Symbols are color-coded according to their litho-stratigraphic group (Figure 1 and www.dinoloket.nl/nomenclator). Squares indicate values derived from logs in hydrocarbon reservoirs, diamonds logs in overburden formations or from exploration wells. a) Static Young’s modulus computed from the dynamic value (Eissa & Kazi, 1988) and b) dynamic Poisson’s ratio for sandstone intervals, and c) static Young’s modulus and d) dynamic Poisson’s ratio for clay-rich intervals (claystone, siltstone, shale), as well anhydrites and carbonates predominantly from the Zechstein Group and to a lesser extent from the Germanic Trias.

The static Young’s modulus E _stat of the sandstone (reservoir) rocks depends primarily on the depth (and hence the porosity), increasing from 15 GPa at 1 km depth to 20–40 GPa at 3 km depth. Note that the scatter in the Young’s moduli is significant, likely due to different mineral content, burial history, etcetera. The dynamic Poisson’s ratio of the sandstone reservoirs decreases with depth. For both E _stat and ν _dyn, no distinct trends with depth can be defined for different stratigraphic ages.

Measured Young’s moduli in other lithologies, including many of the caprocks, show a comparable trend, increasing with depth. Values of E _stat for clay-rich rock (claystone and shale) overlap with those measured in the sandstone intervals, although the Altena Clay has a relatively low stiffness at larger depth (Fig 5c). The values in the Triassic formations appear to fall slightly above the general trend, possibly because many of the Triassic rocks experienced a deeper burial depth than the present-day depth, in for example, the West Netherlands Basin, which results in lower porosities (Van Kempen et al., 2018) and hence stiffer rocks. On average, the E _stat of the Altena clay intervals appears to be lower than that of the other litho-stratigraphic groups at comparable depths. The dynamic Poisson’s ratio of the clay-rich, evaporitic and carbonate rocks is on average higher than that of the sandstones, and the decrease of Poisson’s ratio with depth is much less apparent than for sandstone. Poisson’s ratio appears to be higher in the Rijnland (KN) and Altena (AT) Groups compared to the Triassic, perhaps also related to the larger maximum burial depths seen in some of the Triassic formations.

Pressure and temperature changes

Hydrocarbon production results in net pressure depletion of the reservoir. Most gas and oil reservoirs in the Netherlands are produced by conventional means. For oil fields, the production water is reinjected to maintain pressure drive. In the Schoonebeek field, enhanced recovery occurs through steam injection. Gas field production occurs mainly by primary recovery, where the gas flows freely to the wells, typically recovering 70–95% of initial gas in place (GIIP). Many of the Dutch gas fields show tank-like behaviour with limited aquifer support (e.g. Tichler et al., 2016). Connectivity in most of the Dutch gas fields is good with 80–90% showing little compartmentalisation due to, for example, sealing faults, large fault offsets, or stratigraphic barriers, in particular for the onshore fields (NAM, 2016; Van Hulten, 2010), though compartmentalisation was inferred for some fields, for example, for the Roswinkel field (Fokker et al., 2012). Later in the life of gas fields, as the reservoir pressure has declined significantly, production water may flood the wells. Enhanced Gas Recovery may be applied by for example injecting nitrogen (Goswami et al., 2018), as used for example in De Wijk Field. During production of the gas fields, the initial reservoir pressures are reduced by up to 97% (Roholl et al., 2021). As the initial pressure increases with depth (Fig 6), the magnitude of final pressure decrease in the hydrocarbon fields also increases with depth up to values over -50 MPa. Note that in particularly in the northern onshore and at depth > 2500 m, fields can be significantly overpressured (Verweij & Hegen, 2015), which leads to an even larger pressure decrease with respect to initial conditions (Fig 6a). In terms of spatial characteristics of pressure changes, history-matched pressure distributions within reservoir compartments or the entire reservoirs can be relatively smooth (e.g. in the Groningen field, Fig 6). This implies that a relatively large fault area may experience a relatively uniform pressure decrease over time.

Fig. 6. Pressure changes in hydrocarbon fields. a) Initial pressures (markers) of onshore gas fields. Lines indicate the pressure change down to the final pressure projected at the end of the fields' lifetimes. Pressure data are taken from producation licenses (Roholl et al., 2021). b) example of spatial distribution of history-matched pressure changes for the Groningen field (data from Bourne & Oates, 2017).

A balanced fluid circulation is maintained during geothermal heat production from sedimentary formations. Hot water is produced from the production well, cooled down by a heat exchanger at the surface and reinjected into the injection well. Contrary to the hydrocarbon reservoirs, initial pressures at the geothermal locations are close to the estimated hydrostatic value (Fig 7a). Pressure changes are much smaller than those in hydrocarbon fields. At the production well, drawdown pressures are in the order of 1 MPa. At the injection well, the injection pressures are limited by the regulator to prevent caprock failure and leakage to shallower formations (State Supervision of Mines & TNO-AGE, 2013). According to this protocol, the maximum injection pressure should not exceed a gradient of 13.5 MPa/km. An additional correction is made for the minimum injection temperature when this is lower than ΔT -40°C, with a reduction of 0.1 MPa in pressure per degree below -40°C. Injection pressures and temperatures can however deviate from this value when safe production has been substantiated by extra research. At 2–3 km, the maximum wellhead pressures are typically in the range of 4–7 MPa above the initial reservoir pressure (Fig 7a). Numerical modelling indicates the pressure changes decay logarithmically from the injection well (Buijze et al., 2022; e.g. Willems, 2017; Willems et al., 2017). Corresponding maximum flow rates are in the range 150–400 m³/hr, but average flow rates can be substantially lower as, for example, heat demand is seasonal (Fig 7c). Note that methane gas is co-produced from the formation waters during geothermal operations. In 2021, the total amount of gas co-produced was 22 million Nm³ (Ministry of Economic Affairs & Climate, 2021), comparable to the annual tail end production from some of the smaller gas fields. The gas is dissolved in the formation waters, but returns to the gaseous phase as the pressure in produced fluids decreases and gets below the bubble point as it is pumped up to the surface. This gas is usually captured and burned to generate extra heat, but sometimes the gas is left in solution by maintaining high enough pressures in the surface system. The gas co-production does not lead to a volume change since it was present in dissolved form.

Fig. 7. Pressure and temperature changes in Dutch geothermal doublets. Colors indicate litho-stratigraphic group as in e.g. Fig. 4.  a) Initial pressure (markers) at the top of the targeted reservoir and maximum pressure change at the injector (lines). A hydrostatic gradient of 10.5 MPa/km is drawn for reference. The shaded area indicates the approximate window for injection pressures as set by the regulator (State Supervision of Mines, & TNO-AGE, 2013), note that local injection pressure constraints can deviate depending on a.o. salinity. b) Initial temperature (markers) at the top of the reservoir and the maximum temperature change at the injector (lines). c) Maximum and average hourly flow-rates (data from 2020). d) Estimated radii of the cooled front as a function of production time assuming average flow rates (Equation 1). Data from production licenses: see www.nlog.nl, www.mijnbouwvergunningen.nl, and Supplementary Materials 1.

Contrary to conventional hydrocarbon production, the geothermal operations are accompanied by significant temperature changes ΔT with respect to formation temperature. The injection temperatures are typically 30–35°C, with a minimum ΔT of -40°C set by the regulator to avoid failure of the overlying formations (State Supervision of Mines & TNO-AGE, 2013). Injection temperatures can be lower, in particular now that using a heat pump starts to become more common, but these lower injection temperatures lead to the aforementioned reduction of the maximum allowable injection pressure. The resulting ΔT range from -90 to -30°C for the deeper (>1500 m) doublets (Fig 7b). Heat is transferred from the porous rock volume to the cold injected fluids, and over time the rock volume around the injection well cools down. If the so-called cooling front (boundary of the cooled reservoir volume) reaches the production well (thermal break through), the doublet approaches the end of its productive life. Model results indicate that after an initial period, the cooled rock volume around an injection well expands over time whilst the amount of cooling in the centre of the cooled volume remains relatively constant (Buijze et al., 2022). Under the assumption of lateral convective flow only, the extent of the cooling front can be approximated by (Koning, 1988):

(1) $${r_{cooled}} = {\left( {{{{c_{fluid}}{\rho _{fluid}}q{t_{inj}}} \over {{c_{rock}}{\rho _{rock}}h\pi }}} \right)^{0.5}}$$

where c _fluid and ρ _fluid are the specific heat capacity (Jkg^-1K^-1) and density (kgm^-3) of the injection water respectively, c _rock and ρ _rock are the specific heat capacity and density of the rock formation respectively, q is the flow rate (m³s^-1), t _inj is the injection time (s) and h is the reservoir height (m). Assuming a salinity of 10%, the specific heat capacity of the injected water is c _fluid is 3771 Jkg^-1K^-1 and fluid density is 1052 kgm^-3 (Batzle & Wang, 1992). For an average c _rock of 850 Jkg^-1K^-1, rock density of 2200 kgm^-3, injection time of 35 years, average flow rates and net reservoir thickness, the radii of the cooled volumes for currently producing doublets are 250–1200m. On average, the cooled radii of the doublets in the Rotliegend are smaller, as the net thickness is somewhat larger (Fig 4). The expected areal extents of the cooled volume range from 0.2 to 4 km² which is on the lower end of the areas of some of the smaller gas fields. Several larger gas fields, however, span areas of several 10’s of km², with the Groningen field spanning an area of 900 km².

Summary of observed induced seismicity in relation to reservoir litho-stratigraphy

Since 1986, more than 1800 induced seismic events linked to hydrocarbon production have been recorded by the national network, ranging up to magnitudes of M_L 3.6 (Dost et al., 2017; Muntendam-Bos et al., 2022). The current detection threshold of the network is M 1.5 in many areas of the Netherlands, down to M_L 0.5 in the north of the Netherlands (KNMI Data Platform, Netherlands Seismic Station Magnitude Completeness Map, retrieved: 2022). The distribution of the KNMI catalogue of induced events (www.knmi.nl; KNMI, 2022) is shown in Fig 8a, superimposed on the locations of hydrocarbon fields and geothermal projects. In total, 22 hydrocarbon fields are associated with M_L > 1.5 (Van Thienen-Visser et al., 2012); 11 of these are associated with M > 2.5 events (Fig 8). Magnitudes exceeding M_L 3.0 have occurred in three gas fields: the Groningen gas field with an M_max 3.6, 16 August 2012 (Dost & Kraaijpoel, 2013), the Bergermeer gas field (M_max 3.5) and the Roswinkel gas field (M_max 3.4), all producing from the Upper Rotliegend Group or Triassic. All seismogenic fields related to M > 1.5 events are fields producing from Rotliegend, Zechstein, or the Triassic reservoirs, whereas no seismicity has been observed for younger Jurassic, Cretaceous and Tertiary reservoirs (Fig 8b). Note there may be some ambiquity to which reservoir the events are linked for some stacked fields such as Coevorden and Dalen that produce from both the Limburg and Zechstein Group. The seismogenic reservoirs are mostly located at depths > 2000m. The seismogenic reservoirs are linked to particular tectonic regions in the Netherlands, namely the northern tectonic regions (Rotliegend reservoirs), the Lower Saxony Basin (Zechstein, Triassic, and Rotliegend reservoirs, and potentially the Upper Carboniferous) and the Broad Fourteens Basin & North Holland Platform (Triassic and Rotliegend reservoirs). Although Triassic reservoirs are also found in the West Netherlands Basin, no events with M_L > 1.5 have been reported in region, neither for the Jurassic or Cretaceous oil and gas fields. However, we note that a number of the oil fields in the West Netherlands Basin have ceased production for decennia and monitoring thresholds may have been higher than in the north of the Netherlands. On the other hand, the pressure decrease in oil fields will be much lower than for gas reservoirs. For gas reservoirs, the onset of seismicity typically occurs after significant depletion, mostly after 50% of pressure reduction with respect to the initial reservoir pressure (ΔP/P_ini > 0.5) (Fig 8c). The first recorded events can be relatively large, without any smaller events occurring before. In Bergermeer, for example, the first recorded event was a M_L 3.0 events (M_max 3.5) and in the Roswinkel field the first recorded event had a magnitude of M_L 2.7 (M_max 3.4).

Fig. 8. Occurrence of induced seismicity in relation to operations. a) Map view of induced seismicity with M > 1.5  (black circles). Catalog obtained from www.knmi.nl (December 2022). Hydrocarbon shows and geothermal doublets (triangles) are indicated, colors according to litho-stratigraphic. b) Northing (Rijksdriehoek coordinates, see e.g. Fig. 2) versus depths of onshore hydrocarbon fields and geothermal doublets. Squares: hydrocarbon shows, circles: geothermal projects. Filled symbols indicate fields/projects where M > 1.5 events have been recorded, hollow symbols where no M > 1.5 have been recorded. b) Relative amount of gas depletion vs occurrence of first seismic event with M>1.5.

For geothermal projects, on the other hand, only one M_L > 1.5 event has been reported, a M_L 1.7 event near the doublets in Californië, in the south-east of the Netherlands. These doublets circulated fluid through karstified and faulted carbonates and quartzites of the Lower Carboniferous and are located in the Ruhr Valley Graben, a tectonically active area. The event magnitude is at the border of what could have been felt, but doublet operations have been suspended in fear of larger events (Muntendam-Bos et al., 2022; Vörös & Baisch, 2019). For geothermal operations in sandstone reservoirs, no M > 1.5 events have been reported. In vicinity of a geothermal doublet in the West Netherlands Basin, a M_L 0.0 event has been detected by a local microseismic monitoring network; the relation of the event to operations remains however unclear (Muntendam-Bos et al., 2022).

Discussion

In the previous sections, we have reviewed the similarities and differences between onshore hydrocarbon reservoirs and geothermal reservoirs in the Netherlands. The same litho-stratigraphic reservoirs are targeted by hydrocarbon and geothermal, notably the Rotliegend, Triassic and Jurassic/Cretaceous sandstones, but differences also exist in targeted porosity and permeability, depth range, and the pressure and temperature changes related to both activities. In the following, we discuss how these differences could affect fault reactivation potential and relate this to the observed seismicity presented in the previous section. First, a short review of poro- and thermo-elastic stress changes is given to aid the discussion.

Poro- and thermo-elastic stress changes

Both pressure and temperature changes relating to human operations lead to stress changes in the subsurface and thus on pre-existing faults in the subsurface. These stress changes add to the initial, tectonic state of stress on the fault and can either increase of decrease the fault stress. When fault stresses exceed the fault strength, fault slip can occur, so-called fault reactivation. Pressure and temperature changes result in the following (combinations of) stress changes in porous media:

• Poro-elastic stress changes: the pressure changes cause a volume change of the pressurised medium, for example, a reduction in pressure like occurs for gas production causes the porous reservoir rock to contract, and vice versa. This will affect mainly the total horizontal stress, with depletion leading to a reduction of total horizontal stress. Note that the effective horizontal stress can still increase in that case due to the reduction of pore pressure.

• Thermo-elastic stress changes: similarly, the temperature change also causes a volume change, with cooling resulting in contraction of the reservoir and a reduction in horizontal stress.

• The direct pressure effect: an increase in fault pressure (e.g. near an injection well) will result in a decrease in effective fault normal stresses and destabilisation of faults, whereas a decrease in fault pressure (e.g. near a production well) will stabilise the fault. This effect mainly plays a role for injection into fractured or low permeability formations.

Note that pressure changes, poro- and thermo-elastic stress changes can occur simultaneously, with one dominating over the other depending on location and time (Wassing et al., 2021). In addition, pressure and temperature changes may cause changes in the elastic properties of the reservoir, and chemical changes could influence in fault friction and cohesion, which also affects the reactivation potential (e.g. Hunfeld et al., 2019; Pluymakers et al., 2014). The stress changes induced by hydrocarbon production are dominated by poro-elastic stressing (Buijze et al., 2019b; Roest & Kuilman, 1994; Segall & Fitzgerald, 1998), whereas stress changes in geothermal doublets in porous sandstone reservoir are dominated by hermos-elastic stress changes, as the pressure changes are relatively limited and decrease rapidly with distance from the well (e.g. Buijze et al., 2022; Kivi et al., 2022). Note however that the direct pressure effect plays a role near the well, and for the faulted reservoirs or reservoirs overlying fracture (basement) rocks (Baisch et al., 2010; Hsieh & Bredehoeft, 1981; Keranen et al., 2013). Pressure changes likely played a dominant role for the induced events in the faulted carbonates at Balmatt (Kinscher et al., 2023), although for the events near Californië, Limburg, it is argued that thermo-elasticity and pressure decreases are responsible for generating the observed events (Vörös & Baisch, 2022). For smaller fracture spacing fractured rock masses may behave more as a porous rock mass (Wassing et al., 2021).

First-order magnitudes of poro- and thermo-elastic stress changes

For the comparison between hydrocarbon production and geothermal, it is insightful to consider the differences in magnitudes and directions of poro- and thermo-elastic stress changes. Both poro-elastic and thermo-elastic stress changes depend on the poro- and thermo-elastic parameters of the medium and surroundings, as well as on the shape of the pressurised or cooled volume. For simplified geometries, analytical expressions of the stress changes have been formulated. Poro-elastic horizontal and vertical effective stress changes can be expressed as (Soltanzadeh & Hawkes, 2009):

(2) $$\begin{gathered} \frac{{\Delta {\sigma _h}}}{{\alpha \Delta {\text{P}}}} = \;{\gamma _h};\;\frac{{\Delta {\sigma _H}}}{{\alpha \Delta {\text{P}}}} = \;{\gamma _H};\frac{{\Delta {\sigma _v}}}{{\alpha \Delta {\text{P}}}} = \;{\gamma _v};\;\frac{{\Delta {\sigma _h}^\prime }}{{\alpha \Delta {\text{P}}}} = - \left( {1 - {\gamma _h}} \right);\;\frac{{\Delta {\sigma _H}^\prime }}{{\alpha \Delta {\text{P}}}} \\ \!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\!\! = - \left( {1 - {\gamma _H}} \right);\;\frac{{\Delta {\sigma _v}^\prime }}{{\alpha \Delta {\text{P}}}} = - \left( {1 - {\gamma _v}} \right) \\ \end{gathered}$$

where α is the Biot coefficient and ΔP is the pressure change. γ _h and γ _v are the so-called stress arching ratios or stress path parameters. For thermo-elastic stressing, these are expressed as:

(3) $${{{\rm{\Delta }}{{\rm{\sigma }}_h}} \over {\eta {\rm{\Delta T}}}} = \;{\gamma _{hT}};\;{{{\rm{\Delta }}{{\rm{\sigma }}_H}} \over {\eta {\rm{\Delta T}}}} = \;{\gamma _{HT}};\;{{{\rm{\Delta }}{{\rm{\sigma }}_v}} \over {\eta {\rm{\Delta T}}}} = \;{\gamma _{vT}};\;$$

where η is the linear thermal expansion coefficient. The thermo-elastic ratios relate to the poro-elastic arching ratios as:

(4) $${\gamma _T} = {\rm{\;}}{{\rm{E}} \over {\left( {1 - v} \right){\rm{\;}}}}\gamma ;$$

The arching ratios depend on the shape of the pressurised/cooled volume. For example, for penny-shaped reservoirs with low height/width ratio, the poro-elastic arching ratios are as follows:

(5) $${\gamma _h} = {\gamma _H} = \;{{1 - 2\nu } \over {1 - v}}\left( {1 - {{\pi e} \over 4}} \right);\;{\gamma _v} = {{1 - 2\nu } \over {1 - v}}{{\pi e} \over 2}$$

where ν is Poisson’s ratio and e is the aspect ratio (height/width). For high aspect ratio, the vertical stress changes tends to zero, as is visible in the expression for laterally extensive reservoirs (width >> height):

(6) $${\gamma _h} = {\gamma _H} = \;{{1 - 2\nu } \over {1 - v}};\;{\gamma _v} = 0$$

Note in these formulas, a medium with uniform elastic properties is assumed. Hence, poro-elastic stress changes are primarily dependent on the Poisson’s ratio and Biot coefficient in combination with ΔP, and thermo-elastic stress changes on Young’s modulus and the linear thermal expansion coefficient in combination with ΔT. The resulting horizontal and vertical stresses can be converted to fault effective normal and shear stress components σ _n ’ and τ which determine whether fault slip can occur through the Mohr–Coulomb criterion:

(7) $${\tau _f} \le {\rm{\tau \;where\;\tau }} = {\rm{\;}}C + \;\mu ({\sigma _n} - {\rm{\;P}})$$

An example of depletion-induced and cooling-induced fault stress changes and direct pressure effect are shown in Fig 9 for a 70° dipping fault in a normal faulting environment.

Fig. 9. Comparison of poro-elastic and thermo-elastic stress changes and sensitivity to various geomechanical parameters, for simplified lateral extensive reservoir geometry. Input values: normal faulting regime, depth: 2500 m, vertical stress gradient: 22 MPa/km, minimum horizontal stress gradient Δσ_h /Δy: 16 MPa/km, pore pressure gradient: 10.7 MPa/km, horizontal stress ratio σ_H/σh: 1.07, Young’s modulus E: 15 GPa, Poisson’s ratio ν: 0.15, fault dip θ: 70 degrees fault strike w.r.t. σ_H ϕ: 0, friction coefficient μ: 0.6, and cohesion C = 0 MPa. Poro-elastic stress changes are computed for a -20 MPa pressure drop with pressures within the fault following those of the reservoir, thermo-elastic stress changes for -20 degrees temperature decrease (cf.  Soltanzadeh & Hawkes, 2009). a) poro-elastic stress change b) thermo-elastic stress change, c) direct pressure increase in the fault, with initial stress state in black and final stress state depicted by the dotted line. Black marke: initial fault stress, colored line: stress path, colored maker: final fault stress. Sensitivities of poro-elastic stress changes are shown for  d) Poisson ratio, g) dip, and j) strike, and sensitivities of thermo-elastic stress are shown for e) Young’s modulus, h) dip, and k) strike. Furthermore effect of horizontal stress gradient (f), horizontal stress ratio (i) and effect of fault frictionand cohesion (l) on initial stress are shown.

The largest difference between the stress change magnitudes for depletion and cooling is the sign of the normal stress change. For depletion, the poro-elastic effect leads to an increase in effective normal stress and an increase in shear stress. Depending on the Poisson’s ratio, fault stresses become closer or farther away from failure (Fig 9d). For thermo-elastic stressing on the other hand, shear stress increases but the effective normal stress decreases (Fig 9b). The direction of the resulting stress path is more destabilising than that of the poro-elastic stress changes. The magnitude of the stress change is linearly dependent on Young’s modulus (Fig 9e) and the linear thermal expansion coefficient. For both poro- and thermo-elasticity fault, stress changes also depend on the fault dip (Fig 9g,h) and fault strike (Fig 9j,k). Note that the stress changes can locally become concentrated and deviate from and/or exceed those summarised in equations (3)–(7) as a result form complex geometry or material behaviour (see next section). Furthermore, reactivation potential is not only determined by the stress change but also by the initial stress on the fault, which is a function of σ_h’ and σ_v’ and the fault dip and strike, and σ_H’. The larger the differential stresses σ_v’ – σ_h’ (smaller ratio σ_h’ / σ_v’,) the more unstable the initial state of stress (Fig 9f). Fault strength τ _f on the other hand is influenced by the friction and cohesion (Fig 9l) on the fault. To assess reactivation potential, knowledge of the pressure and temperature changes, poro- and thermo-elastic parameters as well as the initial stress and fault properties are essential. In addition, not only stress change magnitudes but also spatio-temporal evolution of the stress changes on the fault and the interplay with fault properties should be considered to properly assess seismicity potential.

Seismicity in Dutch hydrocarbon fields

Recorded seismic events related to hydrocarbon reservoirs with magnitudes large enough to be felt are limited to Triassic, Zechstein and Rotliegend fields in the northern half of the country at depths > 2 km and occur after substantial decrease in pressure (Fig 8), in line with previous observations (e.g. Muntendam-Bos, 2021; Van Eijs et al., 2006; Van Wees et al., 2014). Several studies have been performed to identify what would cause fields to behave seismically or not. In a statistical analysis of geological and operational parameters in relation to seismicity to assess a priori, the probability of a field becoming seismically active (Deterministic Seismic Hazard Analysis Induced Seismicity, DHAIS) three key parameters were recognised: the aforementioned pressure drop, stiffness ratio between reservoir and overburden, and fault density (Van Eijs et al., 2006). The cut-off values for these three factors were updated several times (Roholl et al., 2021; Van Thienen-Visser et al., 2012), and the analysis is still used to assess seismicity potential in small gas fields (all except Groningen) (State Supervision of Mines, 2016). The correlation with the stiffness ratio might be a causal one, as a stiffer overburden promotes fault reactivation (Mulders, 2003). However, reservoir–seal pairs with stiffness ratio above the cut-off value (predominantly the Rotliegend, Zechstein Carbonate and the Triassic fields) also share other communalities that could influence seismicity, such as the presence of a thick salt layer, depositional environment and related mineralogy, and location in the northern half of the country where, for example, stresses could be different compared to the south (Bakx et al., 2022; Verweij et al., 2016). The causal relation is thus not straightforward. Modeling studies indicate that mechanisms for production-induced fault reactivation include, on top of poro-elastic stressing described in 9.1, stress concentrations along faults offsetting or bounding the reservoir (Buijze et al., 2017; Buijze et al., 2019b; Mulders, 2003; Nagelhout & Roest, 1997; Orlic & Wassing, 2013; Roest & Kuilman, 1994; van den Bogert, 2015; van den Bogert, 2018), salt creep lowering normal stresses on offset faults (Orlic & Wassing, 2013; Wassing et al., 2017) and pressure diffusion in faults (Zbinden et al., 2017). Application of 2D modelling of fault reactivation to the onshore gas field resulted in a match for seismically active fields but overpredicted fault reactivation for fields currently showing no seismicity, in particular those in the south-west (Vörös & Baisch, 2018). The mismatch was explained as a result of higher in situ horizontal stress in this region, but this is contrary to a recent study showing lower horizontal stress gradients in the south-west (Bakx et al., 2022). Several other important model assumptions may have played a role for the (mis)match; fault dip and offset were not available as field-specific parameters, whereas these have a large effect on the stress path (see e.g. Fig 9g), variations in magnitude and orientation of the regional stress field were not considered in detail (Verweij et al., 2016), and pressure in bounding faults did not follow the reservoir pressure, which make the stress path more destabilising (see e.g. Zbinden et al., 2017) and thus will tend to overpredict fault reactivation. In a later study, the observed variations of the in situ stresses were in part accounted for, but the results did not match the observed seismic vs non-seismic signatures of the hydrocarbon fields (Muntendam-Bos, 2021). Instead, the authors proposed the destabilising effect of salt creep on the pre-depletion stress on faults with offset (Orlic & Wassing, 2013; Wassing et al., 2017) as the likely mechanism to explain the occurrence (or lack) of seismicity. Note however that also in this study, field-specific fault dips and offsets were not considered, nor variations in stress orientation or fault properties. In another study, the same regional clustering of seismicity was recognised, but here it was hypothesised that differences in fault strength (cohesion) could explain the lack of seismicity in some fields (de Pater et al., 2020). Note that none of the models above capture the effect of fault zone lithology on the mode of fault slip (seismic vs aseismic fault slip). The occurrence of aseismic fault slip could be another reason for the mismatch between observations and these models, which currently always treat reactivation as being seismic reactivation.

To summarise, despite the fact that a number of studies were devoted to the problem, it is still unclear what exactly causes some fields to cause induced seismicity and others not. Note though, that apart from various geomechanical aspects mentioned above, variability may also arise from the limited sample number, and differences in detection limits in the monitoring network. However, it is generally recognised that there is a correlation with the geological framework, with seismicity occurring in the older reservoir rocks in the northern half of the country (cf. Fig 8a), and a lower potential for induced earthquakes in the south-west and in the younger litho-stratigraphies. It is recommended to update the physics-based modelling with latest insights on the region-, depth- and lithology- specific in situ stresses (Bakx et al., 2022), appropriate fault dips and orientation, and effects of salt creep and perform robust statistical analysis on this. If the seismic vs aseismic behaviour still cannot be explained, this could point to other geological and geomechanical factors that should be taken into account such as the (fault and reservoir) mineralogy, overpressure, elastoplasticity and seismic vs aseismic fault slip.

Differences between hydrocarbon and geothermal operations and effect on seismicity potential

For adequate, region-specific (refinement of) seismic screening methodologies, it is useful to learn from the similarities and differences with the hydrocarbon fields. Also for geothermal projects, an a priori assessment of the seismicity potential is important. Currently, a protocol is in place where penalty scores are given for data quality, proximity to faults and location in the Ruhr Valley Graben (Baisch et al., 2016). The seismic risk analysis protocol will be updated in the near future.

In terms of geological characteristics, localities in sandstone reservoirs currently exploited for geothermal have a relatively high porosity and permeability compared to the hydrocarbon fields (Fig 4). This is expected to translate into a lower Youngs modulus, which is important to consider when assigning elastic parameters in, for example, geomechanical models of geothermal operations. So far, sonic logs in geothermal wells are sparse, but if more well logs are taken in the future the range of elastic parameters for geothermal reservoirs, and their dependence on depth and depositional environment, can be better constrained. Another geological difference is that the operational geothermal projects are situated mostly in regions where salt is absent (West Netherlands Basin) or thin (around the Texel-IJsselmeer High). If salt creep leads to a pre-production state of stress that is closer to failure (see 9.2), this would imply a lower seismic potential for geothermal projects in these regions with little or no salt. This will need to be substantiated by modelling and field measurements. The biggest differences between hydrocarbon and geothermal are however operational, with cooling and an elevated injection pressure around the injection well (Fig 7) leading to markedly different spatio-temporal stress evolution compared to pressure depletion hydrocarbon fields (Fig 6). The direction of stress change is more destabilising for geothermal than for hydrocarbon reservoirs, and stress change magnitudes can be substantial (Fig 9), as is shown in numerical modelling studies (e.g. Candela & Fokker, 2017; Hassanzadegan et al., 2011; Kivi et al., 2022; Van Wees et al., 2020). Even so, no felt seismic events have yet been recorded near geothermal doublets in porous sandstones. This can be due to the gradual growth of the cold front with injection time, with cooled volumes that can still be relatively small as most doublets have not been operational for more than 10 years (Fig 7d). Doublets in the Netherlands are typically placed several hundreds of metres from known faults, and the cold front may not have reached these faults yet. This is different from hydrocarbon fields which are found in structural traps and are therefore often bounded by faults or contain intra-reservoir faults which are affected by the pressure change over a relatively large fault area. Alternatively, the cold front may have reached faults, but as the fault area experiencing a stress change grows gradually generated events might be small. Also, many doublets are situated in the Upper Jurassic/Lower Cretaceous rocks (West Netherlands Basin). Hydrocarbon fields in these rocks have also not been related to any felt seismic events, which may suggest these formations to have a lower seismogenic potential. In addition, the lack of seismicity for geothermal projects may also be a reflection of the initial state of stress which is not considered critical in a large part of the Netherlands (Bakx et al., 2022; Van Wees et al., 2014; Verweij & Hegen, 2015); also the hydrocarbon fields generate felt seismicity only after significant depletion and stress changes. Still, considering the expected large cooling-induced stress change, it is important to keep monitoring the effects of geothermal operations on the subsurface in the coming years.

Geothermal operations in fractured carbonates of the Dinantian are distinctly different and share few similarities with, for example, the carbonate hydrocarbon reservoirs within the Zechstein. Accessible depths for geothermal operations are found in the south-east, in the Ruhr Valley Graben where natural seismicity occurs. Whereas faults are avoided for the sandstone projects, faults, fractures and karsts provide essential permeability for fluid flow in the otherwise low matrix porosity carbonates. Operations near the two doublets in the south-east have been associated with one M 1.7 event. Pressure diffusion likely plays a large role, but likely also poro-elastic stressing and thermo-elastic stressing since a temporal correlation with operation stop was observed (ter Heege et al., 2020; Vörös & Baisch, 2022). In addition, the initial stress is higher and stress changes can propagate to high stress regions through the fracture network, but also the rock types itself are likely to promote larger stress changes and seismic slip (see 9.4.2). Also at several sites, abroad similar geothermal operations through balanced circulation in carbonates have led to M 2–3 events (Broothaers et al., 2021; Seithel et al., 2019).

Other geological factors that can affect seismogenic potential

Besides the various factors presented and discussed in the previous sections, various other geological aspects are important for seismogenic potential.

Fault rock properties

Another aspect to consider is that fault reactivation can occur seismically or aseismically. For seismic slip, a rapid drop in fault strength with progressive fault slip or increased slip rate is required to drive fast fault slip and emit seismic waves. Without such a rapid drop in fault strength, fault slip will occur aseismically. Fault mineralogy and the degree of cementation are important as it affects both fault strength (friction and cohesion) as the strength drop (e.g. Hunfeld et al., 2020; Hunfeld, 2020). For sandstone reservoirs in the Netherlands, not much is known about the mineralogy of in particular the larger-scale faults, as these are typically avoided during drilling. A few findings suggest the large-scale faults in high net-to-gross ratio Rotliegend reservoirs have cores consisting of cataclastic material (fine-grained broken grains) and compaction bands (Van Hulten, 2010). Analysis of smaller scale fractures and faults within the Rotliegend sandstone reservoirs also shows that fractures are often cataclastic (fine-grained broken grain) and/or cemented with quartz or anhydrite (Fisher & Knipe, 1998; Fisher et al., 2000; Ligtenberg et al., 2011). Cataclastic and cemented faults will have (considerable) cohesion, and fault reactivdation on such faults can lead to seismic slip events, in particular if anhydrite is present (Hunfeld et al., 2019; Hunfeld, 2020; Pluymakers et al., 2014; Scuderi et al., 2013). Anhydrite cements derived from the Zechstein are not only present in the underlying Slochteren Formation but could also be present in the overlying Lower Germanic Triassic sandstones (Purvis & Okkerman, 1996). The presence of cohesion or anhydrite cements could be an additional reason why the Triassic & Rotliegend sandstone reservoir plays appear more prone to seismicity, on top of the effect of salt creep which tends to lower the normal stress on faults (e.g. Wassing et al., 2017). Cementation is also expected to decrease with depth as does rock competency, which may set a depth limit from which we could observe seismic slip; this remains to be investigated. Also carbonate rocks such as those in the south-east of the Netherlands are prone to seismic slip. For the Cretaceous and Jurassic sandstones on the other hand, anhydrite cement is less likely due to the marine and fluvial depositional environments. Unfortunately, no studies are present that document fracture or fault mineralogy in these formations. In more clay-rich, lower net-to-gross formations faults are expected to contain high(er) clay content, as seen in a fault drilled in the more clay-rich region of the Groningen field (Van Hulten, 2010). For clay-smeared faults, no cohesion is expected, and such faults are more prone to aseismic sliding (Carpenter et al., 2016; Hunfeld et al., 2020; Hunfeld, 2020; Ikari et al., 2009; Samuelson & Spiers, 2012). It is recommended to systematically investigate the relation between depositional environment, fault kinematics and diagenesis, with fault rock mineralogy, and to perform experimental research into the behaviour of these different fault rocks at the in situ conditions. Additionally, it is important to assess what the effect of pressure and temperature changes on the fault stability are. Note that fault mineralogy is also important for fault permeability and resulting fluid flow and hence the productivity of geothermal doublets.

Effect of lithology-related variations in in situ stress

In 9.2 and 9.3, we have already touched upon the role of the in situ stress preceding operations, which is a key parameter for reactivation potential and seismicity (e.g. Buijze et al., 2021). In situ stress measurements indicate the average horizontal stress gradients to be higher (more stable) in the north and lower (closer to failure) in the south-west (Bakx et al., 2022). This is somewhat contradictory with the observed seismicity which is mainly in the north. However, it is important that this horizontal stress data are combined with vertical stress data and pressure data, as the ratio of the two effective stresses is the key controlling parameter. The study of Bakx et al. (2022) also shows that the horizontal stress gradient is higher in clay-rich and evaporite-rich formations, which is in agreement with other studies (Breckels & Van Eekelen, 1982; Warpinski & Teufel, 1989). For example, mini-frac measurements in Ten Boer vs Slochteren sandstone indicate a 2 MPa higher horizontal stress in the Ten Boer Claystone (Bakx et al., 2022). This has a substantial effect on the initial fault stability (see e.g. Fig 9). It is thus likely that in situ stress is more stable in clay-rich and evaporitic formations, the seals of some of the hydrocarbon or geothermal plays (e.g. Zechstein salt, Ten Boer Claystone, Rodenrijs Claystone, Vlieland Claystone, etc.). This would limit the potential for events being triggered in or propagating into these formations. It is also of interest to study the initial stress as function of depth and look into the stress and stability of the shallower formations. For an improved, region-specific, estimate of fault reactivation potential and seismicity reliable measurements of the in situ stresses through, for example, mini-frac tests, in both reservoir and caprock formations, are highly recommended, in combination with sensitive monitoring of seismicity to build up a dataset with statistical relevance.

Implications for fault reactivation and induced seismicity

Hydrocarbon production and geothermal doublet operations lead to different pressure and temperature changes in the subsurface and different stress paths. The biggest difference in terms of stress change is the effective normal stress change; this component increases for depleting reservoirs, whereas it decreases in the cooled volume around the injection well (Fig 10). Lowering of the effective normal stresses promotes failure but is also related to more stable fault sliding modes (Dieterich, 1992; Mclaskey & Yamashita, 2017; Ruina, 1983). High local injection pressures during, for example, stimulation also lead to low effective normal stresses, which have been linked to the occurrence of aseismic slip in the near-well region (Cappa et al., 2019; De Barros et al., 2018) and lower stress drops (Goertz-Allmann & Wiemer, 2012). Similarly, for geothermal doublets, the low normal stresses resulting from thermo-elastic stressing in the cooled fault area (Wassing et al., 2021) could lead to aseismic slip rather than seismic slip (Im & Avouac, 2021). Note however that the aseismically slipping fault area does lead to stress transfer and could in turn trigger seismic slip. Another related result of the lower normal stress is that even if seismic slip occurs, the stress drop is smaller for a fault governed by slip- or velocity-weakening friction, leading to less motion at the surface. Microseismic monitoring in regions where it is suspected parts of the fault have been experiencing thermal stressing would help to illuminate whether the lowered normal stress (local fault lithologies) indeed leads to aseismic slip, slow fault slip rather than seismic slip. Not only is this important for geothermal but also essential for other sustainable energy technologies exploiting the subsurface.

Fig. 10. Conceptual figure showing possible spatio-temporal signature of stress build-up with production time on a fault in a permeable gas field (a) and near a geothermal doublet injecting in a porous sandstone reservoir (b). Instruction: is referred to in 9.5.

Observations on (the lack of) seismicity near hydrocarbon fields point to certain regions and plays having a higher seismogenic potential, with the Rotliegend, Triassic and Zechstein reservoirs in the Groningen High, Lauwerszee Through, west of the Central Netherlands Basin, and in the Lower Saxony Basin having a higher, and those in the West Netherlands Basin a lower potential. The higher seismogenic potential also correlates with the presence of thick Zechstein salt. This suggests that in a similar manner, geothermal operations in the same regions may have a higher or lower seismogenic potential. However, mechanisms of why these regional differences exist are not yet fully understood. Also, mechanisms behind stress changes in geothermal doublets are different, and geomechanical modelling of these mechanisms requires validation against field data. This underlines the need for detailed monitoring of the effect of geothermal operations on the subsurface, in particular since most doublets have not been operational yet for a very long time. The spatio-temporal evolution of pressures, temperatures and stresses for geothermal operations does suggest that fault reactivation will happen in a more progressive manner than for the hydrocarbon faults and seismicity could be easier to monitor and to mitigate. As a hydrocarbon reservoir is depleted, a large fraction of fault area within the reservoir may experience stress build-up (e.g. van Wees et al., 2019). This could lead to sudden larger faulting events. In fact, some of the gas fields have generated only a few larger earthquakes, without any smaller magnitude events such as the Bergermeer a M3.5 event, which occurred without any preceding events in the range 2–3 that would have been picked up by the network. The spatio-temporal stressing signature in geothermal doublets is very different. As the cooling front reaches a fault, a progressively larger fault area will experience larger stress changes. It may thus be expected that in the case fault reactivation occurs, and if fault slip is seismic, event sizes remain small or become progressively larger as the cooled fault area grows incrementally (Wassing et al., 2021). This suggests that monitoring can be more effective in preventing larger events than in the case of hydrocarbon depletion.

Conclusions

We have reviewed and summarised similarities and differences in the geological characteristics of hydrocarbon reservoirs and geothermal reservoirs in sandstone formations in the Netherlands, looking in particular at the geology, formation depths and thicknesses, porosity/permeability, elastic parameters, pressure and temperature changes, and observed seismicity. Based on this comparison, we have discussed the implications for fault reactivation and induced seismicity, based on geomechanical theory and literature. Our main findings are summarised below.

• The same sandstone formations form primary plays for both hydrocarbon and geothermal, notably the Slochteren Formation, and Jurassic/Cretaceous sandstones, Triassic sandstones and to a lesser extent Tertiary sandstones.

• Geothermal projects are not always developed in the same area, as they do not require a specific source rock – reservoir – caprock sequence and boundary conditions for economic production are different. Most hydrocarbon fields are found in the northern half of the Netherlands, as well as the south-west. Operational doublets are found predominantly in the south-west and centre of the Netherlands, at a more limited depth range (<3 km) than for hydrocarbon fields (<4.5 km), and have higher values of formation-averaged porosity and permeability than hydrocarbon reservoirs. Whereas hydrocarbon fields are mainly found in structural highs, doublets often target structural lows where syn-sedimentary formations tend to be thicker and thus more productive.

• Hydrocarbon depletion results in a pressure decrease in the reservoirs, whereas geothermal operations on the other hand occur through balanced circulation. Pressure changes in hydrocarbon reservoirs are large, decreasing down to < 10% of the initial pressure (several 10’s of MPa). Pressure changes in doublets are limited, but temperature decreases are substantial (up to several 10 of °C of cooling). In both cases, maximum (allowable) pressure change and temperature change increase with depth.

• Twenty-two hydrocarbon shows have been linked to seismic events with M > 1.5. These events occur in Triassic, Zechstein, or Rotliegend reservoirs, at depth > 2 km, located in the centre and northern half of the Netherlands, where substantial Zechstein evaporites are present. No M > 1.5 have been recorded near fields in the Jurassic, Cretaceous and Triassic formations in the south-west of the Netherlands, nor in the Cretaceous reservoirs in the Friesland Platform.

• The cause behind this spatial ‘clustering’ of seismically active and inactive fields remains unclear despite various studies. Besides limited sample size and a potential bias in monitoring, the effect of salt creep on fault stability, variations in in situ stress, variations in fault geometry, and play-dependent variations in fault strength and stability (governing seismic vs aseismic slip behaviour) are important factors of influence, and their effect should be further investigated.

• No seismic events with M > 1.5 have been linked to geothermal operations in sandstone reservoirs in the south-west and centre of the Netherlands.

• Stress change magnitudes due to depletion and cooling can both be significant. The main difference is that effective stress becomes more compressive during depletion but becomes less compressive during cooling. This causes a destabilising stress path but may also lead to more stable fault slip modes and a lower stress drop. The second difference is that for cooling a progressively larger fault area will be stressed when the cooling front reaches a fault, which suggests that if fault slip is seismic, event sizes may remain restricted or increase gradually over time.

• Sensitive monitoring is essential to validate the geomechanical models, understand the reactivation process and mitigate the potential occurrence of larger events by design of adequate mitigation measures. Understanding of the reactivation process can help to define a more detailed region- and depth-specific estimate of seismogenic potential of geothermal projects and other sustainable subsurface technologies.