The interplay of Malm carbonate permeability, gravity-driven groundwater flow, and paleoclimate – implications for the geothermal field and potential in the Molasse Basin (southern Germany), a foreland-basin play

Abstract

The Molasse Basin in Southern Germany is part of the North Alpine Foreland Basin and hosts the largest accumulation of deep geothermal production fields in Central Europe. Despite the vast development of geothermal energy utilization projects especially in the Munich metropolitan region, the evolution of and control factors on the natural geothermal field, more specifically the time-varying recharge and discharge governing groundwater and heat flow, are still debated. Within the Upper Jurassic (Malm) carbonate aquifer as the main geothermal reservoir in the Molasse Basin, temperature anomalies such as the Wasserburg Trough anomaly to the east of Munich and their underlying fluid and heat transport processes are yet poorly understood. We delineate the two end members of thermal–hydraulic regimes in the Molasse Basin by calculating two contrasting permeability scenarios of the heterogeneously karstified Malm carbonate aquifer along a model section through the Wasserburg Trough anomaly by means of two-dimensional numerical thermal-hydraulic modelling. We test the sensitivity of the thermal-hydraulic regime with regard to paleoclimate by computing the two Malm permeability scenarios both with a constant surface temperature of 9 °C and with the impact of paleo-temperature changes during the last 130 ka including the Würm Glaciation. Accordingly, we consider the hydraulic and thermal effects of periglacial conditions like permafrost formation and the impact of the numerous glacial advances onto the Molasse Basin. Thermal-hydraulic modelling reveals the effect of recurrent glacial periods on the subsurface targets of geothermal interest, which is minor compared to the effect of permeability-related, continuous gravity-driven groundwater flow as a major heat transport mechanism.

Introduction

We use the COMSOL Multiphysics® software to perform two-dimensional numerical thermal-hydraulic modelling of the main fluid and heat transport processes in the Upper Jurassic (Malm) carbonate aquifer in the Southern German Molasse Basin (MB) for a better understanding of thermal anomalies induced by the interplay of high aquifer permeability (k), gravity-driven groundwater flow, and paleoclimatic conditions in foreland basins. Specifically, we numerically simulate for the first time the effect of the Würm Glaciation on the hydrogeology and the thermal field of the MB with a focus on the eastern, Bavarian part and particularly the region of Munich and the Wasserburg Trough. In particular, we examine a considerable cold temperature (T) anomaly within the Malm carbonate aquifer to the east of Munich (Agemar et al. 2012; Agemar and Tribbensee 2018), the so-called Wasserburg Trough anomaly, by means of a representative cross section (Fig. 1). The underlying main fluid and heat transport processes are yet poorly understood, which contributes to the geothermal exploration risk in the eastern Munich region, since temperatures in the Malm aquifer are lowered by up to 40 °C (Fig. 2) (Agemar et al. 2012; Agemar and Tribbensee 2018). We test different hypotheses to find the solution to this problem, since the location of the Wasserburg Trough anomaly coincides both with a zone of high Malm carbonate k due to karstification and with the melting zone of the Inn Glacier during the last ice age, i.e. the Würm Glaciation. On the one hand, we focus on two end-member scenarios of Malm k, thereby taking into account that the Malm carbonate aquifer presents a large range of permeabilities as shown in a map by Birner et al. (2012) and Fritzer et al. (2018). On the other hand, we specifically address the possible glacial influence on the present-day thermal and flow regime in the MB by comparison to model runs with a constant surface temperature (T_S) of 9 °C. So far, no numerical thermal-hydraulic model of the subsurface taking into account the effect of paleoclimate, permafrost and glaciers exists for the MB, but subglacial infiltration of meltwater during glaciations has repeatedly been considered as the main recharge mechanism of the Malm aquifer in the deeper MB (Bertleff et al. 1993; Bertleff and Watzel 2002; Heidinger et al. 2019). If this hypothesis is true, we expect that the Wasserburg Trough anomaly occurring in the MB and reaching as deep as the Malm aquifer (Agemar and Tribbensee 2018) is a direct consequence of subglacial recharge by the Inn Glacier during the last glaciations(s). This study focusses on the Würm Glaciation, whose traces and maximum extent are best preserved in the geological and geomorphic record and therefore best studied. We visualize the sensitivity of the model with respect to different combinations of the two Malm k scenarios and the two T_S scenarios notably by T differences between distinct model runs. Despite the existing numerical models, the geostatistical interpolation of the 3D thermal structure based on T measurements in boreholes (Agemar et al. 2012) implemented in the Geothermal Information System (GeotIS) (Agemar et al. 2014; LIAG 2021) to date represents the most reliable prognosis of subsurface T for the entire MB and in particular for the Top Malm (Agemar and Tribbensee 2018). Therefore, we compare and validate our model result of the thermal field with two 2D sections extracted from GeotIS (Fig. 2). In view of the existing dedicated studies of the hydraulic properties of the Malm aquifer by Birner et al. (2012) as well as of the Top Malm T (Agemar and Tribbensee 2018) extracted from the 3D thermal field of the MB (Agemar et al. 2012) to which we refer, a renewed borehole data analysis was not justified and beyond the scope of this work.

Fig. 1

Regional geological map of the Southern German Molasse Basin, modified after Bousquet et al. (2012) and Doben et al. (1996). Maximum extent of the Würm Glaciation according to digital maps by Gibbard et al. (2011). Dark red line AA’ indicates the studied cross section. Inset map shows the location of the map in Germany (MB = Molasse Basin)

Fig. 2

Temperature at top Malm (Upper Jurassic) derived from the geostatistical 3D temperature model of the LIAG (Agemar et al. 2012; Agemar and Tribbensee 2018) and projected on the simplified geological map (Fig. 1). Maximum extent of the Würm Glaciation according to digital maps by Gibbard et al. (2011). Red line AA’ indicates the studied cross section. Dark blue lines (the eastern one being partially covered by red line AA’) indicate the locations of the GeotIS sections presented in Fig. 18

The MB also serves as a reference play for the comparison and understanding of the interplay of high k aquifers, gravity-driven groundwater flow and paleoclimate in other orogenic foreland basins worldwide. As known for instance from the Alberta Basin in Western Canada, which is another foreland basin covered by two large ice sheets (i.e. the Cordilleran and Laurentide ice sheets) during Quaternary glacial periods (Dyke 2004; Seguinot et al. 2014), paleoclimate and the glacial cover can exert a significant influence on subsurface temperatures up to 2 km depth (Demezhko et al. 2018; Gosnold et al. 2011; Gray et al. 2012; Jessop 1971; Majorowicz 2012; Majorowicz et al. 2012).

Geological setting

The MB in Southern Germany is part of the North Alpine Foreland Basin (NAFB) and hosts the largest accumulation of deep geothermal production fields in Central Europe (Goldscheider et al. 2010). Henceforth, we generally refer to the Southern German MB simply as MB. It is up to 130 km wide perpendicular to strike and 400 km long parallel to the Alps (Fig. 1). Three main geological units are distinguished, each separated by a pronounced hiatus or unconformity (Bachmann and Mueller 1992; Bachmann et al. 1987; Bachmann and Muller 1991; GeoMol Team 2015). The oldest unit corresponds to the Variscan crystalline basement including continental clastic sediments in graben structures of Permo-Carboniferous age (Freudenberger and Schwerd 1996; Lemcke 1988). The basement is covered by Triassic to Cretaceous shallow marine sediments that formed as passive-margin shelf deposits of the European plate (Freudenberger and Schwerd 1996; Meyer and Schmidt-Kaler 1990, 1996; Unger and Meyer 1996). The youngest unit corresponds to Late Eocene to Late Miocene shallow marine and terrestrial sediments that accumulated in the NAFB and represent the Molasse sediments (Bachmann and Muller 1991; Kuhlemann and Kempf 2002; Lemcke 1988). Henceforth, we generally refer to the Molasse sediments simply as Molasse. As the only one of seven Alpine foreland basins in Europe (Allen et al. 1986; Moeck et al. 2015b), the NAFB was partially covered by piedmont lobes of glaciers during the Quaternary glacial periods, in particular the last (Würm) glaciation, which lasted from about 115 to 11.7 ka (Seguinot et al. 2018; Van Husen 1987). Figures 1 and 2 show the extent of Würm glaciers in the MB at the Last Glacial Maximum (LGM) (Doppler et al. 2011; Gibbard et al. 2011; Van Husen 1987). Beyond the terrestrial regions covered by glaciers at the LGM, periglacial conditions, notably permafrost, prevailed in vast areas of Europe during each glaciation (Lindgren et al. 2016; Renssen and Vandenberghe 2003; Vandenberghe et al. 2014).

Geothermal use and thermal-hydraulic field

The Malm carbonates currently constitute the main geothermal reservoir. They crop out at the Swabian and Franconian Albs or Jura to the north of the Danube River. To the south, the carbonate platform is dipping beneath the southwards thickening wedge of Tertiary Molasse (Fig. 1). Lying at a depth of about 2000 m below Munich, the top of the carbonates reach a depth of more than 4000 m at the Alpine front (Bachmann and Mueller 1992; Brink et al. 1992; Lemcke 1977; Reinecker et al. 2010). Carbonate aquifers frequently present a high secondary porosity and k due to karstification. The regional pattern of Malm aquifer hydraulic conductivity or k has been studied in detail by Birner (2013) and Birner et al. (2012) for its geothermal potential. Hydraulic conductivities in the Malm are closely linked to karstified zones in relation to the reef carbonate facies, faults and the proximity of the Danube (Birner et al. 2012). The extraordinary high k in the deeper parts of the MB are likely related to the southeasterly emersion and subsequent extensive karstification of the upper 150 to 200 m of the carbonate platform during the Cretaceous (Birner 2013; Birner et al. 2012; Frisch and Huber 2000). By reason of temperatures mostly in the range of 70 to 100 °C in the Munich region and up to 150 °C close to the Alps (Agemar et al. 2012; Agemar and Tribbensee 2018; Dussel et al. 2016; Flechtner et al. 2020), many boreholes tap the Malm carbonate aquifer mainly for district heating and balneological use, but some also for electricity generation (Flechtner et al. 2020; GeoMol Team 2015). In general, zones with hydraulic conductivities higher than 10⁻⁶ m s⁻¹ are suitable for geothermal energy production (Birner et al. 2012). In this context, despite high flow rates, the prominent cold T anomaly of the Wasserburg Trough presents a considerable geothermal exploration risk with respect to the minimum reservoir temperature needed for a specific geothermal use (Stober and Bucher 2012). The Wasserburg Trough anomaly thus clearly restricts the geothermal potential for direct heat applications. Currently, minimum temperatures of about 70 °C are used for direct heat production in the western Munich region (Flechtner et al. 2020). The use of lower reservoir temperatures (< 60 °C) often requires additional heat pumps and thus causes higher costs. Recently, 3D numerical modelling of coupled fluid and heat flow processes in the MB showed that basin-wide fluid flow strongly affects the shallow thermal field (Przybycin 2015; Przybycin et al. 2017). The results from hydrochemical and isotopic studies (e.g. Bertleff et al. 1993; Bertleff and Watzel 2002; Heidinger et al. 2019; Mraz et al. 2019) constitute important evidence for the discussion and validation of our model results.

Recharge of the Malm aquifer occurs on the one hand along the Danube River, and especially to the north of it through the outcropping Malm (e.g. Heine and Einsiedl 2021), and on the other hand by leakage through the overlying Quaternary and Tertiary sediments (Andres and Frisch 1981; Bertleff and Watzel 2002; Birner et al. 2011, 2012; Frisch and Huber 2000; Heidinger et al. 2019; Villinger 1977). The general flow direction in the MB is from the west to the east towards a groundwater depression in the Munich region from where the groundwater flows north and discharges to the Danube between Neustadt and Regensburg (Fig. 2) (Birner et al. 2011). As a hydraulic barrier, the NW–SE-oriented basement structure called Landshut-Neuötting High (Fig. 1), with a throw of up to 1300 m, separates the Wasserburg Trough in the west from the Lower Bavaria Trough in the east (Bachmann et al. 1987; Frisch and Huber 2000). According to Toth (2009), cross-formational gravity-flow of groundwater acts on the geological timescale of 130,000 years considered in the present study. That is why it is obviously very important to quantify its thermal-hydraulic effect in the orogenic foreland-basin geothermal play type, since the driving force of regional groundwater flow is topography, notably the elevation difference between the higher regions close to the mountain belt and the larger streams at lower elevation. According to Moeck (2014), the occurrence of sloping aquifers typical of foreland basins may significantly perturb the geothermal field, especially if faults or high k layers allow advective heat transport from the deeper to the shallower realm of a foreland basin. Permafrost largely hampers groundwater recharge in cold climate, in particular during glaciations (Bertleff et al. 1993; Bertleff and Watzel 2002; Jost et al. 2007; Lebret et al. 1996, 1994; Sloan and van Everdingen 1988; Van Vliet-Lanoë 1989; Woo and Winter 1993). Bertleff et al. (1993) and Bertleff and Watzel (2002) discussed the effect of glaciations using hydraulic, hydrochemical and isotopic data. According to them, the higher hydraulic head produced by foreland glaciers leads to an increased recharge potential in their melting zone, further assuming that larger perennial rivers, lakes and glaciers hamper permafrost formation or preservation in subglacial basins beneath them (Sloan and van Everdingen 1988).

Methods

Modelling approach and model scenarios

We use COMSOL Multiphysics® versions 5.4 to 5.6 (COMSOL 2020a) to simulate coupled heat and fluid flow in the MB along a representative 2D section through the Wasserburg Trough thermal anomaly (Fig. 2) to find an explanation for its occurrence. More specifically, the main objective of the present study centres on the simulation of the thermal-hydraulic processes with regard to the Malm carbonate aquifer as the main geothermal reservoir in the MB. For any details about thermal-hydraulic modelling, we refer to the extensive online documentation of COMSOL Multiphysics® version 5.6 (COMSOL 2020a). The description of the model set-up in the following sub-chapters respects as far as possible the workflow in COMSOL Multiphysics®. The model is set up for multiphysics coupling of the subsurface flow module using Darcy’s Law (COMSOL 2020c) and the heat transfer module (COMSOL 2020b). The fluid flow and thermal processes are fully coupled. The study is time-dependent and simulates the last 130,000 years to capture the effects of surface paleo-T changes during the last glaciation and of gravity-driven groundwater flow on the subsurface thermal field. Our model section represents the geometrical situation to the east of the Munich where a prominent cold anomaly occurs (Fig. 2) (Agemar et al. 2012; Agemar and Tribbensee 2018). The section is oriented NNW–SSE and extends 180 km from the northern margin of the Franconian Alb via Ingolstadt to the Inn valley near Rosenheim close to the Alpine Front (Fig. 1). Model depth amounts to 10 km below mean sea level. We want to emphasize that the focus of our study lies on the significance of Malm k and glacial conditions for the formation of thermal anomalies. Hence, we do not consider 3D groundwater flow which could possibly reproduce more exhaustively the entire flow system of the MB (e.g. Frisch and Huber 2000). In addition, our 2D section is restricted to the MB since only a limited recharge is possible from the Alps into the MB (e.g. Lemcke 1976). Orogenic belts and their adjacent foreland basins usually are hydraulically disconnected from each other by the frontal fault of the foothills (Toth 2009).

The two end members of thermal-hydraulic regimes in the MB are delineated by calculating two contrasting k scenarios of the heterogeneously karstified Malm carbonate aquifer. We test the interplay of both thermal-hydraulic regimes with paleoclimate by computing the Malm k scenarios both with a constant T_S of 9 °C and with the impact of paleoclimatic changes during the last 130 ka including the Würm Glaciation. We consider the hydraulic and thermal effects of periglacial conditions like permafrost formation and the impact of the numerous glacial advances onto the MB. We provide an overview of the model scenarios and the individual model runs in Fig. 3. In addition, a valuable tool to visualize the effect of changes between scenarios or the model sensitivity consists in combining the respective solutions of different model runs to show the difference between the solutions over time or for a specific time. According to the respective columns two and four shown in Fig. 3, model runs A, D and F are based on the low Malm k scenario, whereas model runs B, C, E and G are based on the high Malm k scenario. In the individual model runs, we combined both Malm k scenarios with the two T_S boundary conditions (i.e. paleoclimate and constant T_S of 9 °C, respectively) and modifications thereof to simulate the effect of the Inn Glacier (i.e. thermal and/or hydraulic conditions considered individually), specified in the first column in Fig. 3, to get meaningful and productive results. We show selected results in the specified figures and provide more details about the intentions behind the whole series of model runs in the results chapter.

Fig. 3

Overview sketch on the different Malm permeability scenarios (see also Table 1; scenarios in first row refer to model runs in the corresponding second and fourth columns) and surface temperature scenarios (scenarios in first column refer to model runs in the corresponding rows) including modifications thereof (i.e. thermal and/or hydraulic conditions considered individually). Individual model runs as well as selected results shown in the specified figures are indicated

Table 1 Properties of model units

Model geometry and mesh

For a faster computation of the basic physical, thermal-hydraulic processes over geological times, we built a simplified 2D model section oriented perpendicular to the strike direction of the MB, and therefore containing only the essential geological and structural elements (Fig. 4). With regard to the scale and the geological time scale of the model as well as the necessary simplifications of the stratigraphic and structural complexity, the three-unit model presented hereafter actually reflects the main contrasts in thermal and hydraulic properties. It equally correlates with the main geologic controls on elements of the geothermal-play concept (Moeck et al. 2015a). From bottom to top, the elements consist of the heat source (i.e. the basement), the reservoir (i.e. the Malm aquifer) and the seal or cover (i.e. the Cretaceous sediments and the Tertiary Molasse). We provide further details and explanations about model simplification in the Appendix. In our final model, the Danube and Altmühl Rivers lie at 300 and 305 m altitude, respectively, in the northern part of the model section. The river altitude implemented in our model is lower than the actual altitude of 362 m of the Danube near Ingolstadt, but closer to the 325 m altitude of the discharge zone near Regensburg, because we retained the geometry of the initial generic model. The Inn River lies at 440 m altitude near Rosenheim in the south of the model section. Along the chosen section, due to the occurrence of the Inn valley near Rosenheim, the usual gentle rise of the terrain towards the Alps reaches its maximum altitude of 600 m at 30 km to the north of the Alpine Front.

Fig. 4

Generic model geometry representing the simplified Molasse Basin geology along cross section AA’ shown in Fig. 1. It spans from the Franconian Alb in the north (left) over Ingolstadt to the Inn valley (near Rosenheim close to the Alpine Front) in the south (right). Vertical exaggeration equals five

The first step of model construction concerns the compilation of relevant literature (see above) with thickness maps of the relevant geological units. In step two, we visualized the relevant data and georeferenced the relevant thickness maps using the ArcMap® software. We derived the topographical information from the European Digital Elevation Model (EU-DEM) version 1.1 (European Environment Agency 2016). In step three, based on the trace of the model section, we selected the supporting interpolation points of the model geometry, entered the data in an Excel® spreadsheet and created a first figure in Matlab® (Fig. 4). The model geometry is based on linear interpolation between supporting points. Step four consisted in building the 2D geological model section in COMSOL Multiphysics® by means of the integrated model builder (COMSOL 2020a).

The physics-controlled mesh was set to extremely fine, which produced a mesh consisting of about 8700 elements with an average element area of about 0.2 km². The maximum element size is 1800 m for mesh elements at the base of the model. The minimum element size is 3.6 m in the northern part of the model due to the Altmühl and Danube valleys, a more irregular geometry, and the Molasse unit pinching out. The mesh element area within the Malm aquifer is about 0.17 km², whereas it is smaller than 0.1 km² in the northern shallow part of the model. For verification, we also computed model run G (i.e. a high Malm k scenario) with a maximum element size of 400 m resulting in about 31,300 mesh elements. Only the resolution is increased, but the numerical differences in the results are negligible. Model set-ups with a finer mesh do not converge. We performed all model runs on the identical mesh in order to avoid numerical artefacts when calculating e.g. T differences.

Composition and properties of model units

The lithological composition of the model units, conforming to the properties presented hereafter is summarized based on the available literature where the geological evolution is described in detail (e.g. Bachmann and Mueller 1992; Bachmann et al. 1987; Freudenberger and Schwerd 1996; GeoMol Team 2015; Kuhlemann and Kempf 2002; Lemcke 1988). The Variscan basement is mostly composed of gneiss and intrusions of plutonic rock dominated by granite (GeoMol Team 2015) and is comparable to the rock types exposed in the Black Forest or the Bohemian Massif. The Permo-Carboniferous sediments present a very unequal distribution (Bachmann and Mueller 1992; Bachmann et al. 1987) and are therefore not distinguished from the basement. The thin Triassic sediments with a thickness between 0 and 100 m (Bachmann et al. 1987) are included in the Variscan basement, corresponding to the upper European continental crust. The thin Middle and Lower Jurassic with a thickness between 0 and 150 m (Bachmann et al. 1987) are equally included in the Variscan basement. The contrasts in thermal and hydraulic properties also justify these simplifications. The Malm carbonates consist of various types of limestones, marls and dolomites that developed in a tropical shallow shelf sea including reefs, lagoons and shallow basins (e.g. Meyer and Schmidt-Kaler 1990, 1996). We derive the thickness of the Malm unit ranging between 500 and 600 m from Bachmann et al. (1987) and Meyer and Schmidt-Kaler (1996). Due to the general high degree of karstification of the unconfined Malm along the Danube and the Altmühl rivers, as well as in the Swabian and Franconian Albs to the north of the MB (Birner et al. 2012), we distinguish a highly permeable Malm sub-unit, irrespective of the Malm k scenario applied (Table 1). The Cretaceous period is mainly represented by some hundreds of metres of Middle and Late Cretaceous coastal and shallow marine sediments (Freudenberger and Schwerd 1996) predominantly made up of clay marls (Unger and Meyer 1996), and is here included in the overlying Molasse wedge. The MB itself is composed of a wedge of Tertiary Molasse sediments whose base reaches up to 4000 m below mean sea level at the Alpine Front in the south of the MB. We derive the thickness of the Molasse (depth of Tertiary or Molasse base) from Lemcke (1977), Lemcke (1988), Roeder and Bachmann (1996), Brink et al. (1992), Bachmann et al. (1982) and Bachmann and Mueller (1992). The Molasse sediments are predominantly composed of cyclically deposited clayey and sandy series (Bachmann and Muller 1991; Lemcke 1988), hydraulically corresponding to an equivalent silt composition (Andres and Frisch 1981). According to Schulz et al. (2012), a more detailed subdivision of the Tertiary and Quaternary basin fill is not productive for modelling results with a focus on the Malm aquifer.

Table 1 summarizes the properties used in our model. We refer to the compilation of rock properties and references by Dussel et al. (2016), which they conducted for their thermal-hydraulic 3D model of the Munich region and a multidisciplinary reservoir characterization with a focus on the Malm carbonate aquifer. Birner et al. (2012) studied the hydraulic properties of the Malm aquifer by means of pumping test data and show that the highest hydraulic conductivities of up to 10⁻⁴ m s⁻¹ are found at the northern border of the MB along the Danube. Hydraulic conductivities of 10⁻⁵ m s⁻¹ are characteristic of karstified Malm carbonate to the south of the Danube. The largest contrast in hydraulic conductivities exists in the southern, deeper realm of the MB. In Baden-Württemberg to the west of the Iller River, hydraulic conductivity decreases rapidly to values lower than 10⁻⁸ m s⁻¹ towards the clay-rich Helvetic facies in the south, whereas in Bavaria hydraulic conductivity values of 10⁻⁵ m s⁻¹ are found in the Wasserburg Trough from Ingolstadt down to the Munich metropolitan region. We focus on two Malm k scenarios, thereby taking into account that the Malm carbonate aquifer presents a large range of permeabilities depending on the carbonate facies, the degree of karstification (Birner et al. 2012; Fritzer et al. 2018) and the vertical hydrostratigraphic profile (Bohnsack et al. 2020; Dussel et al. 2016). Since we decided to capture the end members of thermal-hydraulic regimes in the MB and to implement a generalized k structure, we derive the horizontal k variation from the map of hydraulic conductivity values provided by Birner et al. (2012). Because the thermal anomaly imaged by Agemar and Tribbensee (2018) coincides with the prominent zone of high hydraulic conductivity values in the map presented by Birner et al. (2012), our high Malm k scenario specifically reflects the hydraulic conditions in the karstified Malm carbonate in the eastern (Bavarian) part of the MB. We define our low Malm k scenario in a way that it reflects Malm carbonate devoid of karst like in the western part of the MB. We convert the original values of hydraulic conductivity derived from Birner et al. (2012) to k values according to the following equation:

$$\mathrm{k}=\mathrm{K}\cdot \frac{\upmu }{\uprho \cdot \mathrm{g}},$$

(1)

where k is the permeability in m², K is the hydraulic conductivity in m s⁻¹, ρ is the temperature- and pressure-dependent water density in kg m⁻³, µ is the temperature- and pressure-dependent dynamic viscosity of water in Pa s and g is the gravity constant in m s⁻². Except for the Variscan basement, we applied an anisotropic k in agreement with Andres and Frisch (1981) who noted that generally horizontal k is one order of magnitude higher. In the low Malm k scenario, the k of the confined Malm aquifer varies linearly from low anisotropic values of 10⁻¹³ m² horizontally and 10⁻¹⁴ m² vertically at the Danube valley to low anisotropic values of 10⁻¹⁵ m² horizontally and 10⁻¹⁶ m² vertically in the southernmost part of the model. In the high Malm k scenario, the k of the confined Malm aquifer varies linearly from high anisotropic values of 10⁻¹⁰ m² horizontally and 10⁻¹¹ m² vertically at the Danube valley to low anisotropic values of 10⁻¹⁴ m² horizontally and 10⁻¹⁵ m² vertically in the southernmost part of the model. By contrast, the unconfined Malm sub-unit has a constant k of 10⁻¹⁰ m² (Table 1). Rock thermal conductivity is primarily T-dependent and therefore significantly decreases with increasing depth due to the geothermal gradient. For sedimentary rocks, we applied a T correction of thermal conductivity according to Somerton (1992). For igneous and metamorphic rocks, we applied a T correction of thermal conductivity based on Seipold (2001).

Model boundary conditions

Heat transfer boundary conditions

We define the initial model T by applying a thermal gradient of 30 °C km^–1, in agreement with temperatures of 150 °C at more than 4 km depth close to the Alpine front (Dussel et al. 2016). At the vertical model boundaries, we define no-flow conditions for heat transfer. This is a valid assumption since the conductive thermal gradient is largely subvertical between high temperatures at great depth and low temperatures at the surface. Clauser et al. (2005) applied a basal heat flow density (Q) of 80 mW m^–2 for an optimal adaptation of their model. After consideration of radiogenic heat production, we apply a conductive basal Q varying linearly between 60 mW m^–2 in the south and 70 mW m^–2 in the north at the base of our model. In view of the geodynamic setting of the MB and our objective of independently analysing T anomalies, we perceive the published Q values (e.g. Hänel and Hurter 2002) as more reliable input data than the calculated T-depth maps recently published by Przybycin et al. (2015), which are based on own assumptions of the latter authors. At the upper boundary, we apply a paleoclimate forcing (Fig. 5). The paleo-T data are among the most important model input data. We derive climate or paleo-T forcing using time-dependent T offsets from a paleoclimate proxy. For the present study, we chose the EPICA ice core record (European Project for Ice Coring in Antarctica) (Jouzel and Masson-Delmotte 2007; Jouzel et al. 2007) which, in a glaciological study by Seguinot et al. (2018), yielded spatially and temporally consistent ice dynamics. T reconstructions from the EPICA record are provided and can be downloaded e.g. from the PANGAEA data repository (Diepenbroek et al. 1998, 2002). The time resolution of the data decreases with age from 10 to 45 years at 130 ka. For simulation, we use time steps of 50 years. We further amend paleo-T forcing with time-dependent paleo-precipitation reductions and T lapse-rate corrections. Seguinot et al. (2018) reduced precipitation with air T in order to simulate the potential rarefaction of atmospheric moisture in colder climates. The derived relationship results in about 60% less precipitation at the LGM. Accordingly, the linear scaling (f = 1.33) used by Seguinot et al. (2018) was also implemented in this study to calibrate the paleo-T forcing in the light of reduced precipitation. According to a present-day mean annual T_S map (e.g. Hänel and Staroste 1988), and altitude information from the European Digital Elevation Model (EU-DEM) version 1.1 (European Environment Agency 2016), temperatures of 9 °C at 300 m and of 5 °C at 1000 m were determined. The result for Munich at 520 m altitude is 7.7 °C and consequently the T lapse rate in this part of the NAFB corresponds to 0.0057 K m⁻¹. In the Danube and Altmühl valleys, but also in the Franconian Alb, we avoid freezing and thus permafrost formation for numerical stability reasons by disabling the T lapse rate and by restricting the minimum T to 4 °C. In map view, the top Malm negative thermal anomaly (Agemar and Tribbensee 2018) coincides with the confluence of the Isar and Amper valleys. Since permafrost is frequently absent beneath lakes and larger perennial rivers (Boulton et al. 1993; Sloan and van Everdingen 1988), we here apply particular conditions for river valleys, where paleo-T is not allowed to fall below 4 °C. The corresponding average and median paleo-temperatures over the last 130 ka are 3.2 °C and 2.3 °C, respectively. The paleo-T model input data were smoothed by applying a running average of 21 values and adding an additional running average of 11 values on top of the first (Fig. 5). As already mentioned, in order to test the sensitivity of the two Malm k scenarios with regard to different T_S scenarios, we also performed model runs with a constant T_S of 9 °C at 300 m altitude (Fig. 3).

Fig. 5

adapted from the European Project for Ice Coring in Antarctica (EPICA) (Jouzel et al. 2007). Linear scaling for paleo-precipitation forcing is implemented following Seguinot et al. (2018). Red and blue colouring represents mean annual surface temperature above and below zero degrees C, respectively

Paleo-T for a mean annual surface temperature of 9 °C at 300 m altitude

Subsurface fluid flow boundary conditions

We define the initial fluid pressure in the model as a hydraulic head replicating the topography. At the vertical boundaries delimiting the model in the north and the south as well as at the bottom boundary, we define no-flow conditions for groundwater. The Variscan basement is characterized by a very low k and thus is assumed to inhibit groundwater flow from the north and from beneath the studied depth range. Following Toth (2009) as well as previous results from the MB (e.g. Frisch and Huber 2000), only a limited recharge from mountains into foreland basins is possible. Therefore, we assume hydraulic no-flow conditions at the southern boundary. The configuration of the water table largely controls the driving forces of regional gravity-driven groundwater flow systems. In the Altmühl and Danube valleys where groundwater flows towards these major rivers, we define an atmosphere or gauge boundary condition. In the northern part of the model corresponding to the Franconian Alb, we specify a hydraulic head in equilibrium with the Altmühl River. In the southern part, at the surface of the Molasse, we apply a hydraulic head replicating the topography (Andres and Frisch 1981; Clauser et al. 2005). This assumption can be made for low-relief terrain as shown by Freeze and Witherspoon (1967), Forster and Smith (1988), Lopez and Smith (1995), and Toth (2009), which also implies that the Molasse sediments are entirely water-saturated and that the availability of meteoric water is always equal to or higher than or gravity-driven groundwater recharge.

This study aims at studying the effect of periglacial conditions (i.e. permafrost) and numerous glacial advances onto the MB. In order to test the effect of recurrent glacier cover in the southern part of the model (Fig. 4), we include model runs in which the effect of the Inn Glacier is simulated (Fig. 3) by applying an additional hydraulic head fluctuating as a function of paleo-T (Fig. 6). According to Greenwood et al. (2016) and Lemieux et al. (2008), the conduit water level in the foreland part of glaciers (i.e. the piedmont lobes) can reach the glacier surface leading to high meltwater supply for recharge beneath the melting zone of glaciers, a condition which had already been argued for by Bertleff et al. (1993). Since the timing and conditions of actual permafrost occurrence beneath foreland glaciers is unknown, we assume that maximum recharge was possible, and subsequently discuss the significance of our results. We calculate the evolution of Inn Glacier height in the Alpine foreland using the EPICA paleo-T data, and subsequently validate our results with the modelling data presented by Seguinot et al. (2018). Since the evolution of glacier extent modelled by Seguinot et al. (2018) cannot directly be used in our model, we set up an own iterative calibration process for simulating the Inn Glacier in the foreland. Although the approach has no robust scientific basis, the following procedure is the one that reasonably well reproduces the model results for the evolution of the Inn Glacier obtained by Seguinot et al. (2018). The maximum Inn Glacier height of about 500 m in Rosenheim in the south of the model section is in agreement with the maximum extent of the Würm Glaciation at the LGM (Van Husen 1987). As a precondition for foreland glacier formation, paleo-T must fall below a cut-off T of − 9 °C (with respect to the original proxy data). We also implement a linear proportionality of the maximum glacier height with the coldest paleo-T of − 14 °C (with respect to the original proxy data), once the previously mentioned condition is satisfied. We add another condition on top to prevent foreland glacier formation and to force foreland glacier melting if the running average of glacier height is below, respectively, falls below a defined threshold of 200 m over the preceding time interval of about 3000 years.

Fig. 6

Evolution of the Inn Glacier height in the Alpine foreland during the Würm Glaciation (see also Seguinot et al. 2018). Glacier height is applied as time-dependent hydraulic head at the southern model boundary (Inn Glacier in Fig. 4)

Fluid properties

Since glacial and periglacial climate conditions are taken into account, the formation of permafrost by freezing pore water is implemented by means of the phase change material sub-node of the heat transfer module (COMSOL 2020b). Phase change considers specifically the thermal energy needed or freed at the moment of melting, respectively, freezing of water. The transition T is 273.15 K (i.e. 0 °C), the transition interval is 1 K and the latent heat of fusion is 333.5 kJ kg⁻¹. Table 2 presents the contrast in water properties above and below the freezing point.

Table 2 Properties of water and ice above and below the freezing point

To satisfy energy and mass conservation in phase change models, particular attention should be paid to the density in time simulations (COMSOL 2020b). Because volume change and model deformation is expected, we define the transport velocity field and the density to ensure local mass conservation. Therefore, we convert the ice density and ice thermal conductivity by multiplication with the ratio of water and ice densities at the freezing point. Beyond the phase change of water, the simulation of coupled heat and fluid flow requires the following fluid properties: dynamic viscosity, density, heat capacity at constant pressure and thermal conductivity. A particularly important change is the drastic increase in water dynamic viscosity at the freezing point (Fowler 1997). In the model, ice only exists over a relatively small T interval, so that we consider the thermal properties of water below 0 °C to be independent of T and pressure. It is assumed that liquid water is pure, which is a valid assumption since the Malm aquifer has drinking water quality (Frisch and Huber 2000; Stober et al. 2014). Furthermore, we assume that hydrostatic fluid pressure prevails throughout the model domain. This assumption is valid for the Malm aquifer (Andres and Frisch 1981; Lemcke and Tunn 1956) and most of the younger stratigraphic units. We are aware that fluid pressure may tend to lithostatic with depth beneath the Malm aquifer and that the content in dissolved solids certainly increases in the Variscan basement, but no data for this depth realm are available. Water density would remain almost constant with depth and dynamic viscosity increase due to dissolved solids.

Since our model extends to a depth of 10 km, we extend the properties of liquid water beyond T dependency and also take into account the effect of hydrostatic pressure (P) by consulting the international steam tables (Wagner and Kretzschmar 2008). We specifically calculate the individual property values along the chosen T–P gradient. Based on previous model runs, we choose realistic T and hydrostatic pressure gradients of 30 °C km⁻¹ and 95 bar km⁻¹, respectively. The range of T–P conditions covered is 0.01–307 °C and 1–948 bar. Along the selected T–P gradient, all relevant water properties show significant changes, ranging from + 14% for thermal conductivity, + 6% for isobaric heat capacity and − 92.5% for dynamic viscosity (see also Smith and Chapman 1983). The decrease in water density (− 19%) is important because shallower colder water has a tendency to sink, whereas deeper warmer water tends to ascend. The T-driven free convection and the gravity-driven forced convection are united in the expression of the hydraulic head. The latter is a standard solution in COMSOL Multiphysics®. The hydraulic head corresponds to the sum of two components, which are the elevation head z on the one hand and the pressure head φ on the other hand (Toth 2009). The hydraulic head is a direct measure of the fluid’s potential energy in a given point of the flow region. However, since water density is T-dependent in our models, the pressure head and thus the hydraulic head increases with depth. Therefore, in order to be able to compare our results with standardized hydraulic head data published, e.g. by Stober and Villinger (1997) and Frisch and Huber (2000), we calculate the freshwater head according to the following equation:

$${h}_{\mathrm{fw}}=\varphi \frac{\rho }{{\rho }_{\mathrm{fw}}}+z,$$

(2)

where h_fw is the freshwater head in m, ρ_fw is the density of freshwater in kg m⁻³, ρ is the temperature- and pressure-dependent water density in kg m⁻³, φ is the pressure head in m, and z is the elevation head in m. For freshwater, we use a T of 10 °C and a density of 1000.0847 kg m⁻³.

Heat transfer in sedimentary basins is typically controlled by the two main mechanisms of conduction and convection by groundwater. In case the groundwater flow is driven by external forces, i.e. forced convection or advection, then the relative importance of the two heat transfer mechanisms is indicated by the dimensionless geothermal Péclet number (Bachu 1988, 1999) defined as:

$$\mathrm{Pe}=\frac{{\rho }_{w}{{C}_{p}}_{w}q D}{{\lambda }_{m}} \cdot \frac{{\nabla }_{h}T}{{\nabla }_{v}T},$$

(3)

where λ_m is the thermal conductivity of the saturated porous medium in W m⁻¹ K⁻¹, ρ_w is the water density in kg m⁻³, C_pw is the specific heat of water in J kg⁻¹ K⁻¹, q is the Darcy velocity in m s⁻¹, D is the thickness of the respective aquifer in m and ${\nabla }_{h}T$ and ${\nabla }_{v}T$ are the horizontal and vertical components of the thermal gradient in K m⁻¹. For Péclet numbers < 0.1, conduction is the dominant heat transport mechanism (Clauser et al. 2005). With regard to the long geological time scale considered, we consider the groundwater flow as cross-formational and therefore the aquifer thickness used to derive Péclet numbers corresponds to the vertical extent of the entire model.

Results

Model runs A and B: hydraulic head and groundwater flow

Figure 7 shows the hydraulic head, here the freshwater hydraulic head, which develops at different times in the low Malm k scenario as a result of model run A (Fig. 7a, b) and in the high Malm k scenario as a result of model run B (Fig. 7c, d). In particular, Fig. 7a, c visualizes the freshwater hydraulic head which develops in presence of the Inn Glacier. However, at first, we study only the hydraulic effect of the Inn Glacier, i.e. without applying paleoclimatic conditions. The point in time of 26.3 ka (i.e. model time of 103.7 ka) is shown because it corresponds to the maximum height of the Inn Glacier in Fig. 6. At the upper model boundary corresponding to the surface of the Molasse, the freshwater hydraulic head adequately replicates the topography, including the fluctuation height of the Inn Glacier. When Fig. 7a and c, respectively, Fig. 7b and d are compared, the freshwater hydraulic head decrease directly beneath the surface is stronger in the high Malm k scenario. The low Malm k induces a groundwater flow that is essentially limited to the Molasse. Only a minor fraction of groundwater flows towards the Malm aquifer (Fig. 7a, b). By contrast, the high Malm k induces a groundwater flow that is no longer restricted to the Molasse but descends into the Malm aquifer (Fig. 7c, d). A major difference between the LGM, here represented by the hydraulic effect of the Inn Glacier, and the present situation in both low and high Malm k scenarios consists in the contrasting groundwater flow patterns. In our models, this particular situation exists because a region of higher altitude is located between the Danube in the north and the Inn in the south of the section. Presently, two cells of forced convection exist: a southern one draining towards the Inn River, and a northern one draining towards the Danube River (Fig. 7b, d). In the southernmost part of the model, where the surface falls off towards the Inn valley near Rosenheim, groundwater flow descends from the highest freshwater head at 600 m altitude and ascends near the Inn valley (Fig. 7b, d). This groundwater thus forms a separate convection cell, which is decoupled from the flow heading towards the Malm aquifer and to the north. During cold periods, however, each foreland advance of the Inn Glacier onto the MB induces a reorganization of groundwater flow leading to a single cell of forced convection towards the Danube in the north (Fig. 7a, c). This reorganization of groundwater flow as a hydraulic effect of the Inn Glacier is practically immediate (at most a few hundred years) even for the short foreland advance of the Inn Glacier 40 ka ago (Fig. 6).

Fig. 7

Freshwater hydraulic head and streamlines in the low Malm permeability scenario with the sole hydraulic effect of the Inn Glacier (i.e. model run A introduced in Fig. 3) a at the Last Glacial Maximum 26.3 ka ago (time = 1.037E5 a), and b at present (time = 1.3E5 a). By comparison, freshwater hydraulic head and streamlines in the high Malm permeability scenario with the sole hydraulic effect of the Inn Glacier (i.e. model run B) c at the Last Glacial Maximum 26.3 ka ago (time = 1.037E5 a), and d at present (time = 1.3E5 a). Streamlines indicate Darcy’s velocity field resulting from model properties and a constant surface temperature. White lines are isopotential levels at 50 m intervals. Vertical exaggeration equals five

Darcy’s velocity magnitude

Darcy’s velocity magnitude is a measure of how fast groundwater is flowing. We show Darcy’s velocity magnitude resulting from model run B as a primary characteristic of the k structure and groundwater flow in the high Malm k scenario (Fig. 8a). In order to compare Darcy’s velocity magnitude in both k scenarios and to show the related model sensitivity, Fig. 8b shows the difference of Darcy’s velocity magnitude between model runs B and A (Fig. 3). In agreement with the very low k, groundwater flow is extremely slow in the Variscan basement. A contrast of up to three orders of magnitude exists between Darcy’s velocity magnitude in the Malm aquifer and in the Molasse. In both the low and high Malm k scenarios, Darcy’s velocity magnitude within the Malm aquifer is lower than 0.1 m a⁻¹ in the deepest part of the MB, whereas it progressively increases to about 1 m a⁻¹ in the centre of the MB (section km 65) and further increases to about 10 m a⁻¹ near the Danube. Darcy’s velocity magnitude within the Molasse ranges between about 0.01 m a⁻¹ and 0.1 m a⁻¹ (Fig. 8a). In the Franconian Alb to the north of the Danube where the karstified Malm carbonate is exposed, our models show that Darcy’s velocity magnitude exceeds 10 m a⁻¹ with a maximum of more than 1000 m a⁻¹. As shown in Fig. 8b, the higher Malm k causes a significantly higher Darcy’s velocity magnitude within the Molasse and in the direction of the Danube, especially in the northern two-thirds of the MB. The higher Malm k also significantly affects the groundwater regime within the MB by causing an increase and a decrease of Darcy’s velocity magnitude in the deeper and shallower parts, respectively (Fig. 8b).

Fig. 8

a Darcy’s velocity magnitude at present (time = 1.3E5 a) for the high Malm permeability scenario with a constant surface temperature (i.e. model run B), and b difference of Darcy’s velocity magnitude between the high and low Malm permeability scenarios (i.e. difference between model runs B and A). Darcy’s velocity magnitude is expressed as logarithm base 10 and unit m a⁻¹. Vertical exaggeration equals three

Péclet number

Figure 9 shows the Péclet number distribution for the low and high Malm k scenarios, respectively, with a constant T_S condition, as the respective results of model runs A and B (Fig. 3). In both model scenarios, heat conduction characterizes the Variscan basement, which is in agreement with the cataloguing of the MB as conduction-dominated play-type foreland basin. In the low Malm k scenario, the Malm aquifer presents isolated zones dominated by heat conduction as well as zones dominated by heat convection, which shows that heat convection dominates only locally and thus makes large-scale heat transfer much less effective. By contrast, the Malm aquifer in the high Malm k scenario presents a contiguous pattern of high Péclet numbers, which demonstrates that advective heat transfer is significant. Only in the deepest part close to the Alpine Front, where the Malm k equals 10⁻¹⁴ m², does the Péclet number indicate that heat conduction dominates in both Malm k scenarios. The Péclet-number range within the Molasse is also significantly different in the low and high Malm k scenarios (Fig. 9). In the low Malm k scenario, isolated zones of relatively low Péclet numbers consistently characterize the Molasse (Fig. 9a). By contrast, in the high Malm k scenario, the Molasse exhibits higher Péclet numbers especially in its northern part and above the highly permeable Malm aquifer (Fig. 9b). The low Malm k scenario is representative of the overall non- or poorly karstified Malm carbonate, and therefore represents the conduction-dominated benchmark of the MB. We consequently refer to T differences between high and low permeability scenarios as thermal anomaly of the high k, convection-dominated scenario in comparison to the low k, conduction-dominated benchmark.

Fig. 9

Péclet number distribution at present (time = 1.3E5 a) showing the relative importance of heat advection versus conduction with a constant surface temperature in a the low Malm permeability scenario (i.e. model run A introduced in Fig. 3), and b the high Malm permeability scenario (i.e. model run B). The geothermal Péclet number is dimensionless. Vertical exaggeration equals three

Model runs C and E, and difference E–C: hydraulic and thermal effect of the Inn Glacier

In order to test the model sensitivity with respect to the occurrence of the Inn Glacier, we calculate the hydraulic head difference between model runs E and C (Fig. 3). Model run E simulates the high Malm k scenario including the glacial and periglacial features like glaciers and permafrost under paleoclimatic conditions, whereas model run C simulates the high Malm k scenario only with permafrost but without the hydraulic and thermal effects of the Inn Glacier (Fig. 3). Therefore, the hydraulic head difference between model runs E and C reveals the hydraulic effect of the Inn Glacier on subsurface flow in the MB. The results clearly show that the hydraulic head in the southernmost part of the model periodically increases when the Inn Glacier advances into the foreland (Fig. 10a). Similarly, in order to test the model sensitivity with respect to the possible thermal impact of the Inn Glacier, we also calculate the T difference between model runs E and C. A supposedly higher quantity of cold water infiltrating into the basin barely influences the thermal field, except for a small negative thermal anomaly at the southern model boundary (Fig. 10b). In our models, the southern convection cell induces an upflow of groundwater and heat mostly within the Molasse near the Alpine Front and towards the Inn valley near Rosenheim in the south of the section (Fig. 7b, d). In fact, the hydraulic effect of the Inn Glacier produces the small negative thermal anomaly by temporarily precluding this upflow. Flow reversal and cold-water infiltration causes temperatures at the same depth (0–4 km) to decrease by up to 40 °C.

Fig. 10

Hydraulic and thermal effect of the Inn Glacier (i.e. difference between model runs E and C introduced in Fig. 3). a Hydraulic head difference and b temperature difference induced by the Inn Glacier at its maximum extent 26.3 ka ago (time = 1.037E5 a) in the high Malm permeability scenario under paleoclimatic conditions. Vertical exaggeration equals three

Model runs D and E: thermal field as a result of low and high Malm permeabilities

Figure 11a and b shows the thermal field in the MB in the low Malm k scenario with permafrost development under paleoclimatic conditions as a result of model run D (Fig. 3). By contrast, Fig. 11c and d shows the high Malm k scenario with permafrost development under paleoclimatic conditions as a result of model run E. Here, we show the point in time of 20 ka (i.e. model time of 110 ka), because it corresponds to both a considerable Inn Glacier height and a period of long-lasting cold paleo-T towards the end of the Würm Glaciation (Fig. 5). In the low Malm k scenario, groundwater flow limited to the Molasse causes cooling mainly in the Molasse and does not affect the Malm aquifer (Fig. 11a, b). In the high Malm k scenario, groundwater flow descending from the surface into the Malm aquifer causes subsurface cooling extending to significantly greater depth (Fig. 11c, d). In the latter case, the general flow and T pattern results from the downflow or leakage of cold meteoric water through the Molasse of relatively low k, slow heating of this water and finally upflow of the heated groundwater through the karstified Malm aquifer to the surface in the Danube valley. A comparison of Fig. 11a and c also reveals a different maximum permafrost thickness depending on the Malm k scenario. High Malm k causes much more pronounced cooling of the shallow subsurface (Fig. 11c). When the paleo-T starts to decrease, cold-water downflow through the Molasse forms a dark blue curtain moving down towards the Malm aquifer. About 35 ka ago, at the start of the coldest period of the last glaciation (Fig. 5), permafrost develops to greater depth where cold water has previously flown down and thus where the subsurface has been undercooled (Fig. 11c). Permafrost thickness reaches on average up to 600 m in the high Malm k scenario (Fig. 11c). As a comparison, permafrost thickness is only 150–200 m in the low Malm k scenario (Fig. 11a). At the end of the LGM, permafrost has largely melted 15 ka ago, but patches of permafrost remain until about 10 ka ago. Again, after the permafrost has melted, cold-water downflow forms a dark blue curtain. The effect of permafrost formation and thawing is regionalization of cold-water downflow over the entire width of the MB. Indeed, after 10 ka of simulation time, the 25, 50 and 75 °C-isotherms show a maximum downward curvature in the northern part of the MB between section km 40 and 70 where groundwater mass flow, or recharge, is highest due to the very high Malm k and the relatively low thickness of the Molasse cover in the range of 1500 to 2500 m (Fig. 11c, d). After 30 ka of simulation time, the maximum downward curvature of these isotherms has shifted southwards to a zone between section km 55 and 75. After 50 ka of simulation time, the maximum downward curvature of the 20 °C-isotherm is located at section km 70. A second downward curvature begins to develop underneath the highest hydraulic head at profile km 100. After 70 ka of simulation time and the first phase of permafrost formation, the 25 °C-isotherm has become almost horizontal between profile km 60 and 100, i.e. parallel to the base of permafrost. The latter zone extends even further between section km 50 and 110 after 110 ka of simulation time (i.e. at the LGM 20 ka ago; Fig. 11c) with the 25 °C-, 50 °C- and 75 °C-isotherms reaching depths of about 1 km, 1.8 km and 2.4 km, respectively (Fig. 11c, d).

Fig. 11

Thermal field including permafrost development in the low Malm permeability scenario under paleoclimatic conditions (i.e. model run D introduced in Fig. 3) a 20 ka ago (time = 1.1E5 a) and b at present (time = 1.3E5 a), and in the high Malm permeability scenario under paleoclimatic conditions (i.e. model run E) c 20 ka ago (time = 1.1E5 a) and d at present (time = 1.3E5 a). Red lines are isotherms at 25 °C intervals. White line marks base of permafrost (0 °C). Vertical exaggeration equals three

Heat flow density and thermal gradient

As a result of model runs D and E, Fig. 12 shows the vertical component of the conductive Q in the low and high Malm k scenarios, respectively. Figure 13 shows the corresponding thermal gradient in both k scenarios. Overall, a significant perturbation of heat flow and of the thermal field exists in the MB in the high Malm k scenario. In the Franconian Alb to the north of the Danube, the conductive Q increases from 70 mW m⁻² at the base of the model to about 100–105 mW m⁻² at the top of the Variscan basement. The Q rises upwards due to radiogenic heat production within the Variscan basement, which is characteristic of a conductive regime. In the MB part of the model, the conductive Q distribution is in both low and high Malm k scenarios influenced more or less by heat convection. In the low Malm k scenario, the conductive Q ranges relatively uniformly between 65 and 75 mW m⁻² in the central part of the section, not much higher than the Q range of 60 mW m^–2 in the south and 70 mW m^–2 in the north at the base of the model. The conductive Q is extremely low beneath the highest point (about 25 mW m⁻²) and increases to 65 mW m⁻² at 2000 m depth (Fig. 12a). In the high Malm k scenario, the conductive Q ranges between 70 and 90 mW m⁻² within the upper basement in the central part of the section. The Malm stands out with a conductive Q reaching up to 120 mW m⁻². Within the upper part of the Molasse, convective heat redistribution is so intense that the conductive Q is close to zero (5–10 mW m⁻²) at the surface and rapidly increases to 60–65 mW m⁻² at 1500–2000 m depth (Fig. 12b).

Fig. 12

Vertical component of the conductive heat flow density at present (time = 1.3E5 a) under paleoclimatic conditions in a the low Malm permeability scenario (i.e. model run D introduced in Fig. 3), and b the high Malm permeability scenario (i.e. model run E). Vertical exaggeration equals three

Fig. 13

Geothermal gradient at present (time = 1.3E5 a) under paleoclimatic conditions in a the low Malm permeability scenario (i.e. model run D introduced in Fig. 3), and b the high Malm permeability scenario (i.e. model run E). Vertical exaggeration equals three

The coupled fluid and heat flow significantly modifies the initial, homogeneous and conductive geothermal gradient. Three regions with distinct geothermal gradients exist (Fig. 13). In the Franconian Alb to the north of the Danube valley, groundwater mainly cools the shallow subsurface by flowing laterally towards the Altmühl valley as well as towards the Danube valley. In both Malm k scenarios, the relatively undisturbed thermal gradients in the basement mostly range between 0.025 and 0.035 K m⁻¹. The geothermal gradient within the Molasse is relatively undisturbed in the low Malm k scenario, whereas it is significantly reduced in the high Malm k scenario (Fig. 13). Thermal gradients range between about 0.015 K m⁻¹ in the shallow subsurface and 0.045 K m⁻¹ near the Malm aquifer in the low Malm k scenario, and from less than 0.002 K m⁻¹ in the shallow subsurface to more than 0.060 K m⁻¹ just above the Malm aquifer in the high Malm k scenario (Fig. 13). Both the pattern of vertical conductive Q and thermal gradients indicate heat redistribution from the shallow subsurface towards the Malm aquifer at depth, especially in the case of high Malm k. In the region of the Danube valley, the subsurface heats up because of warm water discharge by the Malm aquifer (Fig. 13b).

Model run difference E–D: thermal anomaly with paleoclimatic forcing

In order to compare the different thermal fields and to test the model sensitivity depending on the Malm k scenario under paleoclimatic conditions over time, we calculate the T difference between model runs E and D (Fig. 3). In other words, we subtract the thermal field of the low Malm k scenario from the one of the high Malm k scenario. The T difference between both k scenarios (Fig. 14) is more illustrative than the absolute thermal field shown in Fig. 11. A visible cold anomaly has developed in the northern part of the Molasse wedge after 20 ka of simulation time (110 ka ago) and gets colder with time, reaching a near steady-state minimum T of − 37 °C after 50 ka of simulation time (80 ka ago; Fig. 14a). Cooling is thus very slow and the average cooling rate amounts to about 0.74 and 0.34 K ka⁻¹ within the first 50 and 110 ka, respectively. Later, after 110 and 130 ka of simulation time (i.e. 20 ka ago and at present), it extends through the Malm aquifer and beyond into the Variscan basement (Fig. 14b, c). The maximum T difference generated between high and low Malm k scenarios amounts to about − 38 °C (Fig. 14). The distinct negative thermal anomaly of the high Malm k scenario (i.e. in model run E) continues cooling after the last glaciation. In addition, the Inn Glacier that advances periodically onto the MB produces the cold anomaly observed at the southern boundary of the section after 110 and 130 ka of simulation time (i.e. 20 ka ago and at present; Fig. 14). The Inn Glacier imposes a high hydraulic head and so squeezes a larger volume of cold meltwater into the subsurface in the high Malm k scenario. The yellowish colours in Fig. 14 represent the upflow of warmer water through the Malm aquifer towards the Danube in the high Malm k scenario. This shows again that cooling of the Molasse in the downflow zone results in heat transfer towards the Danube in the upflow zone.

Fig. 14

Temperature difference between both scenarios of high and low Malm permeabilities under paleoclimatic conditions (i.e. difference between model runs E and D introduced in Fig. 3): a 80 ka ago (time = 50,000 a), b 20 ka ago (time = 1.1E5 a), and c at present (time = 1.3E5 a). Vertical exaggeration equals three

Model run difference G–F: thermal anomaly without paleoclimatic forcing

Figure 14 shows that the negative thermal anomaly already exists after 50 ka of simulation time (i.e. 80 ka ago), before surface temperatures fall below the freezing point and permafrost starts to form (Fig. 5), or the first foreland advance of the Inn Glacier (Fig. 6).We therefore tested the model sensitivity with respect to the applied T_S condition. Consequently, in order to evaluate the specific effect of periglacial climate conditions and thus paleo-T forcing on coupled, gravity-driven groundwater and heat flow, we also simulated both k scenarios with a constant T_S of 9 °C at 300 m altitude as a T boundary condition in model runs F and G, respectively. Figure 15 shows the T difference between the corresponding k scenarios without paleoclimate forcing, i.e. the difference between model runs G and F (Fig. 3). As anticipated by the results after 50 ka of simulation time shown in Fig. 14a, the calculated T difference shown in Fig. 15 equally reveals a cold anomaly centred on the MB. Even without paleoclimatic conditions applied, temperatures in the MB are up to 50 °C colder in the high Malm k scenario. Here, the average cooling rate amounts to about 0.82 and 0.44 K ka⁻¹ within the first 50 and 110 ka, respectively. The negative thermal anomaly still shows extremely slow cooling after 110 ka of simulation time. In our 2D model section, without any perpendicular groundwater flow along the Danube, the water ascending through the Malm aquifer in the high Malm k scenario is up to 18 °C warmer near the Danube compared to the low Malm k scenario (Fig. 15a).

Fig. 15

Temperature difference between both scenarios of high and low Malm permeabilities with a constant surface-temperature condition (i.e. difference between model runs G and F introduced in Fig. 3): a 80 ka ago (time = 50,000 a), b 20 ka ago (time = 1.1E5 a), and c at present (time = 1.3E5 a). Vertical exaggeration equals three

Model run differences D–F and E–G: thermal effect of different surface temperature conditions

Based on the low Malm k scenario, we calculate the T difference between model runs D and F (Fig. 3) in order to test the model sensitivity and to compare the thermal fields obtained by applying the paleo-T and the constant T_S, respectively (Fig. 16). Based on the high Malm k scenario, we equally calculate the T difference between model runs E and G (Fig. 3) in order to compare the thermal fields obtained by applying the paleo-T and the constant T_S, respectively (Fig. 17). Both combined sections reveal a negative T difference about 100 °C colder in the southernmost part of the models. Moreover, in both combined sections, the applied paleoclimatic conditions cause a few degrees C colder temperatures in the shallow subsurface at the present day (Figs. 16c and 17c). These are an aftereffect of permafrost formation observable during the LGM about 20 ka ago (Figs. 16b and 17b). In the low Malm k scenario, a relatively small, up to 12.5 °C warm thermal anomaly forms in the southern part of the model (Fig. 16). In the high Malm k scenario (Fig. 17), a prominent positive thermal anomaly forms in the central segment of the MB close to the Malm aquifer. The T difference reaches its maximum of about 25 °C during the LGM about 20 ka ago (Fig. 17b). It is presently weakening and thus cooling very slowly due to resumed groundwater flow. From 20 ka ago to the present-day, the average cooling rate amounts to about 0.5 K ka⁻¹.

Fig. 16

Temperature difference between the paleo-temperature condition and the constant surface-temperature condition in the low Malm permeability scenario (i.e. difference between model runs D and F introduced in Fig. 3): a 80 ka ago (time = 50,000 a), b 20 ka ago (time = 1.1E5 a), and c at present (time = 1.3E5 a). Vertical exaggeration equals three

Fig. 17

Temperature difference between the paleo-temperature condition and the constant surface-temperature condition in the high Malm permeability scenario (i.e. difference between model runs E and G introduced in Fig. 3): a 80 ka ago (time = 50,000 a), b 20 ka ago (time = 1.1E5 a), and c at present (time = 1.3E5 a). Vertical exaggeration equals three

Discussion

Interplay of Malm permeability, gravity-driven groundwater flow, and paleoclimate

Although we did not simulate 3D groundwater flow in the MB, we estimate that the results along our 2D model section through the Wasserburg Trough to the east of Munich are significant. The reason is the existence of the Wasserburg Trough anomaly, which suggests the dominance and a certain independence of the thermal-hydraulic regime caused by the extraordinary high k of the karstified Malm. Specifically, the studied 2D section is appropriate to reproduce maximum hydraulic gradients across the MB, from the Alps in the south towards the Danube in the north. Topography along the section reaches up to 600 m compared to 300 m at the Danube and 440 m at the Inn Rivers, which is why forced convection induces a large downflow zone in the MB. In general, groundwater flow is directed from the recharge to the discharge areas and always in the direction of decreasing hydraulic head (Fig. 7). The particular situation along our section, where a topographic high generates two opposite cells of forced convection, actually illustrates a groundwater divide between the Danube and Inn Rivers (Fig. 7b, d). A similar situation exists between the Danube and the Rhine Rivers in the MB near Lake Constance, where the potentiometric surface provides evidence that the Rhine drains the Malm aquifer in the western part of the MB (Bertleff et al. 1993; Stober and Villinger 1997). Figure 7a, b shows that in the low Malm k scenario the present-day freshwater hydraulic head corresponds to about 400–500 m altitude within large parts of the Malm aquifer. In Fig. 7c, an explanation for the stronger freshwater hydraulic head decrease in the high Malm k scenario is the faster dissipation of the hydraulic pressure induced by topography and the Inn Glacier. When permafrost occurs, notably in model runs C or E, the freshwater hydraulic head also decreases fast beneath the permafrost layer, because it hampers recharge (e.g. Sloan and van Everdingen 1988) while groundwater stored within the Molasse sediments continues flowing towards the Malm aquifer. The water table in the low Malm k scenario rises significantly faster aside the Danube, and is thereby less sub-hydrostatic than in the high Malm k scenario, because the lower Malm k reduces the hydraulic connection with the Danube. By contrast, the freshwater hydraulic head along the highly permeable Malm aquifer is clearly sub-hydrostatic, which is in agreement with Lemcke and Tunn (1956) and Andres and Frisch (1981) who observed a water table up to 200 m below terrain level towards the Alpine front. Accordingly, the high Malm k scenario in Fig. 7c and d indicates that the hydraulic head within the Malm aquifer in the central part of the section can theoretically be as low as the Danube in case the k contrast between Malm and Molasse is very high. Frisch and Huber (2000) indicate a relatively low freshwater head of 370 m in the Munich region, reaching southwards as far as the Starnberger See and being only 40 m higher than the Danube near Regensburg over a distance of more than 100 km. A similar low freshwater head over 80 to 100 km distance results from our high Malm k scenario (Fig. 7d). Therefore, the low hydraulic head values in the Malm aquifer confirm the excellent hydraulic connection with the Danube, which was postulated by Frisch and Huber (2000) and Stober et al. (2014) because of the existence of freshwater with a mineralization below 1 g L⁻¹ in the Malm aquifer in large areas of the MB.

Mraz et al. (2019) also recently provided further evidence that a continuous influence of descending meteoric fluids explains the low salinity and the isotope values observed within the Malm. Mraz et al. (2019) studied the fluid and T evolution of the Malm aquifer by analysing fluid inclusions and performing stable isotope measurements of different cement phases in carbonate rock. Their results reveal that the least saline fluids from the Malm occur relatively close to the Alpine front, namely in the Traunreut well near Lake Chiemsee to the east of the Inn valley. In addition, they conclude that the reservoir T must have decreased more strongly in the eastern part (Traunreut area) than in the central part (Munich area) of the MB. Groundwater dating performed by Heidinger et al. (2019) using Krypton isotopes are in agreement with a relatively dynamic groundwater recharge of the Malm aquifer across the thick Molasse in the Wasserburg Trough region, potentially under subglacial conditions. Strikingly, Tertiary reservoir layers located to the south of the Wasserburg Trough anomaly (Fig. 2) and within the conduction-dominated zone close to the Alpine Front comprise numerous oil and gas fields (Bachmann et al. 1982, 1987; Bachmann and Mueller 1992) and are known for their significant overpressure (Drews et al. 2018). Regarding the impact of groundwater leakage on hydrocarbon reservoirs, as well as brines in the Tertiary Molasse and Cretaceous sediments, we refer to Lemcke (1977), Lemcke and Tunn (1956) and further Andres and Frisch (1981). The latter authors argue that oil and gas have largely been flushed out and destroyed by groundwater flow except for their preservation in sheltered structural traps (see e.g. Bachmann et al. 1982). Similarly, mineralized formation waters have progressively been diluted except for salt accumulations in isolated zones (Andres and Frisch 1981). Hydrocarbon trapping and hydrodynamic regimes are not mutually exclusive, but are frequently coexisting and interrelated in many geological basin types (Wendebourg et al. 2018).

The model results about periodic reorganization of groundwater flow in response to glacier cover (Figs. 7, 10) are in accordance with findings by Boulton et al. (1995) and Piotrowski (1997a) who studied the large-scale patterns of groundwater flow beneath ice sheets. Glacial cycles produce cyclical pulses of pressure and groundwater flow through the subsurface. Thus, glaciations represent the strongest impacts on groundwater during recent geological time. Figure 7 best shows how the model results recreate the observations of Boulton et al. (1995) and Piotrowski (1997a) by exemplifying the reorganization of the relatively small, topographically controlled groundwater catchments typical of interglacial periods, such as the present, to produce basin-scale integration of groundwater flow when glaciers advance into the foreland. Specifically, the hydraulic effect of the Inn Glacier completely suppresses the convection cell of groundwater ascending towards the Inn River and thereby concentrates groundwater flow towards the Danube River (Fig. 7). The consequence is significant groundwater flow in deep aquifers (i.e. in the deep Molasse and in the Malm) and the complete replacement of pre-existing groundwater by colder water in shallow aquifers, where permafrost is absent, such as in the Quaternary and the shallow Molasse beneath larger perennial rivers and lakes, and probably beneath piedmont lobes of glaciers (Fig. 10b). Similar to the present study, Lemieux et al. (2008) numerically simulated the surface–subsurface water exchange flux in Canada and demonstrated that large quantities of meltwater, i.e. on average 43%, episodically infiltrate beneath glaciated areas. However, Piotrowski (1997b) finds evidence that relatively fine-grained substratum like in the MB can drain only about 25% of subglacially produced meltwater. Bertleff et al. (1993) studied the paleogroundwaters that are found in various Tertiary aquifers (i.e. in the Molasse) and in the karstified Malm. According to Bertleff et al. (1993), the interpretation of hydrochemical and hydroisotopical measurements indicate a high altitude origin of the water and a recharge of the intermediate and regional (i.e. deep) circulation systems under colder climatic conditions than today. They outlined a possible hydrodynamic model of groundwater recharge in which groundwater recharge during glacial periods occurred in distinct subglacial environments. Also in the present study, recharge areas under cold conditions in presence of permafrost are essentially confined to deep, relatively permafrost-free glacial basins where glacier cover causes an increased recharge due to a higher water table or hydraulic head (Fig. 10a), which supports similar results obtained by Bertleff et al. (1993). Similar to previous conclusions based on hydraulics, hydrochemistry and isotope signature, our model results generally confirm the subdivision of the Malm aquifer in two main zones (Bertleff and Watzel 2002). The zone along the northern basin margin contains recently regenerated groundwater, whereas the zone in the central and southern parts of the MB contains water of Pleistocene age characterized by a very low regeneration. The Péclet number calculation generally confirms the interpretation by Bertleff et al. (1993) and Bertleff and Watzel (2002) that the southern, deepest part of the Malm is characterized by a relatively stagnant groundwater system dominated by heat conduction. This groundwater can only be renewed by leakage through the Tertiary Molasse (Fig. 7), whereby an ion-exchange water composition formed (Bertleff et al. 1993; Birner et al. 2011).

The region of the Wasserburg Trough to the east of Munich constitutes a remarkable exception. As shown in the present study, the higher k of the karstified Malm regionally induces significant leakage of meteoric water through the Quaternary, Tertiary and Cretaceous overburden, and thereby causes a faster groundwater regeneration and thus an overall younger groundwater age in this more central zone of the MB. Our results thereby support the recent findings of Heidinger et al. (2019) who performed dating of deep Malm groundwater using krypton isotopes and suggested that Malm groundwater recharge occurred during the last glaciation in view of groundwater age. However, in the bulk of the MB, recharge under periglacial conditions certainly was appreciably lower compared to the present temperate climate, mainly because of reduced precipitation (Cohen et al. 2018; Seguinot et al. 2018) and recurrent, persistent permafrost. Therefore, contrary to Bertleff et al. (1993) who put forward that recharge rates under recent climatic conditions are comparatively low, we suggest that recharge in the MB generally was lower during colder climatic conditions. According to our model results for the high Malm k scenario (Fig. 7c, d), it is possible that cold meltwater from the Inn and Achen glaciers (Van Husen 1987) infiltrated down to the Malm aquifer under an increased hydraulic head in the course of several glaciations. However, meltwater infiltration apparently did not produce an appreciable cold anomaly in the deeper MB (Fig. 10b). Equally, the much larger Rhine Glacier (Cohen et al. 2018) apparently does not produce a significant cold anomaly in the deeper MB and the Malm aquifer (Fig. 2) (Agemar et al. 2012; Agemar and Tribbensee 2018). This obvious fact shows that the hydraulic and thermal effect of relatively large, recurrent foreland glaciers is not sufficient and causative for significant cold anomalies such as the Wasserburg Trough anomaly in the relatively deep Malm aquifer.

Controls of Malm permeability on heat transport and thermal field

In both Malm k scenarios and under the paleoclimatic conditions of the last glaciation, groundwater and heat upflow within the Malm aquifer due to forced convection causes the distance between isotherms to increase at medium depth in the upper basement beneath the upflow zone of the Malm aquifer (Fig. 11). This results eventually in a higher geothermal gradient and higher temperatures at shallow depth in and to the south of the Danube valley. Evidence for the regional discharge of the Malm aquifer into the Danube in the region between Neustadt and Regensburg had already been provided by Bertleff et al. (1993), Birner et al. (2011), and Frisch and Huber (2000). Groundwater from the MB circulates deep enough to take up heat and discharge it in the Danube valley. At Bad Gögging, which is located 30 km downstream of Ingolstadt, a 120-m-deep well produces 14 °C warm water in the Malm aquifer (Wirth and Andres 1981), which is significantly warmer than the mean annual T_S of about 8 °C ((Hänel and Staroste 1988). According to our model, the 25 °C-isotherm is located just about 150 m below the surface to the south of the Danube valley (Fig. 11). Temperatures of 25 °C in the shallow subsurface of the MB or Malm aquifer do not exist, mainly because a significant groundwater flow is expected in the Malm parallel to the Danube (Frisch and Huber 2000; Heine and Einsiedl 2021). Nevertheless, the top Malm T presented in GeotIS (Fig. 2) (Agemar et al. 2012, 2014; Agemar and Tribbensee 2018; LIAG 2021) shows that temperatures are higher by up to 20 °C in the realm of Pfaffenhofen an der Ilm and Landshut, a region located about 40 km to the south of the Danube (Fig. 2). Only a 3D model could reveal thermal-hydraulic processes perpendicular to our 2D section.

In order to validate our results, we compare the isotherms of our Malm k scenarios to representative 2D sections (Fig. 18) extracted from the 3D T model of the MB implemented in GeotIS (Agemar et al. 2012, 2014; Agemar and Tribbensee 2018; LIAG 2021). A first section with isotherms and eight projected deep wells, generated in GeotIS just to the west of Munich (Fig. 18a), is comparable our low Malm k scenario (i.e. model run D; Fig. 11b). The subhorizontal and equidistant isotherms in the central part of the MB represent a thermal field typical for a conduction-dominated regime. By comparison, our low Malm k scenario shows a thermal perturbation in the southern part of the model section, which is related to the convection cell induced by deep flow towards the Inn valley in the south (Fig. 11b). The extracted GeotIS section shows a typical ramp towards the Alpine front (Fig. 18a) and therefore the mentioned convection cell is not present. In addition, Molasse k may be lower near the Alpine front due to a higher horizontal and vertical stress and thus a higher compaction. A second section with isotherms and numerous projected deep wells, generated in GeotIS along the trace of our 2D section to the east of Munich (Fig. 18b), is comparable to our high Malm k scenario (i.e. model run E; Fig. 11d). The curved and irregular isotherms are characteristic for a thermal field in a convection-dominated regime. Like our results of the high Malm k scenario, the GeotIS section shows three regions with distinct geothermal gradients (Fig. 18b). Beneath the Franconian Alb to the north of the Danube valley, the isotherms are subhorizontal and show a regular spacing, which is typical for conduction in this northern part of the model dominated by the Variscan basement. In the northern part of the MB, higher temperatures exist at shallow depth and thus produce a higher geothermal gradient. Here, ascending groundwater in the Malm aquifer due to forced convection is responsible for heat upflow. In the southern part of the MB, numerous deep wells mainly drilled by the oil and gas industry confirm a large downflow zone by leakage of relatively cold surface water through the Molasse into the underlying Malm aquifer (Fig. 18b).

Fig. 18:

2D sections extracted from the 3D temperature model (Agemar et al. 2012) implemented in the Geothermal Information System GeotIS (Agemar et al. 2014; LIAG 2021), a being representative of the low Malm permeability known to the west of Munich, and b being representative of the high Malm permeability known in the Wasserburg Trough to the east of Munich. The section locations are indicated in the map slices above each 2D section and in Fig. 2. Full circles in the map slices and schematic drilling rigs with vertical drill paths in the 2D sections indicate boreholes with temperature information. Vertical exaggeration equals five

With regard to model sensitivity to changes in k, the Molasse k is worth mentioning because the Molasse also has a different and complex composition along strike (GeoMol Team 2015; Kuhlemann and Kempf 2002; Schwerd et al. 1996; Zweigel et al. 1998), and presents a varying number of typically more permeable Tertiary sandstone layers (Brink et al. 1992). Bertleff et al. (1993) indicate a very low k range of 10⁻¹⁶ to 10⁻¹⁸ m² for the Lower Freshwater Molasse, which consists mainly of lacustrine and fine-grained fluvial sediments. Comparatively low average permeabilities of the Molasse and a very low recharge by leakage (Bertleff and Watzel 2002) seem to be representative of the western MB in Baden-Württemberg where only a single sandstone reservoir layer is known (Brink et al. 1992). Villinger (1977) specified a k range of 10⁻¹⁵ to 10⁻¹⁶ m². Andres and Frisch (1981) discussed various estimates available in literature and derived a (vertical) Molasse k of 2∙10⁻¹⁶ m² hydraulically representative of an equivalent silt composition of the Molasse sediments. Although the focus of our study lies on the variation of the Malm k, we checked the sensitivity of model runs F and G to a Molasse k of 10⁻¹⁵ m² horizontally and 10⁻¹⁶ m² vertically, which corresponds to a reduction of Molasse k of one order of magnitude (Table 1). In this case, the T difference, in other words the simulated Wasserburg Trough anomaly, between the modified model runs G and F reaches only about − 10 °C. Clearly, in our models, the lower constraint needed to produce the observed negative thermal anomaly corresponds to a minimum Molasse k of 10⁻¹⁴ m² horizontally and 10⁻¹⁵ m² vertically. The latter conclusion suggests that permeable Tertiary sandstone layers within the Molasse are not only a prerequisite for hydrocarbon trapping, but also may be essential for hydrodynamic regimes able to significantly perturb the thermal field. This would support the observation that hydrocarbon trapping and hydrodynamic regimes are frequently coexisting and interrelated in many geological basin types (Wendebourg et al. 2018).

In agreement with the geological setting, the heat input into the MB is characterized by conduction through the Variscan basement, which is characteristic of the foreland-basin play type (Moeck 2014; Moeck and Beardsmore 2014; Moeck et al. 2015a). Our model results confirm the main findings of thermal-hydraulic modelling performed by Przybycin et al. (2017). The Malm and the Tertiary sediment fill of the MB present a sufficiently high formation k for basin-wide gravity-driven fluid flow and convective heat flow. Since the range of Darcy’s velocity magnitude is mostly between 0.01 and 1 m a⁻¹ (Fig. 8), and thus relatively low, Darcy’s law remains valid. Przybycin et al. (2017) suggest that heterogeneities in thermal conductivity contribute fundamentally to the generation of pronounced positive and negative thermal anomalies in the MB. We put forward high formation k and pervading gravity-driven groundwater flow as key factors for significant heat redistribution, because considerable groundwater flow notably in the Wasserburg Trough area is the dominant heat transport mechanism as indicated by the high Péclet numbers (Fig. 9b). In their model of the Upper Austria–Upper Bavaria area, Pfleiderer et al. (2016) concluded that hydraulic and geothermal gradients show similar regional trends caused eventually by the strong effect of groundwater convection on the redistribution of heat in the Malm karst aquifer.

Surface Q data are available from several geothermal atlases (Hänel and Hurter 2002; Hänel and Staroste 1988; Legrand et al. 1979) and updated maps (e.g. Cloetingh et al. 2010). However, the resulting surface Q maps differ largely. A main reason is that surface Q maps generally correspond to interpolated data points and do neither include geological information nor consider physical processes. Data compiled by Hänel and Hurter (2002) yield surface Q values ranging from about 90–95 mW m^–2 in the northwest in the realm of Stuttgart and the Danube River, and values of about 70 mW m⁻² in the realm of the German–Austrian border and the Alpine front. The European Q map by Cloetingh et al. (2010) shows a south to north Q range of 70 to 85 mW m^–2 for the Munich region studied here. The published surface Q values are in the range of values typical for the geodynamic setting of Paleozoic extended crust in Europe (Balling 1995; Norden et al. 2008). The low Malm k model scenario, which we regard as representative of the largest part of the MB, quite confidently reproduces the surface Q pattern presented by Hänel and Hurter (2002). For the high Malm k scenario, our model results clearly demonstrate that the published conductive surface Q values, which are typically derived from T measurements in subvertical wells, are inevitably afflicted with a significant convective Q perturbation due to groundwater flow towards the highly karstified Malm (Fig. 12b). Most importantly, our results provide evidence that Q distribution is heterogeneous and depth-dependent. Therefore, in regions dominated by (free or forced) convection, the determination of background Q, crucial as a constraint for thermal-hydraulic models, is not straightforward. The complexity obviously increases with the magnitude of groundwater flow and convective heat transfer (Fig. 12b). The published Q maps do obviously not reflect the results of the high Malm k scenario. The 2D section extracted from GeotIS (Fig. 18b) shows that a high number of sufficiently deep boreholes with T data mostly from oil and gas exploration exist in the realm of the Wasserburg Trough thermal anomaly. However, for only a minor fraction of these boreholes, Q values have been derived. Therefore, although a wealth of T data are available for a reliable subsurface T prognosis for the entire MB (Agemar et al. 2012; Agemar and Tribbensee 2018), the quantity and quality of Q data are inadequate for a clear representation.

Controls of permafrost and paleoclimate on heat transport and the thermal field in relation to different Malm permeability scenarios

Model run D predicts a maximum permafrost thickness of 150–200 m in the low Malm k scenario (Fig. 11a). These findings are in reasonable agreement with previous modelling results about maximum permafrost thickness in Central and Western Europe. For the last glaciation, numerical modelling satisfactorily shows a permafrost thickness of less than 100 m in the Paris Basin (Lebret et al. 1996), more than hundred metres in north-central Europe (Delisle et al. 2003; Govaerts et al. 2016), and up to 150 m to the north of the Rhine Glacier in Switzerland (Cohen et al. 2018). This is much less than several hundred metres up to some 500 m of permafrost currently subsisting at high latitude in Alaska and Canada (Lebret et al. 1996; Sloan and van Everdingen 1988). The main factors influencing the permafrost thickness are the heat flow density, mean annual ground temperatures, and thermal properties of rocks related to their lithology and their water content (Lebret et al. 1996). Since ground temperatures and rock properties are identical in both Malm k scenarios, the significantly higher maximum permafrost thickness of up to 600 m in the high Malm k scenario in model run E (Fig. 11c) can largely be explained by the drastically reduced conductive heat flow density due to extensive downward groundwater flow (Fig. 12b). Precise and reliable fossil indicators of past deep permafrost are difficult to find and direct evidence often linked to physical markers like indications of frost cracking (Lebret et al. 1996) or disturbed geothermal gradients, like e.g. in the Alberta Basin in Western Canada (Gray et al. 2012; Majorowicz 2012; Majorowicz et al. 2012). In the realm of the Amper and Isar valleys, we implemented a zone where the minimum T_S is limited to 4 °C to prevent permafrost formation. Permafrost is often absent beneath lakes and larger perennial rivers because deep and flowing water locally inhibits freezing of the subsurface (Boulton et al. 1993; Dagher et al. 2014; Sloan and van Everdingen 1988). Besides the general leakage of shallow groundwater into the karstified Malm, unhindered infiltration of cold groundwater in the Amper and Isar valleys even at cold surface temperatures may explain the local correlation of the Amper and Isar valleys with the coldest zone of the underlying thermal anomaly (Fig. 2) (Agemar and Tribbensee 2018). As Bertleff et al. (1993) and Bertleff and Watzel (2002) put forward, our models corroborate the fact that recharge during cold periods can only occur when and where permafrost is absent (compare Figs. 7 and 11). Indeed, a frozen sediment has a hydraulic conductivity many orders of magnitude lower than the same sediment in an unfrozen state (e.g. Piotrowski 1997a, b) because the dynamic viscosity of water drastically increases when it becomes ice, taking a value of some 10¹² Pa s, about 10¹⁵ times that of water (Fowler 1997). Permafrost is absent when T_S is cold but above the freezing point (Fig. 5) and in subglacial basins when periodic glacier cover impedes permafrost formation in the melting zone (Boulton et al. 1993; Piotrowski 1997a). Likewise, our models promote the assumption by Bertleff et al. (1993) that the discharge areas of the Malm between Neustadt and Regensburg remained active even during pleniglacial times, making the Malm aquifer the only possible drainage system in the MB. The present modelling implements a suggestion by Pfleiderer et al. (2016) to set up a transient conductive model for consideration of paleoclimatic signals and subsequently a coupled thermal-hydraulic model to test possible hydraulic concepts.

We regard the simulation of the low Malm k scenario, by contrast to the high Malm k scenario representative of the Wasserburg Trough area (Fig. 11), as benchmark of the primarily conduction-dominated western parts of the MB. Based on model runs D and E, Fig. 14 shows the resulting T difference under paleoclimatic conditions and thus visualizes the strongly modified thermal field. As shown by Birner et al. (2012), the range of hydraulic conductivities or permeabilities spans three to four orders of magnitude between the least permeable, marl-dominated Helvetic facies in the southwest and the most karstified Malm carbonate in the east of the MB. The much higher Malm carbonate k due to karstification drastically increases the gravity-driven coupled groundwater and heat flow, and thus heat redistribution. The higher the Malm k, the stronger and deeper the gravity-driven groundwater flow, which by descending through the Molasse towards the Malm aquifer also significantly cools down the Molasse as well as the Malm aquifer. The geostatistically evaluated thermal field, specifically the top Malm T and the Wasserburg Trough anomaly of about − 40 °C (Fig. 2) (Agemar et al. 2012; Agemar and Tribbensee 2018), validates the thermal anomaly of − 38 °C produced in our study (Fig. 14). Based on plausible permeability values and boundary conditions, our model runs D and E are able to replicate the cold anomaly occurring in the Wasserburg Trough to the east of Munich.

Subglacial infiltration of meltwater during glaciations has repeatedly been considered as the main recharge mechanism of the Malm aquifer in the deeper MB (Bertleff et al. 1993; Bertleff and Watzel 2002; Heidinger et al. 2019). If this hypothesis is true, we expect that the Wasserburg Trough anomaly occurring in the MB and reaching as deep as the Malm aquifer (Agemar and Tribbensee 2018) is a direct consequence of subglacial recharge by the Inn Glacier during the last glaciations(s). Based on model runs D and E, the T difference between both Malm k scenarios under paleoclimatic conditions show that a significant negative thermal anomaly of about − 37 °C forms within 50 ka (Fig. 14a). Strikingly, based on model runs F and G without paleoclimatic conditions applied, temperatures in the Molasse decrease to about − 41 and − 49 °C at 50 ka and 110 ka, respectively, in the high Malm k scenario (Fig. 15). The main explanation for this unexpected and counterintuitive result is that the extreme conditions linked to cold periods such as glacier cover and permafrost formation are relatively sporadic events in comparison to the geological time scale of persistent groundwater flow. Over the 130 ka simulated in the present study, the advances of the Inn Glacier into the foreland and onto the MB are sporadic (Seguinot et al. 2018). The duration of glacier cover makes up only about 26.6 ka (Fig. 6), which represents 20% of the simulated time and 25% of the duration of the Würm Glaciation (Fig. 5). Similarly, the duration of sustained temperatures below the freezing point and thus the duration of continuous permafrost makes up only about 45.2 ka (Fig. 5), which represents 35% of the simulated time and 44% of the duration of the Würm Glaciation (Jouzel et al. 2007). Unlike these relatively sporadic events, the k of geologic units is a permanent, primary geological control factor of the magnitude of forced convection and heat redistribution. Moreover, formation k already controlled coupled groundwater and heat flow probably since the end of MB development in the Late Miocene. The existence and evolution of the Pre-Danube since about 7 Ma (Kuhlemann and Kempf 2002) and of the present-day Danube River (Doppler et al. 2011; Fiebig and Preusser 2003; Jerz and Peters 2002; Schaefer, 1966) is equally important because it strongly influences the discharge of thermal water by the Malm karst between Neustadt and Regensburg. River systems are highly sensitive to environmental changes, including the effects of tectonic, climatic, glacial and anthropogenic forcing (Cordier et al. 2017). In particular, the glaciations in the Alpine realm directly or indirectly control the changes in the drainage pattern of the Danube. It is therefore possible that the thermal anomaly observed to the east of Munich (Agemar et al. 2012; Agemar and Tribbensee 2018) progressively developed over the last 800 ka with the onset of the first major glaciations recorded in the North Alpine Foreland and as a consequence of the dynamic reaction and evolution of the Danube.

Sensitivity of the individual Malm permeability scenarios to paleoclimate and permafrost

Figures 16 and 17 visualize the sole effect of the two contrasting T_S boundary conditions, including permafrost formation and recurrent glacial advances under paleoclimatic conditions, individually on the low and high Malm k scenarios. The explanation of the negative T difference of about 100 °C in the southernmost part of the section is the predicted ascent of warm water towards the depression of the Inn valley near Rosenheim in the scenarios under a constant T_S. By contrast, in the scenarios with paleoclimatic conditions applied, permafrost formation and the hydraulic head induced by the recurrent advance of the Inn Glacier inhibit this ascent of warm water by repeated reorganization of the groundwater flow systems. We directly compare Fig. 17 to Fig. 10 and thereby show, besides the hydraulic and thermal effect of the Inn Glacier, the hydraulic blanketing effect of permafrost, which in turn amplifies the overall effect of the Inn Glacier. In the low Malm k scenario, a relatively small, up to 12.5 °C warmer thermal anomaly forms in the southern part of the model (Fig. 16). The latter likely is an expression of two combined processes leading to a reduced hydraulic conductivity: (1) a primary low Malm k and (2) a higher dynamic viscosity of the infiltrating water due to the lower surface temperatures. First, the relatively low k yields an overall more conductive, warmer model, and second, further reduction in hydraulic conductivity results in even less heat redistribution in the deepest part of the MB. In the high Malm k scenario (Fig. 17), the most surprising result is the prominent positive thermal anomaly in the central segment of the MB and near the Malm aquifer (Fig. 17b). It reaches a maximum T of about 25 °C during the LGM about 20 ka ago when paleo-T is the coldest (about − 5 °C) and permafrost is the thickest. The reason is hydraulic shielding due to the virtually impermeable permafrost, which reduces groundwater flow and heat redistribution. Figure 17 actually captures the interplay of extensive groundwater flow and the hydraulic blanketing effect of permafrost. Up to the present day, despite the on average higher constant T_S of 9 °C compared to an average of 3.2 °C and a median of 2.3 °C for paleo-T (Fig. 5), the thermal anomaly caused by the model with paleo-T forcing is paradoxically about 16 °C warmer than the model with a constant T_S. Actually, the positive thermal anomaly related to the last glaciation is presently weakening and thus cooling very slowly due to resumed groundwater flow.

Clearly, as mentioned e.g. by Bertleff et al. (1993), Lemieux et al. (2008), Piotrowski (1997a) and Piotrowski (1997b), permafrost acts as a barrier for groundwater flow since it drastically reduces the hydraulic conductivity within sediments. Therefore, permafrost formation hampers heat redistribution by forced convection. In a sense, permafrost formation brings the high Malm k scenario closer to the low Malm k scenario. Moreover, in both combined sections in Figs. 16c and 17c, the signature of colder T_S locally remains in the shallow subsurface up to the present day with a T difference of more than − 2 °C. In both k scenarios, the colder zones in the shallow subsurface of the MB are remnants or an aftereffect of permafrost, which retards warming-up. Our results corroborate those obtained by Gray et al. (2012) who investigated the geothermal state of sedimentary basins, specifically the northern Alberta Basin, using oil industry thermal data. Our model shows that cooling of the shallow subsurface and recharge by glaciers essentially during colder climatic conditions is preserved in the low k scenario up to the present, mainly because it reacts slower to warming due to the dominant and relatively low rock thermal conductivity (Fig. 16). As a matter of fact, low k lithologies tend to capture and preserve the long-term climatic evolution of the Pleistocene (e.g. Bertleff et al. 1993), which is clearly dominated by colder climatic conditions in the last 800 ka if not the entire Quaternary (Ehlers and Gibbard 2008). More permeable geologic units, in particular the karstified Malm, and the corresponding groundwater flow systems react much faster or dynamically to warm periods like the Holocene (Bertleff et al. 1993; Bertleff and Watzel 2002). Thereby, the signature of colder climatic conditions is largely overprinted (Fig. 17).

Implications for the geothermal field and potential in the Molasse Basin

The present study aims at providing results that elucidate the main fluid and heat transport processes underlying the formation of the cold anomaly in the Wasserburg Trough by providing an explanation for the existence of 40 °C lower temperatures in the Malm aquifer (Fig. 2) (Agemar et al. 2012; Agemar and Tribbensee 2018). The cataloguing of geothermal resources according to geological criteria and heat transport mechanisms to assess reservoir quality is the central motivation for the recent geothermal play-type concept (Moeck 2014; Moeck and Beardsmore 2014; Moeck et al. 2015a).

Clauser et al. (2005) concluded that no significant heat redistribution due to fluid flow occurs in the Malm in the western part of the MB. This situation refers to the low Malm k scenario. In the eastern part of the MB, in particular in the Wasserburg Trough between Munich and the Landshut-Neuötting High, however, our results demonstrate a substantial heat redistribution. This is best visualized by the differing Péclet number distribution in the low and high Malm k scenarios (Fig. 9). According to Birner et al. (2012), zones in the Malm aquifer with hydraulic conductivities higher than 10⁻⁶ m s⁻¹ are generally suitable for geothermal energy production. However, the present study clearly shows that a sufficiently high Malm k induces a faster groundwater flow and thereby increased cooling of the MB down to the Malm aquifer. The about twice as high Darcy’s velocity magnitude in the high Malm k scenario leads to the heat redistribution shown in Fig. 12b, and eventually to the cold anomaly shown in Figs. 14 and 15. Therefore, since the thermal output of geothermal plants is directly proportional to both production rates and reservoir temperature, we suggest that the thermally and hydraulically most productive area for geothermal energy projects in the MB, and thus the area with the highest geothermal potential, corresponds to the transition zone between the conduction- and convection-dominated areas of the Malm geothermal reservoir. In view of the map published by Birner et al. (2012), the transition zone corresponding to Malm hydraulic conductivities in the range of 10⁻⁵ and 10⁻⁶ m s⁻¹ comprises the Munich metropolitan region.

Conclusions

This study aimed at using numerical thermal-hydraulic modelling in the Southern German Molasse Basin along a representative 2D section through the cold anomaly occurring in the Wasserburg Trough to the east of Munich to assess the magnitude of gravity-driven coupled groundwater and heat flow under the paleoclimatic conditions of the last 130 ka, including the Würm Glaciation. Considering our low Malm permeability scenario as the conduction-dominated benchmark of the Molasse Basin, our convection-dominated high Malm permeability scenario is able to replicate the Wasserburg Trough anomaly. The findings are, however, unexpected concerning the fundamental reason of the observed thermal anomaly. The thermal-hydraulic model scenarios clearly show that neither colder surface temperature, nor permafrost formation, nor glaciers advancing onto the foreland explain the observed cold anomaly. The cause clearly is forced convection in the high Malm permeability scenario. Among the effects of cold periods, permafrost formation has the largest impact on the thermal-hydraulic regime, because it significantly reduces groundwater flow, and thereby forced convection. Permafrost produces not a thermal, but a hydraulic blanketing or shielding effect and effectively reduces heat redistribution. Conversely, models performed without paleoclimatic forcing, but with the constant surface temperature of 9 °C, yield a significantly colder thermal anomaly because groundwater flow reduction due to permafrost is absent. The sensitivity analysis reveals that the development of a negative thermal anomaly comparable to the known one of the Wasserburg Trough is relatively abrupt and forms within a Molasse permeability range of one order of magnitude, which implies a minimum Molasse permeability of 10⁻¹⁴ m² horizontally and 10⁻¹⁵ m² vertically in the high Malm permeability scenario. Our results are in agreement with known recharge of the Malm aquifer along the Danube and through the outcropping Malm north of the Danube, as well as discharge in the Danube valley between Neustadt and Regensburg. Secondly, considerable recharge and leakage through the overlying Quaternary and Tertiary sediments combined with significant drainage by the underlying Malm aquifer are essential to explain the reduction of heat flow density and the related depression of geothermal gradients. Clearly, the main cause and the primary geological control factor for the occurrence and magnitude of the Wasserburg Trough anomaly is the significantly higher Malm permeability due to karstification in the eastern (Bavarian) part of the MB and thus the magnitude of gravity-driven groundwater flow. Unlike previous suggestions that heterogeneities in thermal conductivity contribute to the generation of pronounced positive and negative thermal anomalies in the MB, we put forward that formation permeability and gravity-driven groundwater flow are key factors for heat redistribution, largely overprinting moderate disparities in thermal conductivity. Our results clearly disagree with previous findings of generally increased recharge under cold climate and especially under glacial conditions. Subglacial permafrost-free basins may have promoted increased recharge under a higher hydraulic head. However, in the bulk of the Molasse Basin, recharge under periglacial conditions certainly was appreciably lower compared to the present temperate climate, mainly because of reduced precipitation and recurrent, persistent permafrost. Finally, the comprehensive understanding of the fundamental thermal-hydraulic processes of geological systems over geological time scales requires the consideration not only of the main aquifer or geothermal reservoir, but also of the aquitards like the basement and the overburden.

Appendix: Explanations about model complexity and simplification

On the one hand, the simplification of the actual stratigraphic complexity to only three model units may be justified, if the chosen structural differentiation is representative of the main contrasts in thermal and hydraulic properties observed. On the other hand, the reproduction of the relevant thermal-hydraulic processes and their numerical implementation may pose serious additional constraints on the decisions about model complexity, since the mere model feasibility is at stake. Any model implementation is actually bound to the capabilities of the hardware and the software at hand as well as the kind of physical processes involved. Therefore, in order to test the possibilities of numerical implementation, we primarily built a generic model geometry with the basic three model units, which successfully ran including all thermal and hydraulic processes that we wanted to reproduce and provided the first valuable and meaningful results. We subsequently built a more detailed extended model geometry including six model units. Unfortunately, we encountered numerical problems with poor convergence and persistent model run aborts when we implemented the phase change of water (COMSOL 2020b), which we judge as crucial to simulate the thermal and hydraulic effects of permafrost on the thermal-hydraulic regime in the MB during the Würm Glaciation. Since we were unable to solve these major numerical issues, we deliberately decided to limit the model complexity and resolution. We took up on the simpler generic model, which nevertheless provides essential results about thermal-hydraulic processes for further consideration and future use in larger, more complex and 3D models. Overall, we estimate that the discrepancy between the generic model and the actual geologic situation is negligible, especially with regard to our focus on k scenarios and the main findings.