Abstract
Knowing the thermal and hydraulic conditions below major urban centers is of increasing importance in the context of energy and water supply. With this study, focusing on the major urban center of Berlin, Germany, we aim to gain insights on the coupling of surface water bodies to the subsurface thermal and hydraulic field investigating shallow water to deep groundwater interactions. Therefore, we use a 3D structural model of the subsurface, constrained by all available data and observations, as a base for simulations of the coupled transport of fluid and heat. This model resolves the 3D configuration of the main geological units and thus enables us to account for related heterogeneities in physical properties. Additionally, we resolve surface water body geometries with newly available data. To assess how surface water bodies interact with the deeper groundwater at different depths in the model domain, the influence of different hydraulic boundary conditions is quantified, which indicates that the coupling of surface water bodies and groundwater strongly modifies predicted groundwater circulation. Consequently, changes in subsurface temperatures are also predicted, where lakes may account for temperature differences up to ±5°C and rivers could account for up to ±1°C visible at depths ≤-500 m.a.s.l. These differences are mainly connected to changes in the advective component of heat transport caused by the modifications of the hydraulic boundary condition. Pressure-driven heat transport is most efficient where differences between hydraulic heads of aquifers and surface water bodies are highest. This study therefore illustrates the impact of surface to subsurface water interactions in an urban context.
1. Introduction
In the last decade, the need for big cities to make efficient, environmentally friendly use of the subsurface on which they stand is increasingly getting attention [1]. In order to properly engage in this endeavor, a systematic understanding of the distribution of available natural resources, like groundwater and geothermal energy in the subsurface, is a prerequisite. One of the major limiting factors preventing both an efficient and economically viable planning of georesource utilization derives from a poor understanding of the hydrogeological conditions, including the geological configuration and hydrodynamics below urban areas. This study presents a 3D workflow integrating geological and hydrogeological data, numerical modeling, and available data from monitoring activities into physically based 3D models as a tool to assist in the development of good practices with respect to usage and maintenance of groundwater and geothermal resources. Herein, especially the configuration of groundwater and surface water as boundary conditions (BCs) in these models is the focus of the study, since their realization is often simplified in modeling studies of this type as outlined in the following chapters. We take the example of the city of Berlin (Figure 1(a)), capital of Germany, as a natural laboratory where to develop and test our approach. The model area therefore stands representative of a major urban center with more than four million inhabitants, with largely reshaped surface morphology and shallow hydraulics, also depicting a complex hydrogeological setting.
(a)
(b)
The study area is located in the Northeast German Basin (NEGB), a subbasin of the larger Central European Basin System [2]. Previous 2D and 3D studies focusing on the basin-wide hydrothermal configuration of the NEGB have revealed the presence of a regional hydrothermal regime driven by different physical processes. Those processes include heat diffusion, which is active especially at great depths (below the Zechstein salt), and its coupling to a regional component of pressure-driven groundwater flow within the Mesozoic to Cenozoic aquifers (i.e., [3–5]). This component of fluid-related heat transport has been shown to be mainly controlled by the (1) subsurface porosity and permeability distribution and, of special relevance for this study, (2) the surface hydraulic configuration, imposed in the models as an upper hydraulic BC. In this regard, the major limitation common to previous studies was to make use of a constant hydraulic head at the topographic level, which arguably imposes unrealistically high hydraulic gradients supposedly leading to an overestimation of advective cooling [5–7]. This in turn has a crucial impact on the resulting shallow to deep groundwater circulation, detailed in the deeper infiltration of meteoric water into and uprising of deep-seated water out of the model domain. Therefore, a proper quantification of the urban-wide hydraulic and thermal configuration requires a better, more realistic approximation of the hydraulic BC. In the context of the model area, this is also of crucial interest since most of the fresh water resources for the city of Berlin are produced from shallow aquifers [8] and might be subject to mixing of waters of different ages and salinity content, which might impact water quality (Figure 1(b)). Hence, the investigation of the controlling mechanisms of subsurface hydrodynamics is crucial to first of all reproduce the observed hydraulic and thermal regimes (present-day status) and make reliable predictions of the dynamic response of the system to further anthropogenic forcing (subsurface utilization).
Recent studies [9–11] have demonstrated that the hydraulic and thermal configuration at shallow and intermediate levels (≥-1000 m.a.s.l.) of the study area is controlled to a large degree by the elevation (top and base, Figure 2(a)) and thickness distribution (hereafter called structural configuration, Figure 2(b)) of the Rupelian aquitard. This unit separates fresh water (above) from saline water (below), but it also experienced glacial erosion, resulting in holes through the root of its base, hereafter referred to as hydrogeological windows (Figures 1(b) and 2(b)). Because of its discontinuous nature, former studies focused on quantifying the degree of hydraulic connectivity between the supra- and sub-Rupelian aquifers in the vicinity of these hydraulic windows (i.e., [3, 12]). These studies showed that hydrogeological windows likely act as preferential pathways for deeper groundwater to upcone to the surface and for shallow (meteoric) water to penetrate deeper into the lower aquifers. The depth extent of this mixing zone, as well as the vigor of this process, has been shown to be controlled by the hydraulic configuration at the surface and by the hydromechanic characteristics of the relevant aquifers and aquitards [9, 10]. These studies have highlighted a causal relationship between the vigor and depth extent of the mixing zone and (1) the pressure hydraulic surface forcing (i.e., groundwater heads) and (2) the hydrogeological configuration and spatial configuration of the Rupelian aquitard. For (1), a careful representation of the hydraulic BC is key and the progress in this regard is explained in more detail in Sections 3 and 4.4.1. For (2), a possible leakage through the clay layer in areas where no hydrogeological windows are present is predicted, which stands in contrast to earlier model findings [5, 11]. Another recently published study [13] showcases how an in-depth knowledge of shallow subsurface hydraulic and thermal conditions is of utmost importance for any planned geothermal utilization, thus outlining the importance of the work presented in this paper. Scheck-Wenderoth et al. [13] show that predictions highly rely on the spatial scales considered, stressing that local, highly refined models (like the one presented here) are necessary for robust results.
(a)
(b)
(c)
(d)
(e)
(f)
With this study, we want to investigate deep to shallow interactions as represented in different model scenarios, focusing on the changes in local hydrodynamics resulting from different, increasingly more realistic representations of the upper hydraulic BC. We opt for a most accurate implementation of the hydraulic BC by utilizing a wider database for the groundwater heads compared to previous studies [9] as well as lakes and rivers as new elements, based on newly available data. Hence, we are able to make qualitative and quantitative predictions of changes in hydrodynamics induced by each of these elements. We make use of modeled temperature as a passive tracer which serves as an indicator for the depth and lateral extent of modifications in advective heat transport. Predicted temperatures and respective changes are also used to determine shallow geothermal potentials.
2. Geological Setting
The sedimentary succession in the model area ranges from Permian (Sedimentary Rotliegend) to Neogene (Holocene) in age and consists predominantly of clastics, carbonates, and rock salt [2]. This sedimentary succession is underlain by variably thick volcanics which are Permocarboniferous in age [14]. The Mesozoic evolution of the basin has been largely controlled by halokinetic mobilization of the Permian (Zechstein) salt rock layer, which occurred during several stages starting from the late Middle Triassic onward [15]. This process largely influences the geometry of all overlying units, as illustrated in Figure 1(a). Within the study area, the Zechstein salt layer shows a heterogeneous thickness distribution, with local maxima exceeding 3500 m, envisaged in the NW and E of the model area, and domains of reduced thickness or complete withdrawal, as visible along the southern boundary of the model near the Tempelhof area in the S (Figures 2(e) and 2(f), [15]).
The Mesozoic sedimentary sequences consist predominantly of consolidated clastics or carbonates [16]. The succession also includes the Middle Triassic (Middle Buntsandstein) and Permian Sedimentary Rotliegend, which consist of sandstones consolidated to varying degrees, displaying sufficiently high porosity and permeability distributions, to be considered as target horizons for deep geothermal energy extraction (Table 1, [17]).
The Cenozoic sedimentary succession is predominantly composed of unconsolidated clastics (Figure 1(b)) overlying older strata in a discordant pattern (Figure 1(a)). The basal unconformity is expressed by Cretaceous, Jurassic, or Keuper sediments directly underlying Upper Paleogene or Neogene strata. Moreover, Paleogene sediments display inhomogeneous thickness distributions, as envisioned by the Oligocene Rupelian clay. The initial thickness distribution of this unit was approximately 80 m but has been altered significantly, as manifested in numerous glacial erosional channels and discontinuities (Figures 1(b) and 2(b)).
3. Hydrogeological Setting
The hydrogeology of Berlin has been studied intensively for several decades now, yet a number of open questions remain, owing to the fact that information about the structure and hydraulic parameters is sparse, especially considering the deeper subsurface (i.e., [8, 11, 24–27]).
In this respect, the lower constraint is the crystalline basement which is considered impermeable [5]. The overlying Permian strata are low permeability volcanics (Permocarboniferous), considered as the aquitard, and strongly compacted sediments with comparatively high permeability (Rotliegend), considered as the aquifer, making up the very deep groundwater compartment (Table 1). These are overlain by impermeable evaporites (Zechstein), separating them from the shallower groundwater compartments.
The intermediate groundwater compartment is described by mostly clastic Triassic sediments of varying degrees of sand-, silt-, clay-, and limestone separating them into one lowly permeable unit at the bottom (Lower Buntsandstein), one highly permeable unit in the middle (Middle Buntsandstein), and two lowly permeable units on top (Upper Buntsandstein and Lower Muschelkalk, Table 1). This compartment is separated from above by a layer of tight evaporites (Middle Muschelkalk), considered as a regional aquitard located in depths between 0 and -2400 m.a.s.l. with an average thickness of 90 m [28] (Table 1 and Figures 2(c) and 2(d)). The remaining Mesozoic units show increasing permeabilities from bottom to top, as younger units have experienced less compaction and show increasingly larger proportions of sand and silt in comparison to marl and clay (Table 1).
The hydrogeological configuration of the Cenozoic sedimentary succession is depicted schematically in detail in Figure 1(b). The base is described by a medium permeability unit (Pre-Rupelian) with a brackish to saline pore fluid [29]. This characteristic is also common for all older strata, supposedly connected to the dissolution of rock salt from the Zechstein unit [29–31]. The only exception to this rule is the formation waters of the Rotliegend which are mainly preserved with their original signature of seawater evaporation as the main source [32]. The saline formation waters of the Post-Rupelian are then separated from above to a certain degree by the lowly permeable Rupelian clay, a local aquitard of discontinuous nature, depicting hydrogeological windows (Figures 2(a) and 2(b) and Figure 1(b), [29]). Overlying this strata, five geological units in total are classified as aquifers (Cottbus, Miocene, Elsterian, Saalian, Holocene), only interrupted locally by the Holstein strata (Figure 1(b)) which are classified as an aquitard (Table 1, [33]). Figure 1 clearly shows that these aquifers depict a high level of heterogeneity, deriving mainly from the depositional or erosional character of the glacial and interglacial periods. Here, especially the glacial erosional channels from the Elsterian period are worth mentioning since they locally cut through the base of the Rupelian aquitard (Figure 2(b)). These channels are mainly filled by sand-dominated deposits, thus presenting possible pathways for fluid exchange between the different compartments. We differentiated the Cenozoic strata according to their geometrical relevance, which means that only strata of a certain thickness and spatial extent were resolved in the structural model (Figure 1(a)). However, for each of these strata, lithological distributions derived from the available well logs were analyzed and taken into account (more details in Frick et al. [9]), classifying the units as outlined above and in Table 1, thus resolving the heterogeneity of the units and the main aquifer to a certain extent.
The Cenozoic succession also includes the main aquifers (Cottbus, Miocene, Elsterian, Saalian, and Holocene) for drinking water production (Figure 1(b)). Therefore, the realization of aquifer groundwater levels and surface water body levels is of especial interest for this study. For the model area, earlier studies and newly available datasets suggest a natural direct connection between these two (Figure 1(b)) [34, 35]. Here, depending on the locality, either effluent or influent conditions are observed, depending in turn on local geology as well as anthropogenic overprinting due to groundwater pumping activities. Concerning this specific configuration, predecessor studies focusing on different aspects of the subsurface derived the hydraulic potential solely from measured groundwater data [9, 36] or topographic levels [11]. However, the data provided for this study show that there is a significant difference between interpolated groundwater levels used as a hydraulic BC in the studies mentioned above and surface water levels. In detail, interpolated values over- or undershoot hydraulic potentials in areas where data coverage is sparse, high gradients emerge between two groundwater head measurements, or anthropogenic overprinting leads to a drop in observed heads. As outlined in Section 4.4.1, the scenarios investigated here opt for a more realistic representation taking into account measured surface water heads.
4. Model Development and Scenarios
4.1. Refinement of the Model Surface
To understand the coupling between surface and subsurface water in the presence of major water bearing bodies (i.e., lake and rivers), new data have been integrated into the structural model as described in the previous section [9]. In a first step, we derived a new top structural surface (i.e., Top Holocene) as based on a high-resolution digital elevation model (DEM, grid m, Bundesamt für Kartographie und Geodäsie [37], Figure 3(a)). It represents the water table of all surface water bodies, where present, or the surface of the solid earth.
(a)
(b)
The topography of the model area is typical of a postglacial landscape, which has been formed during the latest glacial and interglacial periods as evidenced by extensive glacial spillways with meltwater valleys, low fluvial terraces, and periglacial basins (i.e., Greater Spree Valley in the center of the model area; see Figure 3(a)). The areas of higher elevation in the NE and S consist of ground and terminal moraines as well as young and old drift plateaus. Overall, erosional and depositional processes have shaped the area in a cyclic pattern associated with the different glacial and interglacial periods [38]. These processes are also well represented in the lithological distribution (Figure 1(b)) and structural configuration (Figure 2(b)) of model units.
To derive the geological surface of the model area, the volume of the water bodies was subtracted from the DEM. This was done by compiling long-term average data for lake water depth and river depth (data provided by MRDEA, SenStadtUm, and WSA; see Acknowledgments). To derive the absolute depth values of the lakes, available contour maps of water depths were interpolated to regularly gridded surfaces within the supplied bounds (Figure 4(b)). These depths were then subtracted from the absolute elevation values of the DEM, thus resulting in the surface horizon (Figure 4(a)). In areas where rivers are present, we derived this surface by combining the newly created lake bottom data (outlet of lakes) with river water depth measurements. Since all investigated rivers are connected to at least one of the lakes, in- and outlet geometries (Figure 4(c)) were used as starting or endpoints for the rivers. To derive the bottom surface of the rivers, an average width buffer was applied to the available centerlines while respective widths have been derived from satellite data (Figure 4(c)). The created outlines of the rivers were then assigned interpolated elevation values between the starting and endpoint, keeping true to an elevation drop of the top surface downriver. The elevation data for the center line were derived by subtracting measured or interpolated water depth data from the absolute elevations of the outline. Finally, the geological surface was constructed by interpolating between the outlines and the centerline (Figure 4(c)).
The resulting geological surface (Figure 3(b)) was then integrated into the structural model (Figure 1(a)) modifying the underlying interfaces accordingly.
4.2. Assignment of Physical Properties
For the numerical simulations carried out in this study, each of the geological layers has been assigned constant values for thermal and hydraulic properties according to Table 1. These physical properties were either derived from measured values or based on literature data according to the dominant lithology. Special attention was given while parameterizing the partly discontinuous Rupelian clay. Assuming all unconformities to be the results of glacial erosion, we have parameterized these areas by giving the same set of parameters as the unit locally overlying the Rupelian. All Cenozoic units were parameterized according to the petrological information supplied in the well database [1]. These lithologies were combined with literature data in order to derive values for the physical properties (Table 1), the exact procedure of which is explained in detail in Frick et al. [9]. Heterogeneities in lithological distributions of the different geological units have therefore been taken into account as part of the parameterization. This approach was chosen since resolving those heterogeneities geometrically would have led to a drastic increase in model elements and therefore simulation time, with little influence on the model result. We chose not to integrate any hydraulic conductivity anisotropies for the study area (), since no data on this property are available and a general analysis of the influence of this parameter is not the scope of this study. However, fundamental work on the topic of geological entropy and the associated field have been published recently and should be considered in future studies of the region, especially in studies opting for solute transport [39]. Since no measured data about the storage coefficient are at our disposal, we used the FEFLOW© default value of 0.0001 (1/m).
4.3. Modeling Method
The structural information described above has been integrated into a 3D structural model after Frick et al. [9]. This model differentiates 18 sedimentary basin fill units (eight of which are Cenozoic in age, Figure 1(a)) and two units for the underlying nonsedimentary units, namely, Permocarboniferous Volcanics and Basement (comprised of Pre-Permian and Upper Crust). The model is part of a series of studies, investigating different aspects of the geological subsurface. In this, models are becoming increasingly more local, thus increasing in the resolution and physical processes considered. These models rely on an extensive database, both in geological information as well as in physical parameters. The study presented here, therefore, makes use of these 3D gravity constrained, calibrated models and builds on them (i.e., [4, 5, 11, 14, 40]; also refer to the supplementary material (available here)). The vertical extent of the model is limited to an elevation of -6000 m.a.s.l., thus integrating all 18 sedimentary units overlying the two nonsedimentary layers, the lowermost one of which has been clipped at -6000 m.a.s.l.
In order to investigate the effects of different setups of the surface hydraulic conditions, we use coupled heat and fluid transport simulations. While the main subject of this study is a quantification of the groundwater dynamics, we also use modeled temperatures as a passive tracer to analyze predicted variations in the hydraulic setting as induced by the model modifications. For this purpose, the 3D geological model was imported into the commercial software package FEFLOW© [41]. The latter solves for (un)saturated groundwater flow in porous media taking into account conductive, advective, and buoyant (fluid-density-related) heat transport processes within a finite element-based computational framework. The mathematical formulation to solve for those processes is based on a system of three coupled equations, expressing the conservation of fluid mass, momentum, and energy, which under the Boussinesq assumption reads as where is the storage coefficient, is the pore pressure, is time, the mass density of the fluid, is the specific discharge (Darcy’s velocity), is the porosity, and is the sink-source mass term.
The specific discharge is given by where is the volumetric flux (Darcy flux), is the permeability tensor, the fluid viscosity, the hydraulic head, the gravitational acceleration, and the reference mass density of the fluid (999.793 kg/m3).
We also use temperatures predicted by the models as a further means to better quantify and visualize the induced changes in the regional hydraulics (e.g., [42, 43]). In this respect, salinity as another tracer was considered to be of secondary relevance for the specific goal of the present study. The paucity of available data which, at the current stage, will not permit an investigation of its impact on the city-wide groundwater and geothermal configuration also influenced this decision. The energy balance equation then reads as where is the specific heat capacity of the system with a fluid and solid phase. is the temperature, is the specific heat capacity of the fluid, the thermal conductivity tensor for the solid and fluid phase, and the heat source function.
The resulting system of partial differential equations (equations (1) to (3)) is closed by setting an Equation Of State for the fluid density (as a function of the system variables) after Blӧcher et al. [44], thus introducing a nonlinear coupling among the equations. More details about the mathematical background and its numerical formulation can be found in Diersch [41].
The horizontal resolution of the model is m (derived from an average mesh triangle area). In detail, the mesh consists of 175856 triangular elements, wherein the refinement of locations where rivers, lakes, and also Rupelian windows are encountered is four times higher compared to the surroundings. The refinement in areas where rivers and lakes are encountered depicts an increase in resolution compared to [9], where rivers and lakes are not yet resolved. To guarantee a good vertical to horizontal element shape ratio, the original model configuration, consisting of 20 geological units (Figure 1(a)) with distinct physical properties, was further refined, so that the final structural model is composed of 56 computational layers. Here, a minimum of two slices per geological unit was realized, whereas new slices were integrated in an equidistant manner between the geological surfaces from the structural model. The resulting numerical model consists of 9847986 elements. Due to the high nonlinearity of the coupled problem, models were run in a transient state for both fluid and heat transport, until reaching quasi-steady-state conditions. The latter were reached after approximately 250k yrs simulation time. The initial time step length was set to 10-3 d, and the maximum time step length was d.
4.4. Boundary and Initial Conditions
Equations (1) to (3) represent an initial and boundary value problem. Therefore, in order for the problem to be well posed, a proper set of initial and boundary conditions must be assigned. We consider all lateral boundaries closed to both fluid and heat transfer. This approach allows for a stable solution but also requires all sources and sinks of heat and fluid to be located within the model boundaries and that no heat or fluid is transferred to neighboring areas beyond the boundaries of the models, the shortcomings of which have already been discussed in Bayer et al. [2]. However, the region under study is located in the middle of the North German plain within an intracratonic area undergoing slow deformation without being surrounded by mountainous regions at its borders. Therefore, cross-boundary flow is likely to be very small and was neglected in this study. Given these conditions, the prescription of the hydraulic load (with prescribed preferential flow areas at the lakes and rivers) as varying across the study area provides enough constraints to capture the groundwater dynamics at the scale of interest. Additionally, the model scenarios described in the following depict fluid and heat budgets within reasonable bounds ().
4.4.1. Hydraulic Boundary Conditions
In order to understand the hydraulic forcing due to rivers and lakes and related changes in the thermal and hydraulic budgets, we set up three different model scenarios. The first of these scenarios (reference model, M1) implements a set of BCs which were parameterized as first order, Dirichlet-type across the top surface of the model domain. We impose fixed hydraulic head values as derived from measured groundwater levels only (Figure 5(a)) for all respective nodes of the uppermost slice of the model. With respect to previous studies [9, 10], we make use of a larger and denser database consisting of more than 2000 wells (data provided by the Department for Urban Development and the Environment of Berlin (SenStadtUm) and Ministry of Rural Development, Environment and Agriculture of the Federal State of Brandenburg (MRDEA), Figure 5(a)). These data represent the recent measured groundwater heads for the main shallow aquifer in the model area (Figure 1(b)). In comparison to the predecessor studies, a more accurate representation of this surface is achieved, whereas differences in elevation are as large as +16.2 m (-12.5 m).
(a)
(b)
(c)
(d)
For this scenario, the hydraulic head BC implemented derives solely from interpolated (convergent interpolation) measured groundwater heads which differ substantially from surface water heads, deriving partly from an overexaggeration of depressions in areas overprinted by anthropogenic effects (e.g., Lake Tegel area, groundwater pumping, Figures 5(a) and 5(b)). Groundwater heads for this scenario range between 24.5 m and 77.2 m, whereas large heads are observed at topographic highs and small heads are observed in the valleys (Figure 3). Comparatively low hydraulic heads are found in the E model domain, where groundwater pumping activities have led to a reduction in groundwater heads (Figure 5(a)). The distribution of the hydraulic gradient displays clear domains of maxima located preferentially in the center of the study area, between the Barnim Plateau and the Greater Spree Valley, consistent with the presence of major gradients in the surface elevation in these areas (Figure 3(a)). This model is used as the reference model for this study and stands representative for a model scenario, where only measured groundwater heads are available since we want to investigate whether such a setup is sufficient to reproduce local flow regimes or if additional data are needed.
The second model scenario aims at quantifying the impact of lakes, represented by the related hydraulic BC, on the resulting subsurface hydraulics and temperature distribution (M2 hereafter). In order to include this additional surface component, the absolute elevations of the water levels of all lakes (derived from the high-resolution DEM, Figure 3(a)) have been implemented as fixed hydraulic heads, which results in the introduction of local gradients in hydraulic potential across the interface of lakes and bounding aquifers (Figure 5(b)). We assume a direct connection between the aquifer and the lakes, depicting either influent or effluent conditions, depending on whether the surface water level is higher or lower than the surrounding groundwater level. The resulting heads range between 24.5 m and 77.2 m again; however, local changes of up to +12.2 m (-7.4 m) can be observed, which correlate spatially with the location of the implemented water bodies.
The last model of this study (M3 hereafter) integrates rivers as connecting elements between the lakes (Figure 4(c)). Therefore, assigned hydraulic heads (Dirichlet BC) display a gradient along the course of the river and more prominently to the surrounding aquifer (Figure 5(b)). Again, differences in head between the different model scenarios are substantial, standing representative for differences to overrepresented depressions deriving from anthropogenic overprinting. The resulting heads range between 24.5 m and 77.2 m again; however, local changes of up to +2.0 m (-8.3 m) can be observed, which correlate spatially with the location of the implemented water bodies.
With these distinct model scenarios, we are able to quantify the impact of surface water bodies (M2 = lakes and M3 = rivers) on the resulting hydraulic and thermal fields.
The pressure and temperature initial conditions for the transient simulations were derived by solving for pressure and temperature separately under steady-state conditions.
4.4.2. Thermal Boundary Conditions
The upper thermal BC was defined by imposing fixed values (1st kind Dirichlet BC) for the temperature across the uppermost four geological surfaces (Holocene, Saalian, Holstein, and Elsterian), the value of which has been derived from interpolated grids based on measured (groundwater) temperatures distributed across the model area (Figure 5(c)) [9, 45]. The top Holocene depicts fixed temperatures of 10.8°C derived from long-term average air temperature data [46]. The distribution of the temperatures of the following three geological surfaces shows that the hottest regions are located below the city centers of Berlin and Potsdam, correlating spatially with the amount of surface sealing (up to 15°C, [47]). In comparison, the lowest temperatures are found below least urbanized regions and in proximity to surface water bodies (7.8°C). This BC was chosen in order to represent the anthropogenic overprinting produced by surface sealing [47] and heat input from industrial activities on the present-day thermal field. By assuming stable temperature conditions at the surface, we were able to investigate the impact of the surface to groundwater dynamics on the thermal configuration of the subsurface. Here, depth and lateral variations of the thermal configuration as induced by predicted changes in fluid dynamics due to each modification of the reference model were systematically investigated.
The lower temperature BC, representing heterogeneously distributed temperatures at -6000 m.a.s.l., was derived from a lithosphere-scale 3D conductive model [9] in order to have a quantifiable heat input coming from the deeper parts (crustal and mantle) into the study domain. Figure 5(d) illustrates the distribution of the imposed temperature, showing a gradual increase from 196.6°C to 220.5°C while moving NW to SE.
5. Model Results
The model results will be presented in terms of pressure, fluid flow, and temperature distribution, predicted by each of the three different model scenarios and their respective comparison. They are shown either as absolute values or relative differences between the model realizations.
5.1. Pressure and Fluid Flow Field Distribution
5.1.1. Reference Model (M1)
The pressure distribution at the uppermost slice of the first model (topography) shows distinct areas of negative (unsaturated conditions) and positive (saturated conditions) pressure (Figure 6). Areas with negative pressure make up most of the model area and reach values of up to -8.11 bar. These areas are located beneath topographical highs (Figure 3(a)). Likewise, a general correlation between the topology of the topography and the pressure is evident (Figures 3 and 6(a)). Areas of high topography display highest negative pressures and vice versa. Accordingly, areas of positive pressure are found in the valleys and lowlands, as particularly visible at Lake Wannsee (Figures 3(a) and 6(a)). The highest positive pressure values are located in the very eastern portions of the model area (8.09 bar, Figure 6) and correlate with the topographical low of the limestone mine in Rüdersdorf. These high pressures also correlate with the location and extent of surface water bodies located in the model area (Figure 5(b)). However, the absolute values display high differences to the present-day situation. This can be seen in water depths at, e.g., Lake m, standing representative for a hydrostatic pressure of 1.6 bar, clearly indicating how the hydraulic potential of these water bodies is only approximated in a very crude manner by such a hydraulic boundary setting ( at bar).
(a)
(b)
(c)
(d)
(e)
(f)
The pressure distribution at -50 m.a.s.l. (chosen as the first continuous surface) shows a similar trend as the one described above. Areas of high pressure are located in the NE, S, and NW, correlating well with the topographic highs as well as areas of high hydraulic potential (Figures 3(a) and 5(a)). Lower pressures are predicted for the central model domain, correlating with areas of low topography, most prominently the Greater Spree Valley (Figure 3). These areas also display a low hydraulic potential (Figure 5(a)). This pressure distribution continues at greater depths, reaching the first regional impermeable layer in the model (Middle Muschelkalk, Figures 2(c) and 2(d)). Predicted fluid flow patterns (Figure 7(a)) follow the general hydraulic gradients as imposed by the upper hydraulic potential. Infiltration mainly occurs beneath areas of high hydraulic potential in the NE, NW, and S almost independently of the geological structure underneath. In comparison, exfiltration is strongly linked to areas of low hydraulic potential, mainly located in the central model domain (Greater Spree Valley) prominently visible at lake sites, the latter not yet resolved in terms of their hydraulic behavior in this model (Figure 7(a)). The example site of Lake Tegel in the east-southeast domain of the study area additionally displays two major domains: (a) influent conditions in the E and (b) effluent conditions in the W (Figure 7(b)).
(a)
(b)
(c)
5.1.2. Lake Model (M2)
The pressure distribution defined by the second model scenario shows some similarities both in trends and magnitude compared to the ones described above for the reference model (M1). Pressure changes between the two model realizations are mainly limited to areas of small lateral extent, directly below or nearby the extent of the surface water bodies implemented (Figure 6(c)). Maximum predicted differences are of the order of ±1.0 bar (Figure 6(e)). Areas of pressure increase mainly relate to lake sites where the hydraulic potential was formerly underestimated (losing lake, influent conditions), and pressure decrease commonly correlates with lake sites where the hydraulic potential was overestimated in M1 (gaining lake, effluent conditions). This hydraulic reversal can be best illustrated beneath the site of Lake Tegel (Figure 6(c)), where pressure differences of +0.5 bar are predicted. This leads to a higher hydraulic potential than that of the surrounding aquifer. In contrast, Lake Stienitz to the very E of the model area displays a drop in pressure of up to -0.6 bar, thereby resulting in a comparatively low hydraulic potential area. The depth extent of these changes reaches as far as the top of the Middle Muschelkalk, considered in all models presented as an aquitard of regional extent, with maximum values of -1920 m.a.s.l. though they are mainly restricted to the shallow model domain down to -300 m.a.s.l. (Figure 6(f)).
The fluid flow field resulting from this pressure distribution is similar to the one described for M1. Major areas of recharge are located in the NE, NW, and S, and major discharge areas are located in the central model domain, namely, the Greater Spree Valley, as well as localized areas of limited lateral extent in the SW and NW. The implementation of the additional hydraulic forcing due to the presence of major lakes leads to complete fluid flow reversal beneath areas that are hydraulically connected to the surface bodies. This is illustrated in the following for the specific Lake Tegel site in the east-southeast domain. While the western part remains a recharge area (Figure 7(c)), though also increasing the vigor of flow (), the eastern part displays a reversal in predicted flow directions. These changes in the hydraulic potential also lead to uprising of groundwater predicted in the direct vicinity of the lake. Here, pressure-driven flow caused by the generally low hydraulic potential of these areas is supported by the buoyant upcoming of heated water which infiltrated at the lake site (see also thermal results described below). This behavior is observed across the whole model area and results in modifications of groundwater circulation mainly located where the hydraulic potential of lakes is significantly higher than that of the surrounding aquifers.
5.1.3. River Model (M3)
The last model scenario predicts pressures which are in the same order of magnitudes as those predicted by M1 and M2 when averaged on a regional scale. Hence, major areas of recharge and discharge show similar distributions. Similarly to M2, pressure changes have a limited spatial extent, located directly below or in vicinity to the implemented river bodies (Figure 6(d)), reaching maximum values of ±0.80 bar (Figure 6(e)). A distinct area of pressure increase is predicted where the River Spree and Havel merge (to the S of Lake Tegel) with maximum values of +0.15 bar. This result is likely caused by the fact that the setup of the previous models (M1 and M2) resulted in effluent conditions (gaining stream) at this site, thus likely underestimating the hydraulic potential carried by these two river bodies (Figure 6(d)). In the very E of the model area, a local drop in pressure can be observed. Here, the hydraulic potential beneath the Mühlenfließ was overestimated in M1 and M2, consequently displaying a maximum drop of -0.4 bar leading to effluent conditions at this location. Pressure differences are observable down to depths of up to -1030 m.a.s.l. (Figure 6(f)) which is approximately half as deep as the changes induced by the lakes and is likely connected to the smaller general volume of water represented by these water bodies. This is also reflected in the majority of changes displaying a vertical extent of ≥ -100 m.a.s.l.
According to the results described above, the hydraulic field as predicted by M3 displays few modifications when compared to the previous model scenarios. These are mainly limited to domains in the direct vicinity of the major river bodies. This is shown clearly in the area where the Spree and Havel merge. At this site, M2 predicts mainly discharge due to low hydraulic potential (Figure 8(a)). In comparison, a complete reversal in flow directions is predicted by M3, changing from an area of discharge to an area where recharge is prevailing due to a local increase in the hydraulic potential (Figure 8(b)). Through this change in local hydrodynamics, discharge is also predicted to the W of the River Havel and S and N of the River Spree. Accordingly, all sites where the pressure dropped in this model scenario show an increase in the vigor of exfiltration (i.e., Mühlenfließ, ).
(a)
(b)
The results presented above show that with implementation of surface water bodies, changes in fluid dynamics lead to a configuration where the gradients between the latter and surrounding aquifers cause the process of bank filtration. Here, infiltration from lake and river sites towards depressions in the hydraulic head is in accordance with fluid dynamics observed in other studies [48–50]. This fit in circulatory patterns increases between the different model scenarios of this study, whereas M3 presents the best fit scenario of this study.
5.2. Temperature
The temperature distributions predicted by the different models are used as a passive marker in order to quantify both the lateral as well as the depth extent of the hydrodynamics caused by the changes in the surface hydraulic configuration. This approach is feasible since the upper model domain is dominated by pressure-driven fluid flow. Thus, any modifications in the pressure configuration and resulting hydrodynamics also lead to changes in advective heat transport and consequently predicted temperatures.
5.2.1. Reference Model (M1)
Shallow temperatures predicted by M1 show a complex distribution of minima and maxima, ranging from 20.8 to 34.2 (Figure 9(a)) at -500 m.a.s.l. Highest temperatures are predicted in the center of the model area and also locally in the E and NW. Negative thermal anomalies at this elevation are mostly confined to the N, E, and NW areas of the model domain. Their locations correlate spatially with areas of high hydraulic potential (Figures 5(a) and 6(b)) driving water to infiltrate from the surface into the deeper model domains (Figure 7(a)). Areas of low hydraulic potential (Figure 5(a)) are characterized by higher temperatures, mainly resulting from the upconing of deeper and warmer water to shallow levels (Figures 7(a) and 9(a)). The temperature distribution at this depth is additionally defined by the background conductive heat flow. Hence, the poorly thermally conductive Rupelian leads to high temperature gradients where maximum thicknesses of the unit are encountered (e.g., NW and SE, Figures 2(b) and 9(a)). The highly thermally conductive Zechstein also displays a positive correlation of its thickness maxima to predicted temperatures at this depth (E, W, and NE, Figures 2(f) and 9(a)).
(a)
(b)
(c)
(d)
(e)
(f)
The temperatures predicted by M1 at -3500 m.a.s.l. range from 108.8°C to 129.6°C (Figure 9(b)). At this elevation, the thermal field is no longer influenced by any active component due to groundwater flow, being only controlled by diffusive processes. The temperature distribution reflects the distribution of the differently conductive layers, with the Zechstein salt playing the most prominent role [2, 6, 11]. Indeed, the highest temperatures are found in domains of minima in the salt thickness (Figure 2(f)) where heat withdrawal and channeling is less effective compared to below salt structures (diapirs). Therefore, the lowest temperatures are predicted, where the highest thicknesses of the Permian unit is present. There is also a minor component of the lower thermal BC present, standing representative for the heat input from the deeper subsurface.
5.2.2. Lake Model (M2)
Despite a similar range in trend and magnitude of modeled temperatures, differences between predicted temperature distributions from M1 and M2 with depth are also found and lie within the range of ±5°C, wherein most are located between -2°C and +1°C with a mean value of -0.2°C (Figure 9(e)). The highest differences at the investigated elevation level (±4.0°C, Figure 9(c)) are located in the direct vicinity to the implemented surface water bodies and have a lateral extent as large as 4 km (e.g., Lake Tegel, = 0.1°C = smallest value where a spatial correlation to surface water bodies is still observable). This is a strong indicator for the changes in fluid dynamics induced by either the modeled increase or decrease in hydraulic loading due to the implementation of lakes as a surface hydraulic parameter. In detail, colder temperatures in M2 are found below surface water bodies, reflecting an increase in the hydraulic load above these areas (Figure 6(c)). This increase in infiltration also leads to higher predicted temperatures surrounding Lake Tegel (Figure 9(c)). This result derives from an increase in uprising of deep groundwater induced by changes of the upper hydraulic head BC in this region (Figures 5(a) and 5(b)). The temperature differences also indicate changes in fluid dynamics with depth, with significant modifications (±1°C) down to depths of more than -1500 m.a.s.l., and differences greater than ±0.1°C are displayed even at elevations up to -4800 m.a.s.l., spatially correlating with pressure changes, thus being likely linked to convective processes (Figure 9(f)). In comparison, pressure differences can only be followed down to -1600 m.a.s.l., indicating that induced changes in fluid dynamics have a more widespread effect than indicated by the pressure change alone.
5.2.3. River Model (M3)
Temperature differences between M3 and M1 amount to ±5°C, approximately at the same domains of maximum difference and lateral extent as described in the previous paragraph. However, comparing M3 and M2, local variations are also visible, though limited in magnitudes of up to ±1.3°C (Figure 9(d)) with most lying within a range of -1°C to +0.5°C (mean: -0.03°C, Figure 9(e)). The distribution of these differences is limited to areas directly beneath or adjacent to the implemented surface water bodies, which is beneath locations where the hydraulic gradients between rivers and adjacent aquifers reach their maximum values. Accordingly, local changes in modeled fluid dynamics result in lower temperatures where the hydraulic potential has been increased (i.e., Figure 9(d) and Figures 2(a) and 2(b)) and higher temperatures are calculated where it has been decreased (Figure 9(d)). These induced changes in fluid dynamics display a small depth range, which is indicated by a maximum elevation of temperature differences of ±0.1°C at -1200 m.a.s.l. (Figure 9(f)). Moreover, these elevations are far greater than the modifications in pressure indicate (-600 m.a.s.l.), thus showing that this tracer has a significantly higher level of sensitivity to capture even small changes more accurately. The horizontal extent of these modifications is also limited compared to M2, reaching its maximum value (~2200 m) at the outlet of Lake Stienitz (Figures 6(d) and 9(d)) and up to 1000 m at the intersection of Havel and Spree (Figure 9(d)).
6. Discussion
The results of the models presented in this study have demonstrated that the shallow hydraulic field below the city of Berlin is mainly controlled by the configuration of the imposed upper hydraulic BCs, specified by the groundwater level distribution and water levels of major surface water bodies. Resulting hydraulic potentials and gradients lead to a fluid flow pattern where infiltration into the deeper parts of the model domain is dominant in areas of high hydraulic heads or hydraulic gradients. Accordingly, exfiltration occurs at sites with low hydraulic potentials (Figures 5(a) and 7(a)).
The first model scenario clearly shows that predicted groundwater circulation is controlled by two factors: (1) the hydraulic head BC and (2) the structure and hydraulic conductivity of model units. To discuss (2) shortly, the results from a previous study [10] are compared with the results derived in this study concerning the feasibility of the parameterization of the Rupelian aquitard. In our previous study, we considered the Rupelian to be completely impermeable, thus only allowing exchange between the different aquifer compartments at local discontinuities (hydrogeological windows). In contrast, a different parameterization of the clay unit has been performed in this study, wherein the Rupelian aquitard still shows a low permeability, though of higher magnitudes than previously done (E-8 m/s, Table 1). To compare the two model realizations, fluid pathways and travel times equivalent to Pleistocene ages were computed, utilizing the base of the Pleistocene as the starting surface. This analysis was carried out, since previous fluid chemistry interpretations indicate mixing between infiltrating Pleistocene and recent waters of different salinity content at m (Tesmer et al. [29] and Möller et al. [30]). These studies indicated a leakage through the Rupelian to depths well beneath its base independent of the location of the windows throughout the model domain lateral extent. Thus, predicted fluid pathways with the time-framing described above should cover the entire lateral model extent. Comparing the different studies then shows that for travel times of Pleistocene origin, an impermeable Rupelian only produces a sparse distribution of lateral mixing of infiltrating water (~40% of model area). In comparison, even a comparatively lowly permeable Rupelian allows for faster infiltration and higher lateral mixing, as indicated by the large extent of infiltration below the Rupelian base (~90% of model area). Therefore, the hydrodynamics predicted by M1 of this study should be considered as more realistically representing real subsurface hydrology compared to the results from Frick [10]. Additionally, isotope ratios point towards an incomplete sealing of Permian and Triassic evaporites towards the overlying sediments [30], which questions the role of the Triassic Middle Muschelkalk as a tight aquitard [28], but was not studied here due to the paucity of data to properly parameterize this unit. It has also be shown by Bianchi and Pedretti [39] that solute transport specifically and fluid pathways in general are influenced strongly by the implementation of heterogeneities in the κ-field, which should be investigated further in future studies.
Even more intensively than the hydraulic conductivity, the results of this study demonstrate how the configuration of the upper hydraulic head BC (1) and any changes to the latter likely lead to strong modifications of the pressure field and the resulting distribution of fluid pathways as outlined in Section 5. In this context, the sensitivity of the model to the implementation of surface water bodies along with the integration of high-resolution topographical data has resulted in significant changes in the predicted hydraulic and thermal fields (Figures 6–9). This is best exemplified by increased pressures in areas below or adjacent to most of the implemented surface water bodies, leading to a complete reversal in the local fluid flow dynamics (Figures 7 and 8). Compared to M1, where the upper hydraulic BC was derived based on measured groundwater heads, differences in the elevation of hydraulic heads up to ±12.2 m ( bar) were implemented in the additional model realizations (M2 and M3). The resulting changes in hydraulic potential and pressure lead to modifications in the predicted local flow field, which is best illustrated at Lake Tegel. Here, a major low in the prescribed hydraulic head representing a major discharge area obtained in M1 features a local potential high in M2 (influent, Figure 5(b)). This in turn induces infiltration of surface water from the lake boundaries into the adjacent aquifers, as is typical for the region [51]. Moreover, the resulting pressure changes also lead to strong uprising in the direct vicinity of the lake shore. These results, visualized by calculating the fluid pathways, show a general downward flow beneath the lakes. When entering the first aquifer, the flow field develops a horizontal component flowing parallel to the surface of the first aquitard encountered at depth and it finally rises to the surface into a major discharge area (Figure 7(c)). The distribution of these flow paths mimics to some extent the observed local fluid dynamics, where intense groundwater pumping activities in these areas have led to considerably lower hydraulic heads (Figures 5(a) and 5(b)). This caused an increase of infiltration of surface water from the lake into the aquifer which is in accordance with the observed hydraulics, where up to 70% of groundwater wells in Berlin is fed by surface water infiltration [50].
From the results of this study, it appears clear that considering a simplified approach in deriving the hydraulic head distribution from measured groundwater heads, as is mostly done in these kinds of simulations [52–56], will only hinder the model to capture the complex dynamics occurring in the subsurface (Figures 7(b) and 7(c)). This can be seen clearly in the thermal signature, indicating changes in the groundwater dynamics as induced by the presence of major lakes reaching depths up to -4000 m.a.s.l. (Figure 9). This is far greater than for variations induced by seasonal fluctuation in surface temperature and/or human activities, illustrating the far-reaching long-term effects of importance for any further study focusing on the aforementioned parameters.
Considering major rivers in addition only leads to local modifications of the modeled hydraulic field. Indeed, even when rivers are hydraulically connected to the underlying aquifers, the overall magnitude of water exchange is comparatively small, reaching its maximum potential beneath areas of highest hydraulic gradients (Figure 5(b) and Figures 8(a) and 8(b)). Similarly to M2, these areas were predicted as discharge areas in M1 but are now featuring infiltration due to a local rise of hydraulic potential, which enables the process of bank filtration in groundwater production which is established in the vicinity of the surface water bodies (Figures 8(a) and 8(b), Massmann et al. [49]). Additionally, modifications of the hydraulic field are identified in regions where rivers connect to lakes. Here, the increase in the hydraulic head at the rivers has the effect of smoothing out local gradients in the hydraulic potential between the different water bodies (e.g., outlet Lake Stienitz: reduction from 8.2° to 0.2°, Figures 5(a) and 5(b)) and therefore results in less vigorous infiltration of surface water in nearby locations. The overall depth influence of these modifications (≥-1200 m.a.s.l., Figure 6(e)) is limited when compared to M2; however, it still encompasses the entirety of the shallow fresh water compartments.
The results obtained in this study also have a more direct and applied merit in the context of a safe fresh water management within the city of Berlin. As an example, M3 showcases the sensitivity of the shallow hydraulic field to major rivers. These results can then be used to analyze the influence of contaminants which are carried into the model domain or might even be coming from inside the model area, thus having a relevance to better quantify the danger of contamination of the fresh water aquifer. For instance, extensive lignite mining in the Lausitz (SE of model area) has led to the exposure of pyrite, which represents a major source of sulfur. This chemical component can be found in river water in the form of (e.g., [57]), which is a major groundwater pollutant. Increased concentrations of this component have already been found in groundwater samples within the model area (e.g., [58]). The source of the latter can be either by river load or through groundwater circulation. This problem and the increased salinity of the shallow groundwater compartments outline the necessity of being able to represent the hydrodynamics near the surface and the related feedbacks to and from deeper groundwater circulation with the highest amount of details. The study presented here therefore depicts a step towards a systematic understanding of river-aquifer interactions in an urban environment. However, to make quantifiable predictions, local studies should be envisaged, depicting even higher resolutions and taking into account contaminant and solute transport.
The thermal field predicted by the models of this study displays heterogeneously distributed temperatures at the investigated elevation levels. At depths greater than -3000 m.a.s.l., the thermal configuration is mainly controlled by conductive heat transport locally shaped by the structural configuration of the units with the highest or lowest values for thermal conductivities, namely, the Permian Zechstein and partly the Paleogene Rupelian. However, in the shallow model domain, this general relationship is overprinted by a component of convective heat transport. At these depth levels (≤-1000 m.a.s.l.), locally colder temperatures are predicted, where major areas of recharge result in infiltration of cold surface water. In contrast, locally warmer temperatures are predicted below major discharge areas, where heated water from the deeper model domain rises to the surface. These findings are in general agreement with a previous study [59] but show distinct differences in predicted temperatures, where focused down- and upward movement of groundwater (Figure 7(a)) is predicted by the models presented here. This is illustrated by predicted temperature differences at -500 m.a.s.l., which are as high as ±7.5°C. These are likely linked to the increase in convective heat transfer, induced by local increases or decreases of hydraulic gradients.
Using the predicted temperatures as a tracer for quantifying the effect of changes in fluid dynamics, the results show that the implementation of lakes produced maximum temperature differences of ±5°C, whereas colder temperatures are more common than warmer (Figures 9(c) and 9(e)). These modifications are traceable down to -4800 m.a.s.l. (mainly above -2000 m.a.s.l., Figure 9(f)). In comparison, thermal perturbations induced by the consideration of rivers are only up to ±1.3°C (Figures 9(d) and 9(e)). The related modifications of the thermal configuration are traceable until -1200 m.a.s.l., but mostly restricted to the very shallow model domain (≥-200 m.a.s.l., Figure 9(f)). Combining the results of the two modified model scenarios, colder temperatures prevail (Figure 9(e)). This likely derives from the general increase in hydraulic potential compared to the reference model. Consequently, most rivers and lakes represent losing conditions in M3, which is in accordance with observations [50] and results partly from reduced hydraulic heads due to groundwater pumping.
The effects of urbanization and climate change on the subsurface thermal field has gained increased focus in the last decade (e.g., [60–63]). The results of these studies show that the shallow subsurface in such areas has been highly modified, both thermally and hydraulically. These modifications are mostly represented by increased temperatures beneath urban centers connected to changes of surficial BCs, like surface sealing, industrial waste heat, and climate change amongst others, commonly referred to as urban heat islands [64]. This excess heat has been proposed for geothermal utilization, providing heating for decades [65]. However, these studies are purely analytical or limited to 1D, thus excluding 3D effects entirely. This study shows that the shallow geothermal field is controlled by numerous factors to which it reacts in a more or less sensitive manner. This in turn entails caution while performing any feasibility study for assessing the utilization of heat resources below major urban areas, where modifications in the local hydrodynamics as related to the aforementioned activities might have far-reaching and unexpected effects. This is supported by the results of this study, where local changes in hydrodynamics lead to temperature differences up to ±5°C, which translates to strong variations in the extractable amount of heat, especially for shallow geothermal utilization. At the same time, the results illustrate that these shallow thermal fingerprints might reach depths up to ≥ -4000 m.a.s.l., depending on the local hydrogeological conditions. All of these factors point towards a careful planning to set up possible city-wide energy provision from geothermal resources, which is especially relevant given that the thermal perturbations are well within the typical thermal breakthrough values during operational activities [66]. Therefore, an accurate description of the shallow hydraulics, represented here by more realistic hydraulic BCs, is of high importance. This is also evident in significant variations in the amount of extractable heat at a given place due to the changes in local hydrodynamics.
The models presented in this study focus on proposing and testing a conceptual physical model of the shallow to deep water dynamics beneath the area of the city of Berlin as based on a detailed geological model. This conceptual model underwent a quantitative and systematic calibration in terms of the gravitational response from the geological structuring, providing constraints on the quality of our geological configuration and range of rock material parameters [67]. Additionally, it has been validated against the two only available deep temperature logs in wells, providing further constraints on the geological reliability of the models as well as their basic parametrization [9]. Moreover, all model parameters entered have been derived from either direct measurements (field and laboratory) or literature research [9, 11]. On top of these data, the study integrates data on the shallow hydraulic conditions, most importantly the hydraulic BC, arriving at the most up-to-date understanding of the geological and hydrogeological configuration of the area beneath the capital city of Berlin. However, these types of studies might at best provide predictions about the system behavior within certain error bounds, which are indeed intrinsic in any study, whether based on extrapolation of measurements or on forward predictive physics-based models. In detail, the models presented here could still be improved and will be, given the appropriate data, in terms of incorporating more information of the geological configuration (well logs and seismics), parameterization of model units (calibration and measurements), or hydraulic BC setup (cross-border flow and denser database). However, we would like to emphasize that at this stage, under consideration of all available data, we are convinced by the robustness of the modeling results in representing the first-order dynamics of the system, which are in accordance with published studies as outlined above.
The results of this study show that lakes and rivers, under the assumption of ideal hydraulic connectivity, might reshape the fluid flow field considerably, especially where surface water bodies with large fluid volumes are involved. However, riverbeds (colmation layer) typically show a lower hydraulic conductivity than the surrounding aquifer, which would likely reduce the exchange of the groundwater and surface water. Studies investigating this aspect are so far mostly restricted to local models (e.g., [68]), and measurements of the physical properties of the colmation layer are rarely available. In the setting of the model area, river gradients and therefore also fluid velocities are rather small, which would suggest a comparatively thick and stable colmation layer as opposed to higher energy settings [69]. Nevertheless, hydraulic gradients between surface water bodies and surrounding aquifers are comparatively high (Figure 5(b)), which has been shown to be sufficient to cause leakage even through low permeability layers (i.e., Rupelian) independent of their thickness. To quantify these processes more accurately, ongoing studies are investigating the sensitivity of 3D thermohydraulic models to the implementation and parameterization of colmation layers. Additionally, the model is forced to present-day observed hydraulic head conditions, and these are prescribed for 250k yrs of simulation time (after Sippel et al. [11]). This implies (1) that we may overestimate advective cooling/warming effects looking at the evolution of temperatures (steady state after approximately 7k yrs) and (2) that we do not consider a “natural” state of the system as present-day hydraulic heads also integrate current pumping activities, which are however not part of the model dynamics. Here, future studies should implement groundwater wells, as effects of the latter on local hydrodynamics, as well as connected thermics are still poorly studied and understood. Especially the impact on predicted shallow temperatures is of importance since current policies opt for a greater share of renewables in the energy mix, also including geothermal.
7. Conclusions
The implementation and consideration of surface water bodies in 3D thermohydraulic modeling is of importance since it significantly influences the predicted hydraulic and temperature fields. First and foremost, changes to the local thermohydraulics are most pronounced where differences in prescribed hydraulic heads are highest between the surface water bodies and adjacent aquifers. Accordingly, we observe more vigorous infiltration of cold surface water into the deeper model domain where hydraulic heads of rivers and lakes are higher than those of neighboring aquifers, leading to forced convective cooling. In contrast, where hydraulic heads of surface water bodies are lower than the surrounding aquifers, a change of the fluid flow direction from downwards to upwards is predicted, leading to higher calculated temperatures since heated waters rise up at these locations. The modifications induced by the implementation of lakes are mostly confined to areas in their direct vicinity since no directional hydraulic gradient is imposed. In contrast, the implementation of rivers might lead to changes in subsurface hydraulics in areas where the absolute hydraulic head remains constant, but a downriver hydraulic gradient is existent.
This study was able to outline the depth influence of surface hydrodynamics on coupled deep and shallow fluid and heat transport. Implemented lakes account for modifications of the hydraulic and thermal fields down to an elevation of -4800 m.a.s.l., and the consideration of rivers modifies these fields down to an elevation of -1200 m.a.s.l. Herein, temperature differences of more than 5°C are predicted down to an elevation of more than -500 m.a.s.l., which translates to major modifications in depths of interest for shallow geothermal utilization. Additionally, these depth ranges encompass the entirety of shallow groundwater compartments, which underlines the importance of implementing surface water bodies, especially concerning the groundwater compartmentalization in connection with groundwater pumping or groundwater contamination from the river load in the model area. Understanding the influence of all physical processes involved serves as basis for future simulations of thermohaline energy and solute transport looking both at the present state and future scenarios.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research is funded as part of the joint initiative of the research field energy of the Helmholtz Association “Energy System 2050.” We would like to thank the Senate Department for Urban Development and the Environment of Berlin (SenStadtUm), Berlin Waterworks (BWB), Ministry of Rural Development, Environment and Agriculture of the Federal State of Brandenburg (MRDEA), and the Waterways and Shipping Office (WSA) for providing the data for this study. The 3D information on the structural configuration of the subsurface of Berlin was compiled, visualized, and interpolated with the commercial software package Petrel (Schlumberger©). All data are available upon request from the listed sources above, and simulation results are available from the first author upon request.
Supplementary Materials
The following files support the conclusions of the manuscript at hand and contain information as follows. Figures 10–18 describe the structural setup of the model in much more detail than was necessary for the readability of the manuscript. We provide an overview of the utilized well database in accordance with Frick et al. [9] (Figure 10). The resulting elevation (Figures 11–14) as well as thickness distributions (Figures 15–18) as derived from the structural modeling done by Frick et al. [9] is shown for all respective geological (model) units. These represent the input grids as produced with Petrel© with a resolution of 100 × 100 m. The only exception is the geological surface of the Holocene (topography) which was integrated with a resolution of 25 × 25 m as necessary for the scope of this study. We also included a map view of the employed upper temperature boundary conditions, which were derived from measured and interpolated groundwater temperatures as described in Frick [10] (Figure 19). Last, we also included map views for the presented model scenarios of recharge and discharge area distributions. These show an increasing alignment with proposed fluid dynamics [50] but were deemed as less illustrative than the stream tracer analysis presented in the main body. Lastly, we include a map of recharge/discharge magnitudes predicted by the models, wherein the most common values are in alignment with reported values [71]. (Supplementary Materials)