Boundary condition control on inter-aquifer flow in the subsurface of Berlin (Germany) – new insights from 3-D numerical modelling

We investigate the degree of hydraulic interconnection between the different (regional to local) groundwater compartments with respect to the choice of boundary conditions and their impact onto the groundwater safety beneath the urban centre of Berlin, capital city of Germany. To this end, we carry out a systematic study based on 3-D hydrothermal models differing in terms of imposed parametric set-ups of the hydrogeology, as well as different surface forcing with respect to their impact on fresh groundwater production.

The study area is part of the Northeast German Basin and consists of a thick sequence (up to 5 km) of differently consolidated sedimentary deposits. This sedimentary succession features a sequence of alternating aquifers (reservoirs) and aquitards which are connected to different degrees, each one depicting a specific composition of its mineralised pore water. The uppermost aquifer system (made up mainly of poorly consolidated siliciclastic rocks) acts as the main freshwater reservoir used for groundwater production by the municipal water supplier. This compartment is incompletely sealed from the brackish to saline aquifers extending at greater depths by a regional clay-enriched aquitard, the Oligocene Rupelian Clay. The latter shows a heterogeneous thickness distribution due to erosion during the latest glacial periods resulting in local discontinuities. This aspect opens to the potential risk of contamination of the drinking water reservoir from mixing with the saline groundwater upconing, locally enhanced by shallow pumping activities.

Based on our results and their correlation with available isotopic and chemical analysis of water samples, we demonstrate how hydraulic connection between the different compartments is indeed likely to occur thus supporting the possibility of a contaminant rise from the saline aquifers below through either natural or anthropogenic (pumping) forcing.

1 Introduction

The aim of this study is to investigate the degree of hydraulic interconnection between the different (regional to local) groundwater compartments and their impact onto the subsurface-groundwater safety beneath the urban centre of Berlin. The basic hydrogeology of the study area is conceptualised in Fig. 1. For this purpose we created two model scenarios featuring different sets of lateral hydraulic boundary conditions.

Figure 1 Hydrogeological setup of model area. (a) Geological cross section and hydrogeological compartmentalisation of the model area (location in upper left corner). Hatched units represent aquitards. Blue arrows represent conceptual flow paths. Vertical exaggeration = 5×. (b) Natural hydraulic head distribution: green dots = well locations. Coordinates [m] in Gauß-Krüger DHDN Zone 4. (c) schematic zoom of shallow-intermediate domain, modified from Limberg and Thierbach (1997). Upper blue line indicates groundwater head distribution modified by pumping well activity. Blue arrows represent conceptual flow paths. Vertical exaggeration = 15×.

For the first model realisation (Model 1, M1 hereafter) we use the classical approach of closing all lateral boundaries to fluid and heat flow, the shortcomings of which have already been discussed in Bayer et al. (1997). Here, all sources and sinks are assumed to be located inside the model area, which hinders a proper representation of the regional scale flow components. In order to quantify the role of the latter component in the overall hydraulic regime beneath the study area, the second model scenario (M2 hereafter) implements open lateral boundaries by applying a fixed pressure head at the lateral boundaries of the model (Fig. 1b, Haacke et al., 2018). In detail, the hydraulic potentials of the outermost bounds of the model have been used to mimic the regional groundwater flow component and to avoid dissimilarities between the different groundwater compartments. Both models also feature the effects of pumping activities on the predicted groundwater dynamics by explicitly implementing production wells with variable pumping rates as active boundary conditions. The initial hydraulic setup for both model scenarios derives from the “natural state”, that is the state before commissioning of the wells, which has been modelled in steady state for both scenarios using the data from Haacke et al. (2018) as upper hydraulic boundary condition.

We analyse the results for both, the regional as well as local characteristics, paying close attention to the interconnectedness between the different compartments at depth. In order to check the feasibility of the regional to local flow characteristics of the models presented in this paper, we compared the predicted groundwater head distributions (extent and depth of depression cones) with measured data, as well as with interpretations from isotopic data and groundwater chemistry. The latter two have been presented by Tesmer et al. (2007), Möller et al. (2008) and Lüders et al. (2010), showing that the different compartments display individual groundwater compositions related to different degrees of mixing. In detail, the authors showed that the shallow, post-Rupelian compartment displays a signature of Recent to Pleistocene precipitation. In contrast, the pre-Rupelian aquifer system represents a mixture of Pleistocene precipitation with older formational water which suggests an infiltration through the main hydrogeological discontinuities within the separating Rupelian aquitard. At the same time, isotope ratios from the aforementioned papers indicate an incomplete sealing of the Permian and Triassic evaporites, which questions the common assumption of the Middle Muschelkalk aquitard as a regional sealing cover (e.g. Pöppelreiter et al., 2005; Kaiser et al., 2011). Here, both, the uprising of dissolved evaporites from below its base, as well as infiltration of Pleistocene precipitation, is indicated in the isotopic data. Generally it can be concluded that the Post-Zechstein compartments were flushed out by meteoric waters displaying only fractions of the older brines formed during and after their deposition. The only exception to this rule is seen in the deep Pre-Zechstein reservoirs, where the original Paleozoic formation waters have been preserved, which points to a complete sealing of the latter (Lüders et al., 2010). With the models presented in this study, we aim to reproduce these characteristics via 3-D hydrothermal modelling. We opted for coupled hydraulic and thermal simulations because we can further use the resulting temperature distribution within the model as a passive tracer, therefore enabling a better visualisation and understanding of the regional and local hydrodynamics. Another focus of interest relates to the possibility of uprising saline water from the deeper compartments into the post-Rupelian successions. Since this compartment is used for the drinking water supply of Berlin (Möller and Burgschweiger, 2008) such uprising, both naturally and anthropogenically driven, is posing a major threat to groundwater safety. Evidence for a natural connection has already been presented in Möller et al. (2008). Therefore, connections between these compartments were also investigated both with respect to their natural occurrence as well as to the over-imposed anthropogenic forcing conditions as represented by shallow water pumping activities.

2 Hydrogeological Setting

The 3-D structural model used in this study has been previously described in Frick et al. (2019) and differentiates 21 geological units, 9 Cenozoic, 8 Mesozoic and 2 Paleozoic, for the sedimentary succession and 2 units for the underlying basement. The sedimentary succession can be subdivided into four major compartments, namely the shallow, intermediate, deep and very deep aquifer systems (Fig. 1a). These consist mainly of consolidated clastics separated to different degrees by variably thick aquitards. The shallowest, fresh water, compartment is separated incompletely from the brackish to saline aquifers below by the Rupelian aquitard (Fig. 1a, c). This aquitard shows erosional unconformities due to glacial erosion, thus opening fluid pathways between the two compartments (Fig. 1c).

The intermediate and deep compartments show higher salinity values of the fluids compared to the shallow compartment mainly deriving from the dissolution of Triassic to Permian evaporites. They are decoupled to a certain degree by the Middle Muschelkalk aquitard, but are not completely sealed off as has been shown by isotopic and chemical data (Tesmer et al., 2007; Möller et al., 2008; Lüders et al., 2010). The very deep compartment is separated from above by a thick layer of Permian Zechstein salt, which shows a highly heterogeneous thickness distribution (Fig. 1a, Frick et al., 2019) with some local discontinuities possibly leading to deep inter-aquifer flow. This compartment also shows a very high salinity of the fluids mainly related to seawater evaporation during depositional times (Lüders et al., 2010).

3 Modelling Method and Scenarios

3.1 Hydrogeological Model

For the problem at hand, we carried out 3-D coupled fluid and heat transport simulations based on the commercial software FEFLOW©. Our hydrogeological model (Haacke, 2018) resolves all relevant geological interfaces (21 geological layers) and the surface hydrology, meaning major lakes and rivers and was created with the commercial software Petrel (© Schlumberger).

The mesh resolution is on average on the order of 20×20 m, whereas the highest resolution is at localities of surface water bodies and groundwater production wells, where highest hydraulic gradients are expected. To ensure a stable vertical to horizontal element size ratio, the model was further subdivided into 57 computational layers (at least 2 per geological layer) summing up to ∼ 18 000 000 elements. The discretisation is therefore higher than the minimum that has been shown to be necessary in these kinds of simulations (Kaiser et al., 2013a).

This structural model was then imported into the commercial software FEFLOW© to solve for coupled fluid flow and heat transport. The mathematical background can be found in Diersch (2013). For the modelling task, we used individual thermal and hydraulic boundary conditions (BCs): (a) the upper thermal BC at the top slice has been derived from long term average annual temperature data (Senate Department for Urban Development and Housing, 2001), displaying values between 7 and 10 ^∘C where maxima are located in the centre of the model area (Frick et al., 2018), (b) the lateral thermal BC at all lateral boundaries of the model corresponds to the results of previously published models (Haacke et al., 2018), (c) the lower thermal BC at slice 58 as derived from Noack et al. (2010), showing values between 196 and 220 ^∘C gradually increasing from Northwest to Southeast. Due to the complexity of the hydraulic boundary conditions we present them in a separate paragraph (Sect. 3.2).

Table 1 Properties of the geological units as used for the calculations. Properties of Cretaceous and older units after Sippel et al. (2013). c^(s) = heat capacity of solid, ϕ = porosity, κ = hydraulic conductivity, Q_r = radiogenic heat production, λ^(b) = bulk thermal conductivity; Values were derived from Otto (2012), VDI (2010), Norden and Förster (2006), Norden et al. (2012), Das (2013), Devlin (2015) and Sippel et al. (2013).

Concerning the basic geometry, all relevant material properties are chosen in a similar fashion to previous studies. We hereafter focus on the parametrisation adopted for the Cenozoic units as listed in Table 1. The reader is referred to Sippel et al. (2013) and Frick et al. (2019) for details on the parametrisation adopted for the deeper model units. Hydraulic conductivity values have been determined following the approach of Devlin (2017), the latter being based on a geometric mean average weighted by the grain size distribution of the respective mineral assemblage. Additionally, an anisotropy ratio of 10:1 (${\mathit{\kappa }}_{xx,yy}:{\mathit{\kappa }}_{zz}$) was taken into account. For all parameters not described in the aforementioned predecessor studies or in this paper (e.g. Dispersivity), the default values of FEFLOW© have been adopted (Diersch, 2013).

3.2 Model Scenarios

To investigate the effect of closed versus open lateral hydraulic boundaries and the influence of pumping activities, two different model scenarios were simulated. Each of these displays a different setup of the hydraulic BCs. The upper hydraulic BC is a Dirichlet type fixed hydraulic head, except for domains which are influenced by groundwater abstraction and those coinciding spatially with locations of production wells (more details in Haacke et al., 2018), here implemented as point-wise sink/source terms. Hence, hydraulic heads can adjust freely in areas influenced by groundwater pumping. The initial condition for the hydraulic head distribution was chosen to represent the natural head distribution (not influenced by any pumping activities). It was obtained by a steady state simulation, without considering the additional (in)output of pumping wells (see Fig. 1b and Haacke et al., 2018). This hydraulic head distribution shows values between 75 m in the Northeast, 50 m in the South and Northwest and 25–30 m in the centre and West. The groundwater abstraction was approximated using available groundwater production data from 2005 to 2014 for the existing 42 well galleries and producing at a constant rate for 100 years.

In the first model realisation (M1), we impose no flow lateral hydraulic BCs. Through this model realisation, we neglect any influence of the regional groundwater circulation in the model, an approach that is often used in this type of simulations (e.g. Kaiser et al., 2013b). Such an assumption is based on the peculiar location of the study area in the middle of the north German plane, not flanked by any substantial mountainous region, therefore presumably characterised by small to negligible cross-boundary flow. Still, it follows that by neglecting the regional flow component, inter-aquifer flow is likely to be overestimated, and therefore the results should be considered as a conservative estimate of the present-day hydrodynamics. Due to this assumption, the influence of pumping activities is likely small, since a general vertical exchange is already favoured.

In order to quantify the effective role of the regional groundwater component on the hydrodynamics of the model area and therefore whether it is negligible or not, we carried out a second set of simulations (M2). Here constant pressure and associated gradients were prescribed along all lateral boundaries as Dirichlet type fixed hydraulic head. Differences between the results obtained from the two model realisations will therefore provide a quantitative measure of the influence of the regional flow on the inter-aquifer hydrodynamics.

Both model scenarios were first solved for coupled fluid and heat transport in transient mode without groundwater pumping, until reaching quasi-steady-state conditions after ∼275 kyr final simulation time. The temperature and pressure distributions derived this way were used as initial condition for the transient simulations of groundwater pumping. These were run for a simulation time of 100 years, which is the time period of large scale groundwater production activities.

To investigate the model results, we first analyse the influence of a regional flow component on the different compartments and their respective degree of hydraulic connection by examining the predicted flow dynamics derived from the final stage of the simulations that neglect pumping activities. In a second step, the impact of pumping activities is quantified by analysing induced changes in the overall dynamics of the groundwater flow, and by making use of the associated changes in temperature as a passive tracer. These results are also used to provide some conclusive remarks on the level of inter-aquifer connection at different depth levels. Therefore, we make use of available interpretations from isotopic and chemical data (for details see Möller et al., 2008; Tesmer et al., 2007; Lüders et al., 2010) to compare to our conclusions. Specifically we look at the role of different forcing, whether natural from regional and local hydrodynamics, or man made from pumping, on the contamination of the shallowest fresh water aquifers by uprising saline water.

4 Results and Discussion

Within the shallow model domain, the groundwater flow dynamics are primarily enforced by the imposed hydraulic surface BC (Fig. 2). Infiltration is predicted at the plateaus in the Northeast, South and Northwest, followed by horizontal flow towards the centre and east, where slow uprising occurs. Overall, the Darcy flux ranges between 10⁻⁵ and 10⁻⁸ m s⁻¹ whereas largest values are predicted at highest elevations in the Northeast. Comparing the results from the two model runs, both have similar flow characteristics. Changes in the local hydrodynamics are mainly limited to domains close to the lateral boundaries, where a gentle drop in flow velocities (magnitudes) is predicted in M2, though the extent and size of these variations can be considered minor. The fluid flow predicted by both simulations is characterised by vigorous infiltration and relative slower uprising and is in agreement with isotopic and chemical data, the latter pointing to mixing of groundwater within the deeper compartments by recent Pleistocene infiltration (Möller et al., 2008).

Figure 2 3-D view of the flow dynamics at different depth levels. (a, b) Shallow domain (above Rupelian), (c, d) intermediate domain (below Rupelian, above Middle Muschelkalk). Coloured flow lines represent preferential flow directions and magnitude of Darcy velocity. Blue arrows represent flow directions. Calculated with Paraview's stream tracer filter (Squillacote et al., 2007); Gray layers from top to bottom: Rupelian, Middle Muschelkalk, Zechstein, Basement. Coordinates [m] in Gauß-Krüger DHDN Zone 4.

Moving to the intermediate compartments, similarities in the hydrodynamics between the two model realisations and to the shallow compartment are also evident (Fig. 2). The Darcy flux at this depth ranges between 10⁻⁷ and 10⁻¹⁰ m s⁻¹ wherein highest fluxes are predicted at similar localities as in the shallow model domain. When compared to the flow regime described above for the shallow domain, it appears that groundwater flow at this depth levels presents a more regional component, being less affected by local systematics in the hydrodynamics. Both models predict extensive infiltration of shallow groundwater to deeper levels in regions where the Rupelian aquitard thins out or is missing, spatially correlating with high hydraulic gradients at the surface (Fig. 1b). In contrast, local uprising of deeper seated groundwater is also observed at discontinuities of the Rupelian located at low hydraulic gradients, as by analysing the forward flow tracers seeded within the intermediate compartment. Based on this information, we could conclude that mixing of shallow fresh waters by uprising saline to brackish waters is likely to occur at these localities. However, flow velocities driving this upward mixing are smaller in magnitudes relative to those responsible for the downward mixing. This last result suggests that, although possible, mixing of saline and fresh waters within the shallow aquifers is a transient phenomenon, the dynamics of which are sensitive to variations in the local hydrodynamics. The influence of the regional flow becomes evident at this depth level, manifested in the general change in flow directions near the boundary, especially in the Northeast and Southeast.

The deep to very deep model domain shows a completely different flow pattern as compared to what is described above (Fig. 3). In M1 we found that major inflow into the deep aquifer is predicted in the Northwest and East near the Zechstein diapirs and nearby major discontinuities in the Middle Muschelkalk aquitard. Along the southern boundary of the domain, focused flow toward a major discontinuity in the Zechstein layer leads to infiltration of groundwater into the very deep compartment. There, the fluids flow radially outwards from the discontinuity with a main northward flow direction. Fluid fluxes drop by up to two orders of magnitude compared to those observed in the aquifers above, ranging in between 10⁻⁹ and 10⁻¹² m s⁻¹. In M2 we found an opposite behaviour, as considering a regional flow component leads to a complete reversal in flow directions at these depth levels. The general flow direction predicted by this model shows a Northeast-Southwest trend for the deep aquifer. This leads to a hydraulic decoupling of this compartment from the very deep one, which means that below the Zechstein salt aquitard considering the regional flow component results in a complete reversal of flow direction in comparison to M1. The downward flow component predicted in M1 is completely hindered by the imposed variations in the overall hydrodynamics as shown in mainly Northeast-Southwest flow.

Figure 3 3-D view of the flow dynamics at different depth levels. (a, b) Deep domain (below Middle Muschelkalk, above Zechstein), (c, d) intermediate domain (below Zechstein, above Basement). Coloured flow lines represent preferential flow directions and magnitude of Darcy velocity. Blue arrows represent flow directions. Calculated with Paraview's stream tracer filter (Squillacote et al., 2007); Gray layers from top to bottom: Middle Muschelkalk, Zechstein, Basement. Coordinates [m] in Gauß-Krüger DHDN Zone 4.

Taken as a whole, the results from this analysis demonstrate the relevance of integrating cross boundary flow of regional characteristics even for hydrogeological settings characterised by relatively small topographic gradients, as in the present study. For the particular case of Berlin, we propose that including this regional flow component leads to a better agreement between model results and available observations and interpretations derived from fluid chemical analysis (Möller et al., 2008; Lüders et al., 2010; Tesmer et al., 2007). This is particularly the case when considering the deeper compartments, where the regional flow introduces a hydraulic decoupling between the respective aquifers as also testified by available isotope and chemical data (Lüders et al., 2010). With respect to the problem related to potential mixing of saline and fresh water within the Cenozoic aquifers, we can conclude that these aquifers show some level of interconnection, as demonstrated by preferential flow exchange by means of Pleistocene infiltration and mixing. However, this dynamic is highly sensitive to any local variation in the hydraulic forcing, an aspect that calls for caution when planning sustainable use of this resource.

In a first attempt to quantify the sensitivity of the inter-aquifer flow dynamics to local variations in the first order forcing conditions, in the following we discuss the results obtained when considering the effect of shallow water pumping in addition to the regional flow component for M2. Given the local character of the pumping activities, we focus on specific areas of interest, where the effects of these variations are at their highest. These include the Lower Havel, Lake Tegel and Spandau, all areas where the major waterworks are located (Fig. 4). A first remark on the results of these simulations is that the models show that there is only very little influence of the regional flow on the well-production induced modifications of the flow field as derived by comparing the final time steps of M1 and M2. Considering the pumping activities therefore leads to a drastic reversal of the fluid flow in both scenarios for the surroundings of the waterworks. In detail, natural discharge zones below the lakes refocus towards the well galleries where accelerated uprising can be observed (Fig. 4a, b). Moreover, a strong flow from the lakes towards the well galleries due to the drop in groundwater heads is observed. This agrees with the data from the local water suppliers (BWB), which approximate the amount of water they produce to be of lake origin (bank filtration) to ∼60 % (Möller and Burgschweiger, 2008).

Figure 4 Impact of well pumping in model 2 (M2). (a) Flow dynamics before pumping (b) flow dynamics after pumping. (a, b) Pathways were derived from the modelled nodal volumetric flux through the Streamtracer filter of Paraview (Squillacote et al., 2007) and plotted on the top of the Rupelian aquitard. The colour-coding represents the Darcy flux in vertically positive (uprising in orange) or negative (infiltrating in blue) direction. Termination points of lines correlate with either sites of recharge (light blue) or discharge (yellow). White outlines = extent of lakes and rivers. (a, b, c) Green spheres = well positions. Left middle panel indicates locations of (a), (b) and (c) with respect to the 3-D model. (c) Thermal tracer for M2: Temperature difference between after and before well pumping. Blue colours = cooling of subsurface, red colours = heating of subsurface. Grey units: upper = Rupelian, lower = Middle Muschelkalk. Integrative view towards North. Coordinates [m] in Gauß-Krüger DHDN Zone 4.

Additionally, extend and depth of the depression cones is in good agreement with the present day situation. However, the most interesting aspect is the change in the inter-aquifer flow, which is induced by groundwater pumping activities. This change can be depicted by the presence of a preferential vertical flow component above discontinuities in the Rupelian clay, but can be visualised more easily by looking at the thermal tracer (Fig. 4c). What we can see is a general cooling of the subsurface for large parts due to an increase in infiltration. This trend should lead to a decrease in the possibility of saltwater upconing. However, some areas show an increase in temperatures induced by groundwater production over 100 years. These variations occur in close proximity of the well bore area, where a focused upflow component is found. The spatial distribution of these areas comprise on the one hand locations where the Rupelian aquitard, here acting as a regional seal to mixing from deeper rooted water, is present and on the other hand, locations where this clay rich aquitard thins out or it is completely absent. Therefore, we could foresee that man made induced acceleration in the intra-aquifer flow might lead to some degree of contamination of the shallow fresh water domain. Looking at the backward fluid traces seeded at the well localities, the models show however, that the source areas are still coming from the freshwater compartments. As these models approximate the status quo under the current climatic conditions, this might change in the future due to decreases in groundwater recharge, surface water runoff and respective lake level drops likely leading to even stronger, far and deep reaching modifications of the fluid flow field. These scenarios are part of ongoing work which will also focus on geogenic mass transport and possible contamination.

5 Conclusions

The models presented in this study illustrate that:

• The modifications to groundwater dynamics caused by open versus closed boundaries increase vertically from marginally influencing to controlling.

• Open boundaries are needed to represent regional trends and might change predicted flow directions intensively leading to a weaker coupling between the groundwater compartments at depth, which was suggested by interpretations from isotopic and chemical measurements in previous studies.

• Groundwater pumping might increase inter-aquifer flow critically, opening the possibility of groundwater contamination which is of major interest for groundwater authorities and suppliers.

Future models should therefore focus on the role of the Rupelian clay in terms of its potential to separate the freshwater compartment from the saline. This can also be extrapolated to the proposed utilisation of the subsurface for geothermal installations, which might also disturb the system and should therefore be planned carefully. Additional components like subsurface pollution from surficial contaminants is also a newly arising problem which could be investigated in a holistic approach.