The origin of deep geothermal anomalies in the German Molasse Basin: results from 3D numerical models of coupled fluid flow and heat transport

The European Molasse Basin is a Tertiary foreland basin at the northern front of the Alps, which is filled with mostly clastic sediments. These Molasse sediments are underlain by Mesozoic sedimentary successions, including the Upper Jurassic aquifer (Malm) which has been used for geothermal energy production since decades. The thermal field of the Molasse Basin area is characterized by prominent thermal anomalies. Since the origin of these anomalies is still an object of debates, especially the negative ones represent a high risk for geothermal energy exploration. With our study, we want to contribute to the understanding of the thermal configuration of the basin area and with that help to reduce the exploration risk for future geothermal projects in the Molasse Basin. For this, we conducted 3D basin-scale coupled fluid and heat transport simulations to reproduce the present-day thermal field of the Molasse Basin by considering conduction, advection, and convection as heat-driving mechanisms. Within this paper, we show how the temperature distribution of the Molasse Basin, including the pronounced thermal anomalies, can be reproduced by coupled fluid flow and heat transport simulations following a multi-scale 3D-modelling approach. We find that the shallow thermal field is strongly affected by basin-wide fluid flow. Furthermore, we show that the temperature distribution at the depth of the Malm aquifer is strongly influenced by the hydraulic conductivity of the Foreland and Folded Molasse Sediments and that hydraulically conductive faults have only a minor influence on the regional temperature distribution. Moreover, we show that the positive and negative thermal anomalies are caused by the superposed effects of conductive and advective heat transport and correlated with the geological structure.

be extracted from different depths by geothermal power plants. The specific temperature pattern at depth of a geothermal area in turn is dependent on local-to-basin-wide geological, tectonic, and hydrogeological conditions and may vary strongly laterally. One area in Europe which is suitable for the production of geothermal energy is the European Molasse Basin.
The European Molasse Basin is a Tertiary foreland basin at the northern front the Alps and extends over parts of France, Switzerland, Germany and Austria (Fig. 1). The typically wedge shape of the basin evolved in consequence of the Euro-Adriatic continental collision and the rise of the Alps since the Cretaceous (Berge and Veal 2005). Today, the basin is filled with mostly clastic sediments originating from erosional processes of the Alps. Those sediments are partly overrun and folded due to the ongoing plate movement of the Euro-Adriatic continental collision (Schmid et al. 2008). The Molasse sediments are underlain by Mesozoic sedimentary successions, which include the Upper Jurassic Malm aquifer (Birner 2013) deposited in the Tethys Ocean. This Malm aquifer shows a high potential for geothermal energy production and is explored by the petroleum and geothermal industries as well as by the scientific community since decades (Sachsenhofer et al. 2006;Büchi et al. 1965).
The exploration efforts have led to a large amount of temperature measurements with which a picture of the subsurface temperature distribution of the basin could be generated by interpolation between the data points (GeotIS: Schulz et al. 2009;Agemar et al. 2014). This picture of the subsurface thermal field shows continuously increasing temperatures from north to south in the basin, as well as areas with deviating temperature patterns. In the latter, large temperature differences within a few kilometres of lateral distance (Agemar et al. 2012) can be observed, so-called positive and negative thermal anomalies. Even though the knowledge about the basin has advanced in response to studies of the industry and the scientific community over the last decades (Berge and  Böhm et al. 2011;Cacace et al. 2013;GeoMol Team 2015;Jodocy and Stober 2009;Kempf et al. 1999;Reischenbacher and Sachsenhofer 2011;Roeder and Bachmann 1996;Pamer and Diepolder 2010), the origin of this temperature pattern and the fluid flow directions in the basin are still not fully understood.
At present, 17 major geothermal power plants are operating in the German Molasse Basin (Agemar et al. 2014), of which the biggest part is located in the area around the city of Munich in Bavaria (Böhm et al. 2012). In this area, two major thermal anomalies are located, a positive and a negative one, with a temperature difference (ΔT) of up to 40 K at a depth of −3,500 m asl (GeotIS, Agemar et al. 2012). Since this strong temperature change over a small lateral distance cannot be explained satisfactorily with the present-day knowledge, it may reduce the willingness for further geothermal drilling projects in the Molasse Basin.
Moreover, a strongly non-uniform distribution of data in the Molasse Basin is additionally impeding the research and exploration: where the extraction temperature or fluid quantity of geothermal drillings has been disappointing in the past , exploration efforts have been reduced leading to a lack of data in these areas and vice versa. To overcome the obstacles of such data gaps, 3D models are useful tools which fill up areas of missing information by applying physically reliable interpolations between points of observation (Cacace et al. 2010). These models consider not only structural heterogeneities, but also different heat transport processes to advance the research in the area of the Molasse Basin.
Different studies have been conducted in the past, based on different physical assumptions and modelling techniques, to understand the temperature and pressure field of the European Molasse Basin, of which Frisch and Huber (2000), Birner (2013), Agemar et al. (2012), Przybycin et al. (2015b), Rühaak (2009Rühaak ( , 2015, and Rühaak et al. (2010) are just some examples. Frisch and Huber (2000) investigated the flow field of the Upper Jurassic Malm aquifer in the German part of the basin using hydraulic potential measurements and thermal water balancing to deduce the general flow direction. They propose a main fluid flow direction from the west to the east inside the Malm aquifer with a significant flattening of the hydraulic potential in the central basin area as well as discharging fluids into parts of the Danube River and in the east of the uplifted crust of the Landshut-Neuöttinger High (LNH). However, they could not explain the fluid flow pattern between the city of Munich and the LNH in the Malm aquifer, an area with prominent thermal anomalies at depth.
Such a fluid flow direction could be mainly confirmed by Birner (2013), who investigated the hydrochemistry and hydraulic regime of the German Molasse Basin with a hydrogeological model of the Upper Jurassic Malm aquifer with the aim to characterize the processes driving the groundwater dynamic and the mass transport. He states that even though the Malm aquifer shows a complex system of pores, joints, and karstrelated cavities, it can be described with lithological units of different but homogenous hydraulic characteristics on a basin scale.
The work of Birner (2013) has been incorporated into the Geothermal Information System (GeotIS, Schulz et al. 2009), a web-based platform integrating all available data about geological structure, lithology and temperature from outcrops, wells, and seismic data for areas of geothermal interest in Germany with the intention to reduce the exploration risk for geothermal drilling projects. Based on 3D interpolation between temperature measurements, maps have been generated, available at the GeotIS platform (Agemar et al. 2012), depicting the temperature distribution at different depths in the Molasse Basin. Unfortunately, the lateral coverage and reliability of these maps are decreasing with increasing depth due to the decreasing data density and with that increasing uncertainty of interpolation. Nonetheless, the thermal model of GeotIS (Agemar et al. 2012;Schulz et al. 2009) is so far the best representation of the thermal field of the Molasse Basin according to the implemented data. Since not much measured temperature data were available for this study, the GeotIS model has been used as a reference for validation.
A different approach to explain the thermal anomalies in the Molasse Basin has been used by Rühaak (2009Rühaak ( , 2015 and Rühaak et al. (2010). They investigated the groundwater flow regime and the thermal field of the western part of the Molasse Basin by using a quasi-steady-state conductive 3D-modelling approach on a local scale with the intention to compare calculated temperatures to measured values. Since they were not able to reproduce the observed thermal anomalies in the western Molasse Basin with a purely conductive approach of heat transport, they proposed E-W striking fault parallel fluid flow to cause lateral temperature differences of more than 10 K in the basin.
Such influence of hydraulically conductive major fault zones on the thermal field has also been studied by Cherubini et al. (2014) for the area of Brandenburg in north-eastern Germany using a basin-scale 3D numerical model of coupled fluid flow and heat transport. Their results, however, indicate only a limited lateral influence of permeable fault zones on the regional thermal field.
Another conductive heat transport approach was followed by Przybycin et al. (2015b), who used a lithospheric-scale 3D model covering the German Molasse Basin as well as the South German Scarpland to the north and parts of the Alps, including the Tauern Body to the south to calculate the present-day 3D thermal field. By choosing a model extent which integrates a large part of the continental collision zone, Przybycin et al. (2015b) investigated the interdependence of the basin with the mountain chain with respect to the thermal field. With their approach, they were able to reproduce the measured thermal field of the Molasse Basin, as shown by GeotIS (Agemar et al. 2012) to a certain extent. The conductively modelled temperature trend shows that the thermal anomalies in the Molasse Basin are generated to some extent by the structural configuration of the crust and the presence of the Tauern Body in the Alps. Przybycin et al. (2015b) propose an additional cooling influence of fluid flow in the sedimentary part of the basin as possible mechanism to explain the observed lower temperatures in the Molasse Basin, an effect already discussed for the Po Plain by Pasquale et al. (2013).
However, being aware of the complex temperature distribution at depth is not enough to reduce the exploration risk as the productivity is another key parameter for geothermal energy production. To enable a reliable temperature prediction for a geothermal reservoir, the significant heat transport mechanisms coupled to the fluid flow processes in the respected area have to be understood on a bigger scale. Thereby, different types of heat transport acting in a sedimentary basin should be considered (Kaiser et al. 2011): (1) conductive heat transport, where the efficiency of heat transfer depends on the thermal conductivity of the lithology of the respective layer and the occurring thermal gradient; (2) advective heat transport, where heat is transported passively by pressuredriven fluid flow; and (3) convective heat transport, where heat is transported by density-driven convective fluid flow. In a fluid influenced geothermal environment, such as the Molasse Basin, a combination of all three heat-driving mechanisms can be expected, as been already shown for other sedimentary basins, e.g. for the North German Basin by Kaiser et al. (2011) and Noack et al. (2013) and for the Po Plain by Pasquale et al. (2013).
In case of the European Molasse basin, such a systematic assessment how strong fluid flow-related heat transport, geological, and hydrogeological heterogeneities or the presence of permeable faults may alter the deep conductive thermal field on a basin-scale in the Molasse Basin is still lacking.
Hence, we have conducted a systematic parameter study by simulating the coupled fluid flow and heat transport in the Bavarian part of the Molasse Basin between the Danube in the north and the Alpine front in the south using a three-dimensional basin-scale numerical FE model. Thereby, we follow a modelling approach, whose advantages have already been demonstrated for different sedimentary basins in the past, e.g. for the North German Basin by Noack et al. (2013) and Scheck-Wenderoth et al. (2014). In this approach, we combine the calculations of the 3D lithospheric-scale conductive thermal field by Przybycin et al. (2015b) with coupled basin-scale simulations of the fluid and heat transport to predict the thermal field in the European Molasse Basin. For this, temperatures from depths dominated by conductive heat transport (e.g. the crust) have been extracted from the lithospheric-scale conductive model of Przybycin et al. (2015b). These temperatures have been prescribed as lower thermal boundary condition to the base of the vertically and horizontally higher resolved basin-scale model of this study which considers coupled transport of heat and fluid. Following this workflow, the obstacle of choosing a proper lower thermal boundary condition for the simulation of coupled fluid flow and heat transport on a basin-scale can be overcome. The resulting thermal and flow field may help to understand the relation between basin-scale fluid dynamics and temperature distribution. Moreover, these results from basin-scale simulations can provide physically reliable thermal and hydraulic boundary conditions for local high-resolution reservoir-scale models as required for areas of exploration interest. To investigate the possible influence of hydraulically conductive faults on the thermal field, we have integrated three large hydraulically conductive faults into the model representing the average strike directions of faults in the Molasse Basin in a second modelling step. However, these faults should only be considered as "testfaults" used to investigate the general qualitative impact of heat transport by fault-related fluid flow on the basin-wide thermal field. For a detailed implantation of realistic 3D faults into the model, more structural data describing the geometry of the faults would be needed.

The basin-scale structural and hydrogeological model
For the simulation of coupled fluid flow and heat transport, a basin-scale 3D structural model has been build based on the lithospheric-scale 3D structural model of Przybycin et al. (2015a). The latter integrates freely available depth and thickness information from wells and seismic lines and has been additionally constrained by 3D gravity modelling. For this study, the upper part of the lithosphere-scale model has been used for the present study. In addition, the vertical resolution of lithostratigraphical units has been increased for the basin-scale model compared to Przybycin et al. (2015a) by additionally distinguishing the Purbeck formation (Lower Cretaceous) from the Cretaceous layer according to StMWIT (2010) and the different Malm layers (Zeta-Alpha, Schulz and Thomas 2012). The base of the structural model of the present study was defined at a constant depth of −7500 m asl, a depth at and below which conduction is the dominant mechanism of heat transport. This basin-scale 3D structural model of the present study covers an area of 180 km in N-S direction and 340 km in E-W direction with a horizontal resolution of ~1 × ~1.7 km 2 corresponding to 185 × 201 grid points. In the vertical direction, the structural model contains 12 lithostratigraphic units from the topography downwards ( Fig. 2).
At the southern border of the model, the Folded Molasse Sediments occur in a thin band of strongly folded and steeply erected sediments along the mountain chain. They consist of conglomerate and sand as well as silt and clay with a moderate average permeability ( Fig. 2a; Table 1). They show apparent thicknesses of up to 7000 m directly in front of the Alps and represent therefore a possible pathway for fluid flow from the surface into larger depths.
The northwards following Foreland Molasse Sediments (Fig. 2b) consist of mostly unconsolidated conglomerate and sand with intercalated silt and clay layers originating from erosional processes of the Alps (Bousquet et al. 2012;Handy et al. 2010). These Molasse sediments are considered as highly permeable for fluid flow (Table 1) and show increasing thicknesses from 0 m in the north to up to 5000 m in the south at the Alpine front.
Below the Molasse Sediments the Cretaceous (Fig. 2c) consists mostly of claystone and limestone and is only preserved in a limited area between the city of Munich and the Landshut-Neuöttinger High with thicknesses of up to 850 m and in a small area east of the Landshut-Neuöttinger High with up to 150 m thickness. This unit has an average hydraulic conductivity even higher than the Molasse Sediments (Table 1).
According to deviating hydraulic properties (higher porosity, but lower hydraulic conductivity) and a different lithology (mainly limestone) compared to the (Upper and Middle) Cretaceous, the Purbeck formation (Lower Cretaceous, Fig. 2d) has been implemented as a separate unit in the hydrogeological model. The Purbeck formation is preserved in a small area between the city of Munich and the Landshut-Neuöttinger High with minor thicknesses between 10 and 40 m.
The Upper Jurassic Malm aquifer consists of mostly limestone and dolomite with small appearance of clay and can be separated into six cycles (overlying Zeta-underlying Alpha, Fig. 2e-j). While the lower two Malm layers (Alpha and Beta) can be regarded as impervious, the upper four Malm layers (Zeta-Gamma) represent a continuous aquifer with a horizontally varying hydraulic conductivity.
For this aquifer, the hydraulic conductivity has been prescribed according to the degree of dolomitization, karstification, and different hydraulic properties determined by Birner (2013;Fig. 3) based on well data: a high hydraulic conductivity was prescribed to the Malm aquifer in the central northern and the most eastern part, with 10 −4 m/s. Towards the south, the aquifer has been characterized by average hydraulic

Table 1 The thermal and hydraulical properties prescribed to the layers of the numerical model
The influence of each property on the predicted thermal field has been tested systematically References: (1) Schulz and Thomas (2012)   Below the Upper Jurassic Malm, the PreMalm layer ( Fig. 2k) summarizes all sediments of the Middle and Lower Jurassic as well as Triassic and consists of mostly claystone, sandstone, and marl. This layer is restricted to the model area north and northwest of the city of Munich and is characterized by an average thickness of 200-500 m with locally higher thickness values.
The lowermost layer of the model is represented by the crystalline crust, which is regarded as impermeable for fluid ( Fig. 2l; Table 1). In response to the constant depth of the model base with −7500 m asl, the thickness of this layer decreases from north to south with a maximum thickness of up to 7500 m directly at the northern border and less than 500 m at the southern border of the model.  (2013) and implemented into the 3D numerical model. Birner (2013) states that these four Malm layers can be regarded as a continuous aquifer on a basin scale and thus described with lithological units of different but homogenous hydraulic characteristics. The purple colour shows the area of the uplifted crust of the Landshut-Neuöttinger High (LNH) above which the Malm aquifer is eroded

The FEM model
The simulations of the coupled fluid flow and heat transport have been carried out using the commercial software package FELOW ® (Version 6.2, Diersch 2009), which is based on a finite-element method (FEM) and allows the consideration of different processes of heat transport in natural porous media, as conductive, advective, and convective heat transport. The coupled calculation of the flow and temperature field in a saturated porous media with FEFLOW ® is done solving three partial differential equations based on Darcy's law, on mass conservation as well as on energy conservation (e.g. Nield and Bejan 2006) which are described in the Appendix.
To transform the geological model into a numerical one for the software FEFLOW ® , all geological layers have been converted into continuous layers and discretized with a three-dimensional finite-element mesh with irregular triangular elements horizontally and prismatic elements vertically. While the horizontal resolution was defined as 600 m in average, the vertical resolution depends on the thickness of the lithostratigraphic units. To avoid a high aspect ratio (horizontal : vertical resolution), layers with large thicknesses (as the Folded Molasse Sediments, the Foreland Molasse Sediments, and the crystalline crust) have been subdivided into sublayers with maximum thicknesses of 1.5 km resulting all in all in 26 numerical layers. The resolution of the horizontal discretization was closely connected to the facies distribution of the Upper Jurassic Malm aquifer, as given by StMWIT (2010) in the Geothermal Atlas of Bavaria: In areas, where no or little fluid flow is expected due to a low permeability of the Malm aquifer, a courser discretization (element size up to 1000 m) has been chosen than in areas of high permeability. Hence, smaller mesh elements had been implemented in the central model area to assure numerical stability (element size down to ~400 m).
To assess the influence of fluid flow along permeable fault zones on the thermal field in the basin (Fig. 4), three faults have been integrated into the numerical model as discrete feature elements-finite elements with lower dimensionality inserted at element faces and edges (Diersch 2009). For those vertical 2D elements, a horizontal buffer zone of 40 m on both sides of the faults has been implemented to which a higher horizontal resolution of up to 20 m element size was set to assure stable calculation during simulations. For these domains (faults and buffer zones), fluid flow according to Darcy's law was assumed.
These testfaults are simplified representations of known fault zones in the Molasse Basin and are oriented in the main fault direction, respectively. Fault 1 has been integrated south of the city of Munich to represent the east-west-oriented fault system in the Malm aquifer (StMWIT 2010). It can be traced vertically from the middle of the Foreland Molasse Sediments down to Malm Epsilon. It represents a normal fault created by the flexural bending of the European plate in response to the load exerted by the Alps. Fault 2 has been integrated at the western side of the Landshut-Neuöttinger High and represents a possible flow pathway along the tectonic border between the Malm and the uplifted crust of the Landshut-Neuöttinger High. Fault 3 has been integrated at the southern border of the model representing the east-west-oriented fault system along the Alpine front in the Folded Molasse Sediments as the tectonical boundary between the Molasse Basin and the Alps. This fault can be traced from the topography to the top of the crystalline crust (Schmid et al. 2008).
All in all, the numerical model including the three faults consists of 2,432,560 prismatic mesh elements and 1,283,391 mesh nodes subdivided into 26 numerical layers and 27 slices and three vertical faults (Fig. 4).
To achieve quasi-steady-state conditions in the numerical simulations and with that predict the present-day thermal and flow field, a simulation time of 100,000 years was chosen with an automatic time step control with a limited maximum time step size of 50,000 days.

Parameterization
All layers have been characterized with thermal (heat capacity, thermal conductivity and radiogenic heat production), mechanical (porosity), as well as hydraulic properties (hydraulic conductivity) according to their dominant lithology (Table 1). Measured values have been favoured for each property. However, in case no measured values were available, average lithology-based values have been assigned. Every layer has been assigned with one uniform value for each property apart from the Upper Jurassic Malm aquifer. For the upper four layers of the Malm (Gamma-Zeta), a variable hydraulic conductivity was adopted following Birner (2013; Fig. 3) as described before. Accordingly, high hydraulic conductivities prevail in the northern and eastern-most model parts, and Like the layers, all faults have been characterized with hydraulic properties following Schulz and Thomas (2012) and Cherubini et al. (2014). Thereby, a hydraulic conductivity much higher than of any other layer of the model was assumed for the faults to open permeable fluid flow paths.
Detailed sensitivity analyses have been carried out to assess the level of influence of all assigned physical properties on the resulting thermal field.

Boundary conditions and initial conditions
To solve the partial differential equations, thermal as well as hydraulic upper and lower boundary conditions are required. The thermal boundary conditions at the base and the top of the model have been defined as fixed temperatures (Dirichlet). At the surface, a mean variable surface temperature measured over 30 years (1960( -1990( DWD 2013Fig. 5a) has been prescribed as upper thermal boundary condition. Thereby, the surface temperature ranges between −2 °C in the south-western corner and +12 °C in the northern most and eastern part of the model area.
A variable temperature distribution, which has been extracted from the 3D lithospheric-scale conductive thermal model of Przybycin et al. (2015b), has been implemented as lower thermal boundary condition at −7500 m asl. The model of Przybycin et al. (2015b) covers the Molasse Basin and the adjoining South German Scarpland and the European Alps and considers the influence of a varying depth of the thermal lithosphere-asthenosphere boundary (LAB), a heterogeneous internal structure of the crystalline crust as well as variations in thermal properties in response to geological structure. Accordingly, the interdependence between the Alps, the Tauern Body and the Molasse Basin with respect to the basin-wide thermal field is captured by the 3D lithospheric-scale conductive model. Hence, lower temperature values in the north (~150 °C) and higher temperature values in the south (~240 °C, Fig. 5b) are prescribed at the base as a thermal boundary condition for the coupled model.
As an upper boundary condition for fluid flow, a constant pressure head (0 Pa, Dirichlet) was assigned to the topography (Fig. 1). Furthermore, fixed hydraulic heads (Dirichlet) were assigned to the river Danube and to the lateral borders of the Malm aquifer in the east and west according to Frisch and Huber (2000) (Fig. 5a).
Furthermore, thermal and pressure initial conditions have been calculated with uncoupled fluid and heat transport simulations for steady-state conditions.

Results and interpretation
In the following chapter, the resulting thermal field predicted by the coupled fluid and heat transport simulations is presented with temperature maps at different depths and prominent thermal effects are highlighted. At first, the results of model 1 without considering faults are presented and compared to an alternative version of the model obtained during the sensitivity study (model 2). Afterwards, the results of the simulation considering hydraulically conductive faults (model 3) are presented.
The basin-scale thermal field (model 1) Figure 6 shows the temperature distribution predicted by Model 1 (without faults) at depths for which temperature maps from GeotIS (StMWIT 2010; Schulz et al. 2009) are available for validation. Compared to the maps of GeotIS (StMWIT 2010), the maps produced in this study show temperature variations changing at a much smaller wavelength at all depths. This effect is related to the assumption of a homogeneous grid in GeotIS, whereas this study considered heterogeneities related to geological structures. In addition, the grid resolution is far higher in this study compared to the GeotIS interpolation grid of measured data (Schulz et al. 2009).
At a depth of −500 m asl (Fig. 6a), temperatures between 8 and 35 °C are predicted for most parts of the model area with the coolest values in the area of the Folded Molasse Sediments. In addition, some positive thermal anomalies of smaller lateral extent with temperatures of up 63 °C are present at the southern border of the model area, in the The lower thermal boundary condition shows the temperature distribution extracted from the 3D lithospheric-scale conductive thermal model of Przybycin et al. (2015b). The upper figure additionally shows the prescribed lateral hydraulic boundary conditions (modified after Frisch and Huber 2000). 20 locations of the geothermal energy production sites are shown for which observed temperatures are available. Our predicted temperatures have been compared to these observations. The residual temperature (calculated-measured) is given. The temperature difference between calculated and measured values is additionally given for the coupled and conductive case in Table 2 western part of the model, at the eastern border of Bavaria as well as in the area of the city of Munich. Of those, the latter one shows the highest temperatures at this depth. A comparable distribution of temperatures and thermal anomalies is shown in the temperature map of GeotIS (StMWIT 2010).
At a depth of −1500 m asl (Fig. 6b), modelled temperatures range between 10 and 90 °C, though most parts of the model area are not warmer than 50 °C. The lowest temperatures are bound to the negative thermal anomalies at the southern border of the model area (~10 °C) and to the negative thermal anomaly predicted in the south-east of Munich (up to 15 °C). The highest temperature values are limited to the positive thermal anomalies, which, compared to the −500 m asl depth level, increased in size and temperatures. At −1500 m asl, three positive thermal anomalies are visible in the western part of the model, one big and two smaller anomalies, with the highest temperatures of up to 90 °C. The positive thermal anomaly in the area of Munich temperatures of up to 75 °C is predicted by the model, which again is comparable to the temperature distribution shown by GeotIS (StMWIT 2010). Furthermore, some smaller positive thermal anomalies are predicted in between the negative thermal anomalies at the southern At −2500 m asl (Fig. 6c), the coupled simulations predict temperatures between 15 and 140 °C. At this depth, the northern model area shows mean temperature values of ~75 °C, while the southern model area shows the lowest temperature values of up to 15 °C within the negative thermal anomalies of the Foreland Molasse Sediments. Compared to shallower levels, these negative thermal anomalies have increased in size as has the negative thermal anomaly south-east of Munich, which shows a temperature of ~25 °C. Likewise, the positive thermal anomalies are larger in size than at shallower levels. While at the depth of −1500 m asl, the three positive thermal anomalies in the west appear as distinct features, and they have merged into one bigger anomaly at the depth of −2500 m asl and show temperatures of up to 140 °C. The small positive thermal anomaly at the eastern border of Bavaria shows temperatures of ~100 °C. The positive thermal anomaly around the city of Munich temperatures between 80 and 100 °C is predicted, which is colder than the temperatures predicted by GeotIS (StMWIT 2010) by ~10 K. However, the trend of the measured temperature distribution could be reproduced even for this depth.
Deeper, at −3500 m asl (Fig. 6d), temperatures between 20 and 160 °C are predicted by the model, with average temperatures of ~90 °C apart from the thermal anomalies. At this depth, the coldest temperatures are limited to the negative thermal anomalies in the Folded Molasse Sediments at the southern border of the model area. The negative thermal anomaly to the south-east of Munich is characterized by temperatures of ~60 °C. The positive thermal anomaly in the west of the model area increased in size compared to −2500 m asl and displays the highest temperature values with up to 160 °C. For the positive thermal anomaly in the area of Munich, the simulations predict a temperature of ~120 °C, and the anomaly is larger in size compared to shallower depths. The temperature of the positive thermal anomaly at the eastern border of Bavaria increased as well to 130 °C. Compared to the measured temperatures shown in the GeotIS maps (StMWIT 2010), the simulations of this study predict temperatures 10 K lower in average for this depth. Nonetheless, the temperature trend of GeotIS is reproduced.
At a depth of −4500 m asl (Fig. 6e), average temperatures of ~100 °C are predicted. The negative thermal anomalies along the southern border of the model area decrease in size and show low temperatures of up to 20 °C. For the negative thermal anomaly in the south-east of Munich, the model predicts temperatures of ~90 °C. To the west, the positive thermal anomaly reaches temperatures to up to 170 °C. The temperature of the warm thermal anomaly around Munich rises to 135 °C. For this depth, no map from GeotIS (StMWIT 2010) is available anymore for comparison.
The last temperature map (Fig. 6f ) shows the temperature distribution at a depth of −5500 m asl. At this depth, an average temperature of 140-150 °C is predicted by the coupled simulations. The negative thermal anomaly at the southern border of the model area decreased strongly in size and predicted temperatures are about 50 °C. The highest temperatures are predicted for the positive thermal anomaly in the most western corner of the model area with up to 190 °C. While the positive thermal anomaly around Munich shows temperatures of up to 170 °C, the negative thermal anomaly south-east of Munich is less pronounced than at shallower depths with values of ~130-140 °C.
All predicted temperatures have been additionally compared to measured bottom hole or extraction temperatures from different depths of 20 measurements (e.g. Agemar et al. 2014). For 12 sites, the modelled temperatures lie in a range of ±10 K of the measured temperatures ( Fig. 5; Table 2). For two sites west of the Landshut-Neuöttinger High, the model predicts temperatures more than 10 K lower and for two sites temperatures higher than the measured ones. For four sites in the area of Munich temperatures, more than 10 K higher, and for one site temperatures, more than 10 K lower than the measured ones are modelled.

Interpretation: influence of deep fluid flow on the thermal field
To compare the temperatures predicted by this study to a purely conductive approach as followed by Przybycin et al. (2015b), the conductive temperature distribution at depth of −3500 m asl is shown in Fig. 7. This temperature map shows increasing temperatures from the north to the south with ~120 °C around Munich. South of Munich three prominent positive thermal anomalies with temperatures of ~140 °C are predicted in the area of the Folded Molasse Sediments by the conductive model with lower temperatures (~130 °C) in between. Compared to this conductive model, the coupled model 1 predicts a more diverse temperature distribution with a variation of much smaller wavelength for the depth of −3500 m asl (Fig. 7a). Where a continuous increase of temperatures from north to south is calculated with the conductive model, an irregular temperature distribution is predicted by the coupled model. Moreover, the positive thermal anomalies predicted by the conductive model in the south of Munich are shifted to the north in the coupled model creating the positive thermal anomalies in the west, around Munich and at the eastern border of Bavaria. In turn, at the position of the positive thermal anomalies in the conductive model at the southern model boundary, negative ones are predicted by the coupled model. This shift of the positive thermal anomalies in the coupled model compared to the conductive case is caused by the downward directed fluid flow through the Folded Molasse Sediments in the coupled case. This fluid inflow from the surface displaces hotter fluid at depth and, thus, relocates the heat coming from deeper parts of the basin by advective heat transport from the south to the north. This relocation can also be seen in profiles through the model area depicting the temperature distribution in the basin (Fig. 8).
The first profile (Fig. 8a) runs east-west through the model and directly cuts the positive thermal anomaly around Munich and the negative thermal anomaly south-east of Munich. In this profile, areas where cold temperatures reach larger depths represent domains of downward flowing fluid (recharge). In contrast, areas of higher temperatures at shallower depth represent domains of upwards flowing fluids (discharge). For the negative thermal anomaly south-east of Munich, this profile illustrates that the fluid turnover creating the negative thermal anomaly is limited to a local area and a comparable shallow depth (~2000 m asl). At this depth, the Purbeck formation hydraulically decouples the Foreland Molasse Sediments from the Upper Jurassic (Fig. 2c, d) due to its low hydraulic conductivity and, thus, prevents a deeper penetration of cold water from the surface. Nevertheless, the overall colder upper 2 km also affect the deeper layers via conductive heat transfer. Profile b (Fig. 8b) shows the temperature distribution in northsouth direction through the negative thermal anomaly south-east of Munich. In addition,  from this perspective, the prevented inflow of cold fluid from the surface to the larger depth is evident. In addition, upward directed flow of warmer fluid from larger depths is prevented as well by the Purbeck layer. This upward directed flow of warm fluid diverts to the west and east and contributes the positive thermal anomalies around Munich and at the western border of Bavaria, as illustrated in profile c (Fig. 8c). In particular, stream tracers show that the positive thermal anomalies are closely related to fluid entering the system vertically through the Folded Molasse Sediments due to a high topographic gradient. This fluid flows to deeper levels until it reaches the impervious crystalline crust. At that depth, the fluid flow direction changes to northward directed fluid flow. This in turn causes a general northward flow pattern through the Malm aquifer and the Molasse Sediments between Munich and the Landshut-Neuöttinger High parallel to the tilted Fig. 7 Temperature maps at a depth of −3500 m asl showing the temperature distribution calculated for the case of coupled transport (above) and the conductive case of heat transport (below, after Przybycin et al. 2015b) top of the crystalline crust and around the Purbeck formation. Thereby, the heat from deeper levels is relocated compared to the conductive case by northward directed flow under advective conditions, leading to the positive thermal anomalies around Munich and at the eastern border of Bavaria and the strong negative thermal anomalies in the Folded Molasse Sediments. In other words: while the negative thermal anomaly in the east of Munich is related to locally restricted downward flow of cold fluids, the positive thermal anomaly around Munich is caused by regional flow of fluid that initially entered the model through the Folded Molasse Sediments. This indicates that the Folded Molasse Sediments have a significant influence on the regional thermal field of the Molasse Basin and the underlying Upper Jurassic Malm aquifer even at larger depths.
Even though the Folded Molasse Sediments are characterized by the lowest hydraulic conductivity of all sediments, their large thickness of up to 7000 m (Fig. 2a) in combination with the high hydraulic gradient caused by the topography in front of the Alps results in an overall flow pathway for cold fluid to larger depths. This combination of factors (hydraulic conductivity, thickness, and high hydraulic gradient) leads to the very low temperatures at large depths. This result is, however, difficult to validate as no temperature measurements are available for comparison for this depths. To assess, how strong the influence of the hydraulic conductivity of the Folded Molasse Sediments on the thermal field in the whole basin is, we decreased the hydraulic conductivity of the Folded Molasse Sediment even further by two orders of magnitude compared to the original model to 2 × 10 −10 m/s. A temperature map at a depth of −3500 m asl (Fig. 9) for this second model is shown. At this depth, a much lower average temperature of ~60 °C is obtained in the second model compared to the model 1 (Fig. 6). While the positive temperature anomalies at the southern model border, which were small in model 1, increase in size and value in the second model, the strong negative thermal anomalies in turn disappear nearly completely. The Folded Molasse Sediments, which were cold in average in model 1, are warm in the second model. The positive thermal anomalies in the west around Munich and at the eastern border of Bavaria are much colder (~90 °C) than in model 1. Likewise, the negative thermal anomaly south-east of Munich also shows lower temperatures in the second model (~30 °C) than in model 1.
Obviously, a lower hydraulic conductivity of the Folded Molasse Sediments leads to warmer temperatures in the Folded Molasse Sediments, but colder temperatures in the rest of the model. Though the general temperature trend remains the same and consistent with GeotIS (StMWIT 2010), the absolute values are better reproduced by model 1 (Fig. 6).

The influence of permeable faults on the thermal field (model 3)
Even though the measured thermal field of the Molasse Basin could be reasonably well reproduced already by coupled fluid flow and heat transport simulations without considering fault-related fluid flow, a possible influence of permeable faults on the temperature distribution cannot be excluded, at least on a local scale. Therefore, three permeable testfaults have been implemented into the model as described earlier to simulate the potential influence of large, permeable faults on the thermal field (model 3). In Fig. 10, the resulting temperature distributions of these simulations are presented with a temperature map at a depth of −3500 m asl. Compared to the temperature distribution at the corresponding depth of model 1 (Fig. 6), only minor changes of the temperature distribution are observable for model 3. Each permeable fault has only a local influence on the thermal field in its direct surroundings (~1 km) by fault-related fluid flow which cools down or heats up the fault-near areas, respectively. These results indicate that permeable faults have no significant influence on the regional flow field and are not particularly necessary to reproduce the basin-wide trend of the temperature distribution in the Molasse Basin area. However, the temperature misfit between measured and predicted values in model 1 of more than 10 K at the three locations in the area of the Landshut-Neuöttinger High could be decreased slightly when considering the permeable faults.

Discussion
As was already shown by Przybycin et al. (2015b), the long wavelength thermal field of the Molasse Basin is dominated by conductive heat transport which in turn is mostly influenced by the structural configuration of the crust and the lithosphere beneath the basin and lateral heterogeneities in the thermal conductivity. Thereby, the depth of the thermal lithosphere-asthenosphere boundary (LAB = 1300 °C isotherm) has a distinct effect on the long wavelength thermal field by creating the principal thermal gradient. A second effect is caused by the upper crystalline crust, which is characterized by a much higher radiogenic heat production than the sediments or the deeper crust and the mantle. In addition, sensitivity analyses of thermal properties have shown that the structure and the high thermal conductivity of the crystalline core of the Tauern Window cause a chimney effect between the thermally more conductive Tauern Body and the thermally less conductive Alpine Body. Moreover, a blanketing effect is caused by the thermally less conductive Molasse Sediments covering the thermally more conductive crystalline crust. These effects may partly explain the positive thermal anomalies around Munich and at the eastern border of Bavaria with a negative thermal anomaly in between. However, the thermal anomalies predicted by the conductive model of Przybycin et al. (2015b) have been too warm compared to measured data, which indicates that other than conductive effects may be relevant.
Using the lithospheric-scale 3D thermal model of Przybycin et al. (2015b) as a base for our basin-scale simulations, we took into account the long wavelength conductive thermal effects by prescribing temperatures extracted from the lithospheric-scale model as lower thermal boundary condition to the coupled model. Following this approach, the fit of the predicted values for the thermal anomalies with observed temperatures could be improved (Figs. 6, 7) compared to the purely conductive approach (Przybycin et al. 2015b). Moreover, the coupled fluid and heat transport simulations were not only able to reproduce the pattern of the temperature distribution, but also the range of the measured temperatures values in the Molasse Basin area.
The remaining average temperature misfit between the measured data (GeotIS, StM-WIT 2010) and the model predictions of ±10 K lies in the range of the standard deviation of bottom hole temperatures and temperature logs (Hermanrud et al. 1990;Förster 2001;Noack et al. 2010;Agemar et al. 2012Agemar et al. , 2014. Moreover, Agemar et al. (2012) recommend a careful treatment of the deeper temperature estimates of GeotIS (StMWIT 2010), since the data base decreased with depth and hence possible local thermal effects may not be considered in the temperature interpolation. In addition, the latter assumes a homogeneous distribution of thermal conductivities in the subsurface which does not account for lithological heterogeneities and related heat refraction. Under this consideration, the reproduction of the observed temperatures in this study can be considered as satisfying for the chosen model size. Albeit, we admit that the limited vertical and horizontal resolution, the assignment of uniform physical properties to most of the layers and the prescribed first kind boundary conditions may be sources of error in our results. Therefore, the sensitivity of the modelling results on the prescribed thermal properties and the boundary conditions has been tested. The results of this analysis show that even if the predicted absolute temperatures are strongly depending on the hydraulic conductivity prescribed to the Molasse Sediments, the temperature trend remains the same for different variations of the model. A comparable though weaker effect could be observed for different boundary conditions. Temperature variation caused by different varying thermal properties is much smaller.
An additional consideration of permeable faults (Fig. 10) had only a minor influence on the regional distribution of thermal anomalies in the Molasse Basin. Ascending warmer fluid in the area of the positive thermal anomalies and descending colder fluid in the area of the negative thermal anomalies driven by density gradients within the faults caused insignificant (1-5 K) local cooling which did not lead to changes in the general temperature trend. Such a limited spatial influence of permeable fault zones on the regional thermal field has been already described by Cherubini et al. (2014). However, an implementation of permeable fault zones into reservoir-scale models might still be of high importance for the local temperature. Such a local significance of permeable faults was implied by the slightly decreased temperature misfit between predicted temperatures and measured values in the area of the Landshut-Neuöttinger High after the implementation of permeable faults into model 3.
Summarizing, our results show that the basin-scale thermal field is caused by a combination of conductive and advective heat transport. In contrast to the work of Pasquale et al. (2013) for the karstified carbonate rocks in the Po Basin, we found no indication for large-scale free thermal convection in the carbonates of the Malm aquifer. This may be explained by the much smaller thickness of the upper four layers of the Upper Jurassic Malm aquifer in Germany of maximum 400 m of permeable thickness compared to the thick carbonate platform of the Po Basin with more than 4000 m. Free thermal convection evolves only if the reservoir is sufficiently thick and lateral gradients in hydraulic head are insignificant (Kaiser et al. 2011) as otherwise pressure-driven advection suppresses density-driven fluid flow. The scenario of the Molasse Basin does not fulfil these criteria for free convection. The permeable layers of the Upper Jurassic Malm (Zeta-Gamma) are not thick enough and the hydraulic head too uniform over the basin area to enable thermal convection. Rühaak (2009Rühaak ( , 2015 and Rühaak et al. (2010) calculated the thermal field of the western Molasse Basin using a local-scale conductive model. Since they were not able to reproduce the observed thermal anomalies in the western Molasse Basin with such a conductive approach, they presumed fluid flow along E-W striking faults to cause temperature differences of more than 10 K in the basin.
We regard the results of this study as consistent with their work, saying that fluid flow in general is needed to reproduce the detected thermal anomalies in the Molasse Basin. However, in our case, the pathways for fluid flow are mainly created by the matrix permeability of the Folded and Foreland Molasse Sediments and the complex permeability variations in the Upper Jurassic Malm aquifer (Birner 2013).
More precisely, sensitivity analyses have shown that the average temperature values in the Molasse Basin are strongly influenced by the hydraulic conductivity of the Folded and Foreland Molasse Sediments (Fig. 9). Hydraulically less conductive Folded Molasse Sediments would lead to higher deep temperatures at the Alpine front, but lower temperatures in the rest of the basin. Hydraulically more conductive Folded Molasse Sediments would lead more cold water to larger depths causing an overall cooling of the system. A similar effect would be caused by hydraulically more conductive Foreland Molasse Sediments, whereas a lower hydraulic conductivity of the Foreland Molasse Sediments would lead to an overall warming of the system. Moreover, our results show that the hydraulic conductivity of the Malm aquifer has only an influence on the flow direction in the Malm aquifer, but not on the overlying layers. Accordingly, the positive thermal anomalies around Munich and at the eastern border of Bavaria appear to be caused by northward directed regional fluid flow (Fig. 8) from deeper parts in the south of the model area. Furthermore, our results indicate locally restricted fluid flow related to the structural and hydraulic configuration of the Purbeck formation causing the negative thermal anomaly in the east of Munich. This combination of basin-wide regional and locally restricted fluid flow as explanation for the development of adjoining positive and negative thermal anomalies in the Molasse Basin has so far not been suggested but can have a significant influence on local and reservoirscale studies. In particular, it may question some results derived from a local-scale 3D structural and thermal model of the city of Munich (Schulz and Thomas 2012) that does not consider different hydraulic boundary conditions at the lateral boundaries.
An investigation of thermal anomalies with local or reservoir-scale models may not capture the hydrothermal dynamics with sufficient detail, since anomalies in the thermal field may be caused by regional fluid flow. In such cases, an adoptable strategy would be to extract pressures and temperatures from basin-scale models as thermal and hydraulic boundary conditions for local and reservoir-scale models to incorporate the appropriate consideration of regional effects.

Conclusions
By following a multi-scale data-based 3D-modelling approach, the thermal field of the German Molasse Basin was investigated with simulations of coupled fluid flow and heat transport. Thereby, existing knowledge about the long wavelength and deep conductive thermal field of the basin was taken into account by prescribing temperatures, extracted from lithospheric-scale 3D conductive calculations, as the lower thermal boundary condition for the basin-scale model. The resulting thermal field reproduces the measured temperature distribution of the basin, including the pronounced negative and positive thermal anomalies satisfactorily for a basin-scale approach. To reduce the remaining misfit between observed and modelled temperatures as well as remaining uncertainties, a higher vertical and horizontal structural resolution and better knowledge of the distribution of physical property are needed.
Considering the results of previous conductive studies, our results confirm that the thermal field of the German Molasse Basin is controlled by conductive heat transport in the first place, especially in the deeper parts, but strongly influenced by advective heat transport within the sediment fill bound to basin-wide fluid flow. The pronounced positive and negative thermal anomalies in the basin are partly triggered by conductive heat transport, but reinforced by a combination of regional and local fluid flow, respectively, and mostly depending on the geological structure and the hydraulic conductivity of the sediments. Faults appear to have only a subordinate, local influence on the thermal field. The results of this basin-wide study of coupled fluid flow and heat transport contribute to a better understanding of the origin of the thermal anomalies in the basin. In addition, the study helps to reduce the exploration risk of geothermal energy project by providing higher resolved predictions of the deep pressure and temperature conditions. Moreover, the model can be used to derive reliable pressure and temperature boundary conditions for high-resolution reservoir-scale models for areas of exploration interest.