Temperature and pressure corrections applied to rock thermal conductivity: impact on subsurface temperature prognosis and heat-flow determination in geothermal exploration

Precise knowledge of the subsurface thermal field plays a key role in the assessment of geothermal targets. Unfortunately, deep underground temperature data is generally scarce and a matter of research. To achieve first estimates for subsurface temperatures, steady-state conductive thermal modeling is commonly applied. Thereby the rock thermal conductivity is an essential parameter, which is usually determined under ambient laboratory conditions. To arrive with in situ thermal conductivity, the ambient values need to be corrected for in situ temperature and pressure. In this paper, we apply different conversion functions for the correction of thermal conductivity and study the impact on the resultant temperature and heat flow prognoses for a synthetic, upper crustal sedimentary and a magmatic scenario along 2-D geological cross sections. Application of the correction functions results in maximum temperature prognosis uncertainties of about 8 °C and 55 °C at 2 km depth and at 8 km depth, respectively. The effect positively correlates with the magnitude of the basal heat flow used in modeling. In contrast to the heat flow determined at depth, the resulting surface heat flow is only minor affected by the different correction functions applied. In addition, the modeled temperature at depth is strongly dependent on the type and sequence of application of the pressure and temperature correction equations.

commonly used as boundary conditions, and the rock-specific thermal properties of the rocks. Under steady-state conductive conditions, one of the most influential parameter is the thermal conductivity (λ, in W m −1 K −1 ) of the rock, which in turn, largely depends on λ of minerals constituting the rock matrix and on λ of the pore fluid (e.g., Blackwell and Steele 1989). Furthermore, λ is pressure (p) and temperature (T) dependent (e.g., Schön 1996). Thus, values of λ determined under ambient conditions in the laboratory differ from values under higher p and T conditions at greater depths. Therefore, for improving the reliability of thermal models, not only a proper mapping and parameterization of geological (and thermal) units is necessary, but also an understanding of the change of thermal rock properties due to different p/T conditions.
Many laboratory experiments were conducted in which λ was measured as a function of either increasing p or T. The general trend from the experiments is often referred in textbooks (e.g., Schön 1996Schön , 2011: λ(Τ) is decreasing and λ(p) is increasing with increasing T or p, respectively. Thus, considering a normal geothermal setting (i.e., where T and p increase with depth), the p and T dependence of thermal properties operate conversely. Usually, p and T dependencies of λ are measured separately. Only few data on limited rock types are available, where λ was measured simultaneously as a function of p and T (Abdulagatov et al. 2006;Emirov et al. 2017). Perhaps this is one reason why corrections to in situ λ are often ignored in thermal modeling (e.g., in the lithosphere thermal models of Noack et al. 2010;Przybycin et al. 2015;Freymark et al. 2017). In other cases, correction to λ is claimed being performed, but the correction functions are not provided (e.g., Hasterok and Chapman 2011). Research that only addresses the T dependence of λ (e.g., Vosteen et al. 2004;Rühaak et al. 2010;Fuchs and Balling 2016;Lemenager et al. 2018) acknowledges the overriding effect of T over p, which results in a decrease of λ with increasing depth (Seipold 1995). Studies in which both, p and T correction functions are deployed, were published, for example, by Norden et al. (2008Norden et al. ( , 2012, Förster et al. (2010Förster et al. ( , 2018, Schütz et al. (2014), and Schintgen et al. (2015a, b).
There is abundant literature on different functions describing individual p and T correction to λ. However, only little effort was directed to the quantification of temperature effects that stem from the diversity of these correction functions. This also holds true for quantifying their relevance for geothermal exploration. In addition to ambiguities inherited from the geological setting being modeled (i.e., the structural model and the considered rock types and their initial (ambient) thermal properties), uncertainty in modeled T is strongly controled by the in situ p and T conditions and the in situ λ. For the central part of the Fennoscandian Shield lithosphere, Kukkonen et al. (1999) discussed variations in calculated temperature and heat-flow density by varying the values of input parameters [λ, heat production (H, in µW m −3 ) and lower boundary condition (constant heat-flow density or constant temperature)] as well as the p/T-correction parameters for λ. This study on parameter variations concluded that the p and T effects on calculated T are minor compared to uncertainties resulting from variations in initial λ and H or the choice of the lower boundary condition of modeling. Lee and Deming (1998) studied the sensitivity of the various T correction equations to λ available for igneous, metamorphic, and sedimentary rocks in the T range < 500 °C. The study resulted in a delineation of absolute and systematic errors in λ, however, did not proceed towards a quantification of the effects on subsurface temperatures. In conclusion, the effects on subsurface temperatures resulting from not applying or only partially applying correction functions remain largely unaccounted for.
The present paper focuses on the thermal impact of p/T corrections at the depth range relevant for the exploration of natural resources and their subsurface utilization (i.e., the uppermost crust). In order to allow an evaluation of different correction schemes for calculating in situ λ, we only consider steady-state and pure conductive conditions. For a specific setting, transient effects may superimpose the conductive regime, thus significantly affect the thermal regime. Here, we want to concentrate on the general effect of p/T corrections on the thermal field and, therefore, disregard any possible transient processes. In such a conductive regime, temperature prognoses for the depth realm of borehole drillings are fundamental, e.g., for drill site selection and project development. For those projects it is essential to know ambiguities in temperature prediction stemming from ambiguities in calculating in situ λ. In addressing this issue, we apply available λ correction functions in modeling the 2-D temperature distribution along two synthetic geological cross sections. While one section (scenario A) considers a crust entirely composed of igneous and metamorphic rocks, the second section (scenario B) includes sedimentary rocks in the upper part of the profile. These 2-D models quantify T by either (I) not using any p/T correction to λ, (II) using either a correction of p or T, and (III) using combined p/T correction algorithms.
The conditions examined in this paper cover crustal temperatures of < 700 °C. At those T, diffusion processes (the lattice or phonon thermal conduction) are typical for polycrystalline rocks exhibiting a decrease in λ with increasing T (Schatz and Simmons 1972;Beck 1988;Seipold 1992). This decrease in λ is variable, depending on rock type (e.g., Sekiguchi 1984;Zoth and Haenel 1988;Somerton 1992;Seipold 1998;Kukkonen et al. 1999;Vosteen and Schellschmidt 2003;Miranda et al. 2019). The p domain applicable to the crust covers values < 1 GPa. In this range, λ increases steeply (in an exponential fashion) from ambient conditions towards a pressure of 50-100 MPa (see Hurtig and Brugger 1970;Emirov et al. 2017, and in the compilation by Fuchs and Förster 2014). At higher pressures, λ increases gently in a linear fashion. Common linear dependencies are those from, e.g., Horai and Susaki (1989) and Seipold (1992Seipold ( , 1995Seipold ( , 2002. Again, different dependencies are observed for different rock types (e.g., Seipold 1998;Seipold and Huenges 1998;Kukkonen et al. 1999). In our study, we focus on consolidated rocks. Changes of λ based on the compaction of sediments during sedimentation are not regarded.
For the reader's convenience, we provide the Additional file 1 to this paper that consists of a compilation of the bulk of the p and T correction functions available in literature, including the conditions under which these relations were derived and the rock type studied. It is apparent that certain functions have limited applicability to the continental crust due to the p and T conditions under which they were constrained. A detailed discussion on the validity of these algorithms, however, is beyond the scope of this paper. Figure 1 shows the variability of common T and p corrections to λ. The corresponding functions are listed in Additional file 1: Tables S1-S3. The graphs are arranged with respect to the investigated rock type in three main groups: (a) sedimentary rocks and (b) mafic (magmatic or metamorphic) rocks and (c) felsic (magmatic or metamorphic) rocks. All graphs show either the p dependency (greyish area, positive change of λ 0 with increasing pressure) or the T dependency (negative change of λ 0 with increasing temperatures). To account for the different initial λ associated with the respective formula, the functions were normalized. So, for each rock type, the shown functions refer to the same ambient laboratory start value (λ 0 ).

Comparison of correction functions
Most of the corrections for T follow a non-linear trend for crystalline as well as sedimentary rocks, while the available p corrections demonstrate a linear shape for most of the p range, except for the parameter range 0-50 MPa, in which a non-linear, strong increase of λ is observed. The combined p-correction equation for sedimentary and magmatic rocks elaborated by Fuchs and Förster (2014;rocks FF in Fig. 1, see also Additional file 1: Table S2) represents a correction which originates from integrating different published data sets within the parameter range of p < 400 MPa and λ 0 -values of 1.5-5 W m −1 K −1 ). For sedimentary rocks, only the p corrections of Fuchs and Förster (2014) and Emirov et al. (2017, for sandstone only) are available, while available T corrections are numerous. The shape of the curve for sediments after Vosteen and Schellschmidt (2003) generally corresponds with the shape of the curve of Zoth and Haenel a b c Fig. 1 Relative change of thermal conductivity (λ) based on different p/T corrections (p corrections with increased λ, grey shaded area; T corrections with decreased λ): a for sedimentary rocks with a start value at ambient conditions (λ 0 ) of 2.8 W m −1 K −1 (except for salt inset with λ 0 = 5.8 W m −1 K −1 ); b for mafic rocks (magmatic or metamorphic) with a λ 0 of 2.4 W m −1 K −1 , and c for felsic rocks (magmatic or metamorphic) with a λ 0 of 3.0 W m −1 K −1 . The main rock type as specified by the respective author is denoted. Authors are abbreviated as follows: BL: Blesch et al. (1983); EM: Emirov et al. (2017); FF: Fuchs and Förster (2014); KK: Kukkonen et al. (1999); MI: Eq. 4 of Miranda et al. (2019); SO: Somerton (1992); SP: Seipold (2001); VS: Vosteen and Schellschmidt (2003); SG: Sekiguchi (1984); ZH: Zoth and Haenel (1988) (1988) derived for limestone, while the correction function supplied by Somerton (1992, cf. also Anand et al. 1973) shows slightly lower correction values running smoothly counter-rotating. The Sekiguchi (1984) correction is based on estimating T-corrected λ values for the rock matrix. In order to calculate bulk λ, a constant porosity of 13% and T-dependent fluid thermal properties, as suggested in Deming and Chapman (1988), are applied in Fig. 1a (TC SG ). The resulting curve is very close to the one based on the Somerton (1992) correction and shows a trend similar to the T correction of Emirov et al. (2017), who measured one type of sandstone. The study of Emirov et al. (2017) considers the p/T dependency of λ in one experimental setup. The authors developed a rock-specific correction function to λ, involving both T and p (see later in text). In Fig. 1, the graphs of the p/T corrections of Emirov et al. (2017) were calculated for an increase of either T or p. Salt rocks display a special thermal behavior: here, a very strong T dependency of λ is observed (shown in the inset graph in Fig. 1a). This result is in contradiction to observations of Durham et al. (1980), indicating an insignificant p-dependent change of λ of salt rocks, at least in the range of 10-50 MPa. The p dependency of λ is better constrained for magmatic and metamorphic rocks. Depending on rock type, the correction functions show different slopes (Fig. 1b, c). All p-correction functions are within a range defined by two generalized functions: the Fuchs and Förster (2014) function for sedimentary plus magmatic rocks at the upper end and the correction function given by Kukkonen et al. (1999) for crystalline crustal rocks at the lower end. For this type of rock, Kukkonen et al. (1999) determined a mean value for the linear p increase. In Fig. 1, their relation was extended to account for the non-linear λ-pressure relation established by Seipold (2001). The T-related change of λ is more diverse for mafic rocks relative to felsic rocks. While most functions for felsic rock types follow more or less the same trend, the application of the function provided by Miranda et al. (2019) results in a notably smaller reduction of λ compared to any other correction function. These authors measured the thermal properties of a granite exhibiting a moderate heat production of about 5 µW m −3 .
The herein conducted analytical study of the different correction functions permits the delineation of a possible minimum and a maximum correction applied to λ (Table 1). Combining the p-correction equation showing the minimum effect on λ 0 with the T-correction equation showing the maximum effect on λ 0 and vice versa enclose a combined maximum and minimum effect on λ 0 , respectively. However, it is still an open question whether calculating the in situ thermal properties by considering the effects of p and Ton λ individually is a valid approach.

Table 1 Applied corrections to estimate the potential impact of p/T-related corrections of thermal conductivity (λ)
For abbreviations see Fig. 1

Rock type
Minimum/maximum change of λ 0 based on p correction (↑λ with pressure) In this study, we calculated λ in situ by three methods: firstly, by applying the T correction before undertaking the p correction of λ [λ(T,p)], secondly, by applying the p correction before performing the T correction of λ [λ(p,T)], and thirdly, by applying a combined correction function of the individual corrections [λ(p + T)] calculated as λ in As long as the respective p/T-correction equation used represents an additive λ 0 system (i.e., λ 0 will not be included as a factor), the third approach will yield the same result as the two other approaches. In this study, most of the considered p-correction equations are respecting λ 0 in a factorized term, while the majority of the T-correction equations do not. For these cases, method three will yield the same results as method two, but differ from those resulting from method one. Furthermore, the T-correction correction formulae for sedimentary rocks of Vosteen and Schellschmidt (2003) and Somerton (1992) are also incorporating λ 0 as a factor. Thus, for their T-correction functions, the third method will show results different from the other two methods. To illustrate the different correction approaches, the following scenario is considered: λ 0 determined under ambient conditions amount to 2.8 W m −1 K −1 and the expected in situ p/T conditions is represented by p = 120 MPa and T = 130 °C, respectively. For the correction of λ 0 to in situ conditions, the equations rocks FF (p correction) and sediments SO (T correction) will be taken into account, respectively (e.g., Fig. 1, Table 1, Additional file 1). The first method [λ(T,p)] will result in a λ in situ of 2.649 W m −1 K −1 , the second method [λ(p,T)] will yield a λ in situ of 2.654 W m −1 K −1 , and λ in situ calculated by the third method [λ(p + T)] will yield a value of 2.767 W m −1 K −1 . Emirov et al. (2017) presented a formulation of the dependency of λ on p and T for one type of sandstone. They observed that p has also a direct influence on the T dependency of λ, implying that a simple arithmetic approach for combining single p/T corrections may represent an unpermitted oversimplification.
To illustrate the impact of the implementation of different correction equations on the temperature field, we studied two synthetic geological models, which will be presented in the next section.

Conceptual crustal models and their parameterization
The 2-D thermal models studied in this paper reflect different geological conditions in the uppermost 7 km of the Earth's crust (Fig. 2). Scenario A covers different types of magmatic and metamorphic rocks. It is inspired by geological models developed for the western Erzgebirge and the adjacent Elbe Zone in the central European Saxothuringian chain of the Variscan orogeny in Germany (Förster et al. 2018), characterized by the emplacement of voluminous igneous rocks in the upper crust. In contrast, scenario B accounts for a sedimentary setting with more or less horizontal layering, only locally deformed by salt tectonics. This scenario encompasses different types of sedimentary rocks and two prominent salt structures (a salt pillow and a salt diapir). This scenario resembles a setting typical for European Permian sedimentary basins in general, and for the North German Basin in particular.
The two conceptual geological 2-D models are subdivided into units/polygons of contrasting thermal properties (λ and H) ( Table 2, scenarios A and B). In the depth range 7-30 km, both upper crustal geological models have the same type of layered crust (Table 2, basement). The structure and composition of these basement units are equivalent to those reported for the Erzgebirge in Germany (Förster and Förster 2000), to account for realistic crustal scenarios. The thermal properties and density values assigned to rock types represent typical ambient textbook values. To investigate the effect of p/T corrections on the temperature field, different simulations are performed with: (1) no p/T correction; (2) only a p correction of λ 0 , to account for the maximum p contribution [λ(maxp)]; (3) only a T correction of λ 0 , to account for the maximum T contribution [λ(maxT)], and (4) combined p/T corrections as discussed previously (see "Comparison of correction functions" section). The correction formulas were selected according to the possible minimum as well as maximum corrections to λ 0 with respect to the main rock type (Table 1, Fig. 1). For scenario B (sedimentary crust), we additionally compare the effect of applying or disregarding the correction formula for the T dependency of the λ of rock salt as well as its possibly p-independent thermal behavior.

Numerical procedure
For 2-D modeling, the equation for the steady-state heat conduction is given by where H represents the internal heat production and λ is assumed to be isotropic. The equation was solved numerically. The temperature distribution T(x, z) within the lithosphere (with (x) the horizontal coordinate and (z) the vertical coordinate) is determined by the temperature-dependent and pressure-dependent thermal conductivity distribution λ(x, z), the heat production distribution H(x, z), and the appropriate thermal boundary conditions. The heat equation is solved by the finite-element approach using the PDE toolbox of the software MATLAB R2018a. The finite-element mesh size is variable according to the different size of the crustal polygons. Thin polygons representing the sedimentary formations or polygons close to the surface have a mesh size of less than 200 × 200 m, and those of larger size, representing the deeper crust, show a mesh size of up to 1 × 1 km. (1) Conceptual 2-D models (uppermost 10 km) a of a magmatic and metamorphic crust (vertical exaggeration approx. fourfold) and b of a sedimentary crust (vertical exaggeration approx. threefold). The model rock units (bold capital letters) and the corresponding thermal conductivity under ambient conditions (in italics, given in W m −1 K −1 ) are denoted (see Table 2)

Table 2 Parametrization of the conceptual 2-D geological models
Thermal conductivity (λ, in W m

Boundary conditions
In the 2-D thermal models, a constant surface temperature of 10 °C is assumed. In addition to the upper temperature boundary, the models involve a lower boundary constraint, representing a constant heat-flow boundary. A value of 30 mW m −2 at a depth of 30 km represents the initial baseline for both scenarios (magmatic and sedimentary). At the site boundaries of the 2-D models, horizontal temperature gradients are assumed to be zero (no horizontal heat transfer). To exclude any significant influence of the site boundary conditions on the modeling results along the main sections, the models were extended by 50 km on each site, respectively. Figure 3 shows the effect of applying the maximum p correction (λ(maxp), considering no T correction) and the maximum T correction (λ(maxT), considering no p correction) on ambient λ values (see Table 2). For comparison, the T calculated with no p/T correction on λ is given. Obviously, the T differences between the different approaches increase strongly with depth. Based on the respective p and T correction applied, the possible in situ T values may deviate from the non-corrected T profile. The maximum difference between λ(maxp) and λ(maxT) could be used as an indicator for the uncertainty of T prognosis using correction functions. It accounts for about 8 °C at 2 km depth and to about 55 °C at 8 km depth, respectively. At the same time, the T dependency of λ is, compared to the p dependency of λ, more pronounced. Although the model scenarios contain different rock types exhibiting different thermophysical parameters, these general observations are the same. However, the rock parameterization gives rise to thermal anomalies in both scenarios. For section A, the radioactive (hot) magmatic monzonite (unit B, Table 2) in the depth interval 2-4 km (Fig. 1) causes a positive thermal anomaly below 3.5 km depth, while the higher λ value of the granitic unit above the monzonite body (unit C) is responsible for a decrease in T in that area (see T profile at 2 km depth, Fig. 3). The high λ value of the basaltic rock complex in the eastern part of section A (unit I) causes the comparably low T values in that section to a depth of about 5 km. For section B, the most prominent feature in the T pattern is due to the salt structures, providing great contrasts in λ values. Due to the chimney effect of the higher conductive salt, T is increased at the salt diapir by about 5 °C close to the surface (section B, 0.5 km). The same chimney effect at the salt diapir is responsible for lower T values (compared to surrounded profile km) at depths of 2.0-8.0 km, with a maximum deviation at about 3.5 km depth. If instead of the correction function sediments VS , the correction formula salt ZH is applied for the salt unit, the chimney effect would slightly increase, however, not changing the overall pattern of the T field. Whereas the maximum T and p corrections account for the maximum possible influence on the simulated thermal field (Fig. 3), the possible influence of a combined correction is shown in Fig. 4.

Results
Due to the differences in thermal properties of the modeled rock units in scenario A, the resulting T distribution for a certain depth level shows a much larger range relative to sedimentary scenario B (Fig. 4). In the calculation, a minimum p/T correction (based on the presented single corrections for p and T) and a maximum p/T correction was considered. Based on the selected correction functions for scenario A (Table 2), different ways of combination of p/T corrections are applied on λ: λ(p,T) with the minimum case λ(maxp,minT) and the maximum case λ(minp,maxT), λ(T,p) with the minimum case λ(minT,maxp) and the maximum case λ(maxT,minp), and λ(p + T) for the minimum case. For the maximum case, λ(p + T) will show the same results as λ(p,T) because the applied T correction is following an additive manner. In 10 km depth, the uncertainty amounts to about 30 °C and increases towards greater depths (not shown in Fig. 4). For scenario B, the same type of combinations of p/T corrections was considered. However, in scenario B, the maximum case (maximum T and minimum p) refers to no p correction (Table 2, scenario B), therefore (although the applied T correction is not additive), λ(minp + maxT) equals λ(minp,maxT) (not shown in Fig. 4b). For both scenarios, calculating T based on λ(maxp + minT) yields a b Fig. 3 Temperature profiles for different depth levels. a For the magmatic-metamorphic scenario and b for the sedimentary scenario considering the respective thermal effect (i) of a maximum pressure (p) correction [λ(maxp)] of the ambient thermal conductivity; (ii) of a maximum temperature (T) correction [λ(maxT)] of the ambient thermal conductivity, and (iii) of no correction of the ambient thermal conductivity on the respective thermal field. For data input see Table 2 slightly larger T differences than based on λ(maxp,minT). In comparison, the bandwidth of the T range for a given depth and correction mode is less pronounced for the sedimentary scenario than for the magmatic-metamorphic cross section (scenario A). This is related to the geometry and the parameterization of scenario A, which is less structurally manifold than scenario B and introduces also lateral heat flow, resulting in a greater T variation. Depending on the configuration of the correction approach, T may differ by about 40 °C in 8-10 km depth.
To shed more light on the possible effects on the T field as function of the geological setting and the linked thermal properties, Figs. 5 and 6 show a selection of T and ΔT profiles of certain positions along the 2-D section for the igneous and the sedimentary scenario, respectively. The uncorrected T profiles (upper part of Figs. 5 and 6) depend on the respective model parameterization and the basal heat flow condition. Thus, the T gradient changes with respect to λ of the involved model polygons. For the igneous scenario, the possible effect of p/T corrections is given by the T difference originating from the different correction modes and for an increased and a reduced basal heat flow. Corrections considering the p dependency of λ may result in a reduction of about 20 °C compared to uncorrected parameters at 10 km depth, while considering only the T dependency may end up with T estimates increased by more than 40 °C compared to T based on uncorrected λ and a basal heat flow of 30 mW m −2 (at 30 km depth). Bordered by the two lines of maximum p and T corrections (given by λ(maxp) and λ(maxT) in Fig. 5) is the branch of possible in situ T conditions. Applying the maximum effect of a combined correction (by correcting firstly for p and subsequently for T or vice versa) has a strong impact on the modeling result. Remarkably, the effect of a λ(maxT) or a λ(T,p) correction on the temperature is on the same order of magnitude as an increase (or decrease) of the basal heat flow in the order of 25% while not applying a correction method of λ (resulting in a ΔT of about 28 °C at 10 km depth, see upper panel of Figs. 5 and 6). If T corrections of λ are applied, a 25% reduction in basal heat flow (22.5 mW m −2 ) results in a less pronounced effect on T, while a 25% increase in heat flow (37.5 mW m −2 ) produces a much greater uncertainty of the calculated T. Figure 6 shows the situation for the sedimentary scenario. If we compare the influence of an increased or reduced basal heat flow on the calculated T correction of λ and the thermal field, the situation is similar to the igneous case: the higher T of an increased heat flow is responsible for an enhanced T correction of λ and, thus, a higher ΔT (Fig. 6, upper and lower panels). In contrast to the general-rock-type correction applied for the igneous scenario (Fig. 5, Table 1), a lithotype-specific scenario (accounting for reasonable rock-specific corrections, see Table 2) is presented in Fig. 6, in addition to applying the extreme corrections [λ(maxp) and λ(maxT)]. While the maximum effect of an applied p correction or T correction is similar to the igneous scenario, resulting in a ΔT of − 20 °C (λ(p)) or + 40 °C (λ(T)) at 10 km depth a1 a2 a3 b1 b2 b3 c1 c2 c3

Fig. 5
Temperature profiles without correction of thermal conductivity (upper part) and temperature differences profiles (subfigures a1-c3) for scenario A (igneous section) at three different locations and for different basal heat flow in 30 km depth. Temperature differences are calculated against the uncorrected temperature profiles (relying on ambient thermal conductivities) based on different thermal conductivity (λ) correction methods: λ(maxp)-maximum of pressure (p) correction on λ only, λ(maxT)-maximum of temperature (T) correction on λ only (both as in Fig. 4); λ(p,T)-combined correction of λ by correcting firstly for p and secondly for T, λ(T,p)-combined correction of λ by correcting firstly for T and secondly for p; λ(p + T)-combined correction of λ by adding the single p/T correction effects. Correction formulas for all ways as shown in Table 2 (with a heat flow of 30 mW m −2 ), the graphs of a combined p/T correction of the lithotype-specific scenario show a much narrower min/max range of about 20 °C (for λ(p,T)) to 25 °C (for λ(T,p)). A distinctive feature of the sedimentary setting in scenario B is the presence of rock salt, which shows a λ higher than most common sedimentary rocks. The rock-specific combined p/T corrections (λ(p,T) and λ(T,p) in Fig. 6) do honor the salt-specific p/T dependency formulations (Table 2). They result in a stronger reduction of the ambient λ and do not correct for p effects, thus lead to comparable higher T values, very close to the maximum T correction [λ(T)] in the salt diapir (a3, b3, c3 in Fig. 6). If, however, a non-salt specific correction formula is applied (e.g., using sediments VS and rocks FF for T and p correction, respectively), the T is very close to the T profile based on the uncorrected ambient λ values for the uppermost 7 km (λ, Fig. 6).  Fig. 4); λ(p,T)-combined correction of λ by correcting firstly for p and secondly for T, λ(T,p)-combined correction of λ by correcting firstly for T and secondly for p; λ(p + T)-combined correction of λ by adding the single p/T correction effects (correction formulas for all ways as denoted in Table 2, rock-specific sedimentary scenario); λ(EM)-using the correction formula of Emirov et al. (2017) for p and T below 400 MPa and 573 K, respectively, otherwise as denoted in Table 2; λcorrection like λ(p,T), but ignoring the salt-specific correction formula for the rock salt layer (unit) and using sediments VS and rocks FF for T and p correction, respectively, instead We additionally applied the Emirov et al. (2017) p/T correction approach to the sedimentary scenario B, assuming that all sedimentary rocks show a similar thermal characteristic as the sandstone investigated in their study (λ(EM), Fig. 6). The resulting T profile (or ΔT in Fig. 6) of the λ(EM) correction shows a remarkable T contour: it follows more or less the reference T profile (based on ambient and uncorrected λ values), indicating also slightly lower values for the uppermost 10 km of the crust. Moreover, it is more or less in line with the no-salt scenario [λ] for the uppermost 5 km indicating that it may be a valid approach for siliciclastic rock types. However, at greater depths, λ(EM) deviates significantly from the other correction modes.

Discussion
The impact of p/T corrections applied to ambient λ values of rocks implemented in thermal field modeling is manifold and depends on the structural-geological setting and the basal heat flow. For the model scenarios treated in this study, the correction equations account for uncertainties in the T simulation on the order of 50-80 °C at 10 km depth, also depending on the basal heat flow (see Figs. 5 and 6, lower panels). With increasing depth, the p effect has generally a minor impact on λ compared to the T effect. Therefore, the p effect remains mostly unaccounted for in thermal modeling (e.g., Lemenager et al. 2018;Miranda et al. 2019). However, disregarding the p correction may give rise to an overestimation of T on the order of 20-30 °C at 7-10 km depth. Although one could argue that such a difference is within the uncertainty of the thermal parametrization of subsurface models, the error reflects a methodological weakness that needs to be addressed in more detail. Our modeling of the synthetic scenarios shows that it is not trivial to account for the in situ λ correction by considering correction functions almost exclusively developed for either a p correction or a T correction to λ. As expected, the correction sequence (first p correction and then T correction, or vice versa) has a direct influence on in situ λ and the calculated T. Normally, λ is corrected first for p and second for T (cf. Kukkonen et al. 1999). From a general point of view, this might be a sound option: p is always increasing with depth and shows normally a much lower order of variability for a certain depth level as T; T may vary strongly with depth due to other heat sources or sinks, such as magmatic intrusions, cooling fluids, etc. De facto, the choice is speculative. In consequence, applying first the p correction and subsequently the T correction will end up with lower T values than the other way round. The reason is that by applying the T correction first, the ambient λ (especially at greater depths and for higher ambient λ values) will be lowered much larger than by the increase of λ in response to the applied p correction (see Fig. 1). Thus, the starting value of λ for the p correction will be significantly lower. This discussion clearly underlines the need of specifically designed experiments, in which the dependence of λ on p and T could be measured simultaneously. The research of Emirov et al. (2017) fills this gap, presenting an attempt for the development of an integrated p/T correction formula. However, their equation relies on too few measurements performed under in situ conditions and is limited to one single rock type (sandstone). Furthermore, the derivation of their equation is not completely transparent, so that supposing a more comprehensive (i.e., independent on rock type) application of the Emirov et al. (2017) equations as exemplarily done in this study (Fig. 6) is not feasible. However, the study of Emirov et al. (2017) may indicate that an individual correction of λ for the respective p and T condition (such as λ(p,T) and λ (T,p) in Figs. 5 and 6) most likely does not meet the real in situ conditions. To sharpen the thermal modeling, more λ measurements considering in situ conditions (accounting for p and T simultaneously) are required to better understand the rock-specific thermal behavior. Until more correction experiments have been conducted, any calculation of the thermal field has to struggle with the discussed ambiguities. In order to get a first approximation of the uncertainty relating to the technical modeling procedure, any of the T models applying (single) p/T corrections to λ should specify in detail, which corrections were applied, and in which sequence.
In connection with the synthetic models the question arises, whether the applied p/T corrections also influence the modeled surface heat flow. In fact, applying or not applying p/T corrections has a minor effect (Fig. 7). Only if strong λ contrasts due to the correction equations are achieved, the surface heat flow is significantly affected. This is the case for salt structures in sedimentary environments. In nearly all other cases, the surface (T−) boundary condition and the similar pressures at shallow depths do not cause any significant change of the surface heat flow. Notably, the thermal properties of the rock types cause a gentle re-distribution of the radiogenic heat. High surface heat-flow values of the magmatic section (Fig. 7a) do correlate with heat refraction processes at rock contacts, while high surface heat-flow values in the sedimentary section are triggered by the chimney effect of the high-conductive rock salt (Fig. 7b).
However, applying or not applying p/T corrections will notably alter the calculated subsurface heat flow. This has consequences for the calculation of heat flow, which generally takes advantage of drillings (drill cores and rock samples) and well-log data (temperature logs). One common way of heat-flow determination is the interval method (e.g., Powell et al. 1988), which involves the calculation of the product of the T gradient and the λ for a certain depth interval according to Fourier's law. If λ is measured on rock samples under ambient conditions, the in situ λ needs to be calculated. Depending on the applied correction functions, the determined heat flow may vary considerably (Fig. 8). In a b general, the higher the starting value of λ 0 and the deeper the considered depth (interval) are, the larger are the deviations. So, just by only considering the available p/T corrections on λ, a heat flow value determined at 4 km depth may have an uncertainty of 10 to > 20%, while the maximum uncertainty at 2-km-depth amounts to only ~ 10%, and would be even less at lower depth. Recent studies investigating the uncertainties of T in basin modeling related to thermal parameterization concentrate on unknowns in the distribution of rock thermal properties in the subsurface (the geological uncertainty of thermal properties, e.g., Elison et al. 2019). However, they did not discuss the influence of the p/T correction on λ and the resulting T field. Elison et al. (2019) compare different modeling approaches for a 2-D generic model, using general T-correction approaches of Sekiguchi (1984) for matrix λ and Zoth and Haenel (1988) for bulk λ, both uncorrelated to rock type. A p correction of λ is not discussed by the authors, but is partly considered by applying Athy's law (Athy 1930), an exponential reduction of porosity with depth that will result in an increasing bulk λ with depth. Noack et al. (2012) investigated the sensitivity of 3-D thermal models to the selection of boundary conditions and thermal properties at a lithospheric scale. They conclude that the values of λ implemented in their models need to be smaller than those determined on rock samples at ambient conditions, to match measured temperature profiles. However, they did not apply any of the p/T corrections on λ in their sensitivity study, which may have solved the problem. Not accounting p/T-dependent thermal properties may lead to wrong model calibrations and false conclusions on the relative importance of conductive versus transient conditions. The ability to assign proper thermal properties in thermal modeling is thus of outmost relevance for interpreting the subsurface thermal field and the thermal processes that define it. In this paper, we discussed heat conduction as the dominant heat transport process at a crustal scale. Transient processes due to the evolution of a sedimentary basin at different time scales (like mechanical compaction and igneous intrusions or changes of the depth of the asthenosphere) a b c Fig. 8 Effect of p/T corrections on thermal conductivity (λ) and determined heat-flow density (HFD), assuming a T gradient of 30 °C km −1 and a p gradient of 26 MPa km −1 for a sedimentary rock with λ = 2.8 W m −1 K −1 , b for a mafic type of rock (λ = 2.4 W m −1 K −1 ), and c for a felsic type of rock (λ = 3 W m −1 K −1 ) applying no correction, only pressure correction, only T correction and combined p/T corrections [λ(p + T)] for the minimum and maximum case (see Table 1) which are undoubtedly in operation in several sedimentary settings, are not considered. In advection-dominated heat transport settings, the conductive thermal field may be totally masked by, e.g., fluid convection. Thus, p/T corrections will have a minor impact on the reservoir simulation of the temperature field of a geothermal target dominated by fluid flow processes. However, also dynamic reservoir simulations require reliable boundary conditions. For instance, constraining the basal (terrestrial) heat flow for a given depth and area requires also a sound understanding of p/T-dependent thermal properties.

Conclusions
Heat flow, thermal conductivity, and radiogenic heat production of rocks are the most influential parameters that affect the temperature distribution in the Earth's interior. The evidence presented in this paper shows that thermal models of the crust heavily depend on the in situ values of λ, which involves the application of p/T correction functions.
To obtain reasonable simulation results of the thermal field at a lithospheric scale, the implementation of both, p and T corrections to λ becomes of increasing relevance with greater depth. Consideration of the entire lithosphere system in thermal models may be advantageous, to better constrain the necessary thermal boundary conditions for small-scale thermal models of the uppermost 5 km of the crust being subject of energy-resources exploration.
Although the modeling performed in this study does not allow to judge on the correctness of a respective function, which is sensitive to rock type and the method of laboratory measurement, the results help to quantify the overall uncertainty of thermal (exploration) models. Advanced laboratory λ data are mandatory, obtained on rock specimens simultaneously exposed to p and T. An improved understanding of the involved thermal processes will support in approaching a new stage of T-field characterization, helping also to evaluate the uncertainty of heat-flow determinations.
Additional file 1. Additional tables with a compilation of p and T correction functions available in literature.