Prediction of Reservoir Parameters of Cambrian Sandstones Using Petrophysical Modelling—Geothermal Potential Study of Polish Mainland Part of the Baltic Basin

Abstract

Important factors controlling the effective utilization of geothermal energy are favorable reservoir properties of rock formations, which determine both the availability and the transfer opportunities of reservoir fluids. Hence, crucial to the successful utilization of a given reservoir is the preliminary recognition of distribution of reservoir parameters as it enables the researchers to select the prospective areas for localization of future geothermal installations and to decide on their characters. The objectives of this paper are analyses and discussion of the properties of quartz sandstones buried down to a depth interval from about 3000 to under 5000 m below surface. These sandstones belong to Ediacaran–Lowery Cambrian Łeba, Kluki and Żarnowiec formations. The source data from the Słupsk IG-1 provided the basis for 1D reconstruction of burial depth and paleothermal conditions as well as enabled the authors to validate of the results of 2D models. Then, porosity distribution within the reservoir formation was determined using the modelings of both the mechanical and chemical compactions along the 70 km-long B’-B part of the A’-A cross-section Bornholm-Słupsk IG-1 well. The results confirmed the low porosities and permeabilities as well as high temperatures of the analyzed rock formations in the Słupsk IG-1 well area. Towards the coast of the Baltic Sea, the porosity increases to more than 5%, while the temperature decreases, but is still relatively high, at about 130 °C. This suggests the application of an enhanced geothermal system or hot dry rocks system as principal methods for using geothermal energy.

1. Introduction

The transformation of energy systems results from the pressure for the elimination of environmental pollution and hazards and reduction in anthropogenic factors stimulating climate change. Such a transformation forces the search for alternative, environmentally neutral energy sources such as the Earth’s internal heat. Exhaustion of hydrocarbon resources provides an opportunity to utilize the closed petroleum wells for supplying geothermal energy to the new installations. Important controlling factors of effective utilization of geothermal energy are favorable reservoir properties of rock formations, which determine both the availability and the transfer potential of reservoir fluids. Hence, crucial to the successful utilization of a given reservoir is the preliminary recognition of distribution of reservoir parameters as it enables the researchers to select prospective areas for localization of future geothermal installations and to decide on their characters [1,2,3,4]. Taking into account the low values for reservoir parameters of the studied rock formations, permeability and effective porosity, not only the “classic”, low-enthalpy geothermal (LEG) installations must be considered, but also the other opportunities: enhanced geothermal systems (EGSs) and even hot dry rocks systems (HDRs).

This paper provides an analysis of geothermal energy utilization in northern Poland. The study area is located in the marginal part of the East European Craton (EEC), in the Polish part of the Baltic Basin (BB). In that area, an extensive hydrocarbon exploration has been carried out for decades in both the Ediacaran–Lower Cambrian sandstones and in the Upper Ordovician–Lower Silurian shales, the latter hopefully hosting the unconventional hydrocarbon accumulations. Numerous exploration projects have provided extensive information, which can be useful in geothermal exploration. The main goal of this publication is to present the results of recognition of geothermal conditions within a particular geological formation using the information acquired from the hydrocarbon exploration projects, thus taking advantage of the hydrocarbon exploration methodologies. From the point of view of evaluation of geothermal conditions, the recognition of temperature distribution within the target rock formation is crucial, as are the reservoir and rock-mechanics parameters responsible for efficient fracking operations. Preliminary analysis of geological data indicates that the Ediacaran–Lower Cambrian sandstones are a promising potential target of geothermal exploration. This quartz sandstones were deposited in the Polish part of the Baltic Basin (BB) and were buried down to the depth interval from about 3000 to more than 5000 m below sea level (b.s.l.). The sandstones belong to the Ediacaran–Lower Cambrian Łeba, Kluki and Żarnowiec sedimentary formations. The analyses are based upon data acquired from the Słupsk IG-1 well, completed in 1974. The wellsite is located in the onshore part of the Darłowo Block, at the line of A’-A cross-section Bornholm-Słupsk IG-1 (Figure 1, Figure 2 and Figure 3). The well truncates the crystalline basement at 5078 m depth and penetrates the Ediacaran–Lower Cambrian quartz sandstones. The bottom portion of the Żarnowiec Formation, located below 5000 m depth, where the sandstones degrade into conglomerates, provides the lower boundary of the studied depth interval.

The porosity and permeability of Cambrian sandstones was recently observed in the study area and very intensive mechanical and chemical compaction values were found for Cambrian sandstones, resulting from both the mechanical and the chemical compactions in the study area. The mechanical compaction was an effect of the load of thick overburden and caused closer packing of grains, forming the rock fabric. The chemical compaction of Cambrian sandstones was represented by quartz cementation. With the increasing burial depth and temperature, more silica was dissolved and reprecipitated onto the free surfaces of closely packed quartz grains. Hence, compaction highly reduced the primary porosity and permeability values to their standard ones.

In the research, the authors used both the 1D and 2D paleothermal modelings together with the modelings of porosity changes. The 1D modelings were performed using the Słupsk IG-1 well (Figure 2 and Figure 3). The 2D modeling was carried out for the B’-B part of the A’-A cross-section (Figure 2 and Figure 3). Crucial to correct modeling of porosity changes in geologic time is the reconstruction of paleothermal conditions. Modeling of diagenetic processes requires not only the knowledge of recent geometry of sediment formations and temperature values, but also the evaluation of subsidence and paleothermal conditions through geologic time. Data acquired from the wells (including the Słupsk IG-1 well) were applied to 1D modelings, from which the estimations were derived, e.g., the mechanical compaction effect, in accordance with Athy’s model [5,6], and the chemical compaction effect after Walderhaug’s model [7]. Both the 2D paleothermal and the 2D porosity distribution models were constructed for the 70 km-long B’-B cross-section, which is a part of the longer A’-A cross-section line (Figure 1, Figure 2 and Figure 3). Both the 1D and 2D modelings of mechanical and chemical compactions were based upon the multistage reconstruction of sediment burial, which included the deciphering of variable density of heat paleoflow determined for the Słupsk IG-1 well. The 1D and 2D models generated with the PetroMod 2016 software, combined with the laboratory measurements of effective porosity and permeability of sandstones, were then applied to the analysis of petrophysical parameters. The 1D reconstruction of paleothermal conditions in relation to burial depth was based upon modeling of vitrinite reflectance curve, temperature and heat flow values. Geothermal conditions of the Polish segment of the BB were compared with those known from the Lithuanian onshore segment of the BB, where a geothermal power plant utilizes the energy of both the Cambrian and the Early Devonian geothermal aquifers [8]. Confirmed low porosity and permeability values as well as high temperatures of sandstone formations enabled the authors to categorize the Cambrian sandstones as a natural geothermal horizon of properties suitable for application of the EGS and/or the HDR technologies for energy tapping.

2. Geological and Geothermal Settings

2.1. Geological Setting

The BB is located in the southwestern part of the East European Craton (EEC) and is bordered by the Teisseyre-Tornquist Zone (TTZ), which, in Poland, coincides with the northeastern boundary of the Trans-European Suture Zone (TESZ) (Figure 1). The TESZ is a major geological boundary in Europe, which extends from the North Sea to the Black Sea, over 2000 km. It separates the mobile Phanerozoic terranes from the stable Precambrian EEC [9].

Figure 1. Location of the Baltic Basin in the East European Craton (after [10]). Dashed purple line—boundary of Baltic Basin, toothed black lines—boundaries of Caledonide, Variscan and Alpine orogens, TTZ—Teisseyre-Tornquist Zone [11], STZ—Sorgenfrei-Tornquist Zone [11]. Red dashed rectangle marks the study area.

Figure 1. Location of the Baltic Basin in the East European Craton (after [10]). Dashed purple line—boundary of Baltic Basin, toothed black lines—boundaries of Caledonide, Variscan and Alpine orogens, TTZ—Teisseyre-Tornquist Zone [11], STZ—Sorgenfrei-Tornquist Zone [11]. Red dashed rectangle marks the study area.

The crystalline basement of the BB comprises Paleo-Mesoproterozoic metaigneous and metasedimentary units originating from the Svekofennian cycle. In its geological history, the EEC was subjected to numerous volcanic episodes. Between 1.55 and 1.45 Ga, the EEC basement was affected not only by intrusions of granitic magma but also by trap volcanism including alkali basalts, trachyandesites and quartz tholeites eruptions, which took place between the Cryogenian and the Ediacaran (650–550 Ma). The last magmatic episode brought local syenite intrusions of the Early Carboniferous age (~350 Ma) [12].

The BB crystalline basement dips from the northeast to the southwest (Figure 2). On the southern side of the EEC margin, the crystalline basement is buried down to about 10,000 m depth [13]. At the EEC margin, the top surface of crystalline basement occurs at about 8000 m depth, as revealed by refraction seismics [14]. In the Polish part of the BB, the greatest depth of Precambrian basement occurs in Słupsk IG-1 well (5078 m).

Figure 2. Geological-tectonic map of sub-Permian surface in the southwestern part of the Baltic Basin (after [15], fragment). Thick green line indicates modeled A’-A cross-section. Thick, black dashed line indicates range of 2D modelling (B’-B).

Figure 2. Geological-tectonic map of sub-Permian surface in the southwestern part of the Baltic Basin (after [15], fragment). Thick green line indicates modeled A’-A cross-section. Thick, black dashed line indicates range of 2D modelling (B’-B).

The crystalline basement is covered by structural units separated by erosional unconformities of regional extension. Thickness of this overburden varies from 2000 m in the central part of the BB to 5000 m in its southwestern periphery located in Poland [16]. In the onshore, southwestern portion of the Polish part of the BB, the Caledonian structural unit includes Ediacaran sediments, which continuously grade into Cambrian, Ordovician and Silurian strata (Figure 3) [17,18,19].

In northern Poland and in the adjacent, offshore part of the Baltic Sea, the following tectonic stages were identified, corresponding to main depositional cycles [13,19]: (i) Ediacaran–Middle Cambrian, (ii) Middle/Late Cambrian–Tremadocian and (iii) Arenigian–Gedinian. Detailed analysis of Lower Paleozoic sediments based upon data from onshore drillings in northern Poland led to identification of depositional sequences separated by regional-scale erosional unconformities. These sequences correspond to the above-mentioned tectonic stages [18].

Figure 3. Geological and tectonic settings of A’-A cross-section Bornholm-Słupsk IG-1 well shown in Figure 2 (after [13]). Explanations: 1—unconformities, 2—stratigraphic boundaries, 3—faults, 4—Caledonian tectonic units, 5—Alpine tectonic units, PP—Paleoproterozoic, Є+NP3—Cambrian+Upper Neoproterozoic, O—Ordovician, S—Silurian, P—Permian, T—Triassic, Cr—Cretaceous, Cz—Cenozoic. Light green rectangle—part of A’-A cross-section included in the 2D modeling (B’-B cross-section).

Figure 3. Geological and tectonic settings of A’-A cross-section Bornholm-Słupsk IG-1 well shown in Figure 2 (after [13]). Explanations: 1—unconformities, 2—stratigraphic boundaries, 3—faults, 4—Caledonian tectonic units, 5—Alpine tectonic units, PP—Paleoproterozoic, Є+NP3—Cambrian+Upper Neoproterozoic, O—Ordovician, S—Silurian, P—Permian, T—Triassic, Cr—Cretaceous, Cz—Cenozoic. Light green rectangle—part of A’-A cross-section included in the 2D modeling (B’-B cross-section).

The Ediacaran–Middle Cambrian structural stage includes sediments of thickness of about 1000 m. The Middle/Late Cambrian–Tremadocian stage comprises over 20 m-thick sediments from the uppermost Middle Cambrian to the Lower Tremadocian. The Arenigian–Gedinian tectonic stage referred to depositional sequences distinguished by Jaworowski [18] can be divided into the two substages: (i) Ordovician, with sediments ranging in age from the Arenigian to the Ashgillian (the latter without the uppermost part), of total thickness up to 140 m, and (ii) Silurian, with very thick (up to 3700 m) sediments of ages from the Llandovery (without the lowermost part) to the Pridolian (without the uppermost part). However, according to Areń and Tomczyk [20], the Caledonian structural unit comprises only the Ordovician (starting from the Arenigian) and the Silurian sediments.

Sequence I was deposited during the extension-related subsidence, which accompanied the break-up of the Precambrian Rodinia supercontinent. It embraced sediments ranging from the Ediacaran to the mid Middle Cambrian, of total thickness of about 1000 m, which directly overlies the crystalline basement surface (Figure 4). Lithology included continental and continental-marine sediments (Ediacaran–Lowermost Cambrian) developed as sandstones and conglomerates deposited in alluvial fans, braided plains, braided and fan deltas grading into Lower and Middle Cambrian marine deposits: transgressive/regressive, tidal and storm-surge sandstones intercalate with shelf sandstone-mudstone heteroliths, which enclose storm deposits and infillings of tidal channels.

Sequence II comprises sediments from the uppermost Middle Cambrian to the Lower Tremadocian, which are more than 20 m thick. These deposits were, transgressive conglomerates and sandstones followed by black shales deposited under anaerobic/dysaerobic conditions of a deep shelf (equivalents of Alun Shales from Scandinavia). Unconformity found at the top of this succession was related to Late Tremadocian erosion.

Figure 4. Chronostratigraphic column of Neoproterozoic–Lower Paleozoic sedimentary cover with indicated depositional sequences from the western part of the East European Craton (after [21,22]; simplified).

Figure 4. Chronostratigraphic column of Neoproterozoic–Lower Paleozoic sedimentary cover with indicated depositional sequences from the western part of the East European Craton (after [21,22]; simplified).

Sequence III comprises sediments from the Arenigian to the Ashgillian (without the uppermost part), with a total thickness of up to 140 m. Two subsequences were evident, separated by unconformity related to submarine erosion in the Middle and the Late Llanvirnian (without the uppermost part). The lower subsequence IIIa comprised both the Arenigian and the Lower Llanvirnian sediments: basal conglomerate, shelf mudstones with glauconite and marly limestones deposited onto a deep carbonate ramp and at the boundary with the shallow ramp. Maximum total thickness of the subsequence IIIa reached up to 50 m. The upper subsequence IIIb included claystones deposited in an anaerobic shelf environment (uppermost Llanvirnian–Caradocian) and alternating mudstones/marls (Ashgillian) sedimented in the border zone between the shelf and the deep carbonate ramp. Maximum total thickness of this subsequence was 90 m. The stratigraphic gap, lasting from the latest Ashgillian to the earliest Llandovery, was related to erosional processes and constituted the upper boundary of Sequence III.

The very thick (over 3700 m) Sequence IV was deposited as an accretionary prism in a foreland basin of the Baltica continental margin which overthrust to the southwest [18]. These were sediments from the Llandovery (without its uppermost part) up to the Pridolian (without the lowermost part), developed as black siltstones deposited in an anaerobic shelf and grey siltstones (with intercalations and lenses of limestones) deposited in a dysaerobic shelf environment. The grey shelf mudstones, sedimented in the Lludlovian, included numerous intercalations of siltstones representing fine-clastic debrites and distal turbidites. The upper boundary of Sequence IV was an erosional surface produced by long-lasting, multistage, subaerial erosion. Sequence IV sediments were covered by Permian and/or Mesozoic strata.

For the majority of the Devonian and Carboniferous periods, the BB was a site of siliciclastic-evaporitic sedimentation. Carbonates were deposited under the open-sea conditions when the BB evolved from the endorheic to the open-sea basin, open to the southwest. The Devonian–Carboniferous deposition ended in the Lower Carboniferous (Tournaisian) when carbonates interbedded with gypsum, claystones and sandstones and were sedimented. Between the end of the Early Tournaisian until the Early Permian, the BB was subjected to erosion, which was responsible for complete removal of Devonian–Carboniferous strata in the vicinity of Słupsk town, whereas in northwestern Lithuania and the southwestern Latvia, the thickness of uneroded sediments reached 1000 m [23].

In the southern, Polish part of the BB, stratigraphic sequence of the Alpine stage was much reduced due to both the initial ranges of consecutive sedimentary basins and the long-lasting denudation, which took place in the Latest Zechstein, the Late Triassic, the Early Jurassic and the Early Cretaceous [24]. The Alpine tectonic stage unconformably covered the eroded surfaces of older stages; however, a significant part of the BB was devoid of both the Permian–Mesozoic and the Cenozoic sediments. This stage comprised the Early Alpine substage with sediments from the Upper Permian to the Cretaceous and the Late Alpine substage represented by Paleogene and Neogene strata developed in the marginal part of the BB (Figure 3).

2.2. Geothermal Setting

The Polish part of the BB generally shows low heat flow values, geothermal gradient and related distribution of temperatures within the rock formations. Heat flow value in the near-surface zone is anomalously low (30–50 mW/m²) in both the central and the southern parts of the Peribaltic Syneclise (PS), in relation to adjacent regions, whereas in both the northwestern and the northeastern parts of the PS, local zones of increased heat flow values appear—up to 70–90 mW/m² [25]. The southwestern part of the PS is presumably influenced by the nearby TESZ, whereas the northwestern one, outside the state border, is supplied by radioactive heat sources in the crystalline basement [8].

Geothermal conditions in the Polish part of the PS were recognized thanks to the analysis and the thermal interpretation of drilling reports [25,26,27] as well as to the results of hydrocarbon exploration projects. Particularly important was the exploration of Ordovician and Silurian shale formations, which are still expected to be unconventional hydrocarbon accumulations of the “shale gas” type [28,29,30,31,32,33,34]. Despite low values of heat flow, the temperatures of potential geothermal horizon of Cambrian sandstones were relatively high and locally reached even up to 180 °C (Figure 2). Relatively high temperatures resulted from significant burial depth of these sediments and from specific distribution of thermal parameters of rock formations (thermal conductivity, heat capacity and location of radiogenic sources) [35]. Additional sources of radiogenic heat were Proterozoic gneisses and Ordovician–Silurian claystones, siltstones and mudstones. Generally, the temperatures increased to the southwest, in accordance with the increasing burial depth of crystalline basement, but were locally modified by basement morphology and anomalies of heat flow values (Figure 2) [25].

In Poland, geothermal energy has not been utilized up to date within the study area. However, geothermal installations supplied by the Cambrian sandstone reservoir were built in Lithuania where both the parameters of geothermal field and the distribution of reservoir parameters were more favorable due to shallower burial depths. Hydrothermal horizons in Lithuania were explored with the seismic survey supported by the analysis of data from producing oil and gas wells, and from geothermal exploration drillings completed in the late 1980s and early 1990s [36].

The Cambrian geothermal reservoir with a thickness varying from 60 to 80 m comprised quartz sandstones interbedded with rare claystones and shales. The results of petrophysical studies of 968 samples collected from deep wells penetrating the Cambrian sandstones revealed that their porosity and permeability decreased in the west, along with the increasing burial depth. Highest porosity values (25–30%) occurred at the eastern margin of the Lithuanian Basin. In its central part, porosity values were distinctly lower (15–18%), whereas in its western part, the lowest values were measured, comparable with those determined in the Polish part of the PS (2–4%). At depths of 400 to 1800 m, porosity changed from 25 to 16% and decreased to 5% below 1800 m depths [8]. Permeability values varied from 0.5 to 1 D at depth interval 1000 to 1800 m, whereas below 1800 m depth, they rapidly decreased down to 10⁻¹–10⁻⁵ D [36]. Decreasing values of porosity and permeability were related to late quartz cementation. Temperatures within the Cambrian geothermal aquifer varied from 15 °C in the east to 70–90 °C in the west, where it was buried at about 2000 m depth. In western-central Lithuania, the Lower Paleozoic succession hosted three hydrothermal reservoir horizons: Middle Cambrian, Lower Devonian and Middle Devonian. High values of geothermal parameters in these horizons resulted from an anomalously high basal heat flow (up to 90 mW/m²) carrying thermal energy from the asthenosphere. Apart from the basal heat flow, important thermal energy sources were also the Mesoproterozoic granitoid intrusions, which invaded the Paleoproterozoic granulites about 1.5 Ga ago. These rocks generated radiogenic heat at the level of 4–18 μW/m³ [8]. The highest heat flow anomalies related to magmatic intrusions in the crystalline basement were recognized in western Lithuania. These were the principal targets of geothermal research projects aimed to find locations suitable for HDR/EGS installations. The most prospective geothermal installation sites were those related to rock formations having temperatures of up to 150 °C at depth intervals of 4–5 km [8].

Geothermal conditions in the Polish part of the PS are less favorable in comparison with those in Lithuania. In Cambrian sandstones, geothermal conditions are controlled by burial depth and by development of geological processes influencing the distribution of reservoir parameters, particularly the burial history and the variability of the geothermal field. Typically, low porosity and permeability values of Cambrian sandstones imply rather poor perspectives for development of low-temperature geothermal installations. The reference site for this research is the Słupsk IG-1 well, where thermal measurements were completed after a 15-day-long well stabilization period [27]. Interpretation of temperature measurements in this well together with the estimations of heat flow values carried out independently by three studies showed remarkable differences: from 54 to 78 mW/m² [25,27,28,37]. These differences will be discussed below, including the results of author’s own thermal modelings.

3. Materials and Methods

The essential parts of the research project were the numerical paleothermal modelings as well as the 1D and 2D modelings of mechanical and chemical compaction processes, all carried out with the PetroMod software. The study target was the succession of Ediacaran–Lower Cambrian quartz sandstones belonging to the Łeba, Kluki and Żarnowiec formations. The reference dataset derived from the Słupsk IG-1 well provided a basis for methodological presumptions and for preliminary 1D modelings whereas 2D studies and analyses were carried outfor a 70 km-long B’-B fragment of the A’-A cross-section Bornholm-Słupsk IG-1 well (Figure 2 and Figure 3).

In the lithostratigraphic column of the Słupsk IG-1 well [38] (Figure 2 and Figure 3), the crystalline basement includes Paleoproterozoic gneisses (81.6 vol.%), which occur at depths from 5120 to 5078 m (

Table 1. Lithostratigraphic column of the Słupsk IG-1 well with thermal parameters to 1D and 2D modelings. Lithology after PertoMod library data.

Stratigraphy	Depth (m)	Thickness (m)	Lithology	Thermal Conductivity of Rock Fabric at 20 °C (W/mK)	Thermal Capacity at 20 °C (cal/gK)
Quaternary	0	47	Sand/silt/mud	2.39	0.209
Miocene	47	103	Siltstone	2.05	0.217
Oligocene	150	12	Siltstone	2.05	0.217
Upper Creataceous	162	400	Siltstone, marl, limestone, sandstone	2.33	0.210
Upper Triassic	562	78	Sandstone clay rich	3.35	0.206
Middle Triassic (Anisian)	640	13	Limestone	3.00	0.200
Lower Triassic	653	373	Siltstone, shale, sandstone	2.14	0.210
Upper Permian (Zechstein)	1026	68	Siltstone, shale, anhydryty, limestone	2.72	0.201
Lower Permian (Rotliegend)	1094	56	Sandstone	3.95	0.204
Middle Silurian (Lludlovian)	1150	2190	Siltstone	2.05	0.217
Lower Silurian (Wenlockian)	3340	1143	Siltstone, shale	1.87	0.213
Lower Silurian (Llandovery)	4483	7	Shale organic rich	1.25	0.215
Ordovician	4490	25	Shale, siltstone	1.87	0.213
Middle Cambrian	4515	97	Shale, siltstone	1.87	0.213
Lower Cambrian	4612	279	Sandstone clay rich	3.35	0.206
Ediacaran–Lower Cambrian	4891	186	Reservoir Formation (Sandstone quartz see text above)	4.70	0.208
Paleoproterozoic	5078	42	Gneiss	2.70	0.191
Crystalline basemenet	5120	-	-	-	-

). The Ediacaran–Lower Cambrian sedimentary succession comprises (from the bottom): (i) a sandstone sequence distinctly dominated by conglomerates in the bottom part, grading into (ii) sandstones intercalated by mudstone laminae rich in glauconite (Żarnowiec Formation) (depth interval from 5078 to 4891.5 m) (Figure 4). This succession corresponds to the Lower Cambrian Platysolenites antiquissimus horizon. Conglomerates encountered at a depth interval from 5000 to 5078 m were rejected from the modelings and analyses. Thickness of the sandstone sequence increased to the southwest and reached a maximum of 186.5 m in the Słupsk IG-1 well. Facies analysis indicates psammitic sediments intercalated with mudstones and heteroliths. The change from psammitic to psefitic fractions is related to the proximity of the source of continental clastics. The Żarnowiec Formation grades into (iii) the Kluki Formation built of sandstones interbedded with claystones (depth interval from 4891.5 to 4825 m). The upper portion of Lower Cambrian succession is occupied by the (iv) Łeba Formation, composed of sandstones, mudstones and claystones (depth interval from 4825 to 4716 m) [39,40].

Analysis of petrophysical parameters includes a depth interval from 4716 to 5000 m, i.e., the middle and the top parts of the Żarnowiec Formation as well as the Kluki Formation and the Łeba Formation (Figure 4). The Żarnowiec Formation sediments contain 72.28 wt.% of quartz, on average, whereas those of the Łeba and the Kluki formations contain average quartz values of 81.11 wt.% each. Additionally, the measured average quartz grains sizes are: 0.35 mm in the Żarnowiec Formation and 0.37 mm in the Łeba and the Kluki formations [41]. Thickness of Ediacaran–Cambrian succession is 563 m.

Significant stratigraphic gaps encountered in the column of the Słupsk IG-1 well, resulted mostly from the post-Caledonian (post-Silurian), post-Variscan and Cretaceous/Paeogene erosion episodes. Paleozoic succession observed in the Słupsk IG-1 well comprised Ediacaran–Cambrian, Ordovician, Silurian and Permian sediments, whereas stratigraphic gaps were the effects of post-Caledonian and post-Variscan erosion. The oldest gap included Upper Cambrian, Lower and Middle Ordovician strata. Along the sequence, a thin layer of the Sasino Claystones occurred together with the Prabuty Marls/Claystones, with a thickness of less than 25 m (Figure 4). These sediments represent a strongly reduced Ordovician succession buried down to depth interval 4515 to 4490 m and strongly reduced in thickness.

Stratigraphic column of Silurian sediments starts with the Pasłęk Formation (Llandovery), which is affected by a fault zone and, hence, its thickness is reduced to 7 m. It is overlain by the Kociewie Formation (claystones, mudstones) and the Puck Formation (calcareous claystones) (Figure 4). Silurian strata reach a maximum thickness over 3000 m and occur at a depth interval from 4490 to 1150 m, as found in the Słupsk IG-1 well. The next stratigraphic gap comprises Silurian strata—more precisely, the younger part of the Upper Ludfordian, the Pridolian as well as the Lower Devonian sediments. This is an erosional surface related to the post-Silurian epigenetic erosion, which removed Lower Lochkovian, Pragian and Emsian sediments. An approximately 350 m-thick gap includes 100 m of the Lludlovian, 200 m of the Pridolian and 50 m of the Lower Lochkovian strata.

The last Paleozoic stratigraphic gap affected sediments from the Middle Devonian to the Lower Carboniferous and was related to the Variscan erosion [30]. It took place presumably in the Late Carboniferous and removed about 1000 m-thick pile of sediments, with 400 m of the Middle Devonian, 400 m of the Frasnian, 500 of the Fammenian and 600 m of the Lower Carboniferous strata.

The Silurian formations are overlain by Permian strata with a total thickness of only 124 m (depth interval 1150 to 1026 m). The Rotliegend is represented by the Darłowo Formation (56.5 m thick) that developed as sandstones with thin interbeds of conglomerates. The Zechstein succession shows anomalously low thickness and includes the lowermost horizons of the PZ1 cyclothem (25 m thick) covered by a terrigenous Zechstein horizon, 42 m in thickness. Between the Early Permian and the beginning of the Late Rotliegend, the primary hiatus occured. The Mesozoic sediments observed in the Słupsk IG-1 well occur at a depth interval from 1026 to 162 m and also have numerous stratigraphic gaps. Lithology comprises Triassic and Upper Cretaceous (Cenomanian–Campanian) marls, limestones and calcareous mudstones. The Mesozoic stratigraphic gap affects sediments of the following units: Rhaetian (Upper Triassic), Toarcian (Lower Jurassic), Middle and Upper Jurassic, and, finally, Maastrichtian (Upper Cretaceous). Moreover, hiatuses are known from Ladinian and Carnian (Middle/Upper Triassic) successions.

Stratigraphic column of the Słupsk IG-1 includes also 162 m-thick sequence of Cenozoic sediments, mostly sandstones with clays and muds belonging to the Paleogene (Oligocene, 12 m thick), Neogene (Miocene, 103 m thick) and Quaternary (Pleistocene–Holocene, 47 m thick). The stratigraphic gap embraces both the Paleocene and Eocene strata, whereas the absence of Pliocene succession is caused by the hiatus.

Measurements of effective porosity and permeability were carried outfor 60 core and cuttings samples [42] obtained during the drilling of the Słupsk IG-1 well. Methodology and measurement devices were provided by the Geoservice Company (Paris, France), which has been working for the Polish petroleum industry for some 70 years [43]. In most samples, permeability values were below the accuracy limit of the device and only few results exceeded 0.1 mD (Figure 5).

Figure 5. Permeability versus effective porosity cross-plot for samples of quartz sandstones from modeled depth interval of 4716 to 5000 m below surface. Laboratory measurements after Modliński [42].

Figure 5. Permeability versus effective porosity cross-plot for samples of quartz sandstones from modeled depth interval of 4716 to 5000 m below surface. Laboratory measurements after Modliński [42].

3.1. Paleothermal Modelings

The modeling of diagenetic processes requires the recognition of not only the recent temperature but also the evolution of paleothermal conditions in geologic time because the latter factor is crucial for diagenetic transformations of reservoir rocks. Variability of temperatures in geological time is evaluated by paleothermal modelings, which practically means the solution of a differential equation of heat transport in rock formation, in 1D/2D variants [44]. For the true conditions, taking into account the diverse geological settings, the finite element calculation method is applied.

The numerical mesh of finite elements was prepared for the B’-B segment of the A’-A cross-section using the relevant procedures available in the PetroMod software (Figure 6). The lithological-structural cross-section was designed for generalized lithologic units. Thermal parameters of rock formations: thermal conductivity and heat capacity were determined using the methodology provided by the software, which related them to the lithology and porosity of analyzed rock formations, with the awareness of rising temperatures (Table 1) [4,45,46]. A strong link between thermal conductivity and porosity of rocks influenced their variability in geologic time, i.e., thermal conductivity increased with decreasing porosity, in accordance with the advancing diagenetic processes. Changes in porosity in geologic time were determined from the solutions of mechanical and chemical compaction equations, as described in the following chapter. Thermal modelings had to also consider the additional radiogenic heat sources, i.e., the Paleoproterozoic gneisses and the thick Silurian siltstones/mudstones. Moreover, the paleothermal modelings also included the influence of high deposition rates of cold sediments having low thermal conductivities and high heat capacities upon reduction in surface heat flow, known in the literature as the thermal blanketing effect [47,48].

Figure 6. Finite element mesh for modeled mechanical and chemical compactions and paleothermal processes, referred to as A-A’ geological cross-section— Figure 3. Colors correspond to lithologic and facies development. Yellow—sandstone facies of Ediacaran–Lower Cambrian succession. Stratigraphic symbols as in Figure 3.

Figure 6. Finite element mesh for modeled mechanical and chemical compactions and paleothermal processes, referred to as A-A’ geological cross-section— Figure 3. Colors correspond to lithologic and facies development. Yellow—sandstone facies of Ediacaran–Lower Cambrian succession. Stratigraphic symbols as in Figure 3.

The solution of the differential equation of heat transport required the determination of boundary conditions for the top, the bottom and the side-edges of the paleothermal model. For the top of the model, the boundary condition was represented by surface temperature, which depended on paleoclimatic conditions. For the bottom of the model, the value of basal heat paleoflow had to be determined. For the vertical side-edges of the model, no heat transfer was presumed. Correct determination of the top and the bottom boundary conditions was crucial for the accurate reconstruction of temperature distribution within the model.

The changes in heat flow values in geologic time were estimated from the calibration of thermal history of rock formations using the vitrinite reflectance. This parameter characterizes the ability of thermally matured organic matter to reflect the visible light [44,49,50]. Taking into account the genetic constrains, the petrographic composition and the forms of occurrence, the organic matter encountered in Paleozoic–Mesozoic succession of the Słupsk IG-1 well is poorly diversified. Moreover, the organic matter from Cambrian and Ordovician sediments differs from that known from the Silurian strata where bituminous matter dominates together with solid bitumens and coalified plant/animal remains [51]. Heat flow was modeled by fitting the depth-controlled trends of thermal maturity of organic matter (measured as changes in vitrinite reflectance) to the calculated trend, determined by the changes in modeled parameters in geologic time. The recent heat flow value was established from the temperature logs and the averaged thermal conductivity values of the given rock formations. Then, the heat flow values were verified by consecutive iteration cycles of paleothermal modeling. However, the values and the spatial distribution of heat flow could be locally modified by distribution of thermal parameters of rock formation: thermal conductivity, heat capacity and the presence of radiogenic heat sources. Moreover, thermal parameters of rock formations were controlled by various geological processes, which modified the heat flow characteristics. One such modifying process of recent heat flow in the study area was the cooling–freezing of the near-surface zone during the Pleistocene glaciations [25,52].

In order to determine the boundary conditions relevant for the study area, the preliminary 1D/2D calibration modeling was executed using the standard model of temperature changes in the near-surface zone, in geologic time, available in PetroMod software [53]. Additionally, a paleoclimatic correction was implemented for temperatures in the near-surface zone, for the glacial and the interglacial epochs, as −3 and +7 °C, respectively [52]. The bottom boundary condition, i.e., the basal heat flow, was established at 7 km depth and this value was successively iterated until the satisfactory consistency was obtained between the modeled and the recent measured temperature distributions as well as the modeled and the measured vitrinite reflectance in the Słupsk IG-1 well.

3.2. Modeling of Mechanical Compaction Process

Mechanical compaction belongs to the principal diagenetic processes and leads to reduction in porosity by closer packing of clastic grains and organic relics. The influence of compaction on porosity depends mostly on effective stress and lithology. Compaction of sandstones results in less porosity reduction than compaction of clayey sediments because clay minerals are much more plastic than quartz or feldspar grains [54]. Hence, the intensity of porosity reduction with the depth in sandstones increases with the increasing content of clay minerals. Moreover, mechanical compaction of sediments is also influenced by cementation. In sandstones affected by early quartz cementation, compaction-induced reduction in porosity is lower. At greater depths, with the increasing temperature, intensive cementation may halt mechanical compaction despite the progressing burial of rock formation [55].

The simplest modeling method of mechanical compaction applies the Athy equation [5], which assumes the hydrostatic conditions and the lithologically homogenous overburden. The Athy equation relates the reduction in total porosity to the maximum burial depth of rock formation using the exponential function determined by initial porosity and the compaction coefficient. However, the application of a “classic” Athy equation to the modeling of porosity reduction with the depth may generate essential errors when, recently or in geologic history, abnormal pore pressures occurred in the studied rock formation or the overburden was untypically lithologically diversified. Laboratory experiments together with geological and drilling data acquired from various sedimentary basins indicate that the driving force of mechanical compaction is the effective stress sensu Terzaghi, which can be calculated from Equation (1) [44,55,56]:

$${{\sigma }^{\prime }}_z=\sigma -p$$

(1)

where:

• ${{\sigma }^{\prime }}_z$—vertical component of effective stress (MPa);
• $\sigma$—stress generated by overburden load (MPa);
• $p$—pore pressure (MPa).

The formula of porosity reduction, corresponding to mechanical compaction process, was proposed by Smith [57] as Equation (2):

$$\varphi ={\varphi }_0{e}^{{-k\sigma {}^{\prime }}_z}$$

(2)

where:

• ${\varphi }_0$—inital porosity;
• $k$—compaction coefficient (GPa⁻¹).

Both the initial porosity and the compaction coefficient depend on lithology of a given sediment. For a sedimentary sequence composed of the Łeba, the Kluki and the Żarnowiec formations, the four-component lithology model was accepted, which applied the standard lithologic determinations taken from the PetroMod software library. It assumed the following percentages of particular lithologic components (in vol.%): quartz sandstone (80.36), clay intercalations (12.25), feldspars (5.79) and carbonate cement (1.60). Corresponding averaged values of parameters ${\varphi }_0$ and $k$ were: 33% and 19.78 GPa⁻¹, respectively. Such an attempt to model the compaction process required the recognition of pore pressure changes in geologic time. In the PetroMod software, this value was determined by the solution of the differential equation, which described fluid flow as a reaction to overburden load [44].

3.3. Modeling of the Chemical Compaction Process

The applied modeling methodology included the solution of the differential equation proposed by Walderhaug [7], in which the precipitation of silica depends on the temperature of and free surface available on quartz grains (3) [7]. Moreover, it was assumed that the transport of silica proceeded due to diffusion in a closed system.

Dissolved silica migrated over short distances and was precipitated at the available surfaces of the closest quartz grains, free of cement and/or clay minerals coatings. Modeling of chemical compaction was carried outfor late cementation (Phase II), which took place at temperatures over 90 °C:

$$\frac{\partial \varnothing }{\partial t}=-\frac{M_q}{{\rho }_q}\frac{6\left(1-f_q\right)f_v}{d_q}\frac{\varphi }{{\varphi }_0}A{e}^{\frac{-E}{RT}}$$

(3)

where:

• $R=8.31447\frac{\mathrm{Ws}}{\frac{\mathrm{mol}}{\mathrm{K}}}$—gas constant;
• $f_q$—quartz grain coating factor (%);
• $f_v$—volume fraction of quartz grains (%);
• $d_q$—average grain size of quartz (mm);
• $A={10}^{-11}\frac{\mathrm{mol}}{\frac{{\mathrm{cm}}^2}{\mathrm{s}}}$—frequency factor;
• $E=61\frac{\mathrm{kJ}}{\mathrm{mol}}$—activation energy;
• $M_q=0.06009\frac{\mathrm{kg}}{\mathrm{mol}}$—quartz mass molar;
• ${\rho }_q=2650\frac{\mathrm{kg}}{{\mathrm{m}}^3}$–—quartz density;
• $T$—absolute temperature (°K);
• $t$—geologic time (Ma).

4. Results and Discussion

4.1. 1D Calibration Modelings

The 1D modeling carried out for the Słupsk IG-1 well aimed to verify the results of 2D modelings, i.e., to determine the time frames of changes in heat paleoflow and the amount of erosion shown by the vitrinite reflectance curve. The 1D modeling demonstrated consistency of vitrinite reflectance values measured in samples from the Słupsk IG-1 with the results of paleothermal modeling, although the trend of R_o% values was complex and its interpretation was ambiguous. The pattern of plotted measurement points for Paleozoic sediments from 1200–5000 m depth interval resembled the rotated lightning (green line in Figure 7C). In the upper part of Silurian succession, down to about 2000 m depth, a distinct linear (in logarithmic scale) increase was observed with regard to vitrinite reflectance (see segment 1 of the green line in Figure 7C). Then, over the 1500 m-long interval, the R_o% values oscillated around 2%, forming a subvertical trend, typical of organic matter maturity profiles known from the other wells drilled in the northwestern margin of the BB (see segment 2 of green line in Figure 7C) [28]. Starting from 3700 m depth, the vertical trend changed again and in both the Lower Silurian and the Cambrian sediments, and the vitrinite reflectance trend can be approximated by a straight line, almost parallel to fragment 1 (compare fragments 1 and 3 of the green line in Figure 7C).

Such a complicated trend of R_o% values with the depth cannot be approximated by a single theoretical R_o% curve calculated for a given heat flow trend through the time and the amount of erosion using the common, standard algorithm after Sweeney and Burnham [58]. In recent decades, several geological models have been proposed (see [28,30]) that aimed to fit the theoretical and the modeled vitrinite reflectance curves. The fitting problem of this subvertical trend in computer models even resulted in the development of hypothesis of organic matter maturation under the overpressured conditions [59].

The iterative 1D modelings enabled the authors to estimate the amount of erosion in the Paleozoic, mostly in the Late Silurian, Devonian and Early Carboniferous to be about 2 km. Exact thickness of eroded sediments, including post-Caledonian epigenetic erosion (350 m) and Late Variscan erosion (1900 m), reached 2250 m (Figure 7). Additionally, the 1D modelings allowed us to estimate 525 m of Mesozoic erosion. These values compared to the thermal model, which assumed a constant value of heat flow, provided the best fitting of the vitrinite reflectance curve to measured R_o % values, except for a single thermal pulse related to the Variscan orogeny. In the area east of Słupsk town, Poprawa et al. [60] presumed a 3000 m-thick erosion of Silurian–Carboniferous sediments in their thermal model with a constant value of heat flow in the geologic time. These authors also noticed that the amount of Paleozoic erosion was anomalously overestimated; hence, the proposed thermal history model had to include the increased heat flow in the Cretaceous. However, the vitrinite reflectance values R_o% measured in samples of Mesozoic sediments, although high, did not justify such a thermal event. Moreover, the measured R_o% values did not support the concept of 1 km-thick erosion during the Laramide orogeny, as quoted in the model after [30] because the Cretaceous sequence described from the Słupsk IG-1 lacked only the Maastrichtian sediments, whereas the total thickness of Cretaceous formation was 400.5 m.

In order to find the most appropriate model of changes in heat flow in time, four thermal history models were tested (Figure 7), presented in the literature concerning the study area by: (i) Botor [30], (ii) Poprawa [28] and (iii) Corrado et al. [31]. In his model, Botor [30] assumed the Variscan thermal pulse, in which the heat flow value increased to 70 mW/m² and then decreased to the recent level of 53 mW/m² (see black line in Figure 7). In the models after Corrado et al. [31] and Poprawa [28], the Cretaceous thermal pulse was considered, which attained the maximum value at the end of the Cretaceous depositional cycle and then decreased down to the recent value. Two variants of Cretaceous paleothermal flow were analyzed, with peak values of 70 and 100 mW/m² (see purple and yellow line in Figure 7). Additionally, a model was considered with a constant value of heat flow: 53 mW/m² during the whole geologic history of the study area (see the blue line in Figure 7). Generally, the black and yellow curves are alternatives but the black curve better fits the vitrinite reflectance measurements for the upper part of the Silurian segment and for the Cambrian one (Figure 7), whereas the violet and the blues curves represent averaged R_o% values for the whole Silurian segment (Figure 7).

It is difficult to decide which proposed model of HF changes through time is the closest to the reality without determination of additional geochemical indicators of paleotemperatures. Hence, the authors accepted, for further calculations, the model with the Variscan thermal event of 70 mW/m² and the basal post-Variscan heat flow value close to the recent one, i.e., 53 mW/m² (see the black line in Figure 7), taking into account the amount of Late Variscan erosion to be about 1900 m. This selection is reasonable from the geological point of view and is supported by well documented magmatic event in the Carboniferous (see Chapter 2.1). Moreover, it explains the high R_o% values found in Cambrian sediments, which are the target of our research.

The authors regard the presumptions of their model to be valid for the marginal zone of the EEC and reject the hypothesis about constant value of heat flow. The 1D modelings reveal that maximum temperature in the Żarnowiec Formation was attained at the end of the Devonian. In the Słupsk IG-1 well, a value of ~300 °C was estimated (Figure 7).

Figure 7. Results of 1D paleothermal modelings: (A) Burial and paleotemperature history of sedimentary formations in the Słupsk IG–1 well. (B) Considered models of basal heat flow variability through geologic time. (C) Fitting of modeled to measured values of vitrinite reflectance for heat flow models presented above. Explanations: red crosses—measured values of vitrinite reflectance [51]. Green line—manual approximation of vitrinite reflectance changes within Paleozoic sediments.

Figure 7. Results of 1D paleothermal modelings: (A) Burial and paleotemperature history of sedimentary formations in the Słupsk IG–1 well. (B) Considered models of basal heat flow variability through geologic time. (C) Fitting of modeled to measured values of vitrinite reflectance for heat flow models presented above. Explanations: red crosses—measured values of vitrinite reflectance [51]. Green line—manual approximation of vitrinite reflectance changes within Paleozoic sediments.

4.2. Paleothermal 2D Modeling

Solution of the heat transport equation provided the time-controlled temperature distribution within the sedimentary cover modeled along the analyzed A’-A cross-section. The modeling indicated two periods of thermal disequilibrium. The first period took place in the Silurian, as disclosed by anomalously low values of vitrinite reflectance measured in samples of Lower and Middle Silurian sediments (Figure 8), and was related to the extreme deposition rate of water-saturated clayey/muddy clastic material. The second period was related to glaciations and the resultant remarkable reduction in near-surface temperature. Finally, after several iterations of calibration modeling, it was demonstrated that the application of a basal heat flow value of 53 mW/m² together with the correction for radiogenic heat and for the influence of glaciation provided a perfect fitting of the synthetic and the true thermal profilings (Figure 8). It must be noticed that the recent heat flow varies with the depth. Generally, heat flow increases at a depth interval from 7 to about 2 km b.s.l. due to the effect of radiogenic heat; then, it decreases towards the Earth’s surface due to cooling of the near-surface zone at the end of the Pleistocene.

Basal heat flow in the Proterozoic rocks (at 7 km depth, consistent with the depth of calculation model) was estimated as 53 mW/m². At shallower depths, this heat flow was supplemented by two additional sources of radiogenic heat originated from: (i) felsic Proterozoic intrusions, which raised heat flow value to 57 mW/m² in the top of the Precambrian, and (ii) thick, Ordovician–Silurian claystone/siltsone/mudstone successions, which gave maximum heat flow values of over 58 mW/m² in the Middle Silurian sediments. In the geological history of the study area, two main periods of thermal disequilibrium took place within the sedimentary cover. The first period was related to very high deposition rate of cold Upper Silurian (Lludlovian) sediments having low thermal conductivity and high heat capacity (Figure 9). The deposition rate in the Lludlovian, calculated after decompaction of sediments, was high (over 1000 m/Ma), which resulted in a significant decrease in the surface heat flow. This effect took place despite the interaction of the opposite process, i.e., the radiogenic heating of thick Silurian mudstones/claystones. The results of modelings for the end of the Lludlovian revealed the HF decrease to do 49 mW/m² in the near-surface zone and its increase with the depth to 55 mW/m² in the top of Paleozoic succession at about 4700 m depth (Figure 9). The thermal equilibrium has been restored at the Silurian/Devonian turn due to the significantly lower deposition rate.

The second important thermal disequilibrum period in the vicinity of the Słupsk IG-1 well was related to an abrupt change in the surface temperature caused by the influence of the Pleistocene glacial epoch. Multiple and long-lasting glaciations covered the study area with a thick ice sheet, which resulted in cooling and freezing of the near-surface zone, thus modifying both the thermal conductivity and heat capacity of rock formations [25,52]. Similarly to that in the Silurian, this event resulted in heat flow decrease to 47 mW/m² but its depth range was much shallower, i.e., below 1000 m (Figure 8). The thermal equilibrium has not been attained, as yet, and the surface heat flow is still influenced by the glaciations, as documented by [52].

As quoted in Chapter 2.2, the geological and geophysical data acquired in the Słupsk IG-1 well were used as the reference dataset in the research project. Moreover, this dataset played a roughly similar role in some earlier studies. Hence, it seems necessary to discuss the differences between the published results and those obtained in this project.

In the Słupsk IG-1 well, advanced thermal measurements were carried out after a 15-day-long well stabilization period. Moreover, three analyses of thermal field parameters were performed by independent authors (Figure 8). However, the results differ remarkably. According to Szewczyk [27], temperatures measured in the well contradicted those corrected after well stabilization due to insufficient stabilization time. Values of heat flow were estimated using the temperature log and the thermal conductivity of rocks interpreted from lithology and porosity logs [37]. Interpretation of thermal parameters in the Słupsk IG-1 well provided: (i) higher temperature, corrected with respect to that observed in the top part of Cambrian succession and at the well bottom, and (ii) anomalous heat flow values. Recent heat flow value measured in the Słupsk IG-1 well was 78 mW/m², whereas in the adjacent areas, this value varied from 55 to 65 mW/m² [27]. However, the recent heat flow value obtained from paleothermal modeling was 54 mW/m² [28]. Additionally, the heat flow value in the area of Słupsk IG-1 well, estimated independently during construction of the map of this parameter (applying the paleoclimatic correction), was about 50 mW/m² (Figure 2) [25]. Such low heat flow values in the near-surface zone confirmed the previous estimations: 54 mW/m² after Poprawa [28] and somewhat over 50 mW/m² after Majorowicz [25]. Hence, the heat flow value of 78 mW/m² published by Szewczyk [27] appeared to be highly overestimated, which presumably resulted from too high a correction of temperature measured in the Słupsk IG-1 well. Such an erroneous correction was possibly caused by thermal disequilibrium still existing after insufficient, 15-day-long stabilization period of that well.

Figure 8. Calibration of the paleothermal model based upon measurements in the Słupsk IG-1 well. Explanations: black line—modeled vitrinite reflectance curve; blue line—recent temperature curve; violet line—recent heat flow curve; red crosses—measured vitrinite reflectance values; black dots—temperature log data. Temperature measurements after Szewczyk [27]; vitrinite reflectance measurements after Grotek [51]. Stratigraphic symbols as in Figure 3.

Figure 8. Calibration of the paleothermal model based upon measurements in the Słupsk IG-1 well. Explanations: black line—modeled vitrinite reflectance curve; blue line—recent temperature curve; violet line—recent heat flow curve; red crosses—measured vitrinite reflectance values; black dots—temperature log data. Temperature measurements after Szewczyk [27]; vitrinite reflectance measurements after Grotek [51]. Stratigraphic symbols as in Figure 3.

Figure 9. Blanketing effect after deposition of Lludlovian sediments.

Figure 9. Blanketing effect after deposition of Lludlovian sediments.

Moreover, the study revealed that the fitting of modeled and true values of vitrinite reflectance measured in samples from the Słupsk IG-1 well required: (i) increase in Paleozoic basal heat flow values to 70 mW/m², (ii) presumption of Variscan uplift of the area followed by (iii) approximately 2 km-deep erosion. Due to the lack of additional reference wells, which would provide data suitable for determination of the changes in heat flow values in geologic time, it was assumed that the bottom boundary condition did not change along the analyzed part of geological cross-section. When referring to the recent value of heat flow, this assumption was validated by the map presented by Majorowicz [25], which documented very low variability of recent heat flow along the cross-section and temperature distribution in the top of Cambrian sandstones succession shown in Figure 2 [26].

Distribution analysis of maximum paleotemperatures along the modeled cross-section (Figure 10) is important for evaluation of the changes in heat flow values and the conditions controlling the development of reservoir properties of studied sedimentary formations. The presented distribution is an effect of 2D paleotemperature modelings, which consider the accepted subsidence/inversion model of the study area and variability of geothermal conditions in time and space (Figure 7). The changes in maximum paleotemperatures in Paleozoic sediments adjust (to some extent) to the relief of the top surface of crystalline basement and are positively correlated with its burial depths. However, such relationships may generate a false impression that the recent burial depth is a dominant factor controlling the distribution of paleotemperatures. Comparison of distributions of maximum paleotemperatures and recent temperatures (Figure 10) reveals obvious differences. The recent distribution of temperatures within the studied rock formations is controlled mostly by their burial depths and is only insignificantly modified by diverse lithology as well as the resultant changes in thermal conductivity and heat capacity of the rocks. The distribution of maximum paleotemperatures shows the presence of a radical boundary located at the bottom surface of Permian formations, which separates both the Lower and Middle Paleozoic formations from Permian, Mesozoic and Cenozoic units of significantly lower maximum temperatures. Moreover, in that younger group of sedimentary units, the maximum temperatures of rocks are almost identical to the recent temperatures (Figure 10). Good representation of lithological variability and related thermal parameters is the distribution of the recent geothermal gradient along the modeled cross-section (Figure 11), which evidences the decreasing trend with the depth. It corresponds to the changes in depth-controlled thermal conductivity of rocks dependent on the gradual reduction in their porosity.

Taking into account the differences between the recent and the maximum values of geothermal parameters as well as the results of maturity modelings of organic matter for various variants of thermal history of the study area, the concept of strong thermal pulse (heat flow = 70–100 mW/m²) affecting the western part of the PS (close to the margin of the EEC) in the Cretaceous seems to be doubtful. However, it does not preclude the possible appearance of Mesozoic thermal pulses in other parts of the BB [34]. The results of Middle and the Late Paleozoic show the maturation of organic matter and stimulation of advanced chemical compaction, with the latter reducing the reservoir parameters of studied rocks almost to zero (see black plot in Figure 7). The direct aim of paleothermal modelings was the recognition of temperature distribution within the Ediacaran Lower–Cambrian sandstones. The 2D model showed the increase in the recent temperature towards the southeast, concordant with the increasing burial depth of these sandstones (Figure 12). The increasing temperature trend is consistent with the decreasing porosity trend. The lowest reservoir temperature—about 100 °C—was found in the uplifted, offshore part of the study area. On the contrary, the temperature in deeply buried reservoir sandstones near the Słupsk IG-1 well reached about 164 °C and, in the deeper parts located to the southeast, it even approached 180 °C.

Figure 10. Models of temperature changes in geologic history of the study area: (A) Maximum temperatures attained by rocks of sedimentary cover in the geologic history. (B) Distribution of recent temperatures within the sedimentary cover. (C) Distribution of differences between the maximum and the recent temperatures within the sedimentary cover.

Figure 10. Models of temperature changes in geologic history of the study area: (A) Maximum temperatures attained by rocks of sedimentary cover in the geologic history. (B) Distribution of recent temperatures within the sedimentary cover. (C) Distribution of differences between the maximum and the recent temperatures within the sedimentary cover.

Figure 11. Recent average vertical temperature gradient showing the decreasing trend with the burial depth.

Figure 11. Recent average vertical temperature gradient showing the decreasing trend with the burial depth.

Figure 12. Model of temperature distribution in Ediacaran–Lower Cambrian sandstones. Stratigraphic symbols as in Figure 3.

Figure 12. Model of temperature distribution in Ediacaran–Lower Cambrian sandstones. Stratigraphic symbols as in Figure 3.

4.3. 2D Compaction Modeling

The results of modeling indicated a more important role of mechanical compaction in porosity loss of Ediacaran–Lower Cambrian sandstones, in comparison to their chemical compaction. The main phase of mechanical compaction has commenced at the beginning of the Silurian due to very high deposition rate of thick siltstone/mudstone succession (Figure 13). Before the end of the Silurian, mechanical compaction reduced porosity of sandstones by a dozen percent. In the Late Carboniferous, the decrease in sandstone porosity was much slower thanks to lower deposition rate and partial hardening of sandstone fabric caused by chemical compaction (Figure 13). After the end of Early Carboniferous sedimentation, mechanical compaction ceased due to the end of burial and the beginning of regional uplift related to the Variscan orogeny (Figure 7). In the Late Carboniferous, the closure of pore space continued, exclusively by chemical compaction; however, this process ended at the Carboniferous/Permian turn, as a result of cooling of reservoir formation caused by the Variscan uplift. Total thickness of Permian–Mesozoic formations did not exceed the thickness of Silurian, Devonian and Lower Carboniferous sediments, neither in the geologic history nor recently. Hence, the load from the overburden and the effective pressure did not exceed the values estimated for the end of Early Carboniferous.

Figure 13. Evolution of diagenetic mechanical and chemical compactions of Ediacaran–Lower Cambrian sandstones. Explanations: red line—porosity changes in Lower Cambrian sandstones, in the vicinity of the Słupsk IG-1 well; black line—porosity changes caused by mechanical compaction; blue line—porosity loss caused by chemical compaction.

Figure 13. Evolution of diagenetic mechanical and chemical compactions of Ediacaran–Lower Cambrian sandstones. Explanations: red line—porosity changes in Lower Cambrian sandstones, in the vicinity of the Słupsk IG-1 well; black line—porosity changes caused by mechanical compaction; blue line—porosity loss caused by chemical compaction.

The spatial image of diagenetic processes resulted from 2D modelings clearly demonstrated that porosity of studied Ediacaran–Lower Cambrian sandstones decreased to the southeast, in accordance with the increasing burial depth of sediments. In the onshore part of the BB, further to the southeast, mechanical compaction reduced the porosity of sandstones to about 7% (Figure 14). Porosity loss due to chemical compaction of sandstones was about 5% in the onshore part of the cross-section, with an increase to the southeast, to over 6% (Figure 14). However, chemical compaction was less advanced in the offshore part of the cross-section, where porosity loss was only 3%–5%. Thus, the role of mechanical compaction in porosity reduction for sandstones was more advanced than that of chemical compaction. The final model of porosity changes in Ediacaran–Lower Cambrian sandstones was generated using the simultaneous modeling of both the mechanical and chemical compactions (Figure 14).

Recently observed total porosity of sandstones decreases towards the southeast, in accordance with the progressing compaction, which is an effect of increasing temperatures and burial depth. The highest and the lowest porosity values of the sandstone reservoir formation are 6.56 and 1.10%, respectively. These values characterize the shallowest (offshore segment) and the deepest (onshore segment, Słupsk IG-1) parts of the geothermal reservoir. Porosity values modeled in the Słupsk IG-1 well are close to the effective porosity measured in the laboratory, i.e., about 1%. Such low porosity of a rock formation demonstrates that in the vicinity of Słupsk, deeply buried Ediacaran–Lower Cambrian quartz sandstones are almost devoid of open pore spaces.

The very low permeability values do not correlate with the effective porosity values [42]. The effective porosity of the Żarnowiec Formation sediments is 1.28%, whereas those of the Łeba and the Kluki formations are 1% each [42]. Hence, the Ediacaran–Lower Cambrian succession shows very poor reservoir properties, closely controlled by mechanical and chemical compactions. Such low values are typical of successions covered by thick Ordovician and Silurian claystones, siltstones and mudstones.

Figure 14. Modeling results of the Ediacaran-Lower Cambrian sandstones mechanical and chemical compaction processes. (A) Model of porosity distribution caused by mechanical compaction. (B) Model of lost porosity distribution caused by chemical compaction. (C) Porosity distribution obtained from simultaneous modelings of mechanical and chemical compactions. Stratigraphic symbols as in Figure 3.

Figure 14. Modeling results of the Ediacaran-Lower Cambrian sandstones mechanical and chemical compaction processes. (A) Model of porosity distribution caused by mechanical compaction. (B) Model of lost porosity distribution caused by chemical compaction. (C) Porosity distribution obtained from simultaneous modelings of mechanical and chemical compactions. Stratigraphic symbols as in Figure 3.

Abrupt porosity reduction took place in the Silurian/Devonian during the burial of Cambrian sediments down to the depths from 2500 to 3000 m and then below 5000 m, which reduced porosity values down to 1% [61]. The most intensive chemical compaction process was silicification. Quartz cementation proceeded in both the early (Phase I) and the late (Phase II) diagenesis (i.e., from eo- to mesodiagenesis). The silica sources were both the clay minerals and quartz, which dissolved under pressure exerted by the overburden.

It cannot be neglected that some additional amounts of silica might have been supplied by illitization of clay minerals, chloritization of biotite or sericitization of feldspars [61,62,63]. Petrographic studies of samples from drillings in the Polish part of the BB enabled Sikorska and Pacześna [61] to demonstrate that the decrease in porosity with the depth was accompanied by the increasing content of quartz cement. Studies on sandstones from the Polish and the Lithuanian parts of the BB indicated that the cementation phase, Phase I, took place at the temperatures of 60–70 °C and Phase II took place at 90–130 °C [61,62,63,64,65].

5. Conclusions

The modeling of diagenetic and paleothermal processes reported above enabled the authors to estimate distribution of porosity and temperature within the Cambrian sandstone succession, regarded as a potential target for supplying the future geothermal installations.

Generally, the intergranular porosity values of studied sandstones are very low, usually close to 1%, and do not exceed 7%, at most. Changes in porosity are negatively correlated with burial depth and paleotemperature of the Ediacaran–Lower Cambrian sandstone formation. Its temperature exceeded 300 °C in the marginal, southeastern part of modeled B’-B cross-section. Insignificant variability ranges of porosity and permeability values together with their statistical dispersion preclude the recognition of relationships between effective porosity and permeability measured on rocks samples in the laboratory (Figure 5). Consequently, prediction of permeability distribution within the sandstone succession related to intergranular porosity cannot be deterministically and numerically modeled. However, the qualitative comparison of modeling results with the laboratory measurements indicates that in the onshore part of the analyzed cross-section, permeability of Cambrian sandstones does not exceed 1 mD. Recent temperatures of sandstone succession are relatively high and pertinent to its burial depths, with some small modifications caused by differences in heat flow values. In the uplifted zones of the northwestern part of the B’-B cross-section, recent temperatures of Cambrian sandstones are up to about 100 °C, whereas in deeply buried zones, these may reach about 180 °C. Generally, the paleothermal modeling led to recognition and interpretation of the changes in heat flow values with the depth and verified the results of previous analyses of thermal field parameters in the vicinity of the Słupsk IG-1 well.

The low values of reservoir parameters reject, in practice, the studied Cambrian sandstones as a potential geothermal reservoir suitable for supplying the future geothermal installations. A practical solution can be the development of some zones showing relatively increased fracture-type permeability by utilization of closed-in oil and gas wells, and/or by application of well stimulation technologies. This first solution seems to be promising taking into account the successive exhaustion of hydrocarbon resources, so that such sites and existing petroleum production equipment can be effectively utilized by geothermal installations. Relatively high temperatures of Cambrian sandstones together with their final quartz cementation, which increased mechanical properties of sediments, favor the application of hydraulic fracturing technologies. Hence, there is quite good potential for effective utilization of these sediments for EGS or HDR installations if their reservoir properties are improved by, e.g., fracking, thus enabling the efficient circulation of a working medium from depth to surface facilities.

The results of studies confirm the applicability of paleothermal modelings to recognition of geothermal aquifers of specific characteristics such as those existing in the marginal zone of the EEC. Important contributions useful for evaluation of geological formations, such as prospective geothermal aquifers, are the results of compaction modelings presented as spatial distributions of porosity in potential geothermal reservoirs. These results enable the researchers to estimate the mechanical parameters of reservoir formations, which are crucial to planning the fracking operations. The methodology presented above may also be useful for petroleum exploration of the BB, particularly its offshore part, which is still regarded as a prospective area where discoveries of new hydrocarbon accumulations could be made, both the conventional and unconventional ones. Similarly, the results of petroleum exploration projects may be adapted to selection of prospective areas for localization of geothermal installations.