Crustal-scale thermal models: revisiting the influence of deep boundary conditions

The societal importance of geothermal energy is significantly increasing because of its low carbon-dioxide footprint. However, geothermal exploration is also subject to high risks. For a better assessment of these risks, extensive parameter studies are required that improve the understanding of the subsurface. This yields computationally demanding analyses. Often, this is compensated by constructing models with a small vertical extent. This paper demonstrates that this leads to entirely boundary-dominated and hence uninformative models. It demonstrates the indispensable requirement to construct models with a large vertical extent to obtain informative models with respect to the model parameters. For this quantitative investigation, global sensitivity studies are essential since they also consider parameter correlations. To compensate for the computationally demanding nature of the analyses, a physics-based machine learning approach is employed, namely the reduced basis method, instead of reducing the physical dimensionality of the model. The reduced basis method yields a significant cost reduction while preserving the physics and a high accuracy, thus providing a more efficient alternative to considering, for instance, a small vertical extent. The reduction of the mathematical instead of physical space leads to less restrictive models and, hence, maintains the model prediction capabilities. The combination of methods is used for a detailed investigation of the influence of model boundary settings in typical regional-scale geothermal simulations and highlights potential problems.


Introduction
Geothermal energy is an important part of the future energy mix on the path to a more sustainable use of resources. Many aspects influence the potential use of a geothermal resource, with one prime parameter being the temperature in the subsurface. To determine expected temperatures on a regional scale, geothermal simulations are often performed (Gelet et al. 2012;Kohl et al. 1995;O'Sullivan et al. 2001;Taron et al. 2009;Watanabe et al. 2010). A common procedure is to start with a geological model, representing the main geological sequences, grouped by similar thermal properties, and to use this information for the parameterization of a geothermal simulation (Cacace et al. 2010;Mottaghy et al. 2011;Sippel et al. 2015). However, the (effective) thermal parameters of subsurface geological units (e.g., thermal conductivity, heat production rate) are generally uncertain and the material parameters are therefore often calibrated on the basis of temperature observations. Extensive parameter studies or full uncertainty quantification studies are non-trivial since basin-scale models tend to be computationally demanding. To overcome this issue, a common approach is to generate models that have a large horizontal extension but a very small vertical extent (often only a few kilometers) that can be up to 40 times smaller than the horizontal extent (Freymark et al. 2019;Kastner et al. 2015;Noack et al. 2013;Pribnow and Clauser 2000;Pujol et al. 2015). The boundary conditions for these models are either based on best estimates or retrieved from larger models (Noack et al. 2013). This work investigates in detail how these typical approaches to treat boundary conditions influence all subsequent analyses, leading partly to fully boundary-dominated models. In this paper, it is demonstrated that they only have very limited capabilities for the analysis and understanding of the physical processes. During the model calibration, a compensation for possible boundary errors through an adjustment of the thermal properties is possible. Consequently, this has no direct impact on the temperature distribution but a significant impact on the physical plausibility of our model. Hence, for scenarios that lay outside of the calibrated regime, any prediction capabilities are lost. This is a major restriction when considering the sparse nature of observations. The models with a small vertical extent are commonly used, although it is well known that diffusion problems are majorly impacted by the boundary conditions. Therefore, this paper illustrates the consequences of this model choice and demonstrates that crustal-scale models are crucial for basin-scale applications.
To investigate the influence of thermal boundaries, full global sensitivity analyses (SA) are employed for several case studies. These types of global SA approaches are usually not performed due to the high associated computational cost. To address these computational challenges, the full finite element solution of the forward solves is replaced with the reduced basis solution. This approach aims to reduce the complexity of the mathematical instead of physical space, yielding fast, accurate, and physics-preserving surrogate models. With these surrogate models, global sensitivity analyses are performed on several model realizations of a regional-scale geothermal basin model in northern Germany (around Berlin and the state of Brandenburg) to demonstrate the influence of the lower boundary condition on the simulation.
Additionally, an automated model calibration is executed to provide an objective and reproducible way to compensate for the errors of both the physical and geological model. Sensitivity analyses for basin-scale models have been performed before in Noack et al. (2012) and also been combined with automated model calibrations (Wellmann and Reid 2014). Also, Fuchs and Balling (2016) consider model calibrations but in their case without sensitivity analyses. Furthermore, local sensitivity studies are presented in Ebigbo et al. (2016). However, none of these can address the computationally demanding nature of the problem. Therefore, they are limited in the number of parameters, sensitivity analyses, and model calibrations they can perform. Using a physics-based machine learning approach instead of the finite element method, the computation time of the forward solve is reduced by several orders of magnitude. It allows, in turn, to perform global sensitivity analysis and full flexibility in the model calibration.
Global sensitivity analyses have been performed for hydrological problems (Baroni and Tarantola 2014;Cloke et al. 2008;Song et al. 2015;Tang et al. 2007;van Griensven et al. 2006;Zhan et al. 2013), for volcanic source modeling (Cannavó 2012), and for geothermal heat exchangers (Fernández et al. 2017). In Degen et al. (2021), the authors have investigated the influence of both local and global sensitivity studies for the Upper Rhine Graben. In this paper, the combination of the global sensitivity study and model calibration, as presented in Degen et al. (2021), is used to investigate the influence of the placement of the boundaries on the model predictions.
The paper is structured as follows: The methodologies and the governing equations are presented in Sect. 2 and in Sect. 3, the problem of the lower boundary condition is conceptually introduced using a simple 1D model. Section 4 presents the impact of the lower boundary conditions, by focusing on a real-case basin-scale application. Therefore, the results of both global sensitivity analyses and model calibrations are presented and discussed.

Materials and methods
In the following, the geothermal conduction problem used for the forward simulations of the temperature is briefly described. Furthermore, the concept of sensitivity analyses is introduced.

Physical model
For the simulation of the temperature field, a geothermal conduction problem with the radiogenic heat production S as the source term (Bayer et al. 1997) is considered: where is the thermal conductivity, and T the temperature. For efficiency reasons and to investigate the relative importance, the nondimensional form of Eq. 1 is considered. Therefore, the nondimensional properties T * = the analysis on the heat diffusion and the boundary condition. Furthermore, for all models Dirichlet boundary conditions are applied at the top and bottom of the model domain.

Sensitivity analysis
Sensitivity analyses aim to determine which model parameters influence the model response to what extent. So, these studies investigate, which thermal conductivities and radiogenic heat productions have a significant impact on the temperature distribution. One distinguishes two types of sensitivity analyses: local and global ones. Local sensitivity analyses consider that all parameters are independent of each other. In contrast, global sensitivity studies investigate also the parameter correlations. A detailed comparison of both methods for hydro-geological problems is presented in Wainwright et al. (2014) and for basin-scale geothermal application in Degen et al. (2021). For the sensitivity analysis (SA), a quantity of interest needs to be defined. Here, the L2-norm of the temperature misfit to the reference model is used as the quantity of interest, for the motivational example. The quantity of interest for the real-case model is the L2-norm of the temperature misfit between the simulated and observed temperature values.
For the global sensitivity analysis, the Sobol method with the Saltelli sampler is used, this is a variance-based sensitivity analysis operating in a probabilistic framework. From the variances, the sensitivity indices are derived as the ratio between the partial and total variance. In this work, the main interest is on the first-and total-order indices. The first-order index is the ratio between the variance of the parameter and the total variance and defines the impact of the parameter itself. In addition, the total-order index captures all parameter correlations. This includes second-order but also any higher-order terms. Second-order terms describe the correlation between two parameters only, whereas higher-order terms define the correlation between multiple parameters. Further information regarding the Sobol method can be found in Sobol (2001), Saltelli (2002) and Saltelli et al. (2010). For the sensitivity analyses the python library SALib (Herman and Usher 2017) is used.

Model calibration
The main aim of this paper is to investigate the influence of the lower boundary condition on the physical interpretation through an evaluation of the temperature distribution. This is the reason why global sensitivity analyses are used. In practical applications, it is often desired to calibrate the model against existing temperature measurements to ensure the correctness of the model. For this, model calibrations are required, which aim to compensate for existing model errors by an adjustment of the model parameters. For deep geothermal applications calibrations are challenging since one usually has only a few shallow data points (Degen et al. 2021). As the real-case study will show, it is possible to adjust a given model to the observed temperatures. However, larger model errors yield unphysical model parameters, imposing the danger of losing the predictability for observation points that have not been included in the calibration. This aspect will be discussed in detail later on.
In this work, a trust region reflective algorithm is employed as the calibration method, which is a suitable choice for constrained problems, meaning that the thermal parameters have defined ranges (Branch et al. 1999). During the calibration, the L1 norm of the misfit between the simulated and observed temperature measurements is minimized. The L1 norm is considered to put less weight on outliers. The analysis is performed through the python library SciPy (Virtanen et al. 2020). For more details regarding the method, refer to Branch et al. (1999) and more details regarding the application to basin-scale models refer to Degen et al. (2021).

Motivational example
This paper investigates the influence of the impact of the lower boundary condition on the temperature distribution. This is an issue concerning geological models in general. For this reason, the problem is first demonstrated using a highly simplified motivational model. The motivational study aims to illustrate the general problems and not to represent a realistic geothermal application. To demonstrate that the issue has a major impact on real-case geothermal applications, the investigation is extended to the real-case study of Berlin-Brandenburg (a sedimentary basin in northeastern Germany which is introduced in Sect. 4).

Forward model
First, the forward problem used for the motivational study is introduced, for which a simplified 1D model is considered.
The 3-layer model, schematically shown in Fig. 1, consists of three layers, where the middle layer is thinner than both adjacent layers. A thermal conductivity of 1.0 is chosen for the top and bottom layer and a thermal conductivity of 1.5 for the thin layer. To recall, throughout the entire section, the dimensionless formulation is used. Consequently, the thermal conductivity has no unit. At the top of the model, a Dirichlet boundary condition of zero is applied for the temperature, and at the bottom a Dirichlet boundary condition of one. The model is solved analytically. Note that the nondimensional form is considered to focus the analysis on the relative difference.
In the following analyses, the influence of the thermal conductivity of the thin middle layer (Layer 2 in Fig. 1) with respect to its distance from the boundary conditions is analyzed. Therefore, the position of the thin layer changes. Three different positions of the thin layer are considered: (i) the thin layer adjacent to the base boundary condition (position P1 in Fig. 2), (ii) the thin layer in the center of the model (position P2 in Fig. 2), and (iii) the thin layer adjacent to the top boundary condition (position P3 in Fig. 2). For the sensitivity analysis, scenario P2 is defined as the reference model, where the thin layer is located around the center (see Fig. 1). Consequently, the reference model represents the case of the lowest possible boundary influence.

Impact of the boundary condition
To determine the influence of the lower boundary condition, a global sensitivity analysis with 100 equally spaced temperature measurements in depth ranging from zero to one is performed. Equally spaced measurements are chosen to avoid any bias induced by the spatial distribution of the measurements in the sensitivity analysis. Furthermore, the thermal conductivities of all three layers have an allowed variation range of ± 50 %.
The results of the global SA are shown in Fig. 2. Before discussing the results for this SA, the terminology needs to be specified. From Fig. 2 first-and total-order terms are obtained. The first-order terms describe the influence from the parameter itself, whereas the total-order term describes the influence from the parameter plus any parameter correlations. Consequently, the correlation is defined as the difference between the total-and first-order contributions. This motivational study investigates the influence of both boundary conditions on the model. Therefore, it needs to take the scenario, where the thin layer is in the center of the model (P2) as the reference case. This means that high influences of the parameters correspond to a high boundary dominance.
For the simple model, all thermal conductivities are dominated by total-order contributions for all three scenarios (P1-P3). This means that the parameters have high correlations. The high correlations are induced by the setup of the model, where the temperature distribution is only determined by the two Dirichlet boundary conditions and by the ratio of the thermal conductivities between adjacent layers. Furthermore, the influence of 2 is at all three positions the lowest, which is an effect of the lower thickness of this layer. Also note that for 2 , nearly no first-order influences are observed.
Focusing on scenario P1, the highest boundary dominance is achieved for 1 , because its corresponding layer is situated at the upper boundary condition. The lowest influence is obtained for 2 , although the layer is at the lower base condition. The reason for the lower influence is the smaller thickness of the layer. Note that the influence of 2 is much higher in P1 than in P2. 3 has a significantly lower influence of the boundary than 1 , which is logical since its associated layer is further away from the boundary. Interesting is that the decrease in the first-order contributions is more pronounced than the decrease in the total-order contributions. This shows that the remaining boundary influences are mainly arising from parameter correlations. By having a detailed look at the SA, one observes that the main correlations are arising from the correlation between 1 and 3 . For scenario P3, the same behavior with reversed roles for The results for all three scenarios are following the expectations since the smallest boundary influences are observed if the layers are further away from the boundaries.
Note that these results can only be returned by a global SA. A local SA would assume that the influence is coming from the parameter itself. As an example, in P1 this would lead to a significant overestimation of the influence of 3 . In the worst case, this yields the misleading conclusion that 3 is still greatly influenced by the boundary.
To conclude, for the motivational example the information about the thin layer is lost when it approaches the boundary condition. Or, as an alternative viewpoint, these two examples highlight the strong influence of boundary conditions on the simulation results. In a typical geothermal simulation setting, the position of the top boundary condition is usually defined as the land surface and cannot be changed. Its impact and possible ways to solve the issue have been discussed in Degen et al. (2020a). In contrast, the position of the lower boundary condition is usually adjustable.

Case study Berlin-Brandenburg
After the demonstration of the general problem of the placement of the boundary for geological models, the consequences for real-case studies are illustrated. Therefore, the simplified 1D example is exchanged with various representations of the Berlin-Brandenburg model, which cover a sedimentary basin in north-eastern Germany (see Fig. 3).

Berlin-Brandenburg models
This paper uses three different versions of the Berlin-Brandenburg (BB) model. The model is located in the southeastern part of the Northeast German Basin, which is part of the Central European Basin System. The formation of the basin started in the Late Carboniferous / Early Permian with a period of extensive volcanism (Benek et al. 1996;Noack et al. 2012). Permian and Cenozoic sediments are deposited above the volcanic rocks (Noack et al. 2012). The model is mainly characterized by mobilized Upper Permian Zechstein salt, which forms salt pillows and diapirs due to halokinetic movements (Noack et al. 2012;Scheck et al. 2003). Also, the deeper crustal domains of the model are further differentiated to account for the different consolidation ages (Noack et al. 2012). For further information regarding the geological background, refer to Noack et al. (2012Noack et al. ( , 2013. The area is of interest for geothermal studies due to a temperature anomaly consisting of high heat flow values. This anomaly stretches from Poland to the river Elbe (Noack et al. 2012).
In the following, the numerical discretizations of the Berlin-Brandenburg models are presented.
The first version of Berlin-Brandenburg, from now on denoted as the Berlin-Brandenburg LAB model (BB-LAB), has already been presented in Noack et al. (2012) and can be seen in Fig. 3a. It has an extension of 250 km in the x-and of 210 km in the y-direction and extends down to the lithosphere-asthenosphere boundary (LAB). The model consists of 15 lithological units and the mesh consists of deformed eight-noded prisms. The grid resolution is one km in the horizontal directions, whereas the vertical length of the layers corresponds to the vertical element length, resulting in a mesh with 840,000 degrees of freedom.
The second model, in the following, referred to as the Berlin-Brandenburg 6 km model, or BB-6 km (Fig. 3b), has the same horizontal extent but extends to a depth of 6 km instead of down to the LAB. It is presented in Noack et al. (2013) and consists of 12 lithological units. The model is discretized into a tetrahedral mesh. In comparison to the Brandenburg LAB model, it is refined in both geological and grid resolution terms. The horizontal element resolution is 0.22 km 2 and vertical resolution is interpolated from the z-evaluations of the geological layers with a minimum thickness of 0.1 m, resulting in a mesh of 1,546,675 degrees of freedom.
Combining the Berlin-Brandenburg 6 km model, the Berlin-Brandenburg LAB model, and removing the minimal vertical thickness of 0.1 m results in the third version of the Brandenburg model, denoted as the Berlin-Brandenburg combined model, or BB-combined (Fig. 3c). Consequently, this model consists of 17 geological layers, where the upper For both the BB-LAB and BB-combined model, a Dirichlet boundary condition of 8 • C, corresponding to the average annual temperature, is applied at the top of the model. Moreover, the Dirichlet boundary condition at the base of the LAB is set to 1300 • C corresponding to the melting temperature of the mantle rocks (Turcotte and Schubert 2002). Additionally, a variation of the temperature at this boundary condition of ± 10 % is allowed to account for errors in the geometrical description of the LAB. The Berlin-Brandenburg 6 km model has the same upper boundary condition, but at the base, various Dirichlet boundary conditions directly taken from the Berlin-Brandenburg LAB model are considered. Furthermore, a lower boundary conditions derived by Kriging is taken into account. For this interpolation, 900 equally spaced temperature observation from the BB-LAB model in a depth of 6 km are considered and the interpolated boundary is derived with a spherical variogram. All thermal properties are summarized in Table 3 in the Supplementary Material. The forward simulations are performed using the DwarfElephant package (Degen et al. 2020b) with a linear and nonlinear solver tolerance of 10 −10 . Due to the nondimensional nature of the problem, no preconditioners are needed for the finite element evaluations.
The reference thermal conductivity ref is equal to the maximum thermal conductivity of the BB-LAB model of 3.95 W m −1 K −1 . For the BB-LAB and the BB-combined model, the maximum temperature of 1300 • C is the reference temperature T ref , whereas for the BB-6 km model a reference temperature of 8 • C is chosen. Homogeneous Dirichlet boundary conditions are used to achieve a better performance of the numerical methods (Degen et al. 2020b). The Berlin-Brandenburg 6 km model has a constant Dirichlet boundary condition at the top. At the base, the model has a Dirichlet boundary condition with a different temperature value for each element. The top boundary condition is normalized to zero by using the value of the top boundary as the reference parameter. The bottom boundary condition is set to zero via a lifting function. In case of the Berlin-Brandenburg LAB and combined model, the models have constant Dirichlet boundary condition values for both upper and lower boundary, and hence one can use both of them as the reference parameter. The value of the lower boundary condition was chosen to better reduce the magnitude of the temperatures, which yields a better performance. The maximum radiogenic heat production of the BB-LAB model of 2.5 W m 3 is the reference radiogenic heat production S ref . The reference length l ref corresponds to the maximum x-extent of all models (250,000 m).
For the validation of the models the temperature measurements presented in Noack et al. (2012Noack et al. ( , 2013 and based on Förster (2001) were used. The observations consist of 81 temperature measurements from 44 wells in the area of Brandenburg. It has been measured at various depth and stratigraphic levels.

Reduced models
The reduced basis (RB) method is a model order reduction technique that aims to significantly reduce the dimensionality of problems resulting from a discretization (e.g., via finite elements) of parameterized partial differential equations (PDE). The method is decomposed into an offline and online stage, where the offline stage, being a one time cost, constructs a reduced basis, and therefore comprises all expensive pre-computations.
The online stage uses this reduced basis to allow very fast forward evaluations, typically in the range of a few milliseconds (Degen et al. 2020b). In contrast to other surrogate models, the RB method has the advantage that the physics is preserved. Other surrogate model techniques build their models upon observations (Miao et al. 2019), without explicitly considering the PDE. The RB method maintains the input-output relationship, meaning that the structure of the original finite element problem (and consequently the PDE) is preserved (Hesthaven et al. 2016). Hence, the method allows an extraction of the entire state vector (e.g. the temperatures at every node of the model). Furthermore, for geothermal conduction problems, it provides an error bound, enabling an objective evaluation of the approximation quality. For further information regarding the RB method, refer to Hesthaven et al. (2016);Prud'homme et al. (2002); Veroy et al. (2003) and for further information in the context of geosciences refer to Degen et al. (2020b).
For using the RB method, the geothermal problem is decomposed into a parameter-dependent and -independent part. In the following, the affine decompositions of the integral formulation of the PDE for the various scenarios of the Brandenburg model are defined. Note that this paper uses the operator representation. Therefore, it presents the bilinear form instead of the stiffness matrix, and the linear form instead of the load vector.
For all Berlin-Brandenburg models, the bilinear form a has the following decomposition: where w ∈ X is the trial function, v ∈ X the test function, "q" denotes the index of the training parameter (for more information see Tab. 3 in the Supplementary Material), X the function space ( H 1 0 (Ω) ⊂ X ⊂ H 1 (Ω) ), Ω the spatial domain in ℝ 3 , ∈ D the parameter, and D the parameter domain in ℝ n . The number of thermal conductivities in the training sample is denoted with n. Consequently, n is equal to thirteen, nine, and fourteen for the BB-LAB, BB-6 km, and BB-combined model, respectively.
For all Berlin-Brandenburg models, except the BB-6 km model with a lower boundary condition derived via Kriging, the linear form f is decomposed in the following way: Here, Γ is the boundary in ℝ 3 , s the scaling parameter for the lower boundary condition, g(x, y, z) the lifting function, T top the temperature at the top of the model, h(x, y, z) the location in the model, z bottom (x, y) the depth of the bottom surface, and d(x, y) the distance between the bottom and top surface.
For the BB-6 km with a Kriging lower boundary condition, the linear form slightly changes to the following: Here, g 1 , g 2 , and g 3 are again the lifting functions, with s 1 being the nugget, s 2 the partial sill, s 3 the scaling parameter for the mean temperature, and a the range.

Parameterization and set-up of the sensitivity analysis
The sensitivity analyses are performed with 13 (BB-LAB model- Fig. 3a), 11 (BB-6 km model- Fig. 3b), 14 parameters (BB-combined model- Fig. 3c) and with 10,000 realizations for each parameter to reduce the statistical error. Note that for the Berlin-Brandenburg 6 km model exemplarily the results using the Kriging lower boundary condition are shown. The results of the sensitivity analyses using the other boundary conditions are analog to the one shown in this manuscript. In this paper, only the thermal conductivities are varied and the radiogenic heat productions are kept constant, to reduce the number of parameters within the (5) reduction and all further analyses. The radiogenic heat productions are fixed and not the thermal conductivities because their influence on the overall temperature distribution is smaller. In Table 1, a list of all rock properties is provided.
A variation of ± 50 % from the initial thermal conductivities is allowed for all thermal conducitivities. Also, for the nugget and the partial sill, a variation of ± 50 % is enabled. For the scaling parameter of the lower boundary of both the Berlin-Brandenburg LAB model and Berlin-Brandenburg combined model a variation ± 10 % and for the scaling parameter of the mean temperature at the lower boundary condition of the BB-6km model ± 20 % is used, to account for the uncertainties related to those boundary conditions.

Results
As for the conceptual study, this work demonstrates the influence of the lower boundary condition. Therefore, first the results from the sensitivity analysis and then the results from the model calibration are presented.

Sensitivity analysis
Before presenting the results of the sensitivity analyses, note that all analyses were performed with the aim to investigate the influence of the lower boundary condition. The paper does not aim to characterize the influences of every single thermal parameter in the model. Nevertheless, some geological impacts can be derived and are presented in the following.
Regarding the sensitivities, the Berlin-Brandenburg LAB (Fig. 4a) is mostly influenced by the Lower Cretaceous/ Jurassic/Buntsandstein layer. The first-order sensitivity index is dominant over the higher-order indices. Furthermore, the model is sensitive to the Quaternary/Tertiary layer and the Lithospheric Mantle. For the Quaternary/Tertiary layer, one again has predominantly first-order influences, whereas the Lithospheric Mantle mostly impacts through higher-order contributions. Less pronounced is the influence from the Zechstein layer. The observed influence has similar firstand higher-order contributions. This is counter-intuitive since one would expect a high influence of the Zechstein layer due to its high thermal conductivity and highly variable thickness resulting in significant property contrast. To explain this discrepancy, a closer look at the set-up of the sensitivity analysis is required. In the analysis, layers with equal thermal conductivities were combined. Therefore, the thermal conductivities of the Lower Cretaceous, Jurassic, and Buntsandstein layer are combined. Consequently, the high influence of this layer is originating from this high combined sediment thickness. Keep in mind that the aim of this analysis is to determine the influence of the boundary condition. For determining which individual thermal conductivity has the highest influence a separate analysis is required. The remaining thermal conductivities have minor influences and are, therefore, disregarded in further analyses.
The Berlin-Brandenburg 6 km model is only influenced by the Basement layer and by the variability of the lower  boundary condition (Fig. 4b). The influence of the scaling parameter of the mean temperature is significantly higher than the one from the Basement layer. Higherorder contributions dominate both parameters. Note that the Basement layer has nearly no first-order contributions, whereas the scaling parameter has non-dominant first-order contributions.
For the Berlin-Brandenburg combined model (Fig. 4c), one observes a similar pattern. The highest influences, dominated by first-order contributions, are arising from the Lower Cretaceous/Jurassic/Buntsandstein layer. The influence of both the Lithospheric Mantle and the scaling parameter of the lower boundary condition increased, but higher-order contributions still dominate both parameters. The Tertiarypre-Rupelian-clay/Upper Cretaceous, and the Zechstein layers are also influencing on the model and comparable first-and higher-order contributions to each other.

Model calibration-temperature distribution
The results from the global sensitivity analysis are taken as an input for the following model calibration. Therefore, only the influencing model parameters are considered as shown in Fig. 4. Hence, four model parameters for the Berlin-Brandenburg LAB, two parameters for the Berlin-Brandenburg 6 km, and five parameters for the Berlin-Brandenburg combined model have to be taken into account for the model calibration. The remaining parameters are kept constant within the calibration since the sensitivity analysis identified them as having no impact on the temperature response. Model calibration is necessary to account for model errors of the Berlin-Brandenburg model.
The calibration of the Berlin-Brandenburg 6 km model is challenging because of the lower boundary condition. The conventional way to define this boundary condition is to extract it from the calibrated BB-LAB model and apply it to the BB-6km model, although it is generally not clear that the calibration for the larger model is also valid for the shallower model. To evaluate the influence of different calibration results, the model calibration for the shallow model using the boundary condition from two uncalibrated Brandenburg LAB model versions and various hierarchical model calibrations are compared. For the hierarchical models, either the boundary condition from the calibrated BC or a boundary condition obtained via Kriging as the lower boundary condition are chosen. Therefore, Fig. 5 compares the model calibrations using various lower boundary conditions of the Berlin-Brandenburg 6 km model. At the top panel, it shows the difference at the observation points. The differences between the various methods are comparably small, which is not surprising since the calibration aims to minimize the difference between the simulated and observed temperatures at these locations. However, if one looks at the three points (P1-P3, positions shown in Fig. 5), one observes differences between the various calibrations that can exceed 50 • C. This means that for temperature prediction for points included inside the calibration data set good fits are obtained (regardless of the chosen boundary conditions). This changes once the points outside the calibration data set (P1-P3) are considered, here significant differences for the different boundary conditions are obtained. This is of great importance for geoscientific applications since many studies face the problem of data sparsity. The model has many regions, where no data is available. Still, these regions might be of major importance. Consequently, it is desired to obtain models that are physically plausible to maintain the predictability of the models. To conclude, one can fit every model to a given temperature data set, with the consequences that the thermal conductivities get partly unphysical. This is less important if the target area coincides with a high data density. However, this is often not the case. Therefore, the need to ensure that the generality of the model is preserved remains. Figure 6 compares the temperature distributions for the interval of the uppermost 6 kilometers of all three versions of the Berlin-Brandenburg model. For the BB-6km model, exemplarily the hierarchical model calibration is shown. The differences for all three points (P1-P3) are comparable among the models. Note that the possible variation range of  Noack et al. (2012Noack et al. ( , 2013 of all layers for and after the automated model calibration The radiogenic heat production is denoted with S, and the initial thermal conductivity with init .

ID
Layer the BB-6km is much larger since the determination of the lower boundary condition is uncertain (see Fig. 5). The BB-LAB and BB-combined model already show the maximum possible variation, whereas the BB-6km model shows only the maximum variation range of the good-fit model. Lastly, Fig. 7 shows the differences in the temperature distributions at the three points (P1-P3) for the entire depth of the BB-LAB and BB-combined models. The major difference between both models is induced by the different treatments of the boundary condition. During the sensitivity analysis of the BB-LAB model, the scaling parameter of the lower boundary condition did not significantly influence the model response, contrary to the analysis of the BB-combined model. Therefore, in the latter model the scaling parameter in the calibration is considered, wheres the value is kept constant for the former model. Although, with a maximum temperature increase of 10 % a great amount of variation is allowed, the possible variations at a depth of 6 km are comparable to those of the Brandenburg 6 km model.

Discussion
In the following, the dangers of constructing models with a small vertical depth are demonstrated. To further illustrate the importance of the placement of the lower boundary condition, first its impact is demonstrated by using the results of the global sensitivity study. Afterwards, the consequences for inverse processes are emphasized, using a deterministic model calibration. Both analyses are presented for the case study of the Berlin-Brandenburg model.

Sensitivity analysis
The impact of the lower boundary condition is apparent by focusing on the difference between the BB-6 km, and the BB-LAB and combined models. For the Berlin-Brandenburg 6 km model, the boundary condition is fixed at 6 km depth, resulting in an entirely boundary dominated model. This is observable due to the enormous sensitivity of the model to the: -Basement layer, -scaling parameter of the respective boundary condition, and -correlation between both parameters.
Consequently, all information that is obtained from the Brandenburg 6 km model is coming from the boundary condition. Hence, the model is uninformative concerning the upper layers. However, these are the layers that are of interest since the target region is within these layers. Loosing the information about the thermal conductivities means that only the boundary is determining the solution. Hence, any errors of the boundary conditions have a possible huge impact on the temperature distribution at the target depth. This demonstrates that generating diffusive models with an extremely small vertical to horizontal length ratio is to be avoided at any cost.
The results of the global sensitivity analysis of the BB-LAB and combined model are matching the expectations. A high sensitivity is observed for the upper layers, which is caused by the shallow measurements (500 m to 6820 m).  First-order contributions of the Lower Cretaceous/Jurassic/ Buntsandstein layers mostly impact the model. That means that the thermal conductivities of these layers are influencing the model themselves and not through a correlation with other layers. For the BB-LAB model, the thermal conductivity of the Quaternary and the Tertiary layer were combined into one training parameter. For the Brandenburg combined model, the thermal conductivities of the Quaternary and Tertiary-post-Rupelian, and the Tertiary-pre-Rupelian-clay and Upper Cretaceous were combined. Comparing the sensitivity analysis of both the BB-LAB model and combined model, one can conclude that the Tertiary-pre-Rupelianclay is the layer that the model is sensitive to. The Quaternary, and the Tertiary-post-Rupelian layer can be ruled out because the Berlin-Brandenburg combined model is insensitive to it. Furthermore, also the Upper Cretaceous can be eliminated because the Berlin-Brandenburg LAB model is insensitive to it. Also, the influence of the thermal conductivity of the Tertiary-pre-Rupelian-clay is mainly originating from the parameter itself and not from interactions between various parameters. Again, the influence of the Tertiary-post-Rupelian-clay seems counter-intuitive due to its low thickness. This influence is a combination of the shallow measurements, which lead to higher influences for the upper layers and the Dirichlet boundary condition at the top. This boundary conditions fixes the temperature for each evaluation to the same value, yielding a reduced influence of the Quaternary and, therefore, a relatively higher influence of the Tertiary layers.
Additionally, for both models a significant influence of the Lithospheric Mantle is retrieved. Higher-order contributions dominate this parameter, and the second-order sensitivity indices show the parameter is correlated to the scaling parameter of the lower boundary condition. The Zechstein layer has similar influences in both model versions and is less significant in comparison to the overall influences.
To conclude, the only meaningful way to construct the model is by inserting the refined model into the original Berlin-Brandenburg LAB model. This results in the BBcombined model, which again shows the expected sensitivity distribution. One needs to keep in mind that this means an increase in degrees of freedom from 1,546,675 to 2,141,550. Nonetheless, both the finite element and the online execution time for both models are comparable since the complexity in these two models remains similar. This demonstrates that a reduction in the mathematical and not in the physical space is advantageous since it is much less restrictive.

Model calibration
At first hierarchical model calibrations seem to be a way to transfer the knowledge from large-scale coarse models to smaller-scale fine discretized models. However, the sensitivities clearly show that the smaller model becomes uninformative towards the upper layers. That is especially dangerous because it is not noticeable looking at the temperature distributions at the observation points only. Hence, at a first glance, one would get to the conclusion that cutting-of the model at 6 km is a valid approach. However, this would only be possible if our sole interests are the temperatures at the measurement points used within the calibration. Naturally, a calibration will match the simulation to the observed temperatures. However, that comes at a cost. For the various model calibrations of the BB-6 km model one obtain thermal conductivities ranging between 1.49 W m −1 K −1 and 2.83 W m −1 K −1 for the Basement layer. Meaning that no longer physical thermal conductivities but effective ones are retrieved. These effective thermal conductivities are tailored to our measurements. However, if a different location (e.g. new drill-hole location) is of interest, one can no longer derive reliable temperatures since the model calibration is not valid for this point and the model lost the information about the physical system. This reveals the next important point. The abovedescribed procedure is valid in a limited application field. However, one should be aware that the model is no longer representative of the physical processes. In contrast, both the BB-LAB and combined model have significant influences from various thermal conductivities. The lower boundary condition is further away from the target area, reducing possible effects from this condition.
In general, one wants to improve through global SA the understanding of the physical model. In this specific case study, a way to determine the most influencing parameters allowing a back correlation to the geoscientific context was demonstrated. Note that both the SA and the calibration focus on the observation locations. Hence, higher influences of shallower layers are observed. A study focusing solely on the temperatures at certain locations is applicable for some geophysical studies but if the interest goes beyond fitting the temperatures it is not advisable to use models that are cut-off at a shallow depth.
Note that the changes for the thermal conductivities were not discussed in detail here. The reason is that the discussion of this paper focuses on the influence of the boundary condition. For further information about the thermal conductivities, refer to the Supplementary Material.

Outlook
Through this study, the path to subsequent tasks is opened. It would be interesting to further investigate the lower boundary condition. For some of the calibrations, very high thermal conductivities of the Lithospheric Mantle were obtained, which might be caused by the geometrical inaccuracies of the LAB. These inaccuracies would impact the lower boundary condition and the calibration would try to compensate for this by adjusting the thermal conductivity of the Lithospheric Mantle. A scaling factor to the temperature value of this boundary to account for these inaccuracies was applied, which slightly improved the results. However, a single parameter is not enough to compensate for the model errors. Therefore, it would be interesting to replace the scaling factor by a function, which could be, for instance, determined through data assimilation. For this reason, a promising next step to take would be to investigate if 3D-Var data assimilation yields improved results. In contrast to classical sequential data assimilation techniques, such as the Ensemble Kalman Filter (Burgers et al. 1998;Evensen 1994), variational data assimilation is a continuous approach, where the entire time frame is considered. Variational data assimilation methods minimize a cost function to obtain an estimate of the state variable. Three dimensional variational data assimilation has been studied intensively in numerical weather forecast by, for instance, Barker et al. (2004) and Lorenc et al. (2000) but is fairly unknown for geothermal simulations. It has been studied in combination with the RB method already by Aretz-Nellesen et al. (2019). However, so far, the study is using benchmark problems only. Therefore, it would be interesting to investigate the performance of the method for complex geophysical problems.

Conclusion
Throughout the entire paper, the high impact of the lower boundary conditions for conductive crustal-scale applications was demonstrated. Using a novel combination of reduced-order modeling techniques and global sensitivity analysis, the paper illustrated that cutting-of models at a shallow depth has severe consequences. For these models, the information content of the geological structures is entirely lost. This is of utmost importance if one aims to derive physical knowledge from the model and or want to perform predictions with the given model. These findings should be well known, still, it is a common procedure to construct models with a small vertical extent. Therefore, this work aims to explicitly show the consequences of this approach. The clear visualization of the boundary problem becomes only apparent through the utilization of a global sensitivity analysis since this method allows also the investigation of parameter correlations. Note that the value of a "too" small vertical extent differs for each model since it is dependent on various factors such as the type of boundary condition, the geological structure, and the governing physical principles. This further highlights the importance of sensitivity analyses to reliably determine whether a model is boundary-dominated. To construct informative models with a smaller vertical extent one could use, for instance, the Moho as the base boundary condition and apply a Neumann boundary condition, which is less restrictive than a Dirichlet boundary condition. Another possibility is to use optimal experimental design techniques to determine a feasible depth of the model.