The Franconian Basin thermal anomaly: testing its origin by conceptual 2-D models of deep-seated heat sources covered by low thermal conductivity sediments

Geological framework

The temperature anomaly of interest is situated within the Franconian Basin [ ²⁶ ] in SE Germany (Fig.  ¹ ). The basin is underlain by Saxothuringian Variscan basement [ ²⁷ ] and comprises Permian, Mesozoic and Cenozoic sedimentary rocks, with sandstones, siltstones, claystones, carbonates and few evaporites forming the prevailing rock types. The sedimentary cover reaches a thickness of > 1600 m in the Staffelstein Graben [ ²⁸ ] NW of well Mürsbach B6 (MueB6) (Fig.  ² b). Upper Triassic sediments of the Keuper unit dominate the surface outcrops in the study area Fig.  ² a).

The basement geology of the Franconian Basin was formed during the amalgamation of Western Europe during the Variscan Orogeny accompanied by the intrusion of numerous voluminous granites [ ²⁷ , ²⁹ , ³⁰ ]. During the Cenozoic, the study area was affected by the development of the European Cenozoic Rift System, which involves the emplacement of Eocene to Pleistocene alkaline and minor tholeiitic intraplate basalt dikes [ ³¹ , ³² ] and adjacent to the study area, the formation of NW and NE striking rift basins.

NW–SE and NE–SW striking faults dominate the regional structural framework, some of which reflecting deep-seated basement faults that experienced multiple reactivation [ ³³ , ³⁴ ]. Permeable and deep-reaching fault zones in the study area are indicated by locally enhanced radon emanations [ ³⁵ ]. Such hydraulic basement connections can be assumed especially for the areas of intersection of NW–SE to NNW–SSE striking faults [ ³⁶ ], which are orientated (sub)parallel to the recent stress field [ ³⁷ ].

Temperature data

All available well data, including temperature readings of deep wells in the case study area are summarised in supplementary Appendix A.1 and A.2. The highest temperature of 66.8 °C at 1204.2 m depth in the centre of the temperature anomaly was recorded by the T -Log1_MueB6 (name given for this study) in the deep well MueB6 [ ³ ] (Fig.  ² a). The temperature log is interpreted to provide undisturbed formation temperatures, as the shut-in time of 18 months (supplementary Appendix A.1) is considered as a sufficient time period to reach complete thermal equilibration after drilling and mud circulation were completed [ ³⁸ , ³⁹ ].

Thermo-hydraulic input data

The presented models were calibrated using geologically and statistically valid thermo-hydraulic data of the Franconian Basin in SE Germany. 509 core samples were collected from the deep wells (> 400 m) Obernsees 1 (OS), Mürsbach B1, Mürsbach B3, Mürsbach B4 (Fig.  ² a), and also from several shallow wells (Fig.  ² a, 163 cores), reaching down to 20–40 m depth. To include stratigraphic units that could not be sampled from cores, additional outcrop samples ( n  = 53) were collected (Fig.  ² a). Because of restricted sample availability of particular formations (e.g., clay- and marlstones of the Upper and Middle Triassic), published thermal–hydraulic data had to be integrated (Table  ¹ ).

Measurements were performed for the main rock types of almost all stratigraphic units of the Permo-Triassic sedimentary cover as well as for selected Variscan basement rocks (Table  ¹ ). The input parameters of the models presented are thermal conductivity ( λ ), thermal diffusivity ( a ), bulk density ( ρ _bulk ) porosity ( Φ ) and specific heat capacity ( C _p ). For convective models, the permeability ( k ) was added as a variable. C _p was calculated from the measured λ , ρ _bulk , and a .

Thermo-hydraulic data of basement rocks could only be measured on core material of well OS (Fig.  ² a). No additional basement core material was available from any other deep wells within the study area due to its restricted exploration status.

Software

For numerical modelling, the heat transfer module and subsurface flow module of the FEM software package COMSOL Multiphysics ^® 5.3a (known as FEMLAB prior to 2005; [ ⁴² ]; Comsol Multiphysics: https://www.comsol.de , accessed 08/02/2019) were used. The software allows modelling of conductive heat transfer in porous media as well as convective Darcy flow of fluids within permeable rocks and fractures.

As the fluid properties exert a major control on the heat transport within a system [ ⁴³ ], the thermo-physical parameters of the fluid ( λ , density, viscosity) had to be adjusted to ambient temperature conditions. The governing equations used in this study are provided in supplementary Appendix B. Specific software settings applied during the numerical simulations are provided in Table  ² .

Tested subsurface configurations

We have tested four conceptual 2-D subsurface configurations, which differ in the type of heat source, the presence or absence of a permeable fault and the type of heat transport (purely conductive or combined with an additional convective transport component). The two different heat sources are:

• Various values for the average heat flow from below (generated by radioactive decay in the mantle and the crust as well as from cooling of the residual heat from the time of Earth’s origin). This is referred as background heat flow (BHT) in this study.

• Variable volumes of a heat-producing granite, representing an additional radioactive heat source (besides the BHT) in the crust.

For applying a pure conductive thermal regime, the two tested subsurface configurations are:

• The basic model (B-model), in which various values for the background heat flow (BHF) of 0.070–0.125 W m ⁻² are applied as the only heat source at the lower boundary of the model (Fig.  ³ a, Table  ⁴ ). A BHF value of 0.070 W m ⁻² generates a temperature–depth profile expectable from a normal continental geothermal gradient of 30 K km ⁻¹ in the B-model setup.

  Fig. 3

  

• The basic model with granite (BG-model), in which a heat-producing granite body provides an additional heat source to the BHF (Fig.  ³ a, Table  ⁴ ). In this model setup, the BHF is not varied and uniformly set to 0.070 W m ⁻² .

The setup of the B-model was calibrated on the stratigraphic data of the deep well MueB6, for which the equilibrated T -Log1_MueB6 is available and where the highest underground temperatures in the study area were measured (Fig.  ² a). The simulations presented do not provide a basin-scale study, but only depict the central region of the temperature anomaly in the Franconian Basin. The stratigraphy of the sedimentary cover in well MueB6 and measured thermo-hydraulic characteristics of the sedimentary rocks (Table  ² ) are valid for this part of the basin.

As the model area in numerical modelling is of finite extension, it is important that it is large enough not to influence the resulting temperature distribution. In our model setup, the total model depth z _model defines the vertical distance between the heat source at the lower boundary (BHF) and the surface of the model. In addition, z _model controls the thickness of host rocks around the pluton in the model setup comprising a heat-producing granite (Fig.  ³ b). Therefore, in a first step, sensitivity analyses of both parameters were performed, in order to determine the minimum model depth and width required not to influence the resulting temperature distribution (“ ^{Model extensions} ”).

The sedimentary cover in the B-model (Fig.  ³ a) is further subdivided into 30 layers of different thicknesses and rock properties, mirroring the main rock types and individual stratigraphic subunit thicknesses of well MueB6 (Table  ¹ ). The lateral thicknesses of the layers are kept constant, considering a simplified “layer cake” stratigraphy, particularly of the Triassic. In the software; each layer was assigned the specific thermo-hydraulic parameters according to Table  ¹ .

As drilling of the adjacent well Mürsbach B1 ceased at 1309.0 m within the Lower Permian (Rotliegend), the exact depth of the basement-cover boundary is not known at the location of well MueB6 and only estimated at c. 1350–1400 m from the nearby DEKORP-3B section [ ⁴⁴ ]. Well OS, located at the E ′ border area of the Franconian Basin, reached Variscan metasedimentary basement rocks at a depth of 1341.7 m [ ⁴⁵ ]. As the deep wells at the location Mürsbach are situated more towards the centre of the sedimentary basin, basement rocks are assumed at a greater depth compared to the location of well OS. In all model configurations, we uniformly set the top of the basement rocks at 1393.0 m assigning the Lower Permian a total thickness of 175 m.

In the B-model, the basement is subdivided into four sublayers (Fig.  ³ b; Table  ¹ ). The subsurface structure of the upper (Variscan) crust is assumed to be composed of 6000 m metapelites (Phyllite) and of 2500 m gneiss (Paragneiss 1) at the base. The upper crust is then underlain by about 7000 m of higher metamorphic grade Cadomian gneisses (Paragneiss 2) [ ⁴⁶ ]. The thermo-hydraulic input parameters for the upper Phyllite layer were calibrated on core measurements of the well OS. Due to the lack of appropriate data for the lower basement layers, the input parameters were taken from the literature (Table  ¹ ).

In order to test the effect of additional convective heat transport in the two basic model setups, we incorporated a permeable fault zone and also allowed convective heat transport in the pore space of the sedimentary cover. The permeability of the basement rock layers was kept at a low level, ensuring only conductive heat transport in this part of the model. The resulting convective model setups are the following:

• The basic model with a fault (BF-model) (Fig.  ³ a).

• The basic model with granite and a fault (BGF-model) (Fig.  ³ a, b).

As it is known that the permeability and the opening width of a fault have a significant impact on the hydrothermal field [ ⁴⁷ ], sensitivity analyses were performed for this two parameters. In addition, the influence of the maximum depth extent of the implemented fault zone was tested.

Initial and boundary conditions

A constant surface temperature of 9.0 °C was applied as upper boundary condition (Fig.  ³ b), which represents the mean annual air temperature in the centre of the case study area (CDC-FTP-Server Deutscher Wetterdienst: ftp://ftpcdc.dwd.de/pub/CDC/ , accessed 14/01/2018). As lower boundary, the BHF was set. So-called infinite element domains were chosen as lateral boundaries (Fig.  ³ b). According to COMSOL Multiphysics ^® , they give the model infinite lateral extension and prevent that the model results are influenced by the lateral model boundaries. The boundary conditions on the outside of the infinite element (thermal insulation) are effectively applied at a very large distance from any region of interest.

The lower and upper model boundaries are set impermeable with respect to flow. Absolute pressure in the model domains is set as pressure gradient, calculated from the specific depth, the density value of the model layer (Table  ¹ ) and gravitational acceleration.

We only applied a T -correction of thermal conductivity ( λ ) data according to Somerton (1992) [ ⁴⁸ ] (Eq.  ¹ ), as the temperature effect dominates the pressure influence on λ values [ ⁴⁸ , ⁴⁹ – ⁵⁰ ]:

Heat production and geometry of the granite

The heat-production rate of granites can vary over a broad scale. Heat-producing granites can be classified as low-, moderate- and high-heat-producing granites [ ⁵¹ ]. In this study we assume a moderate heat production ( A _granite ) of 6 μW m ⁻³ based on the data of Scharfenberg et al. [ ⁵² ] for the Fichtelgebirge and Oberpfalz granites, considered as surface analogues for the deeply buried granites in the study area. There, gravimetric modelling indicates the presence of isolated granite bodies within the buried Saxothuringian basement of the Franconian Basin, located at 2–18 km depth [ ⁵³ ].

Besides the A _granite value, also the geometry and the cross-sectional area of granite bodies exert a control on the subsurface temperatures. Different pluton aspect ratios (length/thickness) are proposed in the literature. Dome and a sheet geometries are the two endmembers discussed [ ⁵⁴ , ⁵⁵ ]. We run our models with these two geometries, resulting in variable granite cross-sectional areas:

• Tabular, sheet-like plutons have a power law geometry according to Eq. ( ² ) [ ⁵⁴ , ⁵⁵ ], where T is the thickness of the pluton and L is its lateral extension. We tested two granite geometries in our 2-D models: G 1 (11.55 km ² ) and G 2 (35.0 km ² ) (Table  ³ ).

  T=0.29L0.6.

  Table 3

  Name	Lateral extension (m)	Thickness (m)	Cross-sectional area (km2)	Reference geometry	Agranite (mW m−3)
  Granite 1 (G1)	10,000	1155	11.55	a	6
  Granite 2 (G2)	20,000	1750	35.0	a	6
  Granite 3 (G3)	21,000	7000	147.0	b	6
  Granite 4 (G4)	30,000	10,000	300.0	b	6

• For the dome-shaped pluton geometry, we chose an aspect ratio of 3:1 (length/thickness). The tested geometries G 3 and G 4 have significantly larger cross-sectional areas (147.0 and 300.0 km ² ) than the tabular geometries G 1 and G 2, resulting in higher heat-producing granite volumes in the models.

To test the impact of the absolute pluton depth, the granite top was set to 2000, 4000 and 6000 m depth, respectively. For all model setups, we applied a time-dependent study with forward modelling of 300 Ma, an average of known granite emplacement ages in the Saxothuringian basement [ ⁵⁶ , ⁵⁷ – ⁵⁸ ]. Heat supply by the granite and the BHT starts at 0 Ma, generating a transient temperature field at the beginning of the simulation. However, at constant heat supply by the heat sources, steady-state conditions are reached after 36 Ma at the latest. Although the tectonic configuration of the modelled part of the basin was not constant over the modelled period, the applied simulation time ensures that steady-state is always reached for all tested underground configurations. Questions regarding the history of the basin and tectonic aspects could not be implemented in the software and would go beyond the scope of the conceptual models presented.

Fault geometries and permeability

The magnitude of fluid flow within a permeable fault is directly related to its hydraulic and geometric aperture [ ⁵⁹ , ⁶⁰ – ⁶¹ ], which in nature can be very complex. However, for numerical modelling a planar fracture model with two smooth and parallel planes is commonly applied for simplification [ ¹⁵ ], representing a reasonable first-order approximation for fractures in crystalline rock environments [ ⁶² ].

For the convective scenarios, a permeable fault was implemented in the basic conductive model setups (B-model + BG-model), resulting in model setups BF and BGF (Fig.  ³ ). In the heat transfer module of COMSOL, this can be done by inserting a line element, which is assigned the property “fracture” with a certain opening width (fracture thickness). The subsurface flow module allows setting gravitational acceleration ( g ) for the entire model domain and (darcy) fracture flow for the line element of the fault. The implemented fault has a dip angle of 70° and crosscuts the investigated model profile P at 970 m depth (Fig.  ³ b). The fault geometry is based on fault interpretations from a 2-D seismic survey in the Mürsbach area [ ¹ , ⁶³ ] and a reported fault zone at 970–977 m depth in the cored section of well MueB6 [ ⁶⁴ ].

The amount of fluid flow along the implemented fault is controlled by the fault permeability ( k _frac ) and the opening width of the fault ( d _frac ). Moreover, the depth of the fault ( z _frac ) was also assumed as controlling factor for the thermal impact of the fault. Therefore, sensitivity analyses have been carried out for this three fault parameters.

We tested three z _frac values of 1500, 2000 and 3000 m, assuming that all fault geometries reach the basement rocks. Since there is no information on possible fault permeabilities in the Franconian Basin, values for potential hydraulically active faults had to be estimated. The hydraulic conductivity of the fault was set much higher than the permeability of the other model layers, in order to generate fluid flow along the implemented fault [ ⁶⁵ , ⁶⁶ ]. Three different values for k _frac have been tested: 1.0 × 10 ⁻¹³  m ² ( k _frac1 ), 5.0 × 10 ⁻¹³  m ² ( k _frac2 ) and 1.0 x10 ⁻¹²  m ² ( k _frac3 ). The tested fault permeabilities represent values which are 10, 50 and 100 times higher than the highest matrix permeability in the sedimentary cover (1.0 × 10 ⁻¹⁴  m ² ; Table  ¹ ). Water was chosen as fluid within the fault zone and within the pore space of the model layers. Further aspects, such as variable flow regimes or salinities of the fluid, were not tested, as the discussion of these parameters would go beyond the scope of this study.

Sensitivity analysis

Model extensions

For the sensitivity analysis of the model extension, we chose the BGF-model with the largest granite geometry G 4 at greatest pluton depth of 6000 m. The incorporated fault has a maximum depth of 2000 m and the lowest fault permeability k _frac1 (Table  ⁴ ). Implications from this scenario are also valid for smaller granite geometries at shallower depths as well as for the conductive model setups. Figure  ⁴ a shows that z _model values of 20, 30 and 40 km result in almost similar temperature curves at the location of the analysed profile P (Fig.  ³ b). By increasing the z _model values, the modelled temperatures slightly decrease, but only to a negligible extent. For computing time and capacity reasons, we, therefore, uniformly set z _model to 20 km.

Table 4

Thermal regime	Model setup (Fig. 3)	BHF (W m−2)	Granite cross-sectional area (Table 3)	kfrac (m2)	Reproduction of T-Log1_MueB6a	Cross-reference
Conductive	B-model	0.070	–	–	No	Figure 6
		0.115	–	–	Yes (> 1000 m)	Figure 14
	BG-model	0.070	G1, G2, G3	–	No	Figure 10
			G4	–	Yes (> 1000 m)	Figure 14
Conductive + convective	BF-model	0.070–0.115 0.070–0.115 0.070–0.105 0.115	–	1 × 10−12 5 × 10−13 1 × 10−13 1 × 10−13	No No No Yes (< 800 m)	Figure 8 – Figure 8 Figure 14
	BGF-model	0.070 0.070	G2, G3 G4	1 × 10−13 1 × 10−13	No Yes (< 800 m)	Figure 12 Figure 14
	BGF-model modified2)	0.070	G4	1 × 10−13	Yes	Figure 14

Fig. 4

In order to eliminate the influence of the lateral model extension ( d _model ), infinite elements were defined as lateral model boundaries (Fig.  ³ b), which prevent boundary effects and ensure that the model results are not influenced by the lateral model boundaries (“ ^{Initial and boundary conditions} ”). The sensitivity analysis of d _model reveals that tested d _model values of 35, 40 and 50 km result in exactly the same temperature curves (not illustrated). The application of infinite elements eliminates the influence of d _model and, therefore, can be regarded as a very appropriate lateral boundary condition for thermo-hydraulic modelling. For all tested scenarios, d _model was set to 40 km, as the largest granite geometry already has a lateral extension of 30 km and isotherms in the surrounding host rocks can be analysed.

Opening width and permeability of the fault

In a first approach for the convective model setups, we carried out sensitivity analyses for the opening width ( d _frac ) of the implemented fault zone, as well as for its maximum depth ( z _frac ) and permeability ( k _frac ). The effect of these model parameters was also tested with the specific BGF-model setup, used for sensitivity analysis of the model extensions. For sensitivity analysis of k _frac , we uniformly set d _frac to 0.3 m.

Figure  ⁴ c shows that, compared to the temperature distribution from a solely conductive heat transfer, the three tested permeability values of the fault ( k _frac1 – k _frac3 ) lead to a considerable temperature rise within the sedimentary cover, but with diverging temperature distributions with increasing depth. This implies an enhanced upward convective heat transport via fluid circulation along the fault. The lowest fault permeability k _frac1 yields a relative continuous temperature rise down to 1600 m depth. In contrast, the higher permeabilities k _frac2 and k _frac3 lead to higher temperatures in the upper part of the investigated profile P (0–800 m). However, in the lower part (800–1200 m), a considerable decrease of the temperature gradient is observed. This is caused by an enhanced upward heat transport from a lower crustal level due to higher fault permeabilities. Within the area of the sedimentary cover, the thermal profile is convexly upward shaped, resulting from convectional heat transfer in this part of the model (Fig.  ⁴ c). In contrast, the temperature profile of the basement rock layers shows a linear increase in temperatures with depth, resulting from the purely conductive heat transport in this part of the model.

The maximum depth of a fault zone ( z _frac ) might also influence the temperature field, caused by convectional fluid circulation along a permeable fault. Therefore, an additional sensitivity analysis of this parameter was carried out. z _frac values of 1500, 2000 and 3000 m were applied in the BGF-model which comprises:

• the largest granite G 4 at 4000 m depth,

• a normal BHF of 0.070 W m ⁻² ,

• the lowest tested fault permeability k _frac1 .

The results of the simulations imply that in our model setup, z _frac does not influence the evolving temperature distribution with depth, resulting in equal temperature-depth profiles from varying z _frac at the location of the profile P (Fig.  ⁴ c). This is also valid for all other scenarios and tested k _frac values and implies that the maximum depth extension of a permeable fault into impermeable basement rocks cannot be predicted from the presented model. As all tested z _frac scenarios reach the impermeable basement layers (where heat is only transferred by conduction), the long simulation period of 300 Ma always led to similar steady-state conditions at the end of the modelling, independent of the applied z _frac value. The isotherms in the area of the fault show steep convex upward shapes (Fig.  ⁵ ), typical for convective heat transport [ ⁴³ ]. A closer look at the temperature distribution in the surrounding of the fault zone reveals that the temperature elevation in the sedimentary cover is restricted to areas very proximal to the fault. The lateral sphere of influence is ~ 500 m and reduces towards the surface (Fig.  ⁵ ). The only localised influence of hydraulically active faults on the temperature field has also been postulated for the German Molasse Basin [ ⁶⁶ ].

Fig. 5

The sensitivity analysis for d _frac was done with three d _frac values: 0.1, 0.3 and 1.0 m. z _frac was held constant at 2000 m and the permeability was chosen to be k _frac1 . Figure  ⁴ d shows that d _frac variance considerable influences the evolving temperature distribution with depth. Temperatures in the upper part of the temperature-depth profile increase with increasing d _frac , whereas a very high d _frac value of 1.0 m leads to lower temperatures below 1400 m, compared to the implementation of smaller d _frac values. As subsurface fault geometries and their hydraulic properties usually are not known very well, they mostly have to be estimated for thermo-hydraulic modelling. In this study, we uniformly set d _frac to 0.3 m as a larger fracture width is assumed to be unrealistic. However, the high impact on the temperature field of both, d _frac and k _frac values demonstrate that the hydraulic parameterization of fault zones is a significant source of error in the quantification of temperature distributions in this study and in thermo-hydraulic simulations in general.

Model setups without granite

Basic model (B-model): Enhanced BHF in a pure conductive regime

For the scenario of an enhanced BHF acting as heat source at the bottom of the model, we established two model setups:

• the basic model (B-model), where only conductive heat transport is allowed and,

• the basic model with fault (BF-model), where additional convectional heat transport is allowed along a permeable fault zone and within the pore space of the sedimentary cover (Fig.  ³ a).

Both model setups were run with BHF values of 0.070–0.125 W m ⁻² . Figure  ⁶ illustrates the resulting temperature–depth profiles for different BHF values in the conductive B-model. The temperature curves for different BHF values imply that they highly impact the temperature distribution in the overburden. At higher BHF values, the temperature in the upper part of the model rises linearly. A BHF of 0.070 W m ⁻² results in a temperature distribution expected from a normal continental grad T of 30 K km ⁻¹ . An increase in BHF of 10 mW m ⁻² then generates about 8% higher temperatures within the sedimentary cover (Fig.  ⁶ ).

Fig. 6

The temperature–depth profiles of the B-model show a linear temperature rise with depth. The temperature log breaks correspond to the thermal conductivity contrasts of the model layers (Table  ¹ ), which is especially obvious at the boundary between the basement rocks and the lower conductive sedimentary cover (Fig.  ⁶ ).

Basic model with fault (BF-model): Enhanced BHF in a convective regime

To test the effect of an enhanced BHF in combination with a deep-reaching permeable fault, the basic model with fault (BF-model) was chosen (Fig.  ³ a). In this setup, convection is allowed in the model layers of the sedimentary cover and along the implemented fault zone. The area of the basement rocks is impermeable and only heat conduction is allowed. We run the simulations for both, the lowest and the highest fault permeabilities ( k _frac1  +  k _frac3 ). As the maximum depth of the fault does not influence the results of the presented model setup (Fig.  ⁴ c), z _frac was uniformly set to 2000 m in all models discussed below.

Compared to the results of the conductive setup of the B-model (Fig.  ⁶ ), the presence of a permeable fault ( k _frac1 ) and convection within the pore space of the sedimentary cover increase the heat accumulation derived from the deep-seated heat source, in the overburden by about 10% (Fig.  ⁷ ).

Fig. 7

However, the temperature distribution derived from the BF-model is strongly controlled by k _frac . The lowest fault permeability k _frac1 yields a relative continuous temperature rise with depth. In contrast, the highest fault permeability k _frac3 leads to significantly higher temperatures in the upper part of the model and causes temperature stagnation in the lower part (800–1200 m), accompanied by significantly lower temperature gradients in this part of the model (Fig.  ⁸ ). Nonetheless, regardless of the set k _frac , a linear increase in heat source power results in a linear temperature increase in the model (Fig.  ⁸ ).

Fig. 8

Model setups with granite

Sensitivity analysis of granite intrusion depth

In a first approach for the model configurations with a heat-producing granite (BG-model + BGF-model, Fig.  ³ a), a sensitivity analysis was carried out concerning the influence of the absolute depth level of the pluton on the resulting temperature distribution. As granites are potential targets for geothermal projects [ ⁶⁹ , ⁷¹ , ⁷² ], a proper estimation of their depth from thermo-hydraulic simulations would be priceless.

However, Fig.  ⁹ implies that the absolute depth of a heat-producing pluton does not exert a considerable influence on the temperature distribution in the sedimentary cover after 300 Ma. This observation is valid for both tested types of heat transfer. Intrusion depths of 2000, 4000 and 6000 m result in almost similar temperature–depth profiles, with only slightly higher temperatures at shallower intrusion depths of the pluton (Fig.  ⁹ ). This is likely related to the fact that the heat source is always situated within the area of the impermeable basement, in which only conductive heat transfer occurs. The time span of the simulation is sufficiently long for reaching similar steady-state conditions independent from the depth position of the pluton. It can be concluded that the absolute depth of a heat-producing pluton—at constant basin configuration and heat production rate—only at shorter simulation times (or younger intrusion ages) would generate a signature in the temperature field.

Fig. 9

Basic model with granite (BG-model): various dimensions of heat-producing granites in a conductive regime

In the conductive basic model with granite (BG-model), we tested four geometries of moderately heat-producing granites, which have various cross-sectional areas of 11.55–300 km ² (Table  ³ ). As the depth of the pluton top does not significantly influence the temperature distribution after 300 Ma (Fig.  ⁹ ), we uniformly set the granite top at 4000 m depth. For all tested underground scenarios with granite, the BHT at the bottom of the model was set to 0.070 W m ⁻² and A _granite is set 6 μW/m ³ (Table  ⁴ ).

By applying pure conduction (BG-model), the small granite volumes G 1 and G 2 with cross-sectional areas of 11.55 and 35.0 km ² only cause a small temperature rise in the overburden, compared to the temperature profile of a normal continental grad T (Fig.  ¹⁰ ). In contrast, the larger granites G 3 and G 4 lead to considerably higher temperatures in the upper part of the model. Figure  ¹⁰ also suggests that the thermal impact of a heat-producing pluton linearly increases with increasing volume.

Fig. 10

A granite geometry G 2 (BG-model) and a slightly enhanced BHF of 0.080 W m ⁻² (B-model) cause a similar temperature rise in the sedimentary cover (Fig.  ¹⁰ ). The temperature effects of the granites G 3 and G 4 (BG-model) are comparable to the heat supply by a BHF of 0.095 W m ⁻² and 0.115 W m ⁻² (B-model), respectively (Fig.  ¹⁰ ).

Figure  ¹¹ illustrates the isotherms of the BG-model applying granite geometry G 4 at 4000 m depth after 300 Ma forward modelling. The heat production of the granite induces a temperature rise below and above the pluton, with its maximum in the central part of the heat source. The isotherms have a broad convex upward shape, which rapidly flattens with increasing distance to the heat source, considered typical for purely conductive cooling [ ⁴³ ].

Fig. 11

Basic model with granite and fault (BGF-model): Various sizes of heat-producing granites in a convective thermal regime

In the BGF-model setup, we tested the effect of additional convective heat transport in the granite model. The simulations were run for the fault permeabilities k _frac1 and k _frac3 .

Compared to the conductive granite model setup (BG-model; Fig.  ¹⁰ ), the granitic heat-source has a higher thermal impact on the overburden, if convection is allowed (Fig.  ¹² ). The effects of different k _frac values are independent of the heat source in the model. As already demonstrated in Fig.  ⁸ , the higher fault permeability k _frac3 leads to an elevated temperature in the upper part of the model, accompanied by a flatter grad T in the lower part of the sedimentary cover. The significantly lower k _frac1 value causes a more continuous temperature rise with depth.

Fig. 12

The temperature field is furthermore controlled by the size of the granite, independent of whether conductive or convective heat transport is considered (Figs.  ¹⁰ , ¹² ). The thermal impact of a deep-seated heat-producing pluton linearly increases with its volume. The simulations confirm that any convectional component generally increases the efficiency of heat transport from the heat source [ ²⁴ ]. However, elevated temperatures due to fluid circulation on a hydraulically active fault zone are limited to areas very proximal to the fault (Fig.  ⁵ ).

The insulation effect of a low conductive sedimentary cover

Another factor contributing to the formation of local temperature anomalies in sedimentary basins might be the presence of thick, low conductive sedimentary rock layers above a heat source. As thermal conductivity exerts a first-order control on the temperature field in the underground [ ²⁴ ], a significant effect on the resulting temperature distribution in the model can be expected by the variation of this parameter. Low conductivity rocks potentially act as insulators below which heat can be accumulated, resulting in high-temperature gradients (e.g., [ ⁵¹ , ⁷³ , ⁷⁴ ], also at economically relevant depths < 5 km [ ⁵¹ ].

To investigate the effect of missing sediment cover with rock layers of low thermal conductivity (e.g., clay- and marlstones) in the modelled part of the sedimentary basin, the thermal conductivity contrast between the sedimentary cover (mean: 2.8 W m ⁻¹  K ⁻¹ ) and the basement (phyllite: 3.7 W m ⁻¹  K ⁻¹ ) was removed. In the BG- and BGF-model setups, the thermo-hydraulic properties of the sedimentary layers were replaced by the input parameters of the uppermost phyllitic basement rock layer (Table  ¹ ). The removal of the thermal conductivity contrast results in a temperature–depth profile (green curve) that provides only slightly higher temperatures than the temperature distribution generated by a normal continental grad T (Fig.  ¹³ ). The absence of an insulating sedimentary cover prohibits the accumulation of heat generated by the heat-producing granite. This demonstrates the important role of a low-conductive, insulating blanket for the emergence of temperature anomalies in sedimentary basins, charged by any deep-seated heat source. Moreover, the effect of decreasing thermal conductivities with increasing temperatures at depth can even reinforce the insulating effect of a sedimentary cover [ ⁷⁵ ].

Fig. 13

The extreme effect of variable thermal conductivity of the sedimentary cover layers in the models presented confirms that this parameter is a first-order control factor for the temperature distribution in the subsurface [ ²⁴ ]. This in turn points out that for valid thermo-hydraulic models, a careful parameterization of the thermal properties of the subsurface lithologies in the target area is essential. The acquisition of statistically valid thermo-hydraulic data by extensive measurement campaigns, therefore, is essential, if there is no local data pool available from the target area.

Implications for the temperature anomaly in the Franconian Basin

Our models allow a first evaluation of the impact of the tested heat sources and geological conditions in generating the Franconian Basin (Mürsbach) temperature anomaly in SE Germany. Because the extrapolation of temperature logs generally bears considerable uncertainty, the modelled thermal profiles are only discussed down to the final depth of T -Log1_MueB6 at 1204.2 m, as the input data of this model section are based on measured data.

The modelled temperature distributions imply that an average BHF (generating a normal continental grad T in the conductive B-model, Fig.  ⁶ ), is not sufficient for achieving the measured temperatures in the case study area for both tested assumptions—a solely conductive or combined conductive/convective heat transport (Fig.  ¹³ ). However, there are several potential scenarios for achieving the elevated temperatures observed in the T -Log1_MueB6. By assuming a purely conductive regime, measured temperatures below 1000 m depth are reproduced by the following:

• an elevated BHF of 0.115 W m ⁻² or

• the presence of a heat-producing granite of large size ( G 4) and a heat production rate of 6 μW m ⁻³ , combined with an average BHF of 0.070 W m ⁻² at the bottom of the model (Fig.  ¹⁴ , Table  ⁴ ).

  Fig. 14

  

However, the measured temperatures in the Triassic section of the well at 300–1000 m depth (Table  ¹ ) could not be reproduced by a purely heat conductive scenario. In this section, the temperature log of well MueB6 is rather convex upward, pointing to a convective heat transport component. The measured temperatures were reached applying additional heat transport along a fault of k _frac1 for both scenarios described above (Fig.  ¹⁴ , Table  ³ ). The temperature curves from fault permeability k _frac1 in combination with a THF of 0.115 W m ⁻² reproduce the log section between at 300 and 1000 m depth. Similar results are achieved by the BGF-model with granite geometry G 4. However, modelled temperatures derived from these scenarios are slightly too high in the well section below 1000 m. This suggests that below 1000 m (Zechstein section; Table  ¹ ), heat may predominantly be transferred by conduction, whereas in the upper part of the sedimentary cover also convection seems to play a key role in heat transport. This assumption is confirmed by the modelled temperature profile of a modified BGF-model setup with granite geometry G 4 (BGF-model mod.), where only conduction is allowed in the Permian section of well MueB6 (Fig.  ¹⁴ ). Such a scenario provides the best reproduction of the measured temperature log, including the 900–1200 m section.

Thus, an elevated BHF in combination with a permeable fault zone of k _frac1 or the scenario of a heat-producing granite of cross-sectional area G 4 in combination with an average BHF of 0.070 W m ⁻² would provide the best explanations for the elevated temperatures in the case study area (Table  ⁴ ). A purely conductive regime in the centre of the investigated temperature anomaly can be excluded, as the measured temperature log cannot be completely reproduced by the purely conductive scenarios tested. Similar conclusions have been drawn in a study of the North Alpine Foreland Molasse Basin, where heat transport in the deeper part is mainly conductive, but strongly influenced by fluid flow within the overlying sediment fill [ ⁶⁶ ].

However, preliminary results of gravity modelling [ ⁵³ ] do not indicate the presence of such a large granite volume in the deeper underground of the well MueB6 location. Therefore, convective heat flow along hydraulically active faults and/or an elevated BHF might provide a more realistic explanation for the local high temperature. As the presence of deep reaching faults is known from local seismic surveys, a key prerequisite for convective heat transport exists. However, a permeable fault alone without additional heat source—may it be an enhanced BHF, a heat-producing granite, a combination of both or other sources—is not sufficient for achieving the measured temperature profiles.

Implications and outlook for future geothermal exploration in the Franconian Basin

The presented conceptual 2-D models show that the thermo-hydraulic parameterization of the implemented fault zone, the basement layers and the sedimentary cover has a very high impact on temperature distributions resulting from thermo-hydraulic simulations. In addition, the volumes and heat production rates of the heat sources tested highly affect the evolving temperature field. With regard to a future geothermal exploration in the Franconian Basin, the lacking knowledge of the structure of the deeper subsurface as well as of the hydraulic properties of fault zones is a major data gap, which must be closed by further research. Currently, an extensive seismic survey program is being conducted in the Franconian Basin aiming to get more detailed information on the basement topography and the structural framework of the deeper subsurface [ ⁷⁶ ].

Relatively simple, conceptual models, as the one presented here, are a valuable tool for the early stage exploration of geothermal reservoirs. For example, they can be used as plausibility tests for the evaluation of the thermal impact of specific heat sources and of the thermo-hydraulic boundary conditions, which are suspected as causing thermal anomalies in comparable geological settings. This is in particular true when dealing with small data sets at different scales. Herein, our model setup can be used as a first approach to verify hypotheses of specific underground scenarios which then can be further developed for a more complex 3-D simulation. They may also be used for an initial assessment, whether the volume of a heat-producing granite would be large enough to provide detectable signatures, e.g. in gravity modelling.