3D basin modeling of the Hils Syncline, Germany: reconstruction of burial and thermal history and implications for petrophysical properties of potential Mesozoic shale host rocks for nuclear waste storage

Abstract

Jurassic sedimentary sequences suitable for nuclear waste storage in northern Germany consist of organic-lean claystone and were uplifted to < 100 m depth in the Hils Syncline area (southern Lower Saxony Basin). This Hils Syncline, showcasing a northwestward increase in thermal maturity, facilitates the study of shale petrophysical properties influenced by burial history. This study introduces a 3D-thermally calibrated numerical model of the Hils Syncline area to analyze its geodynamic evolution and maturity variations. It provides new vitrinite reflectance and sonic velocity data for modeling calibration and erosion estimation. The Hils Syncline area has undergone continuous subsidence, interrupted by a Cretaceous uplift documented by an erosional unconformity. During the latest Early Cretaceous, Jurassic rocks underwent maximum burial reaching up to several thousand meters depth and temperatures up to 160 °C in the northwest. The Late Cretaceous inversion caused stronger erosion towards the northwest removing up to 3300 m of sediment compared to about 1300 m in the south, according to vitrinite reflectance-based estimations. Numerical modeling results along the study area indicate decreasing porosity and permeability northwestward with increasing thermal maturity. Porosity and vertical permeability decreased to 5–14% and 2.8 × 10^–23 to 1.5 × 10^–19 m² [1 mD = 10⁻¹⁵ m²], respectively, while vertical thermal conductivity increased to 1.30–2.12 (W/m/K). These trends of porosity/permeability and thermal conductivity with burial align with sonic velocity and published experimental porosity data, except for the thermally most mature region (Haddessen). This anomaly is tentatively attributed here to localized overpressure generation in the Posidonia Shale during maximum burial, affecting both the underlying Pliensbachian and overlying Doggerian units.

Graphical abstract

3D numerical model of the Hils Syncline and surrounding area revealing that a northwestward increase in maximum burial resulted in higher temperatures and varying maturity levels. While most locations align well with calibration data (i.e. measured vitrinite reflectance and porosity), discrepancies arise in the Haddessen/Bensen area. The mismatch between porosity, vitrinite reflectance, and sonic velocity response indicates local overpressure in the northernmost region mainly during the Cretaceous. It was likely caused by gas generation in the Posidonia Shale affecting nearby Lower and Middle Jurassic units.

Introduction

Based on legally stipulated requirements and criteria (BGE 2022), areas within the Lower Saxony Basin (LSB; Fig. 1) in Northern Germany have been identified as sub-areas where favorable geological conditions can be expected for the final disposal of high-level radioactive waste. Aside from salt diapirs, claystone formations of Jurassic and Cretaceous ages are considered potential host rocks due to their low hydraulic conductivity, high sorption capacity, and self-sealing capacity (BGR 2007; Hoth et al. 2007). The LSB has been subject to strong subsidence during Late Jurassic–Early Cretaceous rifting, and affected by at least one major phase of uplift and erosion (Late Cretaceous inversion) (e.g., Senglaub et al. 2005; Kley and Voigt 2008; Voigt et al. 2008; Doornenbal and Stevenson 2010). This heterogeneous geodynamic evolution resulted in lateral variations of petrophysical and mechanical properties of potential host rock formations for nuclear waste disposal (Fink et al. 2019; Gaus et al. 2022).

Fig. 1

Overview of northwestern Germany where the Lower Saxony Basin is located. The limits of the basin as well as its major faults are illustrated. The location of the 3D model is highlighted with the red rectangle. Within this area, the Hils and Sack Synclines are represented by the dark purple polygons

The Hils Syncline in the southern area of the LSB displays highly variable thermal maturity within the same stratigraphic level (Toarcian Posidonia Shale) over a relatively short distance of 50 km, increasing towards the northwest from 0.48 to 1.45% vitrinite reflectance (VR_r) (Littke et al. 1988, 1991a, b) (Fig. 2). It can, therefore, be considered as a natural laboratory, where effects of burial and thermal maturation on petrophysical, chemical and mechanical properties of sedimentary rocks can be studied. In addition, the thermal evolution in this region, leading to this maturity variation, is still debated and is explained by two different premises. One suggests that the high thermal maturities resulted from a deep-seated mafic intrusion (the Bramsche and Vlotho massifs, located northwest of the Hils Syncline), emplaced during Late Cretaceous times (Bartenstein et al. 1971; Teichmüller et al. 1979, 1985; Rullkötter et al. 1988). A more recent hypothesis assumes instead that deep burial until the Late Cretaceous, followed by strong uplift and erosion is the primary cause of high maturities (Petmecky et al. 1999; Senglaub et al. 2005, 2006; Adriasola Munoz et al. 2007; Brink 2013).

Fig. 2

Study area and map view of the simulated 3D model. Major structural elements in the Hils Syncline and the different stratigraphic units are shown. The boreholes implemented for the calibration of the 3D model and their respective measured vitrinite reflectance values (%VR_r) are illustrated. The boreholes drilled from the 1980s campaign are represented by black dots with a red outline and boreholes from the BGR database are marked in black dots with a yellow outline. Newly drilled boreholes as part of the Maturity Project, which supports the search for a disposal site for highly radioactive waste in Germany, are marked as green dots with black outlines (Main: Mainzholzen, Hun: Hunzen, Har: Harderode and Bens: Bensen). The samples collected from outcrops are symbolized with the black X. The variable thermal maturity in the western flank of the Hils Syncline within the same Lower Jurassic stratigraphic unit can be observed and is derived from %VR_r measurements from Littke et al. (1988, 1991a, b). The measured %VR_r of the new outcrop samples (values in red; Table 2) indicate a thermal maturity increment in a northwest direction that fits basically with the trend shown in previous studies (e.g., Littke et al. 1988, 1991a, b) and the measured vitrinite reflectance of the BGR database

Several 2D and 3D studies on basin modeling of the LSB, investigating burial, temperature and maturity histories, have been published (e.g., Neunzert et al. 1996; Adrisola Muñoz et al. 2007; Uffmann et al. 2010; Uffmann and Littke 2011; Bruns et al. 2013, 2014; Fink et al. 2019). However, they were large-scale, regional models lacking high stratigraphic resolution and/or detailed facies information on the Mesozoic rocks. Furthermore, most of these studies did not include the entire area of the Hils Syncline, but just its northern part. Studies that correlated the petrophysical and mechanical properties of Mesozoic rocks to their paleo-depths and temperatures have been carried out only recently (Fink et al. 2019; Gaus et al. 2022). Nevertheless, they were either large-scale 3D models of more northern areas of the LSB (Fink et al. 2019) or the maximum burial depths were empirically derived from simple correlations that convert VR_r to temperature (Gaus et al. 2022).

In this study, the burial and thermal history of the Hils Syncline is for the first time reconstructed and quantified through 3D numerical modeling by testing different scenarios of basin subsidence and uplift, and thermal fields, seeking to improve the understanding of the geodynamic evolution of the study area and the observed maturity variation. Hence, the main target lies in constraining the maximum temperatures, deepest burial, and erosion amounts that have affected the Mesozoic record of the Hils Syncline area. This article provides new measurements of vitrinite reflectance from outcropping rocks (Fig. 2) collected in the field that are integrated with available vitrinite reflectance data from literature for the reconstruction of maximum paleo-depths and -temperatures.

Accurately assessing the erosion thickness related to the Late Cretaceous inversion event allows to determine the maximum paleo-burial depths and, as a result, the basal heat flow evolution can be better constrained. While vitrinite reflectance data are commonly used to calibrate these models and estimate the erosion thickness, sonic velocity data as an alternative independent parameter can be additionally used (Iyer et al. 2021). Sonic velocity data were acquired from recently drilled boreholes located on the western flank of the Hils Syncline (Fig. 2) as part of the recent research activities supporting the search for a disposal site for highly radioactive waste in Germany. Moreover, this work aims to provide evidence to evaluate the influence of the burial and thermal history on the present-day petrophysical properties of the Cretaceous and Jurassic rocks. Thus, the thermally calibrated 3D model presented here considers the heterogeneity of Cretaceous and Jurassic rocks (facies and lithological information), incorporates direct laboratory measurements on drilled sediment cores from Gaus et al. (2022) as calibration data, and implements these measurements to optimize the mechanical compaction curve of the assigned lithologies to the Jurassic units.

During burial, sedimentary rocks undergo compaction, which can be either mechanically triggered by effective stress, or related to mineral transformation (chemical diagenesis) due to the effects of chemical disequilibria, temperature and time (Bjørlykke 2006). This process results in porosity reduction and an increase in seismic velocity and in density (Storvoll and Brevik 2008). Mechanical compaction involves grain reorientation, pseudo-plastic deformation of ductile grains and brittle grain fracturing; while chemical diagenesis involves dissolution of less stable minerals and precipitation of thermodynamically more stable minerals (e.g., smectite transformation to illite) (Bjørkum et al. 1998; Chuhan et al. 2002; Bjørlykke 2006; Fawad et al. 2011). However, effects of compaction may be counteracted by overpressure generation (Japsen 1998). Overpressure can be generated either by disequilibrium compaction and/or fluid expansion (Mouchet and Mitchell 1989; Swarbrick and Osborne 1998). Disequilibrium compaction occurs when low-permeable rocks are rapidly buried leading to the pore fluids not being completely expelled (Wilson et al. 1977; Mouchet and Mitchell 1989; Osborne and Swarbrick 1997) and is identified by abnormally high porosity at a given depth (Swarbrick and Osborne 1998; Sayers et al. 2002). On the other hand, fluid volume expansion involves processes such as aquathermal expansion, diagenesis of clay minerals and gas generation (Mouchet and Mitchell 1989; Swarbrick and Osborne 1998; Tingay et al. 2007; Lahann and Swarbrick 2011), and produces more likely local overpressure (Swarbrick and Osborne 1998). In this study, anomalous porosity and permeability data of Lower and Middle Jurassic organic-lean mudstones are analyzed and discussed in relation to local overpressure generation.

Geological setting

The LSB is considered an important oil and gas province in Germany and has undergone extensive research with respect to its hydrocarbon potential. The LSB represents an E–W-striking highly differentiated Meso-Cenozoic basin (Betz et al. 1987; Baldschuhn et al. 1991) that is bounded by the Pompeckj Block to the north, Münsterland Basin and Rhenish Massif to the south, the Harz Mountains and Gifhorn Trough to the east, and the Friesland Platform to the west (Fig. 1). This study focuses on the Hils Syncline, which is an NW–SE striking structure located to the south of the LSB and comprises an area of ca. 800 km² (Fig. 1). The Hils Syncline, as part of the LSB, belongs to the Central European Basin system (CEBS) (Littke et al. 2008). Therefore, its geodynamic evolution is linked to the CEBS geological background and is described within a regional context below.

Geodynamic evolution of the LSB

The evolution of the LSB was initiated with the Latest Carboniferous–Early Permian emplacement of the CEBS on the Westphalian foreland basin (Baldschuhn et al. 2001; Maystrenko et al. 2008). During this period, the CEBS experienced faulting and widespread igneous activity (Gast 1988; Bayer et al. 1999; Abramovitz and Thybo 1999; Breitkreuz et al. 2008) followed by thermal subsidence during the Permian (Van Wees et al. 2000). Late Carboniferous and Early Permian times were dominated by terrestrial sediments deposited first in a tropical humid climate, resulting in abundant coal deposits, and later in an arid climate (i.e., Rotliegend formation). Several extensional pulses affected the CEBS during the Triassic, as the result of the beginning of Pangea’s break-up and formation of the North Sea, leading to differentiation into several sub-basins (Maystrenko et al. 2008). Shallow-marine to terrigenous deposition occurred from the latest Permian (Zechstein Formation) to the Late Triassic (Keuper Formation). Terrestrial influence ceased when the onset of Jurassic marine sedimentation started. During the Early Jurassic, transgression from the Tethys Ocean and continued thermal subsidence favored the deposition of shallow marine shales (Liassic unit) (Bachmann et al. 2008). The maximum transgression phase occurred during Toarcian time (Littke et al. 1991a, b; Barth et al. 2018) and was accompanied by anoxic conditions and the formation of organic matter-rich sediments (Posidonia Shale) (Littke et al. 1991a, b; Song et al. 2015, 2017; Ruebsam et al. 2019). Shallowing water depths during the Middle Jurassic are represented by more sandy deposits interbedded within the marine mudstones (Petmecky et al. 1999; Stollhofen et al. 2008).

As Pangea’s break-up continued, Late Jurassic–Early Cretaceous rifting within the CEBS was responsible for the LSB formation (Betz et al. 1987; Kley et al. 2008). Structural blocks neighboring the LSB (i.e., the Rhenish Massif and Pompeckj Block) experienced strong uplift while other parts of the basin were reaching the highest subsidence rates during this period (Bruns et al. 2013; Uffmann and Littke 2013). The Late Jurassic–Early Cretaceous was accompanied by the deposition of paralic to marine Wealden facies (Rippen et al. 2013; Froidl et al. 2021), while marine conditions predominated during the rest of the Cretaceous (Mutterlose 1992; Stollhofen et al. 2008). The convergence between Africa–Arabia and Eurasia resulted in compression and inversion of some parts of the basins along the CEBS in the Late Cretaceous (Bachmann et al. 2008; Kley et al. 2008; Voigt et al. 2008; Kley and Voigt 2008). The LSB underwent significant uplift and erosion due to this compressional deformation (Senglaub et al. 2005, 2006; Adriasola-Munoz et al. 2007; Maystrenko et al. 2008; Voigt et al. 2008), placing highly thermal mature rocks near the surface (Senglaub et al. 2006). The sedimentary record of Upper Cretaceous age is, therefore, scarce in the area. In the Cenozoic, the Pyrenean tectonic pulse might be responsible for another inversion event that affected the LSB (de Jager 2003) while sedimentation continued in other areas such as Pompeckj Block and Münsterland Basin (Bruns et al. 2013). Finally, strong glacial impact characterized the Quaternary (Petmecky et al. 1999).

Jurassic and Cretaceous formations of the Hils Syncline

The Early Jurassic and Middle Jurassic times represent an important change in the depositional conditions within the entire CEBS, the transition from a mainly continental environment (Keuper, Late Triassic) towards a fully marine environment (Early and Middle Jurassic) (Barth et al. 2018). The geological succession deposited during this period constitutes a prominent and widespread sedimentary feature in NW Europe. The Early Jurassic sedimentary succession, also known as Liassic, has been described in Northern Germany based on a core profile from a borehole located towards the east of the Hils Syncline area (Van de Schootbrugge et al. 2019). Hettangian sediments are characterized by coarse sandstone beds, while towards the Sinemurian and the Pliensbachian, the succession is dominated by partly laminated mudstone, partly enriched in organic matter (Van de Schootbrugge et al. 2019). In the Upper Sinemurian and Lower Pliensbachian, sandy limestone beds are locally present (Van de Schootbrugge et al. 2019). Organic matter-rich paper shales occur predominantly in the Toarcian (Posidonia Shale) (Van de Schootbrugge et al. 2019). In addition, detailed studies of the Pliensbachian and Toarcian units have been recently carried out by analyzing two adjacent cores from a borehole located in the western flank of the Hils Syncline area (Burnaz et al. 2024). There, the Upper Pliensbachian comprises organic matter-lean claystone characterized by the presence of abundant siderite concretions (Burnaz et al. 2024) (Fig. 3). The lowermost Lower Toarcian succession, similar to the Upper Pliensbachian interval, contains organic matter-lean claystone, while black shales with high pyrite abundance characterize the rest of the Lower Toarcian and the lowermost Upper Toarcian (Burnaz et al. 2024) (Fig. 3). The Toarcian succession comprises the Posidonia Shale, which is the most important source rock for conventional oil exploration in NW Germany (Leythaeuser et al. 1988; Wehner 1997; Blumenstein et al. 2008) and potentially for unconventional shale gas/oil exploration (Andruleit et al. 2012; Bruns et al. 2013, 2014). The Middle Jurassic succession, also called Doggerian, comprises in its lowermost part the Aalenian Opalinus Clay, which is considered a potential nuclear waste host rock in Switzerland (Nagra 2002; Marschall et al. 2005) (Fig. 3); lithologically, it is very similar to the underlying upper part of the Upper Toarcian. Towards the upper Dogger, the units become sandier and calcareous sediments appear (Weitschaft 1984; Stollhofen et al. 2008). The Upper Jurassic (Malm) in the study area is represented by limestone, sandstone and claystone (Klassen 1984) (Fig. 3).

Fig. 3

Lithostratigraphy of the Lower, Middle and Upper Jurassic, and Lower Cretaceous successions in the Hils Syncline area. Compiled after Mutterlose and Bornemann (2000), Stollhofen et al. (2008), van de Schootbrugge et al. (2019) and Burnaz et al. (2024)

The Lower Cretaceous succession is dominated by non-marine and marginal marine marlstone and claystone that comprise hiatuses towards the south (Mutterlose and Bornemann 2000). Coal seams with interbedded claystone and siltstone beds characterize the Berriasian strata (Wealden Formation; Fig. 3) (Mutterlose and Bornemann 2000). Valanginian sediments have been documented only locally in the central part of the Hils Syncline (Jordan 1994); thus, little information is available. The Hauterivian record is present outside of the Hils Syncline in the area of the Sack syncline and consists on marl- and claystone layers (Mutterlose and Bornemann 2000). Dark grey claystone represent the Barremian as well as Lower Aptian sediments (Mutterlose and Bornemann 2000), which are rich in organic matter (Littke et al. 1998) and overlain by organic matter poor Middle Aptian to Albian sediments. The Albian is subdivided into three members: the Hils Sandstone (lower), Neohibolites minimus Clay (middle) and ‘Flammen-Mergel’ (upper) (Mutterlose and Bornemann 2000). Upper Cretaceous sediments are absent in the study area; they are preserved in adjacent areas towards the northeast (Fig. 2).

Methods and data

Elemental analyses

Total organic carbon (TOC) and total inorganic carbon (TIC) contents were determined using a Liquid TOC II (Elementar) analyzer. About 100 mg of pulverized sample was combusted in the presence of oxygen, first at a temperature of about 550 °C for 600 s, in order to measure TOC, and then during a second heating stage at 1000 °C held for 400 s to quantify TIC. TOC and TIC were calculated based on the amount of CO₂ released during each heating stage, which was assessed by an infrared (IR) detector with a detection limit of 10 ppm. To calibrate the measurements, a standard with a known total carbon (TC) content of 1.6 wt% was used. The total sulfur content (TS) of the samples was measured by implementing a FLASH 2,000 Organic Elemental Analyzer by Thermo Fischer Scientific (detection limit of 20 ppm and accuracy of ± 0.2%); the procedure followed the guidelines outlined in Prinz et al. (2017). Prior to geochemical analysis, samples (n = 30) were ground into powder and then homogenized.

Rock–Eval pyrolysis

Rock–Eval^® Pyrolysis was conducted on 25 outcrop samples utilizing a Rock–Eval^® 7 pyrolyzer (Vinci Technologies) and following the principles outlined by Espitalié et al. (1997) and Grohmann et al. (2023). To begin the process, small crucibles were filled with up to 100 mg of pulverized sample. The next step involved heating the samples in the absence of oxygen. First, the samples were subjected to a constant temperature of 300 °C for three minutes, during which any pre-existing hydrocarbons were evaporated and identified as the S1 peak, expressed in mg HC/g rock. Second, the temperature was gradually raised at a rate of 25 °C per minute until it reached 650 °C, leading kerogen to crack and resulting in the quantification of the generated hydrocarbons, represented by the S2 peak (mg HC/g rock). The temperature of maximum pyrolysis is represented by the Tpeak and depicts the T_max, which is an indicator of maturity. The measurements were calibrated using the IFP 16,000 standard. The analytical procedure’s specifics are described in Behar et al. (2001).

Vitrinite reflectance

Measurements of vitrinite reflectance (VR_r) were performed at a wavelength of 550 nm in oil immersion by employing a Zeiss Axio Imager.M2m microscope equipped with PI 10 × /23 ocular lenses and a 50 × /1.0 Epiplan-NEOFLUAR Oil Pol objective. Prior to the reflector measurements, the microscope was calibrated. For this propose, an Yttrium Aluminum Garnet standard (YAG) with 0.889% reflectance was used. Please refer to Zieger and Littke (2019) for detailed information on the preparation of polished samples and set up of microscopic equipment. For each sample, at least 10 vitrinite particles were measured (when possible) and then the mean random reflectance values were calculated using the DISKUS Fossil software from Hilgers Technisches Büro.

Basin modeling

The numerical basin modeling was performed using Slb PetroMod^® 2022. The simulation works as a finite-element forward modeling approach and allows to reconstruct the basin evolution through geological time intervals, integrating geological events such as sedimentation and erosion with geochemical, geomechanical and petrophysical processes. Detailed principles of basin modeling and the workflow that is required to use the software are described in Hantschel and Kauerauf (2009). The input data required for the model construction include the present-day geological setting and geometry, absolute ages for each stratigraphic unit, erosion events, facies distribution and lithological information, paleo geometries and thermal boundary conditions. After the simulation, the calculated results have to be compared against calibration data in order to validate their geological reliability. For this propose, the thermal and maturity histories were calibrated based on measured vitrinite reflectance and the bottom-hole temperature was used to constrain the present-day temperature.

Model outline and geometry

The 3D basin model is located in the southeastern area of the LSB, and covers the Hils and Sack Synclines, as well as surrounding areas (Figs. 1 and 2). It extends 106 km in an east–west direction and 96 km in a north–south direction, covering a total area of 10,300 km² (Table 1). The resolution of the model is restricted to a maximum cell size of 250 × 250 m.

Table 1 Representative input data for building the 3D model and reconstructing the burial and thermal history of the study area

The 3D model consists of 17 layers from the top of the upper mantle to the present-day topography (Table 1). The present-day stratigraphy and thickness values of the formations present in the study area were established based on different sources. The geometries of Pre-Permian layers (top of the Upper Mantle, Lower Crust, Variscan Crust and Pre-Permian clastics) are based on a lithosphere-scale 3D structural model of the Central European Basin System (Maystrenko et al. 2010, 2020). The Permian (base of the Permian Zechstein), Mesozoic and Cenozoic formations as well as the present-day topography correspond to the recently published TUNB 3D model (BGR et al. 2022), and were built based on the interpretation of available structural geological data and information from maps, boreholes and reflection seismic data. Among the Jurassic and Cretaceous formations, the Upper Jurassic unit is characterized by the highest thickness (up to ca. 1300 m) within the Hils Syncline, while the Lower and Upper Cretaceous rocks are only locally present in the center of the syncline with a maximum thickness of up to 230 m (Fig. 4).

Fig. 4

Present-day thickness maps of the main Mesozoic units (Jurassic and Cretaceous) compiled from the TUNB 3D model (BGR et al. 2022)

Stratigraphy and facies/lithology

The different stratigraphic units, their assigned depositional ages, lithologies and petrophysical properties, are shown in Table 1. The absolute ages were determined based on the Stratigraphic Table of Germany (STD 2016) and the International Chronostratigraphic Chart v.2022/10 (Cohen et al. 2013). The lithological information and facies distribution of Pre-Jurassic units were established according to Littke et al. (2008), Stollhofen et al. (2008) and Doornenbal and Stevenson (2010).

The facies and lithologies of the Jurassic and Cretaceous formations were defined based on Mutterlose and Bornemann (2000), Barth et al. (2018), van de Schootbrugge et al. (2019), Gaus et al. (2022) and the Borehole Database of Lower Saxony (LBEG—NIBIS^® map server). For these units, different user-defined lithologies have been created by mixing the standard PetroMod lithologies to accurately account for the stratigraphic descriptions and depositional conditions. Their lithological characteristics have been automatically assigned according to the defined mixed lithologies (Table 1). Furthermore, lateral variations within the area covered by the Hils Syncline were taken into account, so that the created rock mix can vary from one location (borehole) to another.

Boundary conditions

The upper thermal boundary conditions include the paleowater depth and the sediment–water interface temperature (SWIT) (Peters et al. 2017). The paleowater depth was determined according to the data from TNO-GSN (2020), earlier publications (e.g., Littke et al. 2008) and numerical models previously published for the area (e.g., Rodon and Littke 2005; Adriasola Munoz et al. 2007; Uffmann and Littke 2011; Bruns et al. 2013) (Fig. 5). The SWIT was calculated using PetroMod’s built-in automatic SWIT based on the paleogeographical latitude (52° N) and the paleowater depth evolution (Wygrala 1989) (Fig. 5).

Fig. 5

General boundary condition trends. The boundary conditions displayed in this figure correspond to the simulation perform in the northern boundary of the Hils Syncline (Haddessen)

On the other hand, the lower thermal boundary condition is represented by the basal heat flow (Peters et al. 2017). The heat flow trend assigned to this model was derived from different sources and in general, follows the classical heat flow calculations proposed by McKenzie (1978) (Fig. 5). During the early Permian the first peak of heat flow with values of up to 90 mW/m² (Uffmann and Littke 2011; Bruns et al. 2013) represents a phase of intense volcanic activity. A second peak of high heat flows is set during the Late Jurassic to Early Cretaceous, reaching up to 80 mW/m², as the result of extensional tectonics occurring in the study area (Adriasola Munoz et al. 2007; Bruns et al. 2013). Gradual cooling following the periods of maximum heat flow was assumed, with heat flow values that decreased to 60 mW/m², typical of continental heat flow (Allen and Allen 2013). Present-day basal heat flows were assigned according to heat flow values from Fuchs et al. (2022) and by multi-1D calibration using temperature data from boreholes. Two scenarios were modeled regarding the paleo-heat flow trends. The first one assumes the same heat flow trend over the study area, and the second one takes variations into account based on the geotectonic evolution, maturity level and erosion magnitudes (see chapter “Results”).

Calibration data

In order to validate the thermal and burial histories simulated by the 3D modeling, the model results were calibrated with measured VR_r data. The calibration data consist of unpublished master theses at RWTH Aachen University, the hydrocarbon database of the Federal Institute for Geosciences and Natural Resources (BGR), boreholes from the 1980s campaign and recently drilled boreholes (e.g., Mainzholzen; Burnaz et al. 2024). In addition, analyses conducted on outcrop samples that were taken during a field campaign in 2022 were also integrated (see Fig. 2 for location of the samples and Table 2). From the database, 27 boreholes were selected based on their location (Hils and Sack Synclines and surrounding areas; Fig. 2) and the availability and quality of data. Borehole temperature and vitrinite reflectance data were implemented for the thermal calibration of the 1D models. As a result of the multi-1D calibration approach, erosion magnitudes and heat flow values were determined and subsequently used to build and calibrate the 3D model. Vitrinite reflectance data were calculated using the Easy %R_oDL kinetic algorithm of Burnham et al. (2016) and finally compared to the measured data. The EASY %R_o models incorporate four reactions (generation of water, carbon dioxide, higher hydrocarbons, and methane) in one activation energy distribution that allows the calculation of vitrinite reflectance values between 0.3 and 4.5%VR_r (Sweeney and Burnham 1990; Froidl et al. 2020).

Table 2 Total organic and inorganic carbon (TOC, TIC), Rock–Eval Pyrolysis and VR_r data of outcrop samples collected from the Hils Syncline and adjacent area

Erosion

One major phase of erosion affecting the study area is considered in this model: the Late Cretaceous–Tertiary inversion (ca. 86 to 70 my). This event was responsible for significant uplift and erosion of mainly Lower Cretaceous and Upper Jurassic sediments and even Triassic deposits in the study area.

The erosion map of this event was compiled by applying the multi-1D approach proposed in previous studies (e.g., Bruns et al. 2013; Amberg et al. 2022). Twenty-seven (27) depth extractions were performed from the 3D model at the locations of boreholes where calibration data (temperature and maturity) were available. Every 1D model resulting from the extractions was calibrated using vitrinite reflectance and borehole temperature data. The erosion amounts at the well locations were assigned and then spatially interpolated using the built-in Laplace interpolation method in PetroMod. After each 3D simulation, the amounts of eroded overburden were slightly corrected to reach a suitable match between the calculated maturities and the calibration data. The erosion magnitudes and timing, which resulted once a proper calibration was obtained, are part of the results of this study (see chapter “Results”).

Sonic velocity logs and erosion estimates

Sonic logs, derived from elastic wave travel time (Δt), reveal valuable insights into geological formations, such as velocity, which is influenced by rock properties like matrix, porosity, fluid content, and pressure. These logs aid in calibrating seismic data, identifying lithologies, compaction trends, and stratigraphic correlations (Asquith and Krygowski 2004; Glower 2016; Tixier et al. 1959). A sonic probe measures the time for an elastic wave pulse to travel from transmitter to receiver, undergoing dispersion and attenuation as it traverses the rock. This results in various wave types arriving at different times, including the low-amplitude P-wave followed by the slower, higher amplitude S-wave, alongside other waveforms. The time difference between pulse transmission and the first P-wave arrival is the one-way time, which, when combined with transmitter–receiver distance, allows the calculation of formation slowness (Δt [s/m]) or formation velocity (V = 106/Δt [m/s]).

Sonic logging in this study was conducted in boreholes Bensen 1–2, Harderode, Hunzen 1–2, and Mainzholzen (Fig. 2). Data acquisition and depth matching for the downhole logging data were accomplished using GEOBASE software (Antares, Germany) and, subsequently, WellCAD software was employed for the processing and visualization of sonic logs. The applied sonic tool, known as a borehole-compensated device (BHC), employs two receivers to measure the time difference in the arrival of elastic waves from a transmitter (Tx) at each receiver (Rx1, Rx2). The BHC design enhances measurement reliability by reducing the impact of borehole size variations and tool tilt errors (Kobesh and Blizard 1959; Schlumberger 1972). Calculation of Δt is expressed as

$$\Delta t \left(\frac{s}{m}\right) = \frac{\left(\mathrm{Far }1-\mathrm{Near }1\right)+(\mathrm{ Far }2-\mathrm{Near }2)}{2\times ({\text{Far}}-{\text{Near}})}$$

(1)

where distances Far and Near refer to specific values for the transmitter-to-receiver distances of 3 m and 2.6 m, respectively.

To estimate erosion thickness, sonic velocities of the target formation (Lias) are compared to a reference trend, calculated based on a velocity-depth function by Scherbaum (1982):

$$v=1535 \left(\frac{m}{s}\right)+0.58z$$

(2)

where v and z are formation velocity and depth, respectively. The method involves shifting present-day sonic velocity data for the Lias formation until it intersects the reference trend, determining the paleo-burial depth, and then calculating erosion thickness as the difference between present-day and paleo-burial depths (Lasabuda et al. 2021).

Results

The following chapter describes the results of the analyses performed on the outcrop samples, the analyzed velocity logs and the output of the calibrated 3D model. Elemental analyses and Rock–Eval Pyrolysis were carried out for 25 outcrop samples. Regarding the vitrinite reflectance, not all the samples contained measurable organic particles; therefore, based on the number of pinpoint measurements and the accuracy of the results, 17 samples with sufficient vitrinite particles have been considered for this study and for calibrating the 3D model. Multi-1D models are also presented in this section in order to show in detail the burial and thermal histories at the relevant locations. These locations correspond to the boreholes that present the most reliable thermal maturity data (Haddessen, Harderode, Dohnsen, Dielmissen, Wickensen and Wenzen; Fig. 2). In this chapter, the processed velocity logs from the recently drilled boreholes in the western flank of the Hils Syncline are present as well as the erosion estimations carried out based on the data. In addition, 1D plots of porosity, vertical permeability and vertical thermal conductivity histories are also shown to better understand the evolution of the petrophysical properties of Jurassic units, and their relation with the burial and thermal histories. The description of the maturity will be related to the Lower Jurassic units since the maturity data used in this work was mainly measured at this stratigraphic level at all the locations.

Elemental analyses

The analyzed samples correspond to Mesozoic rocks deposited during the Triassic, Jurassic and Early Cretaceous. The content of organic matter in most of the samples is low to moderate with TOC values between 0.11 and 1.70 wt.% (Table 2). Higher values were measured for two samples located to the north (samples 3 and 5; see Fig. 2 for the location of the samples), with TOC contents of 6.60 wt.% and 24.73 wt.% (a coaly sample), respectively. These samples are probably of Lower Cretaceous age based on the geological map of the area.

TIC contents range from 0.08 to 12.48 wt.% (Table 2) corresponding to calculated calcium carbonate contents of < 1% to 100% (CaCO₃ = TIC * 8.333). There is not a distinctive trend of increasing or decreasing TOC and TIC contents in the analyzed samples and a noteworthy correlation between TOC and TIC contents cannot be observed.

Rock–Eval pyrolysis

Rock–Eval pyrolysis results are listed in Table 2. For all samples, S2 values are reported. T_max values are based on the S2 peak, but only reliable for S2 peaks > 0.5–1 mg HC/g rock. HI values for rocks with lower S2 values were also reported in brackets, but not for rocks with extremely low S2s (less than 0.1 mg/g rock). HI values prove the presence of type III kerogen, which is commonly derived from terrestrial biomass (higher land plants, precursors of vitrinite). This conclusion is confirmed by the very rare occurrence of liptinites, which represent hydrogen-rich phytoclasts. The few reliable T_max values, including those in brackets, range from 430 to 438 °C, which is typical of moderately mature organic matter; they thus fit well to the measured vitrinite reflectance. More information on geochemical maturity parameters in the old (1980s) wells along the western flank of the Hils Syncline is found in Fang et al. (2019) and Stock et al. (2017).

Vitrinite reflectance

VR_r is the main calibration parameter used in this study for reconstructing the thermal history of the sedimentary rocks within the Hils Syncline and adjacent area, serving as a thermal maturity indicator. Vitrinite reflectance values vary between 0.43% and 1.35%VR_r (Table 2), showing an increment in northwest direction (Fig. 2). Sample 6 having a very high TOC content is a coal (or coaly shaly) providing excellent vitrinite reflectance data based on measurements on large collotelinite particles. Samples 4 and 5 contain abundant organic particles in stylolites while the remaining carbonate rocks are devoid of organic particles. In samples 4 and 5, solid bitumen might be present, with a reflectivity somewhat lower than that of “true” autochthonous vitrinite. A widely applied conversion scheme is published by Schoenherr et al. (2007).

Basin modeling

Burial history

The basin evolution is characterized by continuous subsidence along the study area that was interrupted by uplift during the Late Cretaceous resulting in a major erosion event (Fig. 6). During the late Permian (ca. 257 Ma), a phase of stronger subsidence began, represented by steeper curves of the 1D burial plots, when the deposition of the Zechstein and Rotliegend units occurred. The basin continued experiencing subsidence at lower rates during the Triassic and Jurassic. The latest Late Jurassic marks the beginning of very strong tectonic subsidence. Sediments reached the maximum burial depths during the latest Early Cretaceous. The deepest burial of Jurassic and Cretaceous rocks occurred in the northern regions (e.g., Haddessen), while the maximum burial depth decreased towards the south (e.g., Wenzen; Fig. 6). Jurassic rocks were buried to a maximum depth of about 1600 m in the southeast and up to 3600 m in the northwest, whereas Lower Cretaceous sediments reached a depth of 1250 m and 3000 m, respectively.

Fig. 6

Burial and thermal history of the six 1D extractions that correspond to the location of the boreholes in the western flank of the Hils Syncline (Fig. 2). The white lines indicate the temperature isotherms at 100 °C, 200 °C and 300 °C. The red lines represent the Lower Jurassic unit (top and base) that displays the highly variable thermal maturity (Fig. 2). The abbreviations on the right side of the plots indicate the stratigraphic units and the basement. From top to base, they are labeled as follows: UJ Upper Jurassic, MJ Middle Jurassic, LJ Lower Jurassic, BS Upper, Middle and Lower Buntsandstein, Z Zechstein, PPC Pre-Permian Clastics and B Basement

Following the strong tectonic subsidence, major erosion due to uplift took place during the Late Cretaceous inversion of the LSB (Fig. 6). The inversion caused severe uplift and erosion that affected a large part of the stratigraphic column, especially Upper Cretaceous sediments that are only present towards the northeasternmost area and absent within the Hils Syncline. Jurassic sediments were also removed; however, Upper, Middle and Lower Jurassic rocks are partly preserved (drilled) in the Hils Syncline. Based on our calibration and model output, the thickness of erosion increased in a northwest direction, i.e., up to 3300 m of sediment was eroded at the location of Haddessen while ca. 1300 m of sediment was removed at the location of Wenzen (Fig. 7).

Fig. 7

Thickness of eroded sediments during the major erosion event in the Late Cretaceous (72–66 Ma). The erosion map resulted from calibrated multi-1D simulations at borehole locations (black dots) that were then extrapolated to the entire study area. The Hills and Sack Synclines and the main structural elements are represented as black lines

Thermal history

The modeled present-day maturities match the available vitrinite reflectance data (%VR_r) showing good calibration of our model (Fig. 8). In order to achieve this optimal calibration, we assigned paleo-heat flow magnitudes (see chapter “Methods and data”) and evaluated the effect of subsidence and inversion (erosion amount) affecting the stratigraphic record present in the study area. For the Late Jurassic, an elevated basal heat flow was assigned corresponding to the extension occurring in the region during this period. As defined in our model, the heat flow during the Late Jurassic peak reached values of 85 mW/m² towards the north (Haddessen) gradually decreasing in southeast direction to 75 mW/m² (Wenzen) (Fig. 9). Lower basal heat flows were assumed during the Late Cretaceous displaying values from 65 mW/m² in the north to a minimum of 60 mW/m² towards the south. These values were kept constant in our model until the Oligocene. From the Neogene to present day, slightly higher basal heat flows are deduced from present-day borehole temperatures and according to the present-day heat flow map from Fuchs et al. (2022). However, the present-day basal heat flows do not exceed the values reached during the Late Jurassic.

Fig. 8

Present-day maturity calibration plots showing the best-fit model of 1D extractions at the selected boreholes in the western flank of the Hils Syncline (Fig. 2). The red lines represent the calculated maturity using Easy %R_oDL (Burnham et al. 2016) with vitrinite reflectance data. The crosses correspond to the available vitrinite calibration data. The abbreviations on the right side of the plots indicate the stratigraphic units. From top to base, they are labeled as follows: LC Lower Cretaceous, UJ Upper Jurassic, MJ Middle Jurassic, LJ Lower Jurassic, UK Upper Keuper, LK Lower Keuper, UB Upper Buntsandstein, MB Middle Buntsandstein, LB Lower Buntsandstein, Z Zechstein, and PPC Pre-Permian Clastics

Fig. 9

Basal heat flow trend distribution showing heat flow trends through geological time (146 Ma, 90 Ma, and present day). The Hils and Sack Synclines and the main structural elements are indicated as black and white lines

Based on the above-described burial and heat flow histories, the thermal evolution shows constantly increasing temperatures from the late Permian onwards within the sedimentary layers (Fig. 6). The major erosion event starting in the early Late Cretaceous marks the beginning of a steady decrease in the temperature regime, which continues to the present day. The entire stratigraphic succession in the Hils Syncline (Pre-Permian clastics to Lower Cretaceous sediments) experienced maximum temperatures in the latest Early Cretaceous due to deep burial at moderate heat flows. During this period of maximum temperature, the Jurassic rocks reached about 90 °C in the southernmost part of the study area (Wenzen), while the oldest Cretaceous rocks experienced about 75 °C. In the northernmost region (Haddessen), the Jurassic units were subjected to up to 160 °C, and the oldest and deepest Cretaceous rocks to a maximum of 150 °C (Fig. 6).

Figure 10 shows the calculated maturity at the top of the Lower Jurassic succession using the Easy%R_oDL kinetics model of Burnham et al. (2016). During the latest Jurassic (145 Ma), the unit remained immature along almost the entire area with the vitrinite reflectance values ranging between 0.26 and 0.29%VR_r (Fig. 10a); locally, Lower Jurassic rocks entered the oil window towards the northern domain and in the center of the Hils Syncline (Fig. 10a). At the end of the Early Cretaceous and beginning of the Late Cretaceous (99 Ma)—the period of maximum temperatures and deepest burial—the succession reached its highest maturity and a general decreasing trend in southeast direction is observed (Fig. 10b). In the northwesternmost area, about 10 km north of Haddessen, the maximum vitrinite reflectance was up to 2.5%VR_r (wet gas stage). At the location of Haddessen, the Lower Jurassic unit reached 1.48%VR_r, which is very close to the measured vitrinite reflectance of this unit (Fig. 2). On the other hand, the lowest maturities down to 0.45%VR_r (immature with respect to petroleum generation) are calculated for the top of the Lower Jurassic in the southernmost part of the Hils Syncline (e.g., borehole Wenzen) and within most of the Sack Syncline in the southeast (Fig. 10b). The maturation of the succession was halted during and after the Late Cretaceous inversion and the maturity level has remained the same until the present day, i.e., the rocks have not experienced deeper burial and have not been subjected to higher temperatures during the Cenozoic.

Fig. 10

Calculated maturity map at the top of the Lower Jurassic unit at a 145 Ma and b present day. The kinetic model used is Burnham et al. (2016) Easy%R_oDL. Note that at present day, part of the Lower Jurassic sequence was eroded as result of the Late Cretaceous erosion phase

Petrophysical properties evolution

In this section, the petrophysical properties evolution of the Lower and Middle Jurassic units is compared based on 1D extractions at the location of Wenzen (southernmost), Harderode (northwest) and Haddessen (northernmost) (Fig. 11). Wenzen and Haddessen represent endmembers of maximum burial and Cretaceous erosion along the western flank of the Hils Syncline; Harderode is located in between. The lithology assigned to the Lower and Middle Jurassic sequences is considered laterally rather homogeneous and clay-rich. The Cretaceous units are not taken into account for this description due to their absence in most of the study area.

Fig. 11

Burial history plots of wells Wenzen, Harderode and Haddessen extracted from the 3D model (see Fig. 2 for the location of the wells). The overlays correspond to the porosity (a–c), vertical permeability (d–f) and vertical thermal conductivity (g–i) of the Lower and Middle Jurassic units. The labeled values on top indicate the initial values of the properties and the labeled values on the bottom indicate the values during the maximum burial. Differences in the initial values (at deposition) are due to paleowater depth variations (as boundary condition) from one location to another. The strong loss of porosity and vertical permeability is observed, and the increase of vertical thermal conductivity during the deepest burial in the latest Early Cretaceous. The abbreviations on the right side of the plots indicate the stratigraphic units. From top to base, they are labeled as follows: MJ Middle Jurassic, LJ Lower Jurassic, Z Zechstein and PPC Pre-Permian Clastics

The modeled porosity (%), vertical permeability (10^–15 m² [mD]) and vertical thermal conductivity (W/m/K) of Lower and Middle Jurassic rocks show variations through geological time and differ from one location to another. The porosity and vertical permeability of the Lower and Middle Jurassic units decreased by burial due to mechanical compaction, whereas the vertical thermal conductivity increased during burial (Fig. 11).

During the time of maximum burial depth (latest Early Cretaceous), in the southernmost area (e.g., borehole Wenzen), the porosity remained at values of about 14% (Lower Jurassic) and 17% (Middle Jurassic), showing a porosity reduction of about 65% when compared with the initial porosities (Fig. 11a). In Harderode, the porosity of Lower and Middle Jurassic rocks was reduced to 8% and ca. 9%, respectively (Fig. 11b), while in the most deeply buried region (Haddessen), the porosity of the Lower and Middle Jurassic rocks decreased to ca. 5% and 8%, respectively (Fig. 11c). The greater the maximum burial depth, the higher the loss of porosity; therefore, the porosity loss of the Jurassic units within the northernmost region was about 90%.

The vertical permeability behaves similarly to the porosity being reduced the most during the time of deepest burial. Initial values of vertical permeability vary between 0.035 and 0.046 (10^–15 m² [mD]) for the Lower Jurassic unit and between 0.015 and 0.053 (10^–15 m² [mD]) for the Middle Jurassic unit along the study area (Fig. 11d–f). Vertical permeability decreased during maximum burial to a minimum of about 1.5 × 10^–14 mD in the southernmost region (Wenzen; Fig. 11d), and less than 2.8 × 10^–8 mD in the northernmost area (Haddessen; Fig. 11f). The lowest vertical permeability is thus displayed at the location of Haddessen.

On the other hand, the Lower and Middle Jurassic units show higher vertical thermal conductivity at the time of deepest burial (latest Early Cretaceous) than at the initial state. Initial values of vertical thermal conductivity do not display significant variations from one location to another; they range from about 1 to 1.2 (W/m/K) along the analyzed locations (Fig. 11g–i). Towards the south (Wenzen), the vertical thermal conductivity of the analyzed units increased to 1.92 (Lower Jurassic) and 1.33 (W/m/K) (Middle Jurassic) during maximum burial (Fig. 11g), while in the northernmost area (Haddessen), it increased to 2.12 for the Lower Jurassic unit and to 1.53 (W/m/K) for the Middle Jurassic unit (Fig. 11i).

Based on the modeling results, the porosity, permeability and thermal conductivity at present day remain similar to the values displayed during the maximum burial (latest Early Cretaceous) and are not significantly affected by the following uplifting phase during the Late Cretaceous. Thus, based on the simulation performed in this study, a general trend of decreasing porosity and vertical permeability in northwest direction is observed at present day for the Lower and Middle Jurassic sediments. In contrast, the present-day vertical thermal conductivity of Jurassic rocks increases towards the northwest.

Sonic velocity logs and erosion estimates

Sonic logs were measured in recently drilled wells Bensen 1, Bensen 2, Harderode, Hunzen 1, Hunzen 2 and Mainzholzen (Figs. 2 and 12). Besides the above-mentioned sonic logs, sonic velocity data from the borehole Wenzen, drilled in a campaign in the 1980s, were also implemented for the erosion estimation approach (Fig. 13a). The interval analyzed in the logs corresponds to organic-lean shales of Pliensbachian age (Lower Jurassic). The overlying, organic matter-rich Toarcian Posidonia Shale was also drilled in these wells, except for the new Harderode borehole, but not considered for sonic velocity analysis.

Fig. 12

Sonic velocity logs of the P-wave for the boreholes Bensen 1–2, Harderode, Hunzen 1–2, Mainzholzen and Wenzen. Note that two separate boreholes were drilled and measured at close by locations in Bensen and Hunzen

Fig. 13

Velocity trend plots illustrating the erosion estimates from compaction techniques using sonic velocity analysis from boreholes a Wenzen, b Mainzholzen, c Hunzen, d Harderode and e Bensen. f Location of the analyzed boreholes. The velocity well log recorded for the uplifted formation is symbolized by the orange dots and represents the present-day depth of the studied unit (Pliensbachian). The gray dots denote the formation at its maximum burial depth (paleo-depth), which resulted from shifting the present-day depth dataset down until it touches the reference trend (black line). The blue arrow points to the values within the first part of the sonic logs, which are the most representative data. The erosion thickness is the difference between the present-day burial depth and the paleo-depth

In total, the average velocity transmittances of the Pliensbachian shale ranged between 2100 and 3500 m/s across all surveyed locations. A consistent increase in sonic velocity was observed with increasing present-day depth across almost all logs with the exception of the northernmost borehole Bensen (Figs. 12 and 13). The progression of average velocities followed the sequence Harderode (2900 m/s) > Bensen (2620 m/s) > Hunzen (2600 m/s) > Mainzholzen (2280 m/s) > Wenzen (2200 m/s).

It is important to note that the endmember values of sonic velocities in this study were potentially associated with fault zones (relatively small values) and siderite enrichment within the drilled cores (relatively larger values). To address the erosion thickness estimation, values within the initial one-third segment of the sonic logs were considered for downward adjustment until intersecting with the reference trend, as illustrated in Fig. 13. The erosion estimations obtained by applying a straightforward sonic velocity analysis vary from about 1300 m (Wenzen) up to ca. 2100 m (Harderode) of eroded overburden. Overall, the data show the same increasing trend of erosion thickness to the northwest as proposed in the vitrinite reflectance-based estimate (3D model approach) and the magnitudes are very similar as well (Fig. 13).

Discussion

Burial and thermal history

Along the western flank of the Hils Syncline, the thermal maturity of the Lower Toarcian, an about 50 m thick layer including black shales, increases from 0.48 in the southeast (Wenzen) to 1.45% VR_r in the northwest (Haddessen; Fig. 2). Previous studies have been carried out to explain this highly variable maturity (over a short distance of about 50 km) proposing contrasting hypotheses; the reasons are, therefore, still a matter of debate. Figure 14 shows the burial and thermal histories resulting from different thermally calibrated modeled scenarios that represent different theories. Most of the calibration data (%VR_r and temperature) cover up to Jurassic rocks, which are the main focus of this work, resulting in higher uncertainty in deeper domains. The plots compared in Fig. 14 correspond to the 1D extraction at the location of Haddessen (northwest), where the highest thermal maturity was measured (Fig. 2). At locations further to the southeast, where thermal maturity is lower, less elevated paleo-heat flows and lower eroded thicknesses are sufficient to explain the maturity data; therefore, high and low heat flow scenarios would converge. Increased uncertainty in the northern region of the model (outside the Hils Syncline), regarding paleo and present-day heat flows and eroded thicknesses, arises from limited calibration data. Consequently, present-day temperature calibration, in the northern area, was conducted relying on the latest German Heat Flow Database (Fuchs et al. 2022), yet present-day heat flow along this domain exhibits elevated values compared to earlier studies (e.g., Scheck-Wenderoth et al. 2014) possibly caused by hydrothermal activity.

Fig. 14

1D burial and thermal history plots and basal heat flow evolution at the location of Haddessen (see Fig. 2 for the location of the borehole). The white lines indicate the temperature isotherms from 100 to 500 °C. The red lines represent the Lower Jurassic unit (top and base) that displays the highly variable thermal maturity. Four different scenarios are illustrated. a Represent the best-fit model that assumes deep burial followed by strong uplift and erosion as the cause of high maturities. b Alternative scenario that differs from Scenario 1 in higher heat flow values during the extensional event (Late Jurassic–Early Cretaceous). c Deep-seated mafic intrusion scenario. d Basal heat flow trends determined for the tested scenarios. The abbreviations on the right side of the burial plots indicate the stratigraphic units and the basement. From top to base, they are labeled as follows: MJ Middle Jurassic, LJ Lower Jurassic, BS Upper, Middle and Lower Buntsandstein, Z Zechstein, PPC Pre-Permian Clastics and B Basement

Scenario 1 (Fig. 14a) represents the best-fit model described in the previous “Results” chapter and is consistent with the assumption of deep burial until the Late Cretaceous followed by strong uplift and erosion as the cause of high maturities (Schmitz and Wenzlow 1990; Petmecky et al. 1999; Senglaub et al. 2005, 2006; Adriasola-Muñoz et al. 2007; Brink 2013). The calibration of this model was achieved by applying a heat flow peak of 85 mW/m² during the Late Jurassic extension event (Figs. 9 and 14d), and a thickness of eroded sedimentary rocks of 3300 m during the Late Cretaceous inversion phase (Fig. 7). These values of eroded overburden are very similar and follow the same trend as the previously published values (e.g., Senglaub et al. 2006; Bruns et al. 2013, 2014).

An alternative second scenario (Fig. 14b) was tested by assuming higher heat flow values during the extensional event in the Late Jurassic–Early Cretaceous with a peak of maximum heat flow during this period reaching 95 mW/m² in the northwesternmost domain (Fig. 14d). In this scenario, the thermal evolution of the Hils Syncline is similarly controlled by the maximum burial depth reached in the latest Early Cretaceous in comparison to Scenario 1, but with a higher heat flow during Jurassic to Cretaceous extension. The increment in the heat flow peak leads to a higher modeled maximum temperature (approximately 180 °C) experienced by the Lower Jurassic units during the period of maximum burial (Fig. 14b). The thickness of eroded overburden required for proper calibration is lower when compared with Scenario 1 and with the literature, i.e., up to 2900 m of sediment was eroded at the location of Haddessen.

The magnitude of the modeled erosion during the Late Cretaceous inversion in Scenarios 1 and 2 implies the availability of sufficient accommodation space nearby. According to Bruns et al. (2014), the calculated erosion is moderate and consistent with the compressional tectonic regime dominating the study area (Kley et al. 2008; Voigt et al. 2008; Kley and Voigt 2008) and the humid climate conditions during this time (Stollhofen et al. 2008). Increased sedimentation rates in adjacent areas (especially Pompeckj Swell (PS) northern of LSB) during and after erosion in the LSB as well as a shift of the basin axis from south (LSB) to north (PS) are observed (Schwarzer and Littke 2007). These findings indicate a regional transport direction in an N–NW direction during the Late Cretaceous (Bruns et al. 2014), which can also explain the generally greater thickness of the Late Cretaceous/Cenozoic sedimentary units towards the north and the lack of Cenozoic sequences within the Hils Syncline (Fig. 4).

The last scenario (Fig. 14c) was modeled following the premise of a deep-seated mafic intrusion (i.e., the Bramsche and Vlotho massifs, located northwest of the Hils Syncline), emplaced during Late Cretaceous times (Bartenstein et al. 1971; Teichmüller et al. 1979, 1984). The influence of a possible magmatic event on the maturity pattern during the late Early Cretaceous was modeled by locally increasing heat flow values up to 340 mW/m² for a short period of ca. 1 Ma (Fig. 14d) in the northwesternmost domain where the highest thermal maturity has been measured (Haddessen; Fig. 2). In order to accurately calibrate the model, and based on the assigned heat flow trend (Fig. 14d), constant erosion was assumed to have taken place during the Late Cretaceous inversion across the Hils Syncline area. Thus maximum burial depths were low and the same in the entire modeled area for each stratigraphic unit (e.g., the Lower Jurassic unit was buried up to 1600 m deep in the north as well as in the south; Fig. 6). Therefore, Scenario 3 would imply that the eroded amount or differences in the maximum burial depth were not controlling factors in the variation of the thermal maturity. According to Scenario 3, maximum temperatures reached in the Lower Jurassic rocks during and shortly after the magmatic event reached about 200 °C in the north at the location of well Haddessen.

As discussed above, the increasing maturities in the southern LSB and the northern Hils Syncline are clearly the result of high temperatures, but the reason for these high temperatures has been a matter of debate. They may be due to deep Jurassic and Early Cretaceous burial followed by inversion and erosion or due to extremely high heat flows in the context of early Late Cretaceous magmatism (Bruns et al. 2013). Plutonic magmatism was proposed based on magnetic and gravimetric anomalies (e.g., Hahn and Kind 1971; Bachmann and Grosse 1989; “Bramsche Massif”, “Vlotho Massif”). However, magnetotelluric evidence along the Bramsche anomaly shows high electric impedance and low electric conductivities that indicate the absence of a plutonic intrusion in the area (Hoffmann et al. 2008; Brink 2010, 2012). In addition, seismic velocities and high electric conductivities were measured nearby the Bramsche anomaly and are interpreted as Lower Carboniferous/Upper Devonian carbonates rather than a plutonic body (Hoffmann et al. 2008; Doornenbal and Stevenson 2010), below which, at greater depths, an eclogite body of high density and magnetic susceptibility might cause the gravimetric and magnetic anomaly (Brink 2013). Furthermore, intrusive rocks in the western LSB (German–Dutch border) are dated as Permian (van Bergen and Sissingh 2007); a thick Permian volcanic succession is also found in the north and east of the LSB (Maystrenko et al. 2010), while no Cretaceous intrusive or extrusive rocks have been discovered within the Hils Syncline or southern LSB. It is obvious that the high thermal maturities of Mesozoic strata in the Hils Syncline are not related to the Permian igneous activity in the area, which pre-dates the sedimentation and burial.

The presence of hydrothermal fluids ascending through a fault system combined with deep burial during the Cretaceous has also been proposed as the responsible factor for the increasing maturities in the southern LSB that involves the Hils Syncline (Petmecky et al. 1999). This hypothesis is supported by the hydrothermal mineralization in the southern LSB (Stadler 1971; Klassen 1984), hot springs at different locations along the present-day southern LSB (Petmecky et al. 1999), abundant fault systems in the study area (for reference, use the TUNB3D model ref), and zircon fission track data that provide evidence for widespread (hydro)thermal activity during the Late Jurassic–Early Cretaceous (Senglaub et al. 2005). Such hydrothermal systems generally occur in extensional settings (Chaussard and Amelung 2014), which was the tectonic situation of the LSB during the Late Jurassic–Early Cretaceous (e.g., Senglaub et al. 2005; Kley et al. 2008). However, it is assumed that hydrothermal heat transport along fault systems has only very limited lateral effects on nearby rocks, being insufficient to cause the observed high maturities basin-wide in the LSB (Bruns et al. 2013; Cherubini et al. 2014). Even though these hot waters might cause high thermal maturity locally in the northwestern region of the Hils Syncline (e.g., at the location of Haddessen), they would likely not affect the low-permeable shaly Jurassic rocks. Based on the above arguments, it is proposed that deep burial due to rapid sedimentation during the Late Jurassic–Early Cretaceous followed by inversion and erosion is the ultimate cause of high burial temperatures in Mesozoic rocks, and the different magnitudes of burial is the main factor causing the variable thermal maturity that characterizes the Hils Syncline.

Correlating petrophysical properties with geodynamic evolution

The relationship between geodynamic evolution and petrophysical properties is especially significant in the case of Jurassic claystone within the Hils Syncline at a local level and throughout the LSB region on a larger scale due to the extensive burial and subsequent uplift experienced by this claystone. The lack of consideration of poroelastic responses due to pressure relief (e.g., uplift) is a major drawback of the present PetroMod analysis. Mechanical compaction is calculated as an irreversible process, and after effective stress reduction (uplift and erosion), porosities remain unchanged (Fig. 11), resulting in less pressure release during erosion than pressure build-up during burial. Elastic rebound is a known process that facilitates partial porosity recovery (e.g., Hamilton 1976; Neuzil and Pollock 1982) and, depending on the setting, (micro)fracturing during uplift may further enhance permeability (Urai and Wong 1994; Urai 1995). However, porosity and permeability data from Gaus et al. (2022) indicate that the matrix properties measured on uplifted organic-lean mudstones from the Hils area are indeed likely representative of maximum burial conditions. Therefore, alternative models incorporating the porosity data from Gaus et al. (2022) were built as a means to cross-check maximum burial depths modeled from vitrinite reflectance and displayed for well Haddessen in Fig. 15: (1) pure calibration of maximum burial against porosity data assigning the same heat flow as in chapter “Results”—Basin modeling and (2) combined calibration of maximum burial depth against porosity and vitrinite reflectance adjusting the heat flow. Porosity-based calibration in (1) and (2) was achieved in PetroMod by utilization of Athy’s law:

Fig. 15

1D burial and thermal history of an alternative model at the location of Haddessen. The modeled vitrinite reflectance is calculated using Easy %R_oDL Burnham et al. (2016). The white lines indicate the temperature isotherms from 100 to 300 °C. The red lines represent the Lower Jurassic unit (top and base) that displays the highly variable thermal maturity (Fig. 2). The abbreviations on the right side of the plots indicate the stratigraphic units. From top to base, they are labeled as follows: MJ Middle Jurassic, LJ Lower Jurassic, Z Zechstein and PPC Pre-Permian Clastics. a Model calibrated based on the measured porosity in Liassic samples (Pliensbachian–Lower Jurassic). The dashed lines in the modeled vitrinite reflectance and porosity correspond to this alternative model while the continuous lines represent the 3D model further described in the Results chapter. b Model calibrated based on both the vitrinite reflectance and porosity measurements as it is observed in the plots to the right. It is noteworthy that the temperature increased during the Late Jurassic–Early Cretaceous period as a result of the high heat flow peak that was assigned (136 mW/m²)

$$\varnothing \left(Z\right)={\varnothing }_{1}+ {(\varnothing }_{0}{- {\varnothing }_{1})e}^{-kZ}$$

(3)

Here, ∅ (Z) is the present (experimentally determined) porosity, ∅0 is the initial porosity (70%), ∅1 is the minimum porosity (1%), k is Athy’s compaction parameter (0.83) and Z corresponds to the depth. For simplicity, a homogeneous lithology (organic-lean claystone) was assigned.

(1) Pure calibration of maximum burial against porosity with fixed heat flow at the Haddessen location (Fig. 15a) yields lower maximum paleo-depths (3200 m) and -temperatures (147 °C) and, therefore, lower present-day modeled vitrinite reflectance (0.81%VR_r) when compared to the results obtained in the 3D model from this study (Fig. 6, 1.45% VR_r). (2) Combined calibration of maximum burial depth against porosity and vitrinite reflectance adjusting the heat flow yields a heat flow peak of 136 mW/m² during the Late Jurassic–Early Cretaceous period (Fig. 15b). Interestingly, for all other locations analyzed in this study, results of the calibrated 3D regional model show an appropriate fit with both vitrinite reflectance, as well as porosity measurements and no variations in the thermal field (heat flow trend) or erosion thickness, are required. From both models, it is clear that in order to match porosity and vitrinite reflectance increased heat flow is necessary. The modeled heat flow peak (136 mW/m²), however, exceeds the heat flow ranges for extensional regimes (Allen and Allen 2013), which is the most likely tectonic situation of the LSB during the Late Jurassic–Early Cretaceous (e.g., Senglaub et al. 2005; Adriasola Munoz et al. 2007; Bruns et al. 2013) and is only regarded suitable for active ocean ridges and volcanoes (Allen and Allen 2013). However, the absence of Cretaceous intrusive or extrusive rocks in the area negates this possibility (see chapter “Discussion”—“Burial and thermal history”).

Given the imperfect correlation between porosity and vitrinite reflectance for the Haddessen samples, where the Posidonia shale has entered the gas window stage (Leythaeuser et al. 1988; Littke et al. 1988, 1991a, b; Rullkötter et al. 1988; Sundararaman et al. 1993; Bernard et al. 2012), one plausible explanation could involve the development of localized overpressures in the Haddessen region. The generation of gas typically entails volume expansion, which, if the permeability of surrounding rock is suitably low, can result in the formation of overpressure (Osborne and Swarbrick 1997). Considering the very low matrix permeability in the Posidonia Shale and the under and overlying rocks before uplift (Mohnhoff et al. 2016; Gaus et al. 2022), it is reasonable to infer that overpressure within the Posidonia Shale unit occurred during its maximum burial phases, which aligns with the timeframe for gas generation in the study area (Bruns et al. 2014). The relatively high porosities measured in the Posidonia Shale at the Haddessen location (ranging from ca. 10% and 15%) suggest that overpressure, in tandem with secondary cracking (Bernard et al. 2012), might be plausible mechanisms at play. The underlying Pliensbachian and overlying Doggerian units also exhibit relatively high porosities of about 8% (Gaus et al. 2022), at the Haddessen location. Therefore, it is conceivable that overpressure was then conveyed through neighbor-source migration, with microfractures potentially serving as conduits for overpressure propagation (Swarbrick and Osborne 1998; Finkbeiner et al. 2001; Tingay et al. 2007; Fan et al. 2021) into Pliensbachian and Doggerian units. Alternatively, allo-source vertical overpressure transfer could have transpired due to the interconnection between the overpressured unit and neighboring formations facilitated by fault openings (Tingay et al. 2007; Gao and Flemings 2017). Apart from gas generation-related processes, overpressure possibly results due to rapid sedimentation and fluid retention (undercompaction) or clay diagenetic processes, e.g., conversion of smectite into illite. As clay-rich sediments are buried, their porosity reduction can be accompanied by the transformation of smectite, by dehydration, into illite, which releases water into the pore space, resulting in the increase of pore pressure and thereby transferring this pressure to the pore fluid, forming overpressure (Dutta 2002; Lahann 2002; Lahann and Swarbrick 2011). Undoubtedly, the processes outlined above remain speculative in nature. Nonetheless, considering the implausibly high heat flows necessary to reconcile porosity and vitrinite reflectance, coupled with the burial history and a strong likelihood of localized overpressure generation within the Posidonia Shale, this hypothesis emerges as the most plausible explanation. Finally, the overpressure generated would tend to dissipate over time via fluid leakage (Borge 2002) and/or via active faulting (Engelder 1993; Swarbrick and Osborne 1998). The strong erosion occurring during the Late Cretaceous definitely permitted pressures to return to equilibrium condition (normal hydrostatic pressure). However, rocks that have experienced overpressure in the past are more brittle (Nygård et al. 2006) and most likely affected by overpressure-related fracturing, which represents a risk for the feasibility, long-term stability and safety storage capacity of the repository.

Comparison of erosion thickness estimates

The estimation of the magnitude of uplift and erosion is an important aspect with several implications in economic and scientific fields (Iyer et al. 2021). For instance, burial and exhumation histories can affect diagenetic and compaction processes (Iyer et al. 2021) and, therefore, the petrophysical properties of the rocks, ultimately having implications for the suitability of reservoir rocks, e.g., for hydrocarbon accumulations (Ohm et al. 2008) or as nuclear waste hosts (Henriksen et al. 2011; Baig et al. 2016). By accurately determining the erosion thickness associated with the Late Cretaceous inversion, it becomes possible to establish precise maximum paleo-burial depths. This, in turn, may lead to a better understanding of the evolution of the basal heat flow—the main controlling mechanism for the thermal maturity within the study area.

When comparing the eroded thicknesses calculated using vitrinite reflectance data (3D model approach) and sonic velocities of the Pliensbachian (Lower Jurassic) interval, a strong similarity is observed across nearly all of the analyzed boreholes (Table 3). The increasing northwestward trend of eroded thickness is also validated by the different methods, although the estimates show more variability as erosion increases towards the north (Table 3). Accordingly, the overall sonic velocity is increasing from Wenzen towards the Harderode region. However, this trend is not evident in the northernmost location (Haddessen/Bensen), where a reduction in sonic velocity is observed (regional-wide; Figs. 12 and 13). Consequently, a large mismatch is obtained for erosion estimates in the Haddessen/Bensen region; while the vitrinite reflectance-based estimate suggests the highest amount of erosion at Bensen (about 3300 m), the sonic velocity-based estimate indicates smaller erosion (ca. 1800 m; Table 3). Estimating erosion based on porosity according to Athy’s Law gives erosion estimates between the other two methods (2430 m) at Bensen (Table 3). One explanation could be that the analyzed Pliensbachian claystone was overpressured (see chapter “Discussion”—“Correlating physical properties with geodynamic evolution”), which would have preserved the porosity and resulted in a reduction of the sonic velocity (Mouchet and Mitchell 1989; Bowers 1995; Tingay et al. 2013). This, in turn, would lead to an underestimation of erosion thickness (Iyer et al. 2021; Lasabuda et al. 2021).

Table 3 Comparison of erosion thickness among the vitrinite reflectance-based estimate (3D model approach), the straightforward sonic velocity analysis, and Athy’s law porosity-based calculation

Conclusions

The results of this study provide new insights into the maturity variability, as well as thermal and burial evolution of the Hils Syncline and surrounding areas (southern LSB). These factors significantly influence mechanical compaction and petrophysical properties of the sedimentary rocks such as claystones. The implications of the latter are significant when it comes to exploring a nuclear waste repository in claystone formations.

The 3D-thermally calibrated numerical modeling reveals that high thermal maturities of Jurassic rocks in the Hils Syncline, determined by vitrinite reflectance measurements, are neither influenced by a Late Cretaceous magmatic intrusion nor by the Permian igneous activity in the area, which preceded the sedimentation and burial of the Mesozoic units. Rather, the thermal maturities in the study area are the result of the maximum burial depths and temperatures reached during the latest Early Cretaceous, followed by a Late Cretaceous inversion phase that caused severe erosion. Differences in the burial magnitudes increasing northwestward accompanied by higher maximum temperatures are the ultimate cause for the variable maturation levels.

The analyzed sonic velocity data of the Pliensbachian organic-lean claystone show an increase in the velocity from the boreholes located to the southeast towards the northwest. This implies that the sediments experienced deeper burial towards the northwest resulting in higher mechanical loading and, therefore, higher loss of pore space. These observations support the modeling results regarding the burial evolution in the Hils Syncline. Moreover, the erosion estimates resulting from the sonic velocity approach fit very well with those obtained from the vitrinite reflectance analysis performed in the study area—except for one well in the northwestern area. The same trend is obtained from previously published porosity and permeability data if compared to those derived from the model—a very good fit for the Hils area with decreasing values to the northwest due to deeper burial, but a mismatch for the northwesternmost boreholes Haddessen/Bensen (which are already outside of the Hils Syncline sensu strictu).

The association of the mismatch between porosity and vitrinite reflectance for the Haddessen samples, combined with the sonic velocity response, provides compelling evidence for local overpressure in the Haddessen/Bensen region, possibly caused by fluid expansion due to gas generation within the Posidonia Shale which propagated into Pliensbachian and Doggerian units.

The Liassic–Doggerian mudstones, potential host rocks for nuclear waste, exhibit favorable petrophysical characteristics—specifically, low porosity (< 20%) and permeability (< 0.01 mD)—in the studied region, even in the area of maximum thermal maturity (Haddessen/Bensen), where the influence of overpressure is hypothesized. However, their absence at suitable depths (300–1500 m) in the Hils Syncline raises concerns, as does the potential brittleness due to overpressure, posing a risk to long-term nuclear waste storage. It is suggested to reassess the hypothesis of overpressure within the target formations, particularly in the northern LSB, where a deeper basin might allow the host rock formations to reach optimal depths for an effective nuclear waste repository.

Additional porosity and permeability measurements should be performed in the future since the analyzed results in this paper are based on laboratory-measured properties from Gaus et al. (2022), which were measured on a few samples from boreholes drilled 35 years ago. A recently ongoing drilling campaign will provide new and representative petrophysical data, especially on the Pliensbachian interval.