Influence of Quaternary glaciations on subsurface temperatures, pore pressures, rock properties and petroleum systems in the onshore northeastern Netherlands

Abstract

Pleistocene glacial stages were implemented into a 3D basin and petroleum systems model of the northeastern Netherlands to address the influence of low surface temperatures and the mechanical loading of ice sheets on the subsurface. Two ice sheet thickness scenarios were used based on published data. Overall, Quaternary glacial stages have a substantial impact on the temperature and pressure distribution in the subsurface. Subsurface temperatures are significantly reduced during glacial stages, leading to lowered present-day temperatures and a low geothermal gradient in the shallow subsurface. In deeply buried sedimentary formations, pressures build up with every glacial advance resulting in overpressures at the present day. Glacial stages do not directly influence the petroleum generation of petroleum source rocks in the area, but high pressures during loading might have impacted petroleum expulsion of the early mature Coevorden Formation. Hydrocarbon accumulations in the Lower Saxony Basin were simulated to investigate the possible effects of mechanical ice loading and unloading on hydrocarbon migration. A loss of Coevorden Formation-sourced hydrocarbons to the surface was calculated in the Lower Saxony Basin during the glacial stages, indicating an influence of glacial loading on the Mesozoic petroleum system.

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., 1996; Duin et al., 2006; Wong et al., 2007; Littke et al., 2008; Doornenbal & Stevenson, 2010; Nelskamp, 2011; Scheck-Wenderoth & Maystrenko, 2013; Strozyk et al., 2017; Kley, 2018). 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. (2006), Munoz et al. (2007), Uffmann & Littke (2011), Abdul Fattah et al. (2012), Verweij et al. (2012a), Bruns et al. (2013), Bruns et al. (2016), Mohnhoff et al. (2016), Sachse & Littke (2018), and Amberg et al. (2022). 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, 2007). 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., 2011). 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, 2007). Mesozoic source rocks are the Toarcian Posidonia Shale and the Berriassian Coevorden Formation, also referred to as Wealden Shale in Germany (Leythaeuser et al., 1988; Rippen et al., 2013).

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, 1992; Lerche et al., 1997; Lang et al., 2014). 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. (1997), Cavanagh et al. (2006), Grassmann et al. (2010), Rodrigues Duran et al. (2013) and Sachse & Littke (2018). 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., 2006; Grassmann et al., 2010). 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., 2012; Ter Voorde et al., 2014). Ter Voorde et al. (2014) 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., 2022) 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., 2012). In places, Triassic and Permian strata underwent significant erosion during the Jurassic (Duin et al., 2006; Kombrink et al., 2012).

Fig. 1. 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 (2018) 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, 2013). 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, 2007). 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., 2008). A phase of thermal uplift, deep erosion of Carboniferous sediments and local magmatism followed in the Permian (De Jager, 2007). During the Late Permian, evaporitic cycles of the Zechstein group with a decreasing marine influence were deposited in a restricted basin (Strozyk et al., 2017). Continental conditions characterised the basin during Triassic times, with the exception of the Muschelkalk during the Anisian to Ladinian (Geluk, 2007; Stollhofen et al., 2008). 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, 2007). 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., 2011). Marine sedimentation of siliciclastics and carbonates prevailed in the Cretaceous and Cenozoic (De Jager, 2007). 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, 2007). 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, 2007). 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., 2013). 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, 2011; Moreau & Huuse, 2014; Lang et al., 2018). Deposition during interglacial periods was mainly marine, fluvial and aeolian (De Gans, 2007).

Fig. 2. Chronostratigraphic chart of the Neogene according to the general stratigraphic nomenclature of the Netherlands (NAM & RGD, 1980; Van Adrichem Boogaert & Kouwe 1994; Gunnink et al. 2013; TNO-GSN 2020) and Pleistocene ice sheet thicknesses assigned in the modelling (Schokking 1990; Feldmann 2002; Sachse & Littke 2018). 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, 1990; De Gans, 2007; Laban & van der Meer, 2011). Two major glaciations, the Elsterian and the Saalian (Fig. 2), partially covered the area of the northeastern Netherlands (Fig. 2; Ehlers et al., 2018; Lang et al., 2018). Weichselian ice sheets did not advance onto onshore Netherlands (Fig. 2; Laban & van der Meer, 2011). 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., 2018). 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, 1990; Huuse & Lykke-Andersen, 2000; Böse et al., 2012; Moreau et al., 2012; Müther et al., 2012; Lang et al., 2018).

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, 2011; Böse et al., 2012). 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., 2011). 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, 2011; Meijles, 2015). 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, 2007; van der Vegt et al., 2012). 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., 2007).

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., 2005; Laban & van der Meer, 2011; Lang et al., 2018). 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., 2008). 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., 2018). 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, 2007; Gunnink et al., 2013). 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., 1989; Laban & van der Meer, 2011; Böse et al., 2012; Lang et al., 2018). Rappol et al. (1989) proposed three phases of ice advances, which is debated controversially (Laban & van der Meer, 2011). 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, 2011). In contrast to the Elsterian glacial period, incising glacial valleys related to the Saalian glaciation have only been found in offshore Netherlands (De Gans, 2007; Moreau et al., 2012).

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 (2009) and Peters et al. (2017). 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. (2020a) and for the Coevorden Formation from Froidl et al. (2020b) were used for the petroleum generation of Palaeozoic and Mesozoic source rocks in the area (Amberg et al., 2022).

The recently published 3D basin model studying the burial and thermal history of the northeastern Netherlands from Carboniferous times on by Amberg et al. (2022) 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., 2007; Gunnink et al., 2013). Thickness maps of four key Quaternary intervals derived from the BRO DGM 2.2 model are shown in Fig. 4 (Gunnink et al., 2013). Lithologies for Quaternary and Neogene units in PetroMod are set according to Gunnink et al. (2013). 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 (2007), Ehlers et al. (2018) and Sachse & Littke (2018). 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 (1990) near Marum, northern Netherlands. Feldmann (2002) proposed a Saalian ice sheet thickness of 100–600 m in the study area, but thickness gradients near the ice margin vary strongly (Mrugalla, 2020). 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 (1990) and one with a maximum ice sheet thickness of 600 m based on the estimation by Feldmann (2002). 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.

Fig. 3. 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).

Fig. 4. Thickness of key Quaternary sedimentary unit intervals compiled from the BRO-DGM shallow subsurface model (TNO-GSN 2019).

Table 1. 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.

Table 2. 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., 2022). Paleobathymetric data, indicating the submerged amount of an ice sheet below sea level, was taken from Sachse & Littke (2018). The SWIT is calculated from the paleolatitude- and paleowater depth history (Wygrala, 1989) and was strongly influenced by the short term changes in surface and ground temperatures and the presence of ice sheets. Grassmann et al. (2010) and Sachse & Littke (2018) used mean annual ground temperatures based on a temperature model derived from marine proxy data by Delisle et al. (2007) 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., 2001). Paleo-surface temperatures in the Netherlands ranging from −7 to 12°C from Luijendijk et al. (2011) 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., 2010). 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, 2009). 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, 2011). 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., 2001). 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, 2001). 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/m³ (Sachse & Littke, 2018). The basal heat flow, independent of the low temperatures of ice sheets, was not changed from the original model.

Fig. 5. SWIT used in the basin model based on Lujiendijk et al. (2011), Grassmann et al. (2010) and Sachse & Littke (2018). 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, 2009). The extension of the Terzaghi rock stress modelling to a 3D-poroelasticity model by Biot (1941) further enables the calculation and prediction of geomechanical processes such as rock failure (Peters et al., 2017). 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., 2005). Water is assumed to be incompressible, and for the Terzaghi equation, it is assumed that the volume of solids does not change (Neuzil, 2003; Verweij, 2003). The mechanical compaction reduces porosity nearly irreversible (Houseknecht, 1987). 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, 2009). 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, 2003). 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., 2005; Aplin & Macquaker, 2011). 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., 2013). Other possible causes of overpressure generation are fluid flow, sea-level changes, changing groundwater volumes and salinities, oil and gas generation or cementation (Verweij, 2003). 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 (2009).

Heat transfer takes into account heat conduction, convection and radiogenic heat production (Peters et al., 2017). The implementation of permafrost is a purely conductive approach but takes into account latent heat (Grassmann et al., 2010). 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., 2013).

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.

Fig. 6. (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.

Fig. 7. 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.

Fig. 8. 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.

Fig. 9. 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.

Fig. 10. 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.

Table 3. 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).

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. (2020b). 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.

Fig. 12. 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., 2012; Ter Voorde et al., 2014). Luijendijk et al. (2011) 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. (2014) and Gies et al. (2021) 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., 2001; Luijendijk et al., 2011). 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. (2011), Ter Voorde et al. (2014) and Gies et al. (2021) 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., 2010). 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., 2021).

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., 1982), influencing the temperature distribution, particularly in porous sediments (Fjeldskaar & Amantov, 2018). During freezing, latent heat is liberated, while melting leads to an uptake of energy (Delisle, 1998). 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., 2014). Frenzel et al. (2001) 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. (2015), 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. (2021) 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., 2012; Gies et al., 2021) and the CEBS and are referred to as salt chimneys (Rodon & Littke, 2005; Magri et al., 2008). 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., 2010).

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. (2010) and Sachse & Littke (2018). A lower ice sheet thickness, as stated by Schokking (1990) and Feldmann (2002), is assumed in the study area attributed to the more southern location and the closer distance to the ice sheet margin. Winguth et al. (2005) 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, 2018) 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 (1990) and Feldmann (2002).

Fig. 13. 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., 2012b). 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., 2012b) 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., 2012b). 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., 2005; Rodrigues Duran et al., 2013). 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., 1997; Sachse & Littke, 2018). 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, 1995; Person et al., 2007; Rodrigues Duran et al., 2013). Additionally, high pressures from vertical ice loading can result in fault reactivation (Brandes et al., 2011) as well as salt movement (Lang et al., 2014) and affect Pleistocene and Holocene supra-salt sediments as observed in the Eastern Glückstadt Graben in Northern Germany (Huster et al., 2020). 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, 1992; Rodrigues Duran et al., 2013; Daszinnies et al., 2021). Lower subsurface temperatures during ice loading and the reduction of reservoir temperatures to the present day (e.G. Mittelplate oil field; Grassmann et al., 2010) 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., 2005). 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, 2018). 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., 2022) 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, 2007) 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., 2013). 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, 2011). 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, 2018). Hydrocarbon accumulations in the Palaeozoic Rotliegend reservoirs in the northeastern Netherlands contain only natural gas (De Jager & Geluk, 2007); 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, 2018). Microfractures are possible conduits for petroleum migration; this process has been observed in the mature Posidonia Shale (Littke et al., 1988). 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., 2007). 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.