Introduction
The present-day structure and geological evolution of the northeastern Netherlands have been studied by many authors concerning sedimentology, stratigraphy, tectonics, glacial effects, burial and temperature history, as well as petroleum systems (e.g. Rondeel et al., Reference Rondeel, Batjes and Nieuwenhuijs1996; Duin et al., Reference Duin, Doornenbal, Rijkers, Verbeek and Wong2006; Wong et al., Reference Wong, Batjes and de Jager2007; Littke et al., Reference Littke, Bayer, Gajewski and Nelskamp2008; Doornenbal & Stevenson, Reference Doornenbal and Stevenson2010; Nelskamp, Reference Nelskamp2011; Scheck-Wenderoth & Maystrenko, Reference Scheck-Wenderoth and Maystrenko2013; Strozyk et al., Reference Strozyk, Reuning, Scheck-Wenderoth, Tanner, Soto, Flinch and Tari2017; Kley, Reference Kley2018). Large-scale 3D basin models studying the geological and thermal history and the petroleum systems in Central Europe have been conducted by Schroot et al. (Reference Schroot, van Bergen, Abbink, David, van Eijs and Veld2006), Munoz et al. (Reference Munoz, Littke and Brix2007), Uffmann & Littke (Reference Uffmann and Littke2011), Abdul Fattah et al. (Reference Abdul Fattah, Verweij, Witmans and Ten Veen2012), Verweij et al. (Reference Verweij, Echternach, Witmans, Abdul Fattah, Peters, Curry and Kacewicz2012a), Bruns et al. (Reference Bruns, Di Primio, Berner and Littke2013), Bruns et al. (Reference Bruns, Littke, Gasparik, van Wees and Nelskamp2016), Mohnhoff et al. (Reference Mohnhoff, Littke and Sachse2016), Sachse & Littke (Reference Sachse and Littke2018), and Amberg et al. (Reference Amberg, Back, Sachse and Littke2022). Significant hydrocarbon accumulations are the large Schoonebeek oil field and the giant Groningen gas field, both discovered as early as 1943 and 1959, respectively (De Jager & Geluk, Reference De Jager, Geluk, Wong, Batjes and de Jager2007). There are several active petroleum systems, out of which the Palaeozoic petroleum system, which is mainly sourced by thermogenic gas from Carboniferous coal-bearing strata, is by far the most important one (Groetsch et al., Reference Groetsch, Sluijk, Van Ojik, De Keijzer, Graaf, Steenbrink, Grötsch and Gaupp2011). Gas accumulations are mainly found in Rotliegend reservoirs and are sealed by Permian shales and Zechstein salt, while Mesozoic petroleum systems are oil-prone and are found in the Jurassic sub-basins such as the Lower Saxony Basin (De Jager & Geluk, Reference De Jager, Geluk, Wong, Batjes and de Jager2007). Mesozoic source rocks are the Toarcian Posidonia Shale and the Berriassian Coevorden Formation, also referred to as Wealden Shale in Germany (Leythaeuser et al., Reference Leythaeuser, Littke, Radke and Schaefer1988; Rippen et al., Reference Rippen, Littke, Bruns and Mahlstedt2013).
Due to the complex nature of 3D basin and petroleum systems models, most of the previously listed studies did not consider the impact of glacial periods on the subsurface. The imposed mechanical loading and unloading of ice sheets with low surface temperatures and the possible formation of gas hydrates can have effects on the pressure and temperature of the subsurface as well as on fault reactivation and the generation, migration, accumulation and leaking of hydrocarbons (Kjemperud & Fjeldskaar, Reference Kjemperud, Fjeldskaar, Larsen, Brekke, Larsen and Talleraas1992; Lerche et al., Reference Lerche, Yu, Tørudbakken and Thomsen1997; Lang et al., Reference Lang, Hampel, Brandes and Winsemann2014). The influence of glacial loading and unloading on petroleum systems has been studied in areas such as the Barents Sea and Northern Germany by Lerche et al. (Reference Lerche, Yu, Tørudbakken and Thomsen1997), Cavanagh et al. (Reference Cavanagh, Di Primio, Scheck-Wenderoth and Horsfield2006), Grassmann et al. (Reference Grassmann, Cramer, Delisle, Hantschel, Messner and Winsemann2010), Rodrigues Duran et al. (Reference Rodrigues Duran, di Primio, Anka, Stoddart and Horsfield2013) and Sachse & Littke (Reference Sachse and Littke2018). Results of these studies are modelled gas leakages during interglacial and glacial retreats due to changing pressure conditions of the gas phase and an ongoing cooling of reservoir temperatures during glacial times to the present-day (Cavanagh et al., Reference Cavanagh, Di Primio, Scheck-Wenderoth and Horsfield2006; Grassmann et al., Reference Grassmann, Cramer, Delisle, Hantschel, Messner and Winsemann2010). Today, low present-day thermal gradients of 20°C/km are observed down to a depth of 0.5 km in the shallow subsurface of the Netherlands, while higher average thermal gradients of around 30–35°C/km are found in deeper intervals from a depth of 0.5 km down to 3 km (Bonté et al., Reference Bonté, Van Wees and Verweij2012; Ter Voorde et al., Reference Ter Voorde, Van Balen, Luijendijk and Kooi2014). Ter Voorde et al. (Reference Ter Voorde, Van Balen, Luijendijk and Kooi2014) proposed an influence of the cold period during the Pleistocene Weichselian glacial stage on the present-day temperatures in the shallow subsurface.
The recently published 3D basin and petroleum systems model of the northeastern Netherlands (Amberg et al., Reference Amberg, Back, Sachse and Littke2022) was used as the basis for this study. As there is no published 3D basin model with implemented glaciations for the northeastern Netherlands, the 3D model presented in this study integrates the sequential Pleistocene ice sheet advances and retreats leading to mechanical loading and unloading of the underlying sedimentary basin. Thus, this study allows to evaluate ice sheet effects on the thermal and pressure history as well as on rock properties within the subsurface and the impact on petroleum systems. This study aims to assess the possible effect of ice sheets and associated low temperatures on the past and present-day temperature distribution in the deep subsurface. Two scenarios are presented, considering uncertainties regarding the ice sheet thicknesses. Additionally, the possible influence of glacial loading on the pressure distribution and the effect on petroleum systems is discussed.
Geological background
Tectonics and sedimentation
The subsurface of the northeastern Netherlands is part of the large-scale Central European Basin System (CEBS). Several structural elements with different burial and temperature histories are present in the study area (Fig. 1), including: Lower Saxony Basin (LSB), Friesland Platform (FP), Lauwerszee Trough (LT) and Groningen Platform (GP). While Jurassic sediments are solely preserved in a few Jurassic basins (e.g. LSB), most platforms and troughs are characterised by the absence of Jurassic sedimentary units (Kombrink et al., Reference Kombrink, Doornenbal, Duin, Den Dulk, Ten Veen and Witmans2012). In places, Triassic and Permian strata underwent significant erosion during the Jurassic (Duin et al., Reference Duin, Doornenbal, Rijkers, Verbeek and Wong2006; Kombrink et al., Reference Kombrink, Doornenbal, Duin, Den Dulk, Ten Veen and Witmans2012).
Study area and maximum extent of the main Pleistocene glaciations in Central Europe (Elsterian, Saalian, Weichselian); wells with available pressure data are marked in grey, main cities in purple and wells used in the result section are marked in blue and grey with a black frame. The area of the 3D model is outlined with a solid black line. Main structural elements in the study area are marked with dotted black lines, including LSB: Lower Saxony Basin, FP: Friesland Platform, LT: Lauwerszee Trough, GP: Groningen Platform, CNB: Central Netherlands Basin. Location of study area from Sachse & Littke (Reference Sachse and Littke2018) marked as a purple square in the upper right.
The CEBS is filled by Permian to recent sediments up to a thickness of 10 km (Scheck-Wenderoth & Maystrenko, Reference Scheck-Wenderoth and Maystrenko2013). The basin was strongly influenced in terms of tectonics by its basement, which consists of Palaeozoic sediments deformed during the Caledonian and Variscan orogenies (De Jager, Reference De Jager, Wong, Batjes and de Jager2007). A thick succession of marine, coastal- and fluvial sediments was deposited in the foreland basin of the Variscan mountain belt during the Carboniferous, which acts as the primary source rock interval of the Palaeozoic petroleum system (Maystrenko et al., Reference Maystrenko, Bayer, Brink, Littke, Littke, Bayer, Gajewski and Nelskamp2008). A phase of thermal uplift, deep erosion of Carboniferous sediments and local magmatism followed in the Permian (De Jager, Reference De Jager, Wong, Batjes and de Jager2007). During the Late Permian, evaporitic cycles of the Zechstein group with a decreasing marine influence were deposited in a restricted basin (Strozyk et al., Reference Strozyk, Reuning, Scheck-Wenderoth, Tanner, Soto, Flinch and Tari2017). Continental conditions characterised the basin during Triassic times, with the exception of the Muschelkalk during the Anisian to Ladinian (Geluk, Reference Geluk, Wong, Batjes and de Jager2007; Stollhofen et al., Reference Stollhofen, Bachmann, Barnasch, Bayer, Beutler, Franz, Kästner, Legler, Mutterlose, Radies, Littke, Bayer, Gajewski and Nelskamp2008). Pangea broke up in four extensional phases starting in the Triassic and ending in the Late Jurassic (Hardegsen, Early Kimmerian, Mid-Kimmerian, and Late Kimmerian; Wong, Reference Wong, Wong, Batjes and de Jager2007). As a result of this breakup, graben structures formed while adjacent graben shoulders were uplifted and eroded in Jurassic to Early Cretaceous times (Groetsch et al., Reference Groetsch, Sluijk, Van Ojik, De Keijzer, Graaf, Steenbrink, Grötsch and Gaupp2011). Marine sedimentation of siliciclastics and carbonates prevailed in the Cretaceous and Cenozoic (De Jager, Reference De Jager, Wong, Batjes and de Jager2007). Changing stress fields related to the Alpine orogeny and extension in the Atlantic Ocean introduced inversion and erosion during the Late Cretaceous to Palaeogene (De Jager, Reference De Jager, Wong, Batjes and de Jager2007). In the Quaternary, particularly during the Pleistocene, recurring glacial and interglacial periods and shifts in sea levels and shorelines strongly influenced the depositional environment (De Gans, Reference De Gans, Wong, Batjes and de Jager2007). This led to the deposition of various sediments, including glacial and aeolian sediments of the Peelo, Drachten and Drente formations (Fig. 2; Gunnink et al., Reference Gunnink, Maljers, Van Gessel, Menkovic and Hummelman2013). During glacials, ice sheet advances originating in Northern Europe strongly altered the direction and sedimentation of the Rhine and Meuse rivers and shaped the present-day morphology of the northeastern Netherlands (Laban & van der Meer, Reference Laban and van der Meer2011; Moreau & Huuse, Reference Moreau and Huuse2014; Lang et al., Reference Lang, Lauer and Winsemann2018). Deposition during interglacial periods was mainly marine, fluvial and aeolian (De Gans, Reference De Gans, Wong, Batjes and de Jager2007).
Chronostratigraphic chart of the Neogene according to the general stratigraphic nomenclature of the Netherlands (NAM & RGD, 1980; Van Adrichem Boogaert & Kouwe Reference Van Adrichem Boogaert and Kouwe1994; Gunnink et al. Reference Gunnink, Maljers, Van Gessel, Menkovic and Hummelman2013; TNO-GSN 2020) and Pleistocene ice sheet thicknesses assigned in the modelling (Schokking Reference Schokking1990; Feldmann Reference Feldmann2002; Sachse & Littke Reference Sachse and Littke2018). Negative depths of ice sheet thicknesses represent parts of the ice sheet above the surface, positive depths represent the submerged part of an ice sheet.
Ice history
Several glaciations during the Pleistocene partly covered Central Europe. The Elsterian, Saalian and Weichselian are the most well-known glacial periods in Central Europe; yet, indications of Early Pleistocene glaciations in terms of glaciofluvial sediments, deposited near ice margins, of the Tiglian, Eburonian, Menapian and Cromerian stages have been found in the Dutch offshore area as well (Ehlers, Reference Ehlers1990; De Gans, Reference De Gans, Wong, Batjes and de Jager2007; Laban & van der Meer, Reference Laban and van der Meer2011). Two major glaciations, the Elsterian and the Saalian (Fig. 2), partially covered the area of the northeastern Netherlands (Fig. 2; Ehlers et al., Reference Ehlers, Gibbard, Hughes, Menzies and van der Meer2018; Lang et al., Reference Lang, Lauer and Winsemann2018). Weichselian ice sheets did not advance onto onshore Netherlands (Fig. 2; Laban & van der Meer, Reference Laban and van der Meer2011). The glacial periods such as the Elsterian, Saalian and Weichselian should, however, not be viewed as a single glacial advance but instead, be considered as comprising several pulses of advancing and retreating ice sheets (Fig. 2). Therefore, less extensive glaciations tend to be overprinted by more extensive ice advances, and the small-scale stages of advances and retreats are challenging to determine (Ehlers et al., Reference Ehlers, Gibbard, Hughes, Menzies and van der Meer2018). The extent of ice advances in Central Europe has been mainly constrained by large-scale mapping of till sheets, end-moraine landforms, tunnel valleys and recently with methods such as 3D seismic interpretation, luminescence analysis and cosmogenic radionuclide dating (Ehlers, Reference Ehlers1990; Huuse & Lykke-Andersen, Reference Huuse and Lykke-Andersen2000; Böse et al., Reference Böse, Lüthgens, Lee and Rose2012; Moreau et al., Reference Moreau, Huuse, Janszen, van der Vegt, Gibbard and Moscariello2012; Müther et al., Reference Müther, Back, Reuning, Kukla and Lehmkuhl2012; Lang et al., Reference Lang, Lauer and Winsemann2018).
The Elsterian glacial stage comprises the first significant ice advance to partially cover the onshore Netherlands during the Pleistocene and is discussed to have occurred either from 478 to 424 ka or 374 to 337 ka (Laban & van der Meer, Reference Laban and van der Meer2011; Böse et al., Reference Böse, Lüthgens, Lee and Rose2012). The formation of a large ice-dammed lake caused by advancing Scandinavian and British ice in the North Sea and ice sheets in the Baltic Sea altered the drainage of river systems in Central Europe on a large scale (Ehlers et al., Reference Ehlers, Grube, Stephan and Wansa2011). In the early 2000s, sandy tills were identified that indicate an ice advance during the Elsterian onto the onshore Netherlands, as shown in Fig. 1 (Laban & van der Meer, Reference Laban and van der Meer2011; Meijles, Reference Meijles2015). Up to 500 m, deep glacial valleys with meltwater clays and sands assigned to the Peelo Formation (Fig. 2) are found offshore and onshore, indicating the possible presence and extent of the Elsterian ice advances (De Gans, Reference De Gans, Wong, Batjes and de Jager2007; van der Vegt et al., Reference Van der Vegt, Janszen and Moscariello2012). In Northern Germany, two major ice advances are proposed, interrupted by a 100 km ice margin retreat in a short interstadial; yet, the extent in northwestern Germany and the Netherlands stays ambiguous due to strong erosion during the later Saalian stage (Litt et al., Reference Litt, Behre, Meyer, Stephan and Wansa2007).
The Saalian complex, also referred to as the Saalian stage, advanced much further south, covering the entire study area (Fig. 1) and is assigned to an age of 190–130 ka (Beets et al., Reference Beets, Meijer, Beets, Cleveringa, Laban and Van der Spek2005; Laban & van der Meer, Reference Laban and van der Meer2011; Lang et al., Reference Lang, Lauer and Winsemann2018). Three major Saalian ice advances are indicated by till sheets, namely the Drente I and II substages and the more recent Warthe substage, which did not cover the Netherlands (Litt et al., Reference Litt, Schmincke, Frechen, Schluechter and McCann2008). The first Saalian ice sheet advance is discussed to have occurred during an earlier Saalian glacial advance from 290 to around 241 ka (Lang et al., Reference Lang, Lauer and Winsemann2018). In contrast to Northern Germany, tills and associated glacial sediments of the Saalian complex in the Netherlands are solely assigned to the Drente Formation (Fig. 1; De Gans, Reference De Gans, Wong, Batjes and de Jager2007; Gunnink et al., Reference Gunnink, Maljers, Van Gessel, Menkovic and Hummelman2013). Till sheets in the Netherlands comprise two major ice movement directions, indicating an earlier southern directed advance and a southwestern directed later glacial readvance (Rappol et al., Reference Rappol, Haldorsen, Jørgensen, JM van der Meer and MP Stoltenberg1989; Laban & van der Meer, Reference Laban and van der Meer2011; Böse et al., Reference Böse, Lüthgens, Lee and Rose2012; Lang et al., Reference Lang, Lauer and Winsemann2018). Rappol et al. (Reference Rappol, Haldorsen, Jørgensen, JM van der Meer and MP Stoltenberg1989) proposed three phases of ice advances, which is debated controversially (Laban & van der Meer, Reference Laban and van der Meer2011). In the study area, characteristic structures of this Saalian glaciation are NW-SE oriented glacio-tectonic ridges, push moraine deposits and glacial basins (Laban & van der Meer, Reference Laban and van der Meer2011). In contrast to the Elsterian glacial period, incising glacial valleys related to the Saalian glaciation have only been found in offshore Netherlands (De Gans, Reference De Gans, Wong, Batjes and de Jager2007; Moreau et al., Reference Moreau, Huuse, Janszen, van der Vegt, Gibbard and Moscariello2012).
Methods
Model input
Schlumberger’s PetroMod® 2020 was used to create the 3D basin and petroleum systems model. The theoretical background and general workflow of PetroMod and basin and petroleum systems modelling are described in Hantschel & Kauerauf (Reference Hantschel and Kauerauf2009) and Peters et al. (Reference Peters, Schenk, Scheirer, Wygrala, Hantschel, Hsu and Robinson2017). Two models (scenarios 1 and 2) were used comprising the burial history, mechanical loading of ice sheets, low surface temperatures during the Pleistocene and hydrocarbon generation. To investigate hydrocarbon migration, hydrocarbon generation was calculated with open model borders and using the hybrid migration method of PetroMod, in which domain splitting into low (Darcy flow) and high (Flowpath) permeable layers is used to calculate fluid movement. The Zechstein Group acts as the main seal for Palaeozoic hydrocarbons, while the Rotliegend Group acts as the reservoir in the 3D model when simulating hydrocarbon accumulations on the Groningen Platform. Recently published bulk kinetics for Westphalian coals from Froidl et al. (Reference Froidl, Zieger, Mahlstedt and Littke2020a) and for the Coevorden Formation from Froidl et al. (Reference Froidl, Littke, Baniasad, Zheng, Röth, Böcker, Hartkopf-Fröder and Strauss2020b) were used for the petroleum generation of Palaeozoic and Mesozoic source rocks in the area (Amberg et al., Reference Amberg, Back, Sachse and Littke2022).
The recently published 3D basin model studying the burial and thermal history of the northeastern Netherlands from Carboniferous times on by Amberg et al. (Reference Amberg, Back, Sachse and Littke2022) was used and extended to incorporate ice loading and unloading during Pleistocene glacial stages (Fig. 3). The 3D basin model has a high spatial resolution of 250 × 250 m and extends 150 km from north to south and 100 km from west to east. The Upper North Sea Group, comprising Holocene to Miocene sediments, was refined to implement ice loads during glacial periods. Depth maps of BRO DGM 2.2 by TNO (TNO-GSN, 2019), a shallow subsurface model of the Netherlands down to a depth of 500 m, were used to further subdivide the Upper North Sea Group (NU) into four sublayers (Table 1; Holocene – Eemian, Saalian, Elsterian, Rest Pleistocene and Rest NU) to incorporate the glaciations in a valid depositional and time framework (Schokker et al., Reference Schokker, Weerts, Westerhoff, Berendsen and Otter2007; Gunnink et al., Reference Gunnink, Maljers, Van Gessel, Menkovic and Hummelman2013). Thickness maps of four key Quaternary intervals derived from the BRO DGM 2.2 model are shown in Fig. 4 (Gunnink et al., Reference Gunnink, Maljers, Van Gessel, Menkovic and Hummelman2013). Lithologies for Quaternary and Neogene units in PetroMod are set according to Gunnink et al. (Reference Gunnink, Maljers, Van Gessel, Menkovic and Hummelman2013). Lithologies of all units used in this 3D model and their initial geomechanical properties are listed in Table 2. Four additional layers, representing two major Elsterian and two major Saalian glacial advances during the Pleistocene (Fig. 2), were incorporated in the stratigraphic framework described above. The extent of ice sheets was set according to De Gans (Reference De Gans, Wong, Batjes and de Jager2007), Ehlers et al. (Reference Ehlers, Gibbard, Hughes, Menzies and van der Meer2018) and Sachse & Littke (Reference Sachse and Littke2018). Assumptions regarding the ice sheet thickness are generally more uncertain. The minimum thickness of the Drenthe ice advance during the Saalian was estimated to be around 200 m by Schokking (Reference Schokking1990) near Marum, northern Netherlands. Feldmann (Reference Feldmann2002) proposed a Saalian ice sheet thickness of 100–600 m in the study area, but thickness gradients near the ice margin vary strongly (Mrugalla, Reference Mrugalla2020). For the Elsterian ice advances, a smaller ice sheet thickness is proposed due to its closer location to the ice margin (Fig. 2). To account for the general uncertainty regarding ice sheet thicknesses, two scenarios were used (Fig. 2), one with the minimum thickness estimated by Schokking (Reference Schokking1990) and one with a maximum ice sheet thickness of 600 m based on the estimation by Feldmann (Reference Feldmann2002). For the implementation of ice loading, a workaround was constructed in which pseudo-erosion maps represent periods of ice loading. The ice sheet advances of one glacial stage are separately implemented in the 3D basin model (Fig. 2; Table 1) and differ in their ice sheet thickness and associated boundary conditions. Supplementary pressure data from various stratigraphic intervals of the Germanic Triassic Group down to the Carboniferous Limburg Group were accessed from TNO’s NLOG database (www.nlog.nl) to study the effect of ice loading on pressure trends across the study area.
Present-day geometry and extent of the 3D basin and petroleum system model with main sedimentary units down to the crystalline basement and imposed ice sheet during the second Saalian glacial advance. Subdivison of the sedimentary units according to the general stratigraphic nomenclature of the Netherlands (NAM & RGD 1980; TNO-GSN 2020).
Thickness of key Quaternary sedimentary unit intervals compiled from the BRO-DGM shallow subsurface model (TNO-GSN 2019).
Input data of the Upper North Sea Group used in the modelling process. Glacial Weichselian and interglacial Eemian deposits are summarised in the Holocene & Post Saalian sedimentary layer.
Lithologies and geomechanical properties of sedimentary layers used in the 3D model.
Boundary conditions
The upper thermal boundary conditions, the paleowater depth and the sediment-water interface temperature (SWIT) were adjusted from the original model for Quaternary times (Amberg et al., Reference Amberg, Back, Sachse and Littke2022). Paleobathymetric data, indicating the submerged amount of an ice sheet below sea level, was taken from Sachse & Littke (Reference Sachse and Littke2018). The SWIT is calculated from the paleolatitude- and paleowater depth history (Wygrala, Reference Wygrala1989) and was strongly influenced by the short term changes in surface and ground temperatures and the presence of ice sheets. Grassmann et al. (Reference Grassmann, Cramer, Delisle, Hantschel, Messner and Winsemann2010) and Sachse & Littke (Reference Sachse and Littke2018) used mean annual ground temperatures based on a temperature model derived from marine proxy data by Delisle et al. (Reference Delisle, Grassmann, Cramer, Messner, Winsemann, Hambrey, Christoffersen, Glasser and Hubbard2007) in numerical modelling studies in Northern Germany. Due to the proximity to the study areas above, ground temperature data ranging from −6 to 10°C were also used in this study for interglacial times and ice-free areas during the Elsterian stages as one part of the upper thermal boundary conditions. Temperatures during the Eemian interglacial (∼120 ka) were up to 2°C higher than present-day temperatures in the study area (Frenzel et al., Reference Frenzel, Pécsi and Velichko2001). Paleo-surface temperatures in the Netherlands ranging from −7 to 12°C from Luijendijk et al. (Reference Luijendijk, Ter Voorde, Van Balen, Verweij and Simmelink2011) were used for the Eemian, Weichselian and Holocene times. At ground temperatures below 0°C, the pore water in the shallowest subsurface is frozen and is present as permafrost impacting the migration of fluids due to reduced rock porosities and permeabilities and the distribution of heat in a sedimentary basin (Grassmann et al., Reference Grassmann, Cramer, Delisle, Hantschel, Messner and Winsemann2010). Permafrost was calculated using the built-in simulator in PetroMod, in which permafrost lithologies replace lithologies with ice-filled pores below a temperature of 0.7°C (Hantschel & Kauerauf, Reference Hantschel and Kauerauf2009). Even though ice sheets during the most recent glacial Weichselian stage did not advance to the onshore Netherlands, the presence of permafrost, as a consequence of low ground temperatures and changes in sea level, did have an influence on the subsurface and sedimentation patterns within the study area (Laban & Van der Meer, Reference Laban and van der Meer2011). During the Weichselian stage, the southernmost extent of permafrost was at a latitude of around 50°N, therefore comprising the study area, with a permafrost thickness of around 200 m and an increased thickness towards northern latitudes (Frenzel et al., Reference Frenzel, Pécsi and Velichko2001). There are permafrost indications such as ice-wedge casts and periglacial features in western Germany and the Netherlands during the Elsterian and Saalian stages (Vandenberghe, Reference Vandenberghe, Paepe and Melnikov Vladimir2001). Ice sheets and glaciers can differ in their thermal structures (e.g. warm- or cold-based). In this study, a warm-based ice sheet with water and temperatures of around 0°C at the sediment to ice contact is assumed that is most common in low latitude glaciers. The presence of meltwater at the base of the glaciers is expected, but meltwater flow at the base of the ice sheet is not modelled with the software used. In comparison to the base of a temperate ice sheet, temperatures at the surface of a temperate ice sheet are significantly lower and are set to values ranging from −5 to −20°C. This results in an ice sheet base temperature of around 0°C. In areas not covered by ice sheets during Elsterian times, a SWIT temperature of −2°C was set. Paleowater depth and SWIT maps with the properties mentioned above were created and assigned to the different glacial and interglacial stages. The evolution of the SWIT is illustrated in Fig. 5. The ice lithology was set to be impermeable, without porosity, and is characterised by thermal conductivities of 2.39 W/mK at −20°C and 2.22 W/mK at 0°C and a density of 0.917 g/m3 (Sachse & Littke, Reference Sachse and Littke2018). The basal heat flow, independent of the low temperatures of ice sheets, was not changed from the original model.
SWIT used in the basin model based on Lujiendijk et al. (Reference Luijendijk, Ter Voorde, Van Balen, Verweij and Simmelink2011), Grassmann et al. (Reference Grassmann, Cramer, Delisle, Hantschel, Messner and Winsemann2010) and Sachse & Littke (Reference Sachse and Littke2018). Weichselian, Saalian and Elsterian glacial periods are indicated with different colors.
Geomechanical modelling and fluid flow in sedimentary systems
External loading, such as sedimentation or glacial loading, acting on the subsurface can severely influence sedimentary rocks and the fluids within. In basin and petroleum system applications, the external loading of rocks, expressed as rock stress, is modelled using the Terzaghi rock stress model, assuming only vertical stress as the maximum principal stress identical to the lithostatic stress/pressure (Hantschel & Kauerauf, Reference Hantschel and Kauerauf2009). The extension of the Terzaghi rock stress modelling to a 3D-poroelasticity model by Biot (Reference Biot1941) further enables the calculation and prediction of geomechanical processes such as rock failure (Peters et al., Reference Peters, Schenk, Scheirer, Wygrala, Hantschel, Hsu and Robinson2017). The compaction of rocks refers to a reduction of the sediment bulk volume and depends on the outflow of fluids and the capacity of a rock to compact, the compressibility (Broichhausen et al., Reference Broichhausen, Littke and Hantschel2005). Water is assumed to be incompressible, and for the Terzaghi equation, it is assumed that the volume of solids does not change (Neuzil, Reference Neuzil2003; Verweij, Reference Verweij2003). The mechanical compaction reduces porosity nearly irreversible (Houseknecht, Reference Houseknecht1987). A lower compressibility during decreasing effective stress, caused by e.g. erosion or uplift, results in a certain degree of porosity preservation and lower than initial porosities (Hantschel & Kauerauf, Reference Hantschel and Kauerauf2009). Rock permeability describes the ability of fluids to pass through a porous material and governs pore pressure fields by influencing fluid pathways and fluid flow rates. Overpressure is generated by several processes, such as rapidly increased vertical stress or an enhanced fluid volume due to gas generation, clay dehydration or fluid movement (Verweij, Reference Verweij2003). In this study, pore pressure and overpressure refer to water pressure not taking into account increased fluid pressures from hydrocarbon generation. The most important process of overpressure generation is compaction disequilibrium caused by overburden load and incomplete sediment compaction (Broichhausen et al., Reference Broichhausen, Littke and Hantschel2005; Aplin & Macquaker, Reference Aplin and Macquaker2011). An under compaction is caused by incompressible pore water in fine-grained sediments inhibited from escaping due to very low permeabilities (Rodrigues Duran et al., Reference Rodrigues Duran, di Primio, Anka, Stoddart and Horsfield2013). Other possible causes of overpressure generation are fluid flow, sea-level changes, changing groundwater volumes and salinities, oil and gas generation or cementation (Verweij, Reference Verweij2003). Erosion or the withdrawal of an ice sheet results in the slow decrease of overpressure mainly due to the reduction of vertical stress coupled with the preservation of porosity and smaller compressibility. In impermeable facies, such as salt, the pressure is equal to the lithostatic pressure resulting in a strong pressure difference at the interface of a permeable to an impermeable layer. In-detail insight into how geological processes are calculated with traditional basin modelling software such as PetroMod is found in Hantschel & Kauerauf (Reference Hantschel and Kauerauf2009).
Heat transfer takes into account heat conduction, convection and radiogenic heat production (Peters et al., Reference Peters, Schenk, Scheirer, Wygrala, Hantschel, Hsu and Robinson2017). The implementation of permafrost is a purely conductive approach but takes into account latent heat (Grassmann et al., Reference Grassmann, Cramer, Delisle, Hantschel, Messner and Winsemann2010). When pore water is replaced by ice (permafrost), the thermal conductivity and permeability of the respective layer are changed. Ice sheets are assumed to consist of pure ice with its typical thermal and mechanical properties (Table 1). A warm glacier base is assumed with ice at its pressure melting point. Meltwater flow underneath the ice sheet is not considered (Rodrigues Duran et al., Reference Rodrigues Duran, di Primio, Anka, Stoddart and Horsfield2013).
Modelling results
Influence on temperature
By introducing Pleistocene glaciations, temperatures in the subsurface are lowered during and after the glacial stages. Calculated effects on the temperature distribution and history are approximately similar for both ice thickness scenarios due to using the same upper boundary conditions and assuming a temperature of around 0°C at the base of the ice sheet (SWIT); therefore only results of the ice thickness scenario two are shown. Figure 6 depicts the Pleistocene to recent temperature history in scenario two of three sedimentary units in well USQ-01 in the northeast, influenced by both the Elsterian and Saalian glaciations, and in well CLD-01 in the south, where solely the Saalian glaciation was present. In general, temperatures show the most substantial decrease in the shallow subsurface and this effect declines with further depth. In areas that experienced ice advances during both glacial stages, the lowest temperatures are observed before the Elsterian I and Saalian I glacial advances (Fig. 6a). In places where the Elsterian glaciations are absent, low temperatures are calculated throughout the Elsterian stage due to permafrost and low ground temperatures and at the beginning of the first Saalian glaciation (Fig. 6b). During maximum ice sheet coverage, temperatures are not decreased but elevated (Fig. 6a, b). This effect weakens with smaller ice thicknesses used in scenario one. This results from two processes counteracting each other: firstly, increasing temperatures are caused by a subsidence of the underlying sedimentological succession in response to the load of the overlying ice sheet. This effect is strengthened with an increasing ice sheet thickness. Secondly, low ice sheet base temperatures propagate down into the subsurface. This interplay limits times of minimum temperatures to pre-and post-maximum ice sheet thicknesses. During the Eemian interglacial (around 126 ka), temperatures increased to a pre-glacial state (Figs. 6 and 7). Even though ice masses during the Weichselian glacial stage did not advance to the onshore Netherlands, low associated ground temperatures influence the temperature distribution in the subsurface. Varying degrees of permafrost are present from Elsterian glacial times until the Holocene. A maximum permafrost thickness of 180 m is calculated during the Last Glacial Maximum (LGM) at 20 ka. Permafrost thicknesses of up to 150 m are also calculated prior and in between glacial advances during the Saalian and Elsterian stages. Additionally, the highly variable thickness of Zechstein salt strongly influences the temperature distribution in the subsurface. Above salt domes, temperatures in uppermost sedimentary layers are increased up to 10°C due to the high conductivity of salt impeding and slowing down the permafrost formation. Throughout the study area, low temperatures in the shallow subsurface prevail until the present day (Fig. 8), indicating that no thermal equilibrium has been reached in the shallow subsurface. Temperature through sedimentary formations in the shallow subsurface is lowered down to a depth of around 3000–4000 m resulting in a difference of around 3°C at a depth of 3 km (Fig. 8). At greater depths, the influence of low temperatures from the top becomes insignificant.
(a) Influence of glacial loading on temperatures in Lower Pleistocene sediments (depth: ~50 m), the Middle North Sea Group (depth: ~500 m), the Rijnland Group (depth: ~2100 m) from the Pleistocene to the present day at well location USQ-01 in scenario two; (b) Influence of glacial loading on temperatures in Lower Pleistocene sediments (depth: ~40 m), the Middle North Sea Group (depth: ~200 m), the Rijnland Group (depth: ~2000 m) from the Pleistocene to the present day at well location CLD-01 in scenario two; Locations of wells USQ-01 and CLD-01 are marked in Fig. 1.
Temperature maps of the the Lower North Sea Group (NL) before, during and after ice loading at six points in time (380 ka; 360 ka; 190 ka; 175 ka; 170 ka; 130 ka; 20 ka; present day) using scenario two. The unit is buried to a depth of up to 1000 meters in the north and down to 200 meters in the south. The dotted black line represents the maximum extent of Elsterian glacial advances.
Present-day temperatures (solid lines) and local geothermal gradients (dotted lines) illustrating the difference of implemented Pleistocene glaciations to a no ice scenario on the subsurface in the vicinity of the USQ-01 well location on the Groningen Platform and the NSL-01 well on the Friesland Platform; Blue lines show temperatures and geothermal gradient with implemented glaciations, red lines show the no ice scenario.
Compared to the non-ice scenario, a low local geothermal gradient is present in the uppermost part of the subsurface in the ice scenario (Fig. 8). Based on the dominance of sandy and coarse-grained lithologies in the Miocene to Holocene Upper North Sea Group, this relates to a geothermal gradient of around 15°C/km for the Upper North Sea Group to a depth of around 460 m in the vicinity of the USQ-01 well (Fig. 1). This trend is observable across the study area and seems to be mainly related to low ground temperatures during the Weichselian (Fig. 8).
Influence on pressure
Glacial loading has the potential to alter pressure conditions in the subsurface on a large scale. Lithostatic, hydrostatic and pore pressures build up during glacials by extra loading and then decrease to a pre-loading state with time, depending on the strength of the loading before (Fig. 9). In areas with both Elsterian and Saalian ice coverage, overpressure is observed during both glacial periods (Figs 9 and 10). Pressures build up during the Elsterian, and in deeper buried formations, pressures do not equilibrate fast enough to a pre-glacial state and further increase during the Saalian glaciation. Pore pressures during the Elsterian stage do not increase as strongly as during the Saalian stage due to the compaction of the shallow buried Middle North Sea Group (NM) at the well location USQ-01 during the Elsterian stage (Fig. 9). In the southern part of the study area, where ice loading in the Elsterian stage is absent (Fig. 1), pressure changes are only observed from the Saalian stage (Fig. 10). A further distinction must be made regarding the depths and the dominant lithology of the observed formations. Within relatively shallow formations of the North Sea Group in areas with both glaciations, pore pressures decrease strongly after the retreat of an ice mass and go back to the initial values relatively fast (Fig. 9). It takes longer to equilibrate back to the initial state in the sub-salt domain and Cretaceous to Triassic formations because of the burial depth and the dominance of fine-grained lithologies. Maximum overpressures, induced by mechanical ice sheet loading, of up to 5 and 11 MPa are observed in scenarios one and two in the LSB and the northern part of the study area in the Germanic Trias Group (Fig. 10). In the deeply buried sub-salt Carboniferous Limburg Group, maximum overpressures of up to 14 and 22 MPa in the north, decreasing to 10 and 14 MPa in the south are calculated for scenarios one and two. In scenario one, the Cretaceous Rijnland group, the Niedersachsen Group and shallower sedimentary units are not or only slightly overpressured at the present day. In scenario two, overpressures of up to 7 MPa are calculated in the Rijnland to Altena Groups at the present day. In the Coevorden Formation that is part of the Niedersachsen Group, overpressures of up to 4 MPa in scenario one and 10 MPa in scenario two are calculated during the Saalian glacial stage, but pressure decreased strongly after the retreat of the ice sheets. In the southern part, overpressures are only observed in the sub-salt interval due to the missing Elsterian glacial coverage, an absence of Mesozoic formations and the resulting shallow depth of these sedimentary units. Geomechanical properties of selected layers discussed above are listed in Table 3. Prior to glacial loading, no or just slight overpressures of up to 1 MPa are observed in the Niedersachsen, Altena, Germanic Trias and the Limburg Groups. Overpressures of around 10 MPa are calculated in the Carboniferous Limestone Group. A strong influence of salt walls and diapirs is observed, slowing down the pressure decrease due to the very low permeability of salt. This is attributed to the salts rock properties and its missing porosity and permeability, highlighting the lithologies importance in such a model.
Computed hydrostatic, lithostatic and pore pressures for the no ice model (solid lines) and ice model (scenario one, dotted lines) for the shallow Middle North Sea Group at the well location USQ-01 on the Groningen Platform. Times of overpressure are marked with red line patterns. E1= Elsterian glacial advance 1; E2 = Elsterian glacial advance 2; S1 = Saalian glacial advance 1; S2 = Saalian glacial advance 2.
Effect of glacial loading on overpressures of the deeply buried Lower Germanic Trias Group (RB) in scenario one. The formation shallows to the south and is absent in parts of the study area. The dotted black line represents the maximum extent of Elsterian glacial advances.
Geomechanical properties (porosity, permeability and compressibility) of selected sedimentary layers prior to glacial loading at 400 ka.
Influence on porosity
Porosities in the shallow subsurface are decreased during the first Elsterian ice advance and the first Saalian ice advance in the southern area. The Upper North Sea Group consists of mainly sandy sediments where pore space is permanently reduced. This permanent reduction due to mechanical compaction is observable in both scenarios down to a depth of around 500 m by up to 4% depending on the used lithology and its compressibility. A porosity reduction in deeper intervals and sub-salt layers is not observable.
Influence on source rocks and hydrocarbon migration
Elsterian and Saalian glaciations in the Pleistocene do not severely influence the temperature or the thermal maturity of sub-salt source rock intervals (Fig. 11). The temperature of supra-salt source rocks such as the Posidonia and Coevorden Formation in the Jurassic basin is effectively lowered, but an influence on maturities is not calculated (Fig. 11).
Calculated burial depths, temperatures and maturities of the sub-salt Caumer Subgroup (DCC) and the supra-salt Altena Group (AT) and Niedersachsen Group (SK) at the well location EMM-07 in the Lower Saxony Basin in scenario one.
One of the most significant hydrocarbon accumulations of Europe, the Groningen gas field, is located in the study area. The hydrocarbon accumulations on the Groningen Platform were simulated, and the effects of ice loading and unloading were studied. A small volume change in hydrocarbons is calculated at the beginning of ice loading and unloading in scenario two. During loading, hydrocarbon volume decreases and after unloading, it increases reaching the original volume. Changes in hydrocarbon mass seem to be restricted to horizontal and vertical redistribution within the Rotliegend reservoir related to changes in volume due to compression and extension of the hydrocarbons. The migration of hydrocarbons sourced from the Berriasian Coevorden Formation in the LSB during the Pleistocene was simulated using a kinetic from Froidl et al. (Reference Froidl, Littke, Baniasad, Zheng, Röth, Böcker, Hartkopf-Fröder and Strauss2020b). Fig. 12 shows a top loss of hydrocarbons through the overburden in scenario two. Using scenario one, this effect is lowered. In PetroMod, the top loss or top outflow, is the amount of hydrocarbons that reached the upper boundary of the model, the surface-water-interface. This loss coincides with the retreat of the first Elsterian advance at around 370 ka and the decrease of lithostatic pressure. Top losses also occur during the retreat of the second Elsterian and both Saalian ice sheets (Fig. 12). This strongly indicates the influence of glacial loading and unloading on the migration pathways of hydrocarbons in the subsurface.
Top loss of calculated hydrocarbons of the Wealden Shale to the surface during the Pleistocene and a general burial history of the Niedersachsen Group in the LSB indicating times of glacial loading and unloading in scenario two. Top loss masses are approximations due to kinetics used and simplifications of the 3D basin model.
Discussion
Temperature
Low present-day geothermal gradients of around 20°C/km are found in the subsurface across the Netherlands down to a depth of 500 m, while a geothermal gradient of around 33°C/km is present in deeper intervals (Bonté et al., Reference Bonté, Van Wees and Verweij2012; Ter Voorde et al., Reference Ter Voorde, Van Balen, Luijendijk and Kooi2014). Luijendijk et al. (Reference Luijendijk, Ter Voorde, Van Balen, Verweij and Simmelink2011) proposed an effect of low paleosurface temperatures and associated groundwater flow in permeable layers during the Weichselian glacial stage on the present-day temperature distribution in the subsurface. Recent studies by Ter Voorde et al. (Reference Ter Voorde, Van Balen, Luijendijk and Kooi2014) and Gies et al. (Reference Gies, Struijk, Békési, Veldkamp and van Wees2021) only implemented low ground temperatures during the Weichselian stage, when ice sheets did not advance onto the area of the Netherlands. Results of our model show that temperatures recovered after the Elsterian and Saalian glacial stages during the warm Eemian interglacial stage to a pre-glacial state (Figs. 6 and 8). During the Eemian interglacial at around 125 ka, temperatures in Central Europe were about 0–2°C higher than at the present day (Frenzel et al., Reference Frenzel, Pécsi and Velichko2001; Luijendijk et al., Reference Luijendijk, Ter Voorde, Van Balen, Verweij and Simmelink2011). The study results support the hypothesis that mainly low ground temperatures related to the Weichselian glacial stage and not to the Elsterian and Saalian stages are the driving factor of the present-day thermal disequilibrium in the shallow subsurface of the Netherlands. In contrast to the transient state of the geothermal gradient in the shallow subsurface, the deep subsurface does not seem to be severely influenced by Pleistocene glacial stages. Luijendijk et al. (Reference Luijendijk, Ter Voorde, Van Balen, Verweij and Simmelink2011), Ter Voorde et al. (Reference Ter Voorde, Van Balen, Luijendijk and Kooi2014) and Gies et al. (Reference Gies, Struijk, Békési, Veldkamp and van Wees2021) calculated a temperature difference of up to 3°C down to a depth of up to 3-4 km using low Weichselian surface temperatures. In a basin modelling study of the Mittelplate oilfield in Northern Germany, a temperature difference of 4°C was calculated in the Mid-Jurassic reservoir at a depth of around 3 km from pre-Pleistocene times to the present day (Grassmann et al., Reference Grassmann, Cramer, Delisle, Hantschel, Messner and Winsemann2010). This is in accordance with our results showing diminishing influence of temperatures at intervals below 4 km and the steady state temperature of the sub-salt sedimentary layers. Locally, other factors such as hydrothermal fluid flow can strongly alter the temperature field even at greater depths (Gies et al., Reference Gies, Struijk, Békési, Veldkamp and van Wees2021).
Permafrost is formed in the shallow subsurface due to low ground temperatures. Frozen pore water increases the thermal conductivity of rocks significantly (Lachenbruch et al., Reference Lachenbruch, Sass, Marshall and Moses1982), influencing the temperature distribution, particularly in porous sediments (Fjeldskaar & Amantov, Reference Fjeldskaar and Amantov2018). During freezing, latent heat is liberated, while melting leads to an uptake of energy (Delisle, Reference Delisle1998). This mechanism mitigates the effect of temperature variations caused by low ground temperatures during glacials and elevated temperatures during interglacials in the subsurface (Ter Voorde et al., Reference Ter Voorde, Van Balen, Luijendijk and Kooi2014). Frenzel et al. (Reference Frenzel, Pécsi and Velichko2001) proposed a permafrost thickness of 0–200 m in the study area with permafrost temperatures ranging from −3°C to −5°C during the Last Glacial Maximum in the Weichselian glacial stage. In a more recent study by Govaerts et al. (Reference Govaerts, Beerten and Ten Veen2015), a permafrost depth of 120–180 m was calculated for the Last Glacial Maximum in the Weichselian. This range is in accordance with permafrost thicknesses calculated with the model presented. Implementing permafrost with the method in PetroMod results in a maximum difference of 0.1–0.4°C on calculated present-day temperatures in different depths in the subsurface. Gies et al. (Reference Gies, Struijk, Békési, Veldkamp and van Wees2021) used a different method to incorporate permafrost but concluded that the effect of an implementation of permafrost and latent heat only has a minor influence on present-day temperatures in the Netherlands and is of much higher importance in areas of more recent permafrost presence.
Salt effectively transports heat in the subsurface due to its high thermal conductivity. Increased temperatures above salt structures are found across onshore and offshore Netherlands (e.g. Bonté et al., Reference Bonté, Van Wees and Verweij2012; Gies et al., Reference Gies, Struijk, Békési, Veldkamp and van Wees2021) and the CEBS and are referred to as salt chimneys (Rodon & Littke, Reference Rodon and Littke2005; Magri et al., Reference Magri, Littke, Rodon, Bayer, Urai, Littke, Bayer, Gajewski and Nelskamp2008). This increased heat transport from the basement can impede the downward propagation of permafrost and reduce the maximum depth of permafrost by up to 100 m above salt domes compared to surrounding areas with more minor salt thicknesses (Grassmann et al., Reference Grassmann, Cramer, Delisle, Hantschel, Messner and Winsemann2010).
Ice sheet thicknesses
The ice sheet’s thickness, which remains a parameter of uncertainty, and the duration of the ice sheet coverage, strongly control the pressure distribution and evolution in the subsurface. In Northern Germany, a greater ice sheet thickness was used in basin modelling studies by Grassmann et al. (Reference Grassmann, Cramer, Delisle, Hantschel, Messner and Winsemann2010) and Sachse & Littke (Reference Sachse and Littke2018). A lower ice sheet thickness, as stated by Schokking (Reference Schokking1990) and Feldmann (Reference Feldmann2002), is assumed in the study area attributed to the more southern location and the closer distance to the ice sheet margin. Winguth et al. (Reference Winguth, Mickelson, Larsen, Darter, Moeller and Stalsberg2005) used ice sheet thickness gradients of around 1 km thickness per 100–150 km distance along the flowline of the ice sheet during the Late Weichselian, depending on the hard-rock structure at its base in western Norway. Comparing both scenarios with ice thicknesses used in Northern Germany (Sachse & Littke, Reference Sachse and Littke2018) with available pore pressure data from the study area, the impact of different ice thicknesses is visible (Fig. 13). Ice thicknesses from Northern Germany result in an excess pore pressure at all locations at the present day, in some wells near the lithostatic pressure. Scenario two results in pore pressures close to measured pore pressure data, while scenario one and the no-ice scenario do not greatly increase the pore pressures (Fig. 13). By neglecting other possible factors of overpressure generation and assuming a sole influence of glacial loading on the pressure distribution of the subsurface, both scenarios one and two could be reasonable scenarios for the ice thicknesses during the Pleistocene. Ice sheet thicknesses greater than scenario two would result in excess pore pressures at the present day (Fig. 13) and make such scenarios unlikely. This is in agreement with the estimation of Schokking (Reference Schokking1990) and Feldmann (Reference Feldmann2002).
Comparison of measured fluid pressure data with computed hydrostatic, lithostatic and pore pressures using a non-ice scenario, the scenarios depicted in Fig. 2 and one scenario with ice thicknesses from Northern Germany on different well locations (WRM-02, VLW-02, N07-01). The location of wells is shown in Fig. 1.
Pressure and porosity
In the presented 3D model with implemented mechanical ice loading, high overpressures are calculated in the Germanic Trias and Carboniferous Limburg Group, while the Rijnland Group, the Altena Group and the Niedersachsen Group are only slightly overpressured. High subsurface pressures indicated from fluid and leak-off pressure data in the sedimentary groups listed above are present in the northeastern Netherlands with increasing pressures to the northwest (Verweij et al., Reference Verweij, Simmelink, Underschultz and Witmans2012b). Even though simplifications regarding the facies distribution are made in such a large-scale model, calculated pressures show trends similar to those as observed based on measured pore pressures in the area. Salt in the Upper Germanic Trias Group and the Zechstein Group, in which a simplified pure salt lithology with only 1% porosity and no permeability is used, acts as a permeability barrier (Verweij et al., Reference Verweij, Simmelink, Underschultz and Witmans2012b) slowing down the pressure dissipation in the layers below after the retreat of ice shields. In the northern part of the study area, overpressures are also related to the relatively rapid sediment accumulation since Cenozoic times (Verweij et al., Reference Verweij, Simmelink, Underschultz and Witmans2012b). Calculated pressures in the southern part of the study area cannot be compared to measured data due to lacking subsurface pressure data. In deeply buried sedimentary units such as the Germanic Trias Group, pressures did not equilibrate back to the pressure conditions prior to ice loading (Fig. 10). The primary mechanism of abnormal pressures with implemented ice advances is a compaction disequilibrium caused by a fast reduction in porosity and permeability from high overburden and the impeded escape from pore water due to permeability barriers (Broichhausen et al., Reference Broichhausen, Littke and Hantschel2005; Rodrigues Duran et al., Reference Rodrigues Duran, di Primio, Anka, Stoddart and Horsfield2013). Therefore, advances of ice sheets during the Pleistocene possibly contributed to the overpressures in rock volumes that are sensitive for overpressure generation. An influence of gas accumulations in Rotliegend and Triassic reservoirs on subsurface pore pressures is expected but not considered in this study. Hydrostatic pressures are increased during glacial loading (Fig. 9) by the increased overburden and glacial meltwater in a warm-based glacier setting (van Weert et al., Reference Van Weert, Van Gijssel, Leijnse and Boulton1997; Sachse & Littke, Reference Sachse and Littke2018). Melt- and groundwater flow, not included in this modelling study, and the presence of aquifers in the subsurface connected to the ice sheet base can be essential factors in the distribution of subsurface pressures and the dissipation of overpressures (Boulton and Caban, Reference Boulton and Caban1995; Person et al., Reference Person, McIntosh, Bense and Remenda2007; Rodrigues Duran et al., Reference Rodrigues Duran, di Primio, Anka, Stoddart and Horsfield2013). Additionally, high pressures from vertical ice loading can result in fault reactivation (Brandes et al., Reference Brandes, Polom and Winsemann2011) as well as salt movement (Lang et al., Reference Lang, Hampel, Brandes and Winsemann2014) and affect Pleistocene and Holocene supra-salt sediments as observed in the Eastern Glückstadt Graben in Northern Germany (Huster et al., Reference Huster, Hübscher and Seidel2020). A similar setting with thick salt sequences piercing to a shallow depth is present in the Netherlands. Therefore, reactivation of faults and induced salt movement is another possible scenario in the study area.
Sedimentary layers in the upper 500 m of the subsurface show a decrease of porosity of up to 4% caused by advancing ice sheets. Porosities are decreased more strongly in scenario two due to the higher overburden. A comparison to porosity data, available in reservoir intervals of the Rijnland, Triassic and Rotliegend Group, is not possible due to the vertical resolution of the model and its layers.
Source rocks and hydrocarbon migration
Low surface temperatures and sequential ice sheet loading and unloading during the Pleistocene can have effects on different parts of a petroleum system, such as lowering source rock and reservoir temperatures (Fig. 11) as well as leakage of hydrocarbons from the reservoir (Kjemperud & Fjeldskaar, Reference Kjemperud, Fjeldskaar, Larsen, Brekke, Larsen and Talleraas1992; Rodrigues Duran et al., Reference Rodrigues Duran, di Primio, Anka, Stoddart and Horsfield2013; Daszinnies et al., Reference Daszinnies, Plaza-Faverola, Sylta, Bünz, Mattingsdal, Tømmerås and Knies2021). Lower subsurface temperatures during ice loading and the reduction of reservoir temperatures to the present day (e.G. Mittelplate oil field; Grassmann et al., Reference Grassmann, Cramer, Delisle, Hantschel, Messner and Winsemann2010) could potentially negatively affect hydrocarbon quality by allowing microbial growth and facilitating biodegradation when reservoir temperatures are below 60–80°C and when, additionally, hydraulic connectivity to shallower layers exists (Peters et al., Reference Peters, Walters and Moldowan2005). Even though temperatures below 80°C prevailed during ice loading and unloading in the Mittelplate oil field, no signs of biodegradation are visible which is attributed to the short duration of the low temperature events (Sachse & Littke, Reference Sachse and Littke2018). In the Palaeozoic petroleum system of northeastern Netherlands, an effect on hydrocarbon quality is not probable because reservoir temperatures stayed above 80°C during glacial stages. In the LSB, however, Mesozoic reservoir intervals potentially sourced by the Coevorden Formation (Amberg et al., Reference Amberg, Back, Sachse and Littke2022) encountered lower temperatures that might impact the oil quality of fields in the area. The shallow (800 m depth) Schoonebeek oil field in the LSB contains indeed oil of relatively low API gravity of 25°. However, this fact is assumed to result from early oil expulsion out of the source rock (De Jager & Geluk, Reference De Jager, Geluk, Wong, Batjes and de Jager2007) and not to biodegradation based on unpublished geochemical data.
A leakage of hydrocarbons through the top seal due to capillary seal failure following the expansion of gas volumes during ice retreat was modelled in the Hammerfest Basin in the Norwegian Barents Sea (Rodrigues Duran et al., Reference Rodrigues Duran, di Primio, Anka, Stoddart and Horsfield2013). Ice loading during phases of maximum ice sheet thickness could induce overpressures and finally fracture failure in brittle top seals and thus a hydrocarbon loss from reservoirs. In the study area, natural gas accumulations in the Rotliegend reservoirs are sealed by the Zechstein Group containing ductile salt. Due to this fact, fracture failure in the top seal is unlikely. Additionally, it should be noted that the software treats salt as a perfect seal, so no matter how thin the salt is, no migration through the salt is possible. In reality, however, some gas loss through the Zechstein seal has occurred, very probably in areas, where Zechstein salt is thin or missing; otherwise much more gas would be present in the Rotliegend and Carboniferous (Uffmann & Littke, Reference Uffmann and Littke2011). In Northern Germany, a decrease of the gas–oil ratio was calculated after ice loads, which can be expected due to the higher solubility of gas in oil at higher pressure (Sachse & Littke, Reference Sachse and Littke2018). Hydrocarbon accumulations in the Palaeozoic Rotliegend reservoirs in the northeastern Netherlands contain only natural gas (De Jager & Geluk, Reference De Jager, Geluk, Wong, Batjes and de Jager2007); therefore, an additional solution in oil is not possible.
High pore pressures and overpressures are calculated in some supra-salt sedimentary formations such as the Germanic Trias Group at the present day due to compaction disequilibrium and the inhibition of an escape of fluids. Overpressures have been calculated in the Cretaceous Coevorden Formation during the entire Elsterian and Saalian glacial stages. Such overpressures can lead to the formation of microfractures and enhanced migration of petroleum in and out of source rocks, referred to as “glacial pumping” (Sachse & Littke, Reference Sachse and Littke2018). Microfractures are possible conduits for petroleum migration; this process has been observed in the mature Posidonia Shale (Littke et al., Reference Littke, Baker and Leythaeuser1988). This mechanism might lead to an increased amount of hydrocarbons migrating and accumulating in the supra-salt petroleum systems sourced by the Coevorden Formation and Posidonia Shale in the Dutch LSB during glacial times. However, no pressure measurements in the Niedersachsen Group (Coevorden Formation) and Altena Group (Posidonia Shale) are available in the study area to further calibrate the model (Fig. 1). Particularly in places of lower overburden thickness, e.g. near salt domes, hydrocarbons might leak through the relatively thin top seal and release hydrocarbons in shallower units or into the atmosphere during glacial stages. Furthermore, the alteration of groundwater flow patterns and an influx of meltwater might lead to increased microbial activity and subsequent generation of hydrocarbons in organic-rich sediments (Person et al., Reference Person, McIntosh, Bense and Remenda2007). Even though this large-scale model does not represent small scale changes in the lithology, a top loss of hydrocarbons has been calculated, especially during the first Elsterian glacial stage retreat. This indicates the influence of glacial loading and unloading on the migration patterns in the subsurface and might impact the release of hydrocarbons into the atmosphere during the Quaternary glacial and interglacial stages.
Conclusions
The effects of Pleistocene glaciations on the temperature and rock properties of the subsurface were simulated using a 3D basin and petroleum system model of northeastern Netherlands in which (1) Quaternary sedimentary units and (2) sequential glacial loading and unloading were implemented. Integrating a 3D basin and petroleum systems model with a shallow subsurface model allowed the implementation of Pleistocene glaciations in a valid sedimentary framework. Simulated glacial loading and unloading revealed significant changes in the paleo- and present-day subsurface temperature and pressure distribution, rock properties and migration of hydrocarbons.
Two thickness scenarios were used based on published literature to account for the general uncertainty of the ice sheet thickness. Low ground temperatures during glacial stages effectively lowered temperatures in the Dutch subsurface until the present day. In the shallow subsurface, temperatures equilibrated back to a pre-glacial state during the warm Eemian interglacial period in which temperatures were slightly increased compared to the present day. During the Weichselian stage, low surface temperatures decreased subsurface temperatures, resulting in lower than expected shallow subsurface temperatures at the present day. The presence of salt diapirs strongly influences the temperature distribution and the propagation of low temperatures during glacial stages. Permafrost calculations revealed maximum permafrost thicknesses during the Weichselian and in between glacial advances in the Elsterian and Saalian glacial stages. In places with great salt thickness, the formation of permafrost is impeded due to the high thermal conductivity of salt. While source rock temperatures are lowered during and after glacial stages, the maturities of Palaeozoic and Mesozoic source rocks are not affected. This is the case in places of maximum burial today and in places where maximum burial occurred in the past.
Subsurface pressures are significantly altered by the loading and unloading of thick ice sheets during the Elsterian and Saalian glacial stages. In fine-grained Mesozoic formations and deeper intervals, subsurface pressures increased with every ice sheet advance and did not equilibrate back to a pre-glacial state at the present day. Therefore, glacial loading seems to be a potential contributor to elevated subsurface pressures in areas covered by Pleistocene ice sheets. Overpressure in supra-salt formations might contribute to enhanced petroleum expulsion and migration in Mesozoic petroleum systems. Porosities and permeabilities in the subsurface were permanently reduced during glacial loading in the Elsterian stage.
Hydrocarbon accumulations of the northeastern Palaeozoic and the Lower Saxony Basin Mesozoic petroleum systems were simulated using recently published kinetics. Sequential mechanical loading and unloading with ice sheets did result in compression and expansion of hydrocarbons in sub-salt layers. In the supra-salt layers, a top loss of hydrocarbons is calculated during the retreat of the Elsterian and Saalian ice sheets, indicating a possible trigger for migration of hydrocarbons and a loss of hydrocarbons to the surface.
Acknowledgements
This study was funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation), Project SA 3094/4-1. We thank TNO for providing the public databases NLOG (www.nlog.nl) and Schlumberger for providing the software Petromod under an academic licence agreement. Finally we thank reviewer Victor Bense, an anonymous reviewer and NJG associate editor Hanneke Verweij for their constructive comments which helped to improve the manuscript.
Open access
