Geothermal Resources and ATES Potential of Mesozoic Reservoirs in the North German Basin

Abstract

Mesozoic sandstone aquifers in the North German Basin offer significant potential to provide green and sustainable geothermal heat as well as large-scale storage of heat or chill. The determination of geothermal and subsurface heat storage potentials is still afflicted with obstacles due to sparse and partly uncertain subsurface data. Relevant data include the structural and depositional architecture of the underground and the detailed petrophysical properties of the constituting rocks; both are required for a detailed physics-based integrated modeling and a potential assessment of the subsurface. For the present study, we combine recently published basin-wide structural interpretations of depth horizons of the main stratigraphic formations, with temperature data from geological and geostatistical 3D models (i.e., CEBS, GeotIS). Based on available reservoir sandstone facies data, additional well-log-based reservoir lithology identification, and by providing technical boundary conditions, we calculated the geothermal heat in place and the heat storage potential for virtual well doublet systems in Mesozoic reservoirs. This analysis reveals a large potential for both geothermal heating and aquifer thermal energy storage in geologically favorable regions, and in many areas with a high population density or a high heat demand. Given the uncertainties in the input data, the applied methods and the combination of data from different sources are most powerful in identifying promising regions for economically feasible subsurface utilization, and will help decrease exploration risks when combined with detailed geological site analysis beforehand.

1. Introduction

1.1. General Background

In the context of the energy transition, the focus has been on the provision of electricity from renewable energy sources. However, in many regions of the world, especially in the Northern Hemisphere, more energy is consumed for heating and cooling purposes than for electricity. In Germany, for example, 56% of all primary energy consumption is in the heating-and-cooling sector, while only 25% are needed for electricity provision [1]. This proportion of energy consumption by the different sectors is in direct contrast to the energy provided from renewable sources: While the installation of renewable energy for electricity has reached levels of more than 50% in Germany by 2021, only 14% of the heat is generated from renewable sources [1]. Hence, an efficient reduction of greenhouse gas emissions from the energy sector needs to address the heating and cooling market.

One of the challenges in the provision of heat based on renewable and sustainable energy sources is the large seasonal variability on the demand and supply site. Electricity consumption fluctuates in much narrower bands than the consumer demand of heat. In addition, demand and supply from (fluctuating) renewable sources are highly asynchronous, with a surplus in heat supply in the summer months from solar installations, but also from excess heat in buildings requiring cooling. An efficient heating and cooling system without the flexibility provided by conventional solutions like the burning of natural gas will require seasonal heat storage in large volumes, especially in densely populated urban areas, where the heating and cooling demand often exceeds 75% of the primary energy consumed (e.g., [2]). The large volumes required for this purpose could be found in the subsurface, primarily in natural aquifers [3].

While aquifer thermal energy storage (ATES) systems are well established in some parts of the world, primarily in the Netherlands, where shallow aquifers are widely used for seasonal heat storage (e.g., [3]), they are rather scarce in most other countries. In Germany, one of the reasons for their limited distribution is the potential interference with the use of shallow freshwater aquifers for drinking water supply. Therefore, deeper and saline aquifers have recently received more attention from licensing authorities and developers. Such deeper ATES systems have the additional advantage of higher temperatures, making the heat storage even more energetically favorable. On the other hand, porosity usually decreases with depth, lowering the available potential storage volume and the overall hydraulic connectivity. In addition, several non-technical issues hamper the utilization of the deeper aquifers currently. In many countries, a distinction is made between deep and shallow geothermal applications, which imply different licensing criteria and procedures. The general understanding of the geothermal resource classification varies from country to country. In Germany, for example, geosystems with depths greater than 100 m are regulated by the mining law, whereas regulations according to the water management apply to any depth, but focus on shallow formations with water of drinking water quality. Another crucial limitation is the commonly poor data availability at greater depths. In regions with hydrocarbon exploitation or deep geothermal operations, information at depths below 1500 m is often available, while the shallower depths have not received the same attention for subsurface characterization in the past. However, even if the respective data do exist, they are not necessarily publicly accessible. While countries like the Netherlands or Norway made all subsurface data accessible for the public a few years after acquisition, other countries are not as transparent and open in allowing the use of available subsurface data.

In Germany, recent developments in legislation regarding the usage of subsurface data now allow access to primary geological data. In this study, we focus on the screening of the available data to map Mesozoic target reservoirs and to determine the geothermal and thermal storage potential. Five of the most promising and widely spread Mesozoic reservoirs in Northern Germany were identified for this purpose, and their dimensions as well as their thermal and hydraulic properties were evaluated. The calculated thermal storage capacity was compared to the heat demand in densely populated areas at the surface. This exercise was performed as a contribution to a Roadmap for Net-Zero emissions (www.netto-null.org (accessed on 20 January 2021)) in Germany, developed in the framework of the Helmholtz Climate Initiative HI-CAM (www.helmholtz-klima.de/en (accessed on 21 January 2021)).

1.2. Geothermal and Heat Storage Potentials

The evaluation of the geothermal and heat storage potential of a given area can be approached in numerous ways of varying complexity. These approaches range from analyzing the subsurface temperature distribution (e.g., [4]), over investigating the heat stored in a certain reservoir (e.g., [5]), to estimates of the producible electricity of a geothermal power plant (e.g., [6]). In the study presented here, we focus on the estimation of the heat in place and the heat storage potential of a given area as indicated below.

1.2.1. Geothermal Potential—Heat in Place

The estimation of the heat in place (HIP) has a long history (e.g., [7]) with investigations by Grant [8] and, most recently, Ciriaco and coauthors [9]. For the estimation of the geothermal potential in this paper, we focus on the heat in place method after [7,10]. This methodology relies on the principle that the temperature in the subsurface is higher than that at the surface, and the difference between those two temperature levels can be translated into a potential energy resource. To keep the uncertainty of such an estimate at a low level, making as few assumptions as possible, the heat in place after the USGS method requires the lowest amount of input parameters [7,10]. In detail, the input parameters are limited to the reservoir volume, the reservoir temperature, the reservoir volumetric heat capacity, and the reinjection temperature. Optionally, conversion efficiency, plant lifetime, and plant load factor could be considered but were omitted since these are only increasing the uncertainty in the output data. Another important factor when looking at the HIP of a certain reservoir is the consideration of the over- and under-burden rock, effectively increasing the estimated potential due to conductive heat transport [11,12]. We do not consider this parameter, since it introduces uncertainty again, and the method presented in this paper only serves to outline areas of higher and lower potential, which should then undergo further investigation for their appropriateness for actual site development.

1.2.2. Heat Storage Potential—Aquifer Thermal Energy Storage (ATES)

The ATES principle utilizes the geological subsurface to store large quantities of hot or cold water. It requires suitable geological formations that act as aquifers and offer the potential to store a fluid, e.g., in porous (sandstone) or karstic (carbonate) rocks. Both the formation rock and the containing fluid each have distinctive characteristics making them less or more favorable for the intended use. In detail, highly permeable rocks of high thickness have been proven to represent good targets. Here, commonly, two wells are needed, where one well serves as the “cold well” and the other one as the “warm well”. The warm well then represents the storage medium, while the cold well, including the surrounding aquifer, is the fluid reservoir. This enables a pendulum operation consisting of two phases: charging and discharging of the storage. Similar to the heat in place, the main parameters required for the estimation of the heat storage potential are the reservoir volume, the reservoir temperature, the reservoir volumetric heat capacity, and the injection temperature. Optionally, a recovery factor and a packing correction factor can be considered to provide a more conservative estimate.

1.2.3. Parameters Considered in Geothermal and Storage Potential Studies

As outlined above, in order to evaluate the geothermal and heat storage potential of a region, several parameters have to be considered. The first and most important here is the subsurface temperature, which can be deduced in numerous ways. The most straightforward and reliable determination is done by direct measurements in boreholes. Alternatively, less cost-intensive, but more indirect and imprecise, estimations of the subsurface temperature are derived from other observations. These observations include seismics (e.g., [13]), satellite data (e.g., [14]) and magnetotellurics (e.g., [15]). Another means of deriving temperatures is numerical modeling. Here, 3D models considering both the geological subsurface in high detail as well as the physical processes involved in heat transport have proven to be a powerful tool (e.g., [16,17,18,19]).

The second important parameter is the geometry of the targeted geological reservoir. Here, both the identification of its volume and the precise knowledge of its extent pose a challenging task. For this task, borehole data are of major importance which could be jointly interpreted with seismic data and result in integrated structural modeling to provide the best estimate of the reservoir geometry.

For the potential analysis, further parameters need to be considered like the volumetric heat capacity of the reservoir and its constituents, namely the rock and fluid density, the rock and fluid specific heat capacity, and the coupling parameter porosity have to be determined or at least confined within certain bounds. Depending on the given focus on a geothermal (doublet) production or ATES system, the temperatures at which water is reinjected or the temperature at which heat is stored have to be determined as well. Another important factor to be taken into consideration is the reservoir depth and the associated drilling depth. Here, target reservoirs that are either too shallow (<100 m.b.s.), possibly in conflict with groundwater utilization, or too deep (>2500 m.b.s.), requiring too expensive boreholes, are excluded as potential targets.

For this paper, we want to provide an estimate of the geothermal and storage potential based on the parameters mentioned above. These estimates represent potentials on a basin-scale and could not consider the detailed geological situation at a certain location or the respective fluid (brine) composition. They cannot replace local feasibility studies.

1.3. Previously Published Estimates of Heat in Place and Heat Storage Potentials

Over the last decades, many studies on geothermal potential were carried out worldwide. The derived potentials can be distinguished by local, regional, and global estimates. Depending on the aim of the study, authors took more to fewer parameters into consideration and relied on a varying number of assumptions. For the global scale, a good overview is given by Limberger and coauthors [20] who showed the generalized technical potential [21], which is the theoretical HIP reduced by assumptions of reservoir lifetime and a recovery factor. They calculated this potential on all continents based on database compilations for all main parameters described above and generalized assumptions of the others. The worldwide computed potential reaches up to 200 GJ/m², while the calculated potential in the North German Basin averages in a range between 75 and 150 GJ/m². However, the resolution of this dataset is coarse (resolution of 0.1°~10 km) and the underlying structural model allows no distinction of the different reservoirs, which limits its informational value.

For the regional level of national studies, we want to outline just a few examples in the neighborhood of Germany. The Polish geological survey published a compilation of maps [22] showcasing the static resources (HIP) of several stratigraphic units, where potentials of up to 250 GJ/m² are predicted. Under consideration of certain limiting criteria, this drops down to a maximum of 75 GJ/m². In the Netherlands, advanced models are used to estimate the geothermal potential in several different aspects. For our study, we want to determine the value of HIP for comparative purposes which reaches up to 220 GJ/m², with a value of 1.1 GJ/m² when considering the potentially recoverable heat [23,24]. In Germany, a common approach focuses on the availability of the two main parameters: temperatures above a certain threshold (i.e., above 40 °C, 60 °C, 100 °C, [25]) and a sedimentary system with suitable sandstone reservoirs [25]. In their 2018 paper [26], Agemar and co-authors conclude that a revised quantitative geothermal resource assessment is overdue since the last one is almost 20 years old [26,27]. Consequently, a new interest for more accurate resource assessments is emerging [28]. There is no published methodology for estimating the (HSP) storage potential that considers the general design of ATES systems. Therefore, the method presented in this paper was developed by the authors considering the basic design criteria of ATES systems (e.g., as described in [29]).

2. Methodology and Basic Data

For the analysis of geothermal and heat storage potentials of Northern Germany, we are following a workflow outlined in Figure 1. Both heat in place and heat storage potentials are calculated for five Mesozoic stratigraphic units: the Lower Cretaceous, the Middle and Lower Jurassic, the Upper Keuper, and the Middle Buntsandstein. Each of these units contain sandstone successions of sufficient thickness to be considered for geothermal utilization [30]. Within these stratigraphic units, 9 formations are included that have been analyzed separately: Upper Bajocian, Upper Aalian, Upper and Lower Toarcian, Upper Exter, Lower Exter II and III, and Upper and Lower Schilfsandstein. Tertiary sand successions could also be considered for geothermal storage utilization, but are not in the focus of the present study.

2.1. Methodology

The heat in place (HIP) of a reservoir is given by the aquifer volume that contains the heat, the volumetric heat capacity of the aquifer, the reservoir temperature and the reinjection temperature. Considering the aquifer as a porous medium, for example, sandstone, the volumetric heat capacity $c_A$ of the aquifer is given by the specific heat of the solid $c_S$ and the liquid phase $c_W$, and the porosity $\varphi$ of the aquifer:

$$c_A=\varphi {\rho }_W{c}_W+\left(1-\varphi \right){\rho }_S{c}_S$$

(1)

Based on the definition of the volumetric heat capacity of an aquifer (Equation (1)), the heat in place HIP of the target reservoirs can be calculated as follows:

$$HIP=Q_{th}=V_A{c}_A\left(T_I-T_R\right)$$

(2)

The description of the constituents and equation parameters is listed in

Table 1. Input parameters for the different calculations. HIP = heat in place, HSP = storage potential.

Parameter	Symbol/Unit	HIP	HSP	Comment
Specific heat capacity (water/solid)	$c_{w/s}$/$\frac{J}{kgK}$	4180/723	4180/723	1, 2
Reservoir area	$A/{\mathrm{m}}^2$	1000 × 1000	1000 × 1000	3
Reservoir thickness	$H/\mathrm{m}$	>10, <20, if >20 =>20	>10, <20, if >20 =>20	4
Effective reservoir volume	$V_A/{\mathrm{m}}^3$	$AH$	$0.5AH{F}_P$	5
Injection temperature	$T_I/°\mathrm{C}$	8	90	6
Reservoir temperature	$T_R/°\mathrm{C}$	38 < $T_R$	$T_R$< 90	7
Porosity	$\Phi /--$	0.1	0.1	8
Volumetric heat capacity (bulk)	$c_A$/$\frac{MJ}{m^{3}K}$	2.125	2.125
Density (water/solid)	${\rho }_{w/s}$/$\frac{kg}{m^3}$	1000/2644	1000/2644	9
Packing correction factor	$F_p/--$	–	0.8	10
Recovery factor	$F_r/--$	–	0.7	11
Minimum depth	$d_{min}/\mathrm{m}$	100	100	12
Maximum depth	$d_{max}/\mathrm{m}$	2500	750	13

.

The heat storage potential (HSP) is given by the specific heat of the aquifer, the aquifer volume that can be used for thermal energy storage, the storage temperature (the temperature at which the thermal energy is stored in the aquifer or the temperature of the hot water that is injected into the aquifer) and the natural aquifer temperature. Including Equation (1), the storage potential $HSP$ of the target reservoirs can be calculated as follows:

$$HSP=Qth=V_A{c}_A\left(T_R-T_I\right)F_r{F}_p$$

(3)

The description of the constituents is listed in Table 1.

In order to outline the parameter choices made, we elaborate on each one individually, relating to the comment number in Table 1. (1) The value of specific heat capacity of water for ambient conditions was taken from [33]. The dependence on temperature and pressure is known, but for the targeted depth range, these effects are smaller than 2% [33] and thus ignored. (2) For the calculations (following the procedure outlined in [34]), we considered fixed heat capacities and densities for the reservoirs, which are representative for a sandstone of the composition 70% quartz, 20% feldspar, and 10% illite. The parameter values for the mineral constituents were taken from [35,36,37,38,39]. (3) We chose a grid resolution of 1000 × 1000 m for the structural model and therefore also for the calculations of HIP and HSP. (4) We chose an (accumulated) minimum thickness threshold of 10 m as anything less is assumed to be economically unfeasible, and a maximum threshold of 20 m, as larger thicknesses of consecutive sandstone are rarely encountered, and we want to provide a conservative estimate. (5) The effective reservoir volume for the storage potential is assumed to be half the total, since we only consider the heat storage side of the system. The volume is additionally reduced by a packing factor representing the areal coverage of the thermally influenced cone around the wellbore which should not overlap. See Appendix A for full description. (6) The injection temperature for the HIP calculations was considered as the surface average temperature (ca. 8 °C, [19,40,41]). For the HSP, we chose 90 °C, as this is considered as a general target storage temperature for HT (high temperature) ATES systems [42]. (7) To ensure a sufficiently large temperature difference between the extraction and injection temperature, the reservoir temperature should be more than 30 °C above the reinjection temperature for the calculation of the HIP. For the HSP, the reservoir temperature has to be below 90 °C as higher temperatures lead to negative potentials. (8) Little is known about regional reservoir-specific porosity distributions. Therefore, we chose a conservative estimate for the porosity, which lies at the lower end for a sandstone of the composition described above [43]. (9) The density of water was chosen for 8 °C at atmospheric pressure. The density of the solid was derived as described in (2). (10) The packing correction factor for the HSP is derived by circle packing density. For uniform circles, the highest circle packing density would be around 90% when arranging them hexagonally. Since we consider a regular grid (squares) and we want to provide a conservative estimate, we consider the squared variant of circular packing, representing a density of around 80%. (11) As with almost all energy storage systems, we have to consider energy losses. Although the thermal energy stored in aquifers is not lost, since the thermal energy remains in the subsurface, not all of the stored energy can be recovered, as the heat in the aquifer is present at a lower temperature level. For estimating the heat that can be extracted from the aquifer, we consider a heat recovery factor (HRF) of 0.7. This value is based on [44] for higher temperature ATES and represents a good practice average value for well-designed ATES systems. (12) We chose a minimum depth of more than 100 m below surface as the uppermost 100 m are frequently used for drinking water purposes, especially in urbanized areas, and we want to avoid any interference here. (13) The maximum depth for HIP was chosen to be 2500 m below surface as drilling costs and risks increase exponentially with depth, and this is a feasible cutoff point [45]. For the HSP, 750 m below surface was chosen in order to reduce drill costs and production costs, as the potential is higher, the shallower the storage site is located.

The temperature data necessary to calculate the HIP and HSP were taken from two different datasets. These are from the 3D numerical process model of the Central European Basin System (CEBS) [32] and from the geostatistical 3D temperature interpolation of the collected borehole temperature database of Northern Germany (mapping: GeotIS, data: Fachinformationssystem Geophysik FIS GP; cf. [46]). For more details, see Section 2.2.3.

Since calculating the available HIP and HSP alone only has limited informative value, we included published data on population density and heating demand. These two measures are used to restrict the calculated values to areas where there is a high need for utilization. For more information see Section 2.2.4.

2.2. Basic Data

2.2.1. Structural Data and 3D Geological Model

The structural data used for the calculations were taken from a recently published 3D geological model of the North German Basin (TUNB model, https://gst.bgr.de/ (accessed on 20 March 2021) [46]. The TUNB model was separately constructed by the geological state agencies of the Northern German federal states, and later on, merged into a consistent regional model for Northern Germany. The model resolves the base of 13 stratigraphic formations and their respective lateral extent. Base layers are available for the Rupelian (within the Tertiary), the Tertiary, the Upper Cretaceous, the Lower Cretaceous, the Upper Jurassic, the Middle Jurassic, the Lower Jurassic, the Upper Keuper, the Lower Keuper, the Upper Buntsandstein, the Middle Buntsandstein, the Lower Buntsandstein, and the Zechstein. The horizontal resolution varies as the structural model is meshed more coarsely in some states and more refined at state boundaries and fault traces. Therefore, the data were translated into the mesh nodes first and regularized using the commercial software Petrel^©. For the units’ surfaces, regular grids were created with a resolution of 1000 × 1000 m and clipped where the unit is not present. As outlined in the introduction, the data needed from the structural model are the reservoir thickness as well as the reservoir depth. The reservoir depth is needed for extracting the temperature from the respective datasets (see Section 2.2.3). A common practice is taking the reservoir mid as being representative for the calculation of the potential later on (e.g., [47]). Based on the regular grids, we calculated the thickness of the geological units and the respective middle depth of each stratigraphic unit, and then exported the results to csv files for further processing. For the nine formations investigated, we extracted the elevation data from the stratigraphic unit which they are part of. As an example, we show the depth of the middle of the Upper Keuper formation which shows an elevation distribution of the reservoir mid between 0 m.b.s. and 5500 m.b.s. (Figure 2b), where it is shallowest in the NE, SE, and S, and maximum depths are reached in the NW towards the North Sea. Generally, for the central model domain, the formation shows larger depths compared to the margin model domain.

2.2.2. Data for Assessing Potential Reservoirs

The general structural-geological model allows us to consider the main reservoir units associated with the different formations. To identify which units of the aforementioned large-scale structural model are, in principle, suitable for geothermal energy extraction and storage, we approached the task at hand in two different ways. First, we analyzed available bore logs from 15 wells to derive the sandstone fraction of each geological formation considered as a model unit (Figure 2b). The spatially distributed sandstone fraction of each unit is used as a reducing factor for the total formation thickness according to Equation (4):

$$thic{k}_{ss}=thic{k}_{tot}fra{c}_{ss}$$

(4)

As a result of this analysis, six stratigraphic formations were identified as suitable for either heat extraction or heat storage (sandstone fraction ≧ 0.2,

Table 2. Sandstone fraction for all analyzed units for the 15 selected wells. Unit index as indicated in Figure 2a, well locations in Figure 2b. Color-coding from red = low, over yellow = medium to red = high. Av. = average.

Well	1 + 2	3	4		5	6	7	8	9	10	11	12	13	14	15
1	0.00	–	0.10		0.04	0.25	0.01	0.44	0.11	0.00	0.54	0.00	–	–	–
2	0.23	0.00	–		–	–	–	–	–	–	–	–	–	–	–
3	0.18	0.00	0.02		–	–	–	0.00	0.04	0.04	0.32	0.00	0.00	–	–
4	0.10	0.00	0.00		–	–	0.00	0.29	–	–	–	–	–	–	–
5	0.25	–	–		–	–	–	0.12	–	–	–	–	–	–	–
6	0.45	0.00	0.88		–	0.40	0.35	0.15	0.05	0.00	0.10	0.00	0.03	0.30	–
7	0.36	0.00	0.88		0.00	0.58	0.35	0.85	–	–	–	–	–	–	–
8	0.00	0.00	0.65		–	–	0.38	0.17	0.27	0.05	0.71	–	–	–	–
9	0.36	0.00	0.00		–	–	0.84	0.00	0.00	0.04	1.00	1.00	0.00	0.04	0.12
10	0.40	0.00	–		–	0.50	0.34	0.13	0.07	–	–	–	–	–	–
11	0.19	0.00	0.00		–	0.52	0.18	0.14	0.18	0.00	0.17	0.00	0.00	–	–
12	0.75	–	–		–	–	0.46	0.21	0.00	0.04	0.29	0.08	0.00	0.00	0.17
13	0.50	0.00	0.58		–	0.70	0.53	0.20	0.00	0.00	0.37	0.04	0.00	–	–
14	0.33	0.00	0.00		–	–	0.62	0.03	–	–	–	–	–	–	–
15	0.45	0.00	0.73		–	0.30	0.56	0.09	0.14	0.17	0.66	–	–	–	–
Av.	0.30	0.00	0.35		0.02	0.46	0.39	0.20	0.08	0.04	0.46	0.16	0.01	0.11	0.15

). We exclude the uppermost (Cenozoic) due to potential usage conflicts with groundwater production, especially in urbanized areas. The remaining units are the Lower Cretaceous, the Middle and Lower Jurassic, the Upper Keuper, and the Middle Buntsandstein.

In addition, the thickness distributions from regional facies studies were considered (Sandsteinfazies, GeoPoNDD [30,48]). The authors of the latter studied the facies of several sandy reservoirs of the Bajocian and Upper Aalenian (both Middle Jurassic), Upper Toarc and the Lower Toarc (both Lower Jurassic), Upper Exter, Lower Exter III, and Lower Exter II (all Upper Keuper), and finally, the Upper Schilfsandstein, and the Lower Schilfsandstein (both Lower Keuper, Figure 2). In the following, we explain and apply the methodology for two units: (1) the Upper Keuper stratigraphic unit and (2) the Upper Exter Formation as one subsection of the Upper Keuper. We focus on this stratigraphic unit and this formation to cover a depth range both suitable for heat extraction and heat storage as well as to have a high coverage of the model area (Figure 2b,c). Arguably, the Middle Buntsandstein unit has a large lateral extent but it is located at greater depths making it less accessible for the targeted utilizations.

The facies distribution of the Upper Exter Formation, in turn, shows an intricate pattern of fluvial channels and their associated fans and bars (Figure 2c). These were analyzed by Franz et al. [30] for the thickness distribution of sandstone facies which reduces the extent and thickness significantly. The facies data were combined with the elevation data of the Upper Keuper unit for later processing (Figure 2b,c). The favorable areas are mainly located in the N, NE, and center, tracking from the central S towards the Baltic Sea in the NE.

2.2.3. Temperature Data

Temperature data from two different sources were used for the present study: (1) from a physics-based 3D numerical process model of the Central European Basin System (CEBS, named M1 hereafter, [32]) and (2) from the geostatistical computed 3D temperature interpolation of the collected borehole temperature database of Northern Germany (mapping: GeotIS, data: Fachinformationssystem Geophysik FIS GP, named M2 hereafter, [25]). The CEBS model represents a transient, that is time varying, 3D thermohydraulic model which resolves 16 geological units. It has been run in transient considering hydraulic loading (also including glacial loading for glaciated times) between 26 ka before present and 310 ka after present. It has a rather coarse resolution of 4 × 4 km but is the most up-to-date 3D thermohydraulic model that covers the entire model region. The CEBS model is in general agreement with available temperature measurements [32], but of unknown accuracy on the local scale. For the GeotIS data [25], thousands of borehole measurements were interpolated to a resolution of 2000 × 2000 × 100 m. We sampled both datasets at the elevation derived in the step before (see Section 2.2.1). As an example, we show the extracted temperatures for the Upper Keuper and Upper Exter formation (Figure 3). For the Upper Keuper, temperatures range from 8 °C to 154 °C for the CEBS dataset (Figure 3c) and range from 14 °C to 175 °C (Figure 3d) for the GeotIS dataset [25], whereby for both datasets, the lowest temperatures are found in the SE and NE, and highest temperatures are located in the NW. In comparison, temperatures for the Upper Exter Formation range from 13 °C to 140 °C (Figure 3e) for the CEBS dataset and between 20 °C and 141 °C (Figure 3f), with the lowest temperatures in the N and highest in the NW and SW. In general, the CEBS dataset has lower temperatures compared to the GeotIS dataset which possibly derives from the forced convective cooling component in the CEBS model. Compared with measured data, the CEBS model shows small deviations but is generally in good agreement [32].

2.2.4. Heat Demand and Population Density

For the heating demand, we used the data from the Pan-European Thermal Atlas 5.1 [50]. This atlas resolves the heating demand of Europe in a 100 × 100 m resolved grid. In the model area, the values range from 0 to 68 GJ/m²/a for the grid cells. For our analysis, we clipped the dataset to regions where more than 0.02 GJ/m²/a are needed and combined them into a 1000 × 1000 m area resulting in a minimum heating grid capacity of ~5 GWh/a, as this is assumed to be a suitable size for modern district heating systems. As can be seen in Figure 4, the remaining areas are numerous but amount to only ~40% of the model area, clustering around the major urban centers.

The population density data are acquired from a different dataset [51]. The latter has a resolution of 30 × 18 m and ranges from 0 to 16 people per grid cell. We aggregated those data to a regular 1 km × 1 km grid whose range is then 0 to 11,440 ppl/km². We then clipped the raster to areas with more than 500 ppl/km² and added a buffer of 1000 m again. The resulting area is only ~13% of the model area, clustering around the large urban centers.

3. Results

3.1. Hydrothermal Resources in the North German Basin

We present the results for the hydrothermal resources for the Upper Keuper stratigraphic unit and for the Upper Exter formation contained therein. Here, we compare the calculated values for the two different temperature models. Figure 5 shows that the predicted HIP for the Upper Keuper ranges between 0 and 3.7 GJ/m² for M1 (Median = 1.3 GJ/m², Table 2) and between 0 and 4.8 for M2 with distinct areas of high and low HIP. The identified sum for M1 amounts to 109,115 GJ (

Table 3. Statistical HIP values of the reservoirs. ¹ = Excluding areas of no potential, ² = percent of total area with potential greater 0. Underlined values represent maximum and italics the minimum value of investigated units for the respective parameter.

HIP (GJ/m2)	Lower Cretaceous	Middle Jurassic	Lower Jurassic	Upper Keuper	Middle Buntsandstein
Maximum	5.963	3.914	3.949	3.688	3.970
Median	0.0	1.169	1.180	1.296	0.0
Median 1	1.844	1.933	2.069	2.060	2.591
Mean	0.590	1.026	1.076	1.135	0.847
Mean 1	1.990	2.024	2.108	2.102	2.456
Area (km2)	1.20 × 103	4.78 × 104	9.17 × 104	9.93 × 104	1.48 × 105
Thickmax (m)	1973.0	1621.0	1827.2	820.8	3111.6
Sum (GJ)	69,170.9	47,078.1	95,002.7	109,115.2	122,327.5
Area% 2	29.626	50.712	51.056	54.012	34.485
HD Area% 2	11.888	25.452	26.363	27.732	12.801
PD Area% 2	3.662	7.485	7.353	7.342	4.039
HD Sum (GJ)	27,064.3	23,850.7	49,409.5	56,073.9	42,474.5
PD Sum (GJ)	8675.4	7,320.4	14,156.4	14,860.2	12,900.3

). In the central model domain, extensive areas of high potential are found wherein maximum values coincide with the position of salt rim synclines where increased thicknesses, depths and therefore, also temperatures are encountered. This trend is mainly dominant in the eastern model domain. In the western model domain, where salt walls are more frequently encountered than diapirs, elongated areas of high potential are situated in between the walls. Here, strong gradients in potentials are encountered due to the rapid spatial change in thickness, depth and temperature of the unit. There are also extensive areas where the unit is not present at all, denoted by the gaps. Towards the basin rims at the model area’s outer bounds, large areas of low to no potential are found due to the low depths and hence low temperatures. Comparing the results of the two temperature datasets, we can see that the area and values of the predicted HIP are larger for the M2. This derives from the higher average temperatures in M2 as compared to M1. While this results in higher average predicted HIP, the areas where there is any significant potential (>= 0.3 GJ/m²) also increase, most strikingly visible in the NE and SE.

In a next step, we superimposed the predicted HIP with available data on population density (PD) and heat demand (HD, Figure 5). Consequently, areas that are both favorable geologically and close to densely populated areas or areas of high heating demand are much smaller in size than the geologically favorable regions alone. Looking at the population density superimposed maps, only smaller regions in the model area remain which are located where the major urban centers are (e.g., Hamburg, Berlin). For Hamburg, while large areas fulfill both requirements, only the E and the S show significant potential (up to 3.6 GJ/m², Figure 5). The same is true for Berlin, where only locally in the S and W, some potential is predicted (up to 2.5 GJ/m², Figure 5). The maps focusing on the heating demand paint another picture (Figure 5). Here, larger regions, that is, many individual localities, fulfill the requirements. There is a clustering around the major urban centers, but it is not as clear-cut as for the population density. Looking at Berlin and Hamburg, comparable regions are identified, like for the population density. The areas are somewhat smaller and more individual regions could be identified. The maximum values remain the same as for the comparison before. When looking at the comparison of M1 and M2, the population density maps show similar trends for both models but the magnitude of the predicted values is again higher for M2. There are several localities where both parameters are high for M2 (PD, HIP) where only one is met for M1 (PD). The same is true for the heating demand maps. Generally, identified areas match closely, but there are numerous localities where both parameters are high for M2 (HD, HIP) but only one for M1 (HD).

We show the final results for the Upper Exter formation in the lower panels of Figure 5 to illustrate exemplarily the outcome for the sandstone facies data workflow, also because this formation is part of the Upper Keuper stratigraphic unit. Here we can see that the lateral extent of the identified reservoirs (compare Figure 2) is limited to the center and NE of the model domain, which is then further reduced when considering heating demand. In the center of the identified reservoirs where the channels reach their highest thickness (Figure 2), highest HIP values are encountered (Figure 5g,h). However, these are rarely located also in areas with high HD. Remaining localities where both a high HIP and HD are encountered, are located S of Schwerin and S of Hamburg, both in the focus of current geothermal drilling campaigns. Comparing M1 and M2, more areas have both high HD and high HIP, most prominently between Berlin and Schwerin, and NE of Hannover.

For the other stratigraphic units, we present the statistical values in the paper text (Table 3) and add additional map views in the Supplementary Materials. Here, the highest potentials of the HIP are found in the Lower Cretaceous unit with 5.963 GJ/m², while the lowest potentials refer to the Upper Keuper when looking at the maximum values. For the median and mean values of the HIP, maximum values of 1.296 GJ/m² and 1.135 GJ/m² are calculated respectively when taking the whole reservoir into consideration. In areas with proven potential, the highest potential is calculated for the Middle Buntsandstein for both parameters. Even though only 34.5% of the reservoirs’ area fulfills our criteria, the unit still has the highest cumulative potential with ~122.3 TJ. The highest percentage of the total area to qualified area is calculated for the Upper Keuper unit with 54%, which drops significantly when taking the heating demand into consideration (27.7%), and even more so when looking at the population density (Upper Keuper: 7.3%, Middle Jurassic: 7.5%). With these areas, the highest cumulative potentials are still located in the Upper Keuper with ~56 TJ for the HD and ~15 TJ for the PD. Looking at the overall values, all units exhibit high potential.

3.2. ATES Potential in the North German Basin

Similar to the hydrothermal resources, we present the results for the ATES potential, that is, the heat storage potential, for the Upper Keuper stratigraphic unit and for the Upper Exter formation contained therein. Here, we compare the calculated values for the two different temperature models. Figure 6 shows that the regions with high potentials are less present compared to the areas of enhanced potential in the HIP. In detail, areas of high HSP are located in the NE, SE, and S. Subordinate, small areas are located in the central model domain. Their location correlates with the basin edges as the stratigraphic unit is shallower there. This is also true for the areas of high potential which are not located on the basin edges but rather in regions where the salt diapirism lifted the Upper Keuper unit, making it more accessible, both from the depth level and the reservoir temperature. Values range from 0 to 0.97 GJ/m² for M1 (Median 0 GJ/m², Figure 6) and from 0 to 0.91 GJ/m² (Figure 6). As indicated by the median, only less than 10% of the unit’s area has been identified to be favorable for storage of heat. Comparing the two temperature datasets, the results show that M2 predicts lower storage potentials than M1. The overall area suitable for heat storage is also smaller for M2, meaning some areas identified to have storage potential in M1 show no potential in M2. This derives from the fact that the average temperature of M2 is higher, which results in more areas not fulfilling the requirements for heat storage as compared to M1.

When overlaying the HSP of the Upper Keuper with the population density and heating demand, the remaining favorable area is significantly reduced. Similar to the HIP, large parts of high potential areas do not qualify (e.g., NE and SE area of high potential), especially the local highs associated with salt diapirism in the central model domain. Looking at Hamburg, no favorable areas remain for this geological unit (0 GJ/m² max, Figure 6), while in Berlin, high potentials are predicted for large parts of the city (center, NE, 0.9 GJ/m² max, Figure 6). In comparison, when looking at the HSP superimposed with the heating demand, a significant reduction in the potential areas is observable again, however, numerous areas remain. For Hamburg, no potential areas were identified again (0 GJ/m² max, Figure 6), while Berlin shows promising regions for large parts of the city’s area (center, NE, 0.9 GJ/m² max, Figure 6), and also in its direct vicinity. The comparison between the different temperature datasets reveals a similar trend as before. While there is generally good agreement in the results, there are several areas where both HSP and HD or HSP and PD are high for M1 while only HD or PD are high for M2.

Similar to the HIP, we only show the final results for the Upper Exter formation HSP calculations in the lower panels of Figure 6 to exemplarily illustrate the outcome for the sandstone facies data workflow, because this formation is part of the Upper Keuper stratigraphic unit. Again, large areas are unfavorable, as indicated by the range from 0 to 0.81 GJ/m² and a median of 0 GJ/m² (Figure 6). Here, the highest potentials are located in a larger area in the N, while there are some local highs in several locations associated with local highs in the reservoir’s elevation above Zechstein diapirs. As for the results above, M1 predicts higher values in larger areas for the HSP as compared to M2.

Similar to the HIP, we describe the heat storage potential (HSP) of the other reservoirs only briefly, focusing on the statistical values (

Table 4. Statistical HSP values of the reservoirs. ¹ = Excluding areas of no potential, ² = percent of total area with potential greater than 0. Underlined values represent maximum and italics represent the minimum value of investigated units for the respective parameter.

HSP (GJ/m2)	Lower Cretaceous	Middle Jurassic	Lower Jurassic	Upper Keuper	Middle Buntsandstein
Median	0.000	0.000	0.000	0.000	0.000
Median 1	0.788	0.786	0.821	0.754	0.822
Max	0.988	0.979	0.977	0.961	0.988
Mean	0.094	0.132	0.147	0.072	0.034
Mean1	0.761	0.775	0.818	0.719	0.813
Sum	11,079.1	6057.3	12,951.1	6878.7	4878.1
Area% 2	12.4	17.0	17.9	9.9	4.1
Area HD% 2	6.4	8.7	9.1	5.4	2.1
Area PD% 2	2.3	3.1	3.2	2.0	0.7
Sum HD	5789.0	3126.5	6587.9	3728.6	2499.0
Sum PD	2177.9	1123.1	2302.3	1398.6	870.1

). Here, all units show a median value of 0 GJ/m² due to the fact that a large proportion of the units’ area is unfavorable due to it being either too hot, too shallow or too deep. When looking at the median of those areas with a potential higher than 0, the Middle Buntsandstein unit has the highest value (Table 3, 0.822 GJ/m²) which derives from the fact that in areas where the unit is comparatively shallow, it also displays low temperatures. The median storage potential of all other units is only slightly lower due to the strict restrictions on valid areas. The same is true if we look at the mean predicted HSP for areas with a potential greater 0 GJ/m². Looking at the overall average, the highest average HSP of 0.147 GJ/m² is predicted for the Lower Cretaceous unit, while the lowest average HSP of 0.034 GJ/m² is predicted for the Middle Buntsandstein. This was derived from the former displaying a valid area percentage of 18% and the latter only showing a validity percentage of 4%. Hence, the highest accumulated potential of ~13 TJ is also predicted for the Middle Jurassic. Comparing this with the HIP, we predict a maximum value which is one order of magnitude lower. Looking at the HSP for areas of high heating demand and population density, the same trends are predicted. The Middle Jurassic displays the highest percentages of valid areas and cumulative HSP, while the Middle Buntsandstein has the lowest values for both parameters. Comparing the HIP and HSP predictions, the deeper reservoirs are more favorable for HIP and the shallower units are more favorable for HSP.

4. Discussion

The heat in place (HIP) and heat storage potentials (HSP) outlined here describe geothermal reservoirs and the associated depth levels that are most promising for a technical utilization. For the HIP, a clear correlation between the reservoir depth, reservoir temperature, and the predicted potential is evident, meaning that the deeper, hotter reservoirs hold the highest theoretical potential. Here, the combination of temperature data and detailed facies analysis of the reservoir sandstones are powerful for further in-depth exploration. For the HSP, an inverse trend is predicted, where the shallower and thus colder the reservoir is, the higher the storage potential. Both of these relationships derive from the fact that the temperature difference between the reservoir and the injected fluid has a major impact on the potential. For the HIP, the hotter the reservoir fluid, the more energy can be extracted since we assume a fixed reinjection temperature. For the HSP, the colder the reservoir fluid, the higher the contrast to the injected fluid, which we assume to have a constant temperature of 90 °C, and the higher the predicted potential.

To illustrate which areas hold at least one stratigraphic unit that exceeds 1 GJ/m² for the HIP and 0.3 GJ/m² for the HSP, we plotted the results for M1. Figure 7 then illustrates that most parts of the study area are covered by at least one studied unit exceeding the HIP threshold. In the central model domain, high HIP is predicted in the Lower Cretaceous already, being the youngest and shallowest unit of this study, also underlain by the other stratigraphic units. Moving further towards the basin margins, the number of units exceeding the threshold decrease, resulting in the only unit with considerable HIP at the basin margin being the Middle Buntsandstein. The NW of the model domain shows a different trend where the Upper Cretaceous is the only unit with high HIP in the very NW followed by the Middle Buntsandstein. For the HSP, only the SE over SW to NE hold any unit of high potential. Here, the E shows predominantly high potentials in the Jurassic and Cretaceous, while the S also shows some areas where the Middle Buntsandstein is the only unit suitable for HSP.

4.1. Improvement and Comparison to Previous Work

Compared to previously published results for Germany [26,46], we can demonstrate that regions with high geothermal potential are a lot more widespread (Figure 7 and Figure 8). The general area identified for any level of geothermal utilization, low or high enthalpy, has also significantly increased. We want to stress, however, that the datasets we compare here, while being comparable in nature, have different input and output data. Identifying geothermally promising areas based on the aquifer temperature alone only has limited quantitative and qualitative value. The predictions from this study not only use an updated background structural model [31] but also updated temperatures [32]—also on a pan-national level. We also make conservative calculations of the energy contained in the reservoir units. With this approach, a more evenly spread distribution of areas with high potential is calculated. Future updates should be based on a physics-based simulation of the NGB, implementing the new structural data [31] refined by further reservoir facies analysis (e.g., Mesotherm project, [52]), and a spatial parameter model adjusted to the local in situ conditions at reservoir depths [53].

As a further step, we calculated the storage potential of all reservoirs, which is done here for the first time. The areas with high potential identified this way partly overlap with those from the heat in place calculations (Figure 8), but also new areas emerge where the heat in place is low, but the storage potential is high. This identification is a first step in differentiating regions that can either be used for storage alone, for both storage and heat extraction, or only for heat extraction. Comparing promising heat storage areas to each other, for example, reveals that these are almost exclusively located where earlier studies have failed to identify any potential (Figure 8, [46]).

In a global comparison, the North German Basin shows extensive areas with high geothermal theoretical potential. The study by Limberger et al. [20] already showed that the North German Basin displays HIP values of up to 150 GJ/m², which is approximately 30 times more than our predictions. This derives from the fact that the study makes some rough assumptions about the thickness of the reservoir under study, which we consider very conservatively, probably more realistic. Additionally, we consider conservative cutoff values for the depth under consideration. For the HIP, we consider everything at greater depth than 2500 m.b.s. to be of increased cost and risk [34], and therefore disregard it. We apply a similar threshold for the HSP of 750 m.b.s. These cutoff values were chosen exemplarily but do not represent any technical boundary. We merely adopted these values to consider conservative depth ranges which likely represent lower risk at the exploration stage and also come with lower costs. In regions of exceptionally high geothermal gradients and respective temperatures at depths higher than 2500 m.b.s., these increased costs and risks might, however, be feasible and should be investigated further.

4.2. Data Availability and Quality Limitations

In order to calculate the heat in place and storage potential, a lot of different factors have to be considered. The work presented in this paper tries to include as many of those as possible and necessary, while also making some assumptions. With this approach, we want to show what is possible with limited data availability while still showing the general trends in both output parameters.

The underlying structural model is based on an extensive database, effectively including available seismic, gravitational, geological, and other observational data [31]. It is the most detailed geological-structural representation of the NGB up to date. It is, however, limited by its vertical stratigraphical differentiation, the data availability and, therefore, comes with some uncertainty in the elevation distribution of the mapped horizons, and thus also the derived thickness distributions. The thickness data acquired from the sandstone facies projects, while also based on a large number of borehole data points, sedimentation concepts, and related interpretations, have some uncertainty in the extent of the mapped horizons as well as their vertical position. This is evident by sandstone deposits with several meters of thickness in this dataset, where there is no mapped thickness of the containing unit in the TUNB model. An example of this would be the Upper Exter Formation vs. the Upper Keuper unit containing the latter (Figure 9). A matching of both approaches to solve these obvious discrepancies between structural and facies interpretation would be required to reduce uncertainties of any future resource assessment based on these data (Appendix B, Figure A2). A first attempt at this has been published in Franz et al. [48] but an update with the newly available structural model from BGR et al. [31] was necessary and partly carried out for this paper.

Moreover, the temperature data used for the calculations are based on geostatistical interpolations of borehole temperature data [46] and on numerical process modeling at a pan-national scale [32] which introduces inherent uncertainties. The choice of the temperature input data becomes apparent when comparing the results of the different datasets. Here, an emphasis on relying on the most up-to-date temperature predictions and measurements should be given, and also the predictions of the geothermal potentials should be updated as soon as new data and temperature models are available. Since the results presented in this paper are most powerful in identifying trends and regions, any local study should undergo a crucial examination of the temperature predictions in this area. As stated in the Introduction, this paper does not consider over- or under-burden rocks and their effects on the HIP [11,12]. Here, the estimates given in the Results section of this paper would likely be higher on average since a reheating of the reinjected fluid would increase the potential energy retained in the reservoir and should be considered in more detailed studies of the regions identified in this paper.

We consider fixed values for the heat capacity of the fluid and solid, the density of the fluid and solid, and a fixed porosity. These assumptions were made to keep the models as simple as possible and not to introduce unnecessary uncertainties. There are similar potential studies using depth-dependent values for the estimation of these parameters [6,23], however, due to the vastly different lithologies of the different reservoir units [30] and the complex nature of fluid pathways [17], a simple depth relationship of these parameters might lead to erroneous results. Another possibility to investigate the effects of different realizations for the input physical parameters would be to run a Monte Carlo simulation (e.g., [54]) where parameter ranges could be tested. This statistical approach would improve the lower threshold predictions even further (low likelihood members could be even more conservative), however, such an approach applied to the size of our study area would only be of limited informational value. Here we want to emphasize, that the implementation of the regional parameter variations, whose significant effect for temperature predictions is known from [55,56,57], would require a far more detailed knowledge of the basin filling and its parameter variation, and remains for future work. For the regions identified with the tool presented in this paper, these parameters should therefore be included when calculating the geothermal potential, given detailed lithological and geochemical information.

Due to the nature of this paper to identify promising geothermal targets, we do not consider parameters such as plant design and other technical aspects of the final geothermal installation. For the HSP, we merely introduce a recovery factor of 0.7 based on [44] for higher temperature ATES, which represents a good practice average value for well-designed ATES systems.

5. Summary and Conclusions

This study shows for a first time:

• Basin-wide heat storage potentials (HSP) and heat in place (HIP) estimates for the North German Basin;
• Quantitative evaluation of HSP and HIP in regions of high heating demand and high population density in Northern Germany;
• Recommendations for conservative geothermal potential assessment.

The major conclusions are as follows:

• The potential for geothermal heat production from depths lower than 2500 m covers the societal demand for Germany by orders of magnitudes;
• The potential for subsurface heat storage covers the demand for mid- to long-term large-scale storage solutions for heat or cold;
• The selected workflow and database reveal more regions with high geothermal potential than previous studies;
• Future improvements in the geothermal assessment of the North German Basin requires a detailed physics-based numerical simulation, based on a profound understanding of the subsurface parameter fields.