Global sensitivity analysis and uncertainty quantification for design parameters of shallow geothermal systems

.

is available at all times and almost everywhere (Stober and Bucher 2014).The thermal potential of shallow geological space can be exploited using heat pumps.The overall system is referred to as ground source heat pump (GSHP) system.An optimized mode of operation is crucial for low energy consumption as well as the economic use of geothermal resources.
The presented outcomes are results of investigations within the joint research project EASyQuart (Bucher et al. 2024).This research project, which was funded by the German Federal Ministry for Economic Affairs and Climate Protection, aimed to improve the planning of shallow geothermal systems and had a strong application focus.For the planning and design of shallow geothermal systems, it is important to know the conditions at the site as precisely as possible.The performance of the geothermal system and whether applicable regulations are actually complied with depends on this design.However, information about the geological subsurface is usually only available to a limited extent or in insufficient detail.In addition, exploration is often time-consuming and cannot be carried out comprehensively for all parameters at the site due to cost and time constraints.To find out which data are particularly relevant, global sensitivity analyses can be performed.They are used to identify relevant parameters and compare their influences.Furthermore the design of geothermal systems is subject to various uncertainties.These are mainly due to the geological subsurface heterogeneity and to the availability and quality of the related data.Investigations of the subsurface, e.g. by thermal response tests (TRTs), usually provide only selective, averaged information about material parameters at a specific location but not for its surroundings.These uncertainties are amplified if assumptions have to be made using literature data due to a lack of measurement data.As a consequence, over-or undersizing and the possibility of not complying with design regulations may occur.Uncertainty analyses, usually known as uncertainty quantifications, can help to account for these uncertainties in the design process and thus help to improve the design of the geothermal system.In difference to the sensitivity analysis, which considers the influence of design parameters of BHEs in general, the uncertainty quantification deals with uncertainties of the design of a geothermal system at a specific site.Uncertainty analyses are state of the art for investigations in the field of deep geotechnical applications like geological storage of energy carriers or waste or deep geothermal energy (Neuman 1973;Yeh 1986;Ewing and Lin 1991).The investigations within the EASyQuart project dealt with the applicability of this kind of analysis to shallow geothermal systems.In this context, the long-term results were of particular interest to the researchers.
The present study is considering borehole heat exchangers (BHE).Current research projects consider various aspects of BHEs and the influence of heat extraction from the subsurface.This includes sensitivity analyses on load profiles (Gao et al. 2022) or on existing as well as newly developed grout materials to analyze their positive influence on heat transport (Chicco and Mandrone 2022;Badenes et al. 2020).Furthermore, Focaccia (2013) performed a sensitivity analysis of the effects of varying grout properties on the results of TRTs.For deep geothermal systems, sensitivity analyses were carried out by Doran et al. (2021) and Pan et al. (2020).These mainly considered the conditions of the geothermal system itself (e.g. with parameters describing the pipe and grout materials, the borehole depth and diameter or the flow rate).In contrast, a study of Ebigbo et al. (2016) analyzed variations in depth and subsurface temperature, as well as the influence of different types of heat sources.Simulating various scenarios for shallow geothermal systems, Gunawardhana et al. (2015) investigated influences of parameters such as the groundwater flow direction, the surface temperature, and the injection rate.Additional scenarios for subsurface and BHE properties were considered by Casasso and Sethi (2014), where discrete value ranges were specified for the eight parameters considered there.
Geostatistical methods are often used to quantify uncertainties in the geologic subsurface (Chilès and Delfiner 2012;Pyrcz and Deutsch 2014).Furthermore, for thermohydro-mechanical coupled processes, Watanabe et al. (2010) presents another approach to uncertainty quantification using Monte Carlo methods and provides additional parallelization approaches by means of domain decomposition.In the field of deep geothermal energy, model uncertainties were analyzed for the Bavarian Molasse Basin by Ziegler and Heidbach (2020).For shallow geothermal systems, uncertainty quantification have been performed on topics such as borehole resistance (Choi et al. 2022) or the actual position of BHEs in the geological subsurface (Steinbach et al. 2021).
In this paper, a global sensitivity analysis is performed for shallow geothermal systems that use BHEs.To investigate the influence of the design parameters, the performance of one BHE is considered.Initially, 20 parameters are included, covering both the characteristics of the BHE (such as flow rate, dimensions and material properties of pipes, refrigerant and grout) and subsurface properties (thermal and hydraulic).For each of these parameters, value ranges with associated distribution functions were defined to represent reality as closely as possible.The sensitivity analysis is global in the sense, that the whole value range of each input parameter is considered.In contrast, local sensitivity analyses are based on derivations and focus on the effects of small deviations in the input data (Saltelli et al. 2008).
The uncertainty quantification that is described afterward in this paper, is based on the results of the global sensitivity analysis.It is performed for a fictitious example site.This is intended to give an impression for the potential of uncertainty quantification as a planning tool for shallow geothermal systems.

Methods and workflows
This section presents the methods and workflows of the sensitivity analysis and the uncertainty quantification.

Sensitivity analysis
In the sensitivity analysis, the influence of parameter variations on a quantity of interest is examined.For this purpose, the values of the parameters are varied within a previously defined value range.The change of the quantity of interest due to this variation is examined to determine the parameter influences.The workflow of the presented analysis is summarized in Fig. 1.
The first step of the workflow is parameter definition.The sensitivity analyses presented are based on numerical simulations with the open-source software OpenGeoSys (OGS) (Kolditz et al. 2012;Naumov et al. 2022) with a 3D finite element (FE) model.The model is described in more detail in "Computational Model" section.
A software written in Python is used to conduct the sensitivity analysis.It is based on a development by the Helmholtz Centre for Environmental Research (UFZ) in Leipzig for sensitivity and uncertainty analyses in radioactive waste repository research by Buchwald et al. (2020) and was adapted and further developed for application to BHEs.The Python software is coupled with the simulation software OGS via the application programming interface (API) ogs6py (Buchwald et al. 2021).This enables the automated generation, calculation and evaluation of a variety of simulation models, as required in sensitivity analysis.
In the second step of the workflow presented in Fig. 1, parameters of low relevance are excluded in the parameter screening.For this purpose, the screening techniques One-Variable-At-a-Time (OVAT) and Morris method are used.These methods provide information about parameter influences based on the smallest possible representative number of parameter variations (Campolongo et al. 2007).The OVAT method is used in this study as described by Buchwald et al. (2020).Starting from a base setting of the parameters, two calculations are carried out for each parameter with its maximum and minimum values.The other parameters remain at the base value.The deviations of the calculated results to the one with the base values can be considered as a sensitivity measure.With the sampling strategy of the Morris method, more values are calculated for the parameters from their value range.Furthermore, in contrast to the OVAT method, the remaining parameters are not fixed at a reference value with Morris.This creates a more differentiated impression of the parameter influences, which are evaluated with variables µ , µ * and σ .The value of µ describes the mean influence of a parameter.The variable µ * , in contrast to µ , refers to the influence in terms of absolute values, where no distinction is made between decreasing and increasing influences on the result quantity.The larger the value of µ or µ * , the greater the influence of the parameter.A large value for σ , on the other hand, indicates an interaction of the considered parameter with others or a non-linear relationship between them and the result quantity.In this work, the improved version of the Morris method by Ruano et al. (2012) was applied, using an optimized sampling strategy.Further information on the Morris method can also be found in Morris (1991) and Campolongo et al. (2007).Parameters of lower relevance can be excluded after analyzing the results of the OVAT and Morris methods.Since a very large number of sampling points or computationally expensive numerical simulations are necessary for the global sensitivity analysis, a machine learning approximated model (proxy model) was used.This enables the calculation of parameter variations in a fraction of the time needed for numerical simulations.For this purpose, numerical models were calculated in OGS and a proxy model was trained from the results using the technique of Gaussian process regression (Stone 2011;GPy 2014).A detailed description of this method can be found, e.g. in Lubashevsky (2022).In the Python software used, a link to statistical methods such as Gaussian process regression (also known as a variant of the geostatistical method Kriging) is already implemented in such a way that they can be directly incorporated into the sensitivity and uncertainty studies.The parameter constellations for the proxy building were created using a Latin Hypercube Sampling (LHS).The Latin hypercube method was developed by Mckay et al. (2000) and, in comparison to purely random Monte Carlo samplings, offers the advantage that the selection of parameter values is only partially random.By splitting the parameter value ranges into subranges a more uniform coverage of the entire value range is ensured.The values are then randomly selected within the subranges, resulting in a smaller sample size required.Therefore one also speaks of a quasi-random method or quasi-Monte Carlo method.
The global sensitivity analysis is carried out for the remaining parameters using the variance-based method of Sobol' (2001).As results of this method, so-called Sobol' indices are calculated as a sensitivity measure.Sobol' indices are dimensionless sensitivity indices and describe the contribution of a parameter to the total variance of the quantity of interest for the considered value ranges of the input parameters.A distinction is made between Sobol' indices of different orders.First-order Sobol' indices describe the influence of a single parameter.Sobol' indices of second or higher order characterize the influence or the contribution of the interaction of two or more parameters to the total variance of the quantity of interest (Sobol' 2001).The greater the Sobol' index of a parameter or parameter combination, the greater the influence of the parameter or parameters on the quantity of interest.This procedure requires a larger sample size of parameter variations and in return provides more detailed results (Saltelli et al. 2008).So-called Sobol' sequences are used here to create suitable samples.One of the properties of these sequences is that, with the help of scaling, they can represent any distribution of the parameter value ranges and offer good coverage with a comparatively small number of points.Using Sobol' sequences, sampling points are not chosen randomly, even if they imitate randomness.Therefore the method is called a quasi-Monte-Carlo method, as in the case of the LHS presented.More details on the properties of such sequences and their construction rules can be found in Dick et al. (2013).

Uncertainty quantification
Uncertainty quantification is related to the sensitivity analyses.Here, a similar approach is used as in the sensitivity analysis presented.In this approach, the common technique of Monte Carlo simulations is used, for which numerical calculations in OGS were utilized.Figure 2 summarizes the workflow used for the uncertainty quantification.
As in the sensitivity analysis, different variations of parameter values are calculated.Each variation of the parameter values corresponds to a possible constellation that may arise due to the uncertainty of input parameters.The value ranges the parameters vary within for uncertainty quantification are the uncertainty ranges.For this purpose, in a first step, the uncertainty ranges of the input parameters must be determined with a parameter definition.
Since the uncertainty quantification also requires the calculation of a large number of parameter variations, a proxy model was again created with results from numerical simulations in OGS.The proxy building was performed after the parameter definition in the same way as described above for the sensitivity analysis.
The uncertainty quantification is evaluated in two ways.On the one hand, an assessment is made based on the statistical distributions of the model output resulting from the variation of the parameter values.For this purpose, the values of the model output are represented in a cumulative distribution function.The distribution function can be used, for example, to determine the probability that the model output is above or below a certain value.On the other hand, an evaluation is again carried out using Sobol' indices.In the uncertainty quantification, Sobol' indices provide information on the contribution of each parameter to the overall uncertainty in the result quantity.Thus, sources of uncertainty can be distinguished depending on their influence on the result quantity.

Computational model
As mentioned earlier, both parts of this research were based on numerical simulations.For this reason, an FE model was created following the dual-continuum approach of OGS.As described in Al-Khoury et al. (2010) the model consists of 3D-prism elements representing the soil and 1D-line elements representing the BHE.This methodology was further developed by Diersch et al. (2011aDiersch et al. ( , 2011b) ) and implemented in OGS.Furthermore, the numerical calculation is based on capacitor-resistor models (Bauer et al. 2011) defining the interactions between the different parts of a BHE, which enables a detailed analysis of the heat transport in the subsurface.
A basic model was first developed for the purpose of sensitivity analysis and later also used for uncertainty quantification with different parameter settings.To determine fundamental relations between parameters in the sensitivity analysis and to limit the computational effort, the investigation was carried out on a model with one BHE.A double U-tube was selected as BHE type, since this is currently the most commonly used type in shallow geothermal energy in Germany.The subsurface was assumed to be homogeneous for the model, which enables isolated examination of individual parameter influences.The model also includes a layer with groundwater flow, whose depth and thickness are varied in the analyses.The FE model has a cuboid structure with a larger extent of the base surface in the direction of the groundwater flow.The subsurface model is discretized with prismatic elements.Consistent with the dual-continuum approach of OGS, the double-U BHE is represented as the second continuum by a chain of line elements (Shao 2016).The structure of the model is schematically shown in Fig. 3.
To optimize the performance of the model, convergence analyses were carried out with respect to the element sizes and the time-stepping in the transient calculation.In the area around the BHE and at the transitions between areas with and without groundwater flow, a finer meshing was selected than in the rest of the model.The time steps in both studies were set as t = 86400 s = 1 d .To exclude influences by the edges of the model, appropriate distances to the BHE were also maintained.The base surface of the models is 160 m × 230 m .Since the BHE length also varies in the investiga- tions, the model depth was varied accordingly so that a corresponding distance to the bottom of the model is also maintained below the BHE.For meshing, the Python tool bheEASyMesh 1 was developed to be able to automatically create models with varying parameters within the framework of the previously mentioned Python software.

Sensitivity analysis
As mentioned before, a sensitivity analysis is used to investigate the influence of parameter variations within certain value ranges on a quantity of interest.The quantity of interest in this study is the mean temperature of the BHE's refrigerant, hereinafter also referred to as fluid temperature, which results from the mean value of the in-and outflow temperatures of the BHE.The fluid temperature of the BHE is suitable as a quantity for measuring the sensitivities with regard to the design of BHE parameters, as it is the result of the designing process and reflects its quality.To take seasonal variations into account, the simulation period was set to 1 year.The simulation delivers a time-dependent curve of the mean temperature for the simulation period of 1 year.However, instead of a time-dependent curve, a single value is needed for the evaluation of sensitivity.Within this study, the mean fluid temperature was evaluated at different times for this purpose.The results presented here refer to the mean fluid temperature at the time after 1 year of simulation T 365 (simulation start on January 1st).This temperature is a kind of final result of the whole simulation year and due to the fact that it is in the heating period of the year, it gives a good insight how well the BHE performs.In addition, prestudies have shown that an integral measure of the annual variation in temperature is less representative (Richter et al. 2023a).

Parameters and value ranges
For the sensitivity analysis, 20 parameters were examined.These refer to the BHE properties and dimensions, the type of use of the BHE and geological as well as climatic conditions, as shown in Table 1.
Regardless of the used method, the quality of the results of sensitivity analyses is directly related to the definition of the value ranges.Accordingly, careful preparation and definition of the value ranges is required for meaningful results.The value ranges of this sensitivity analysis are intended to represent parameter values in Germany.Various sources were used for their definition.The data are based on manufacturer information, literature references and experience as well as surveys from practice (Stober and Bucher 2014;Hölting and Coldewey 2013;Sanner et al. 2005;VDI 4640 2020;Casasso and Sethi 2014;Ingenieure Verband Beratender 2012;Kretzschmar et al. 2011).For the sensitivity

Geothermal gradient T geo
Variation of the undisturbed subsurface temperature** T geo analysis with Sobol' indices, statistical distributions are additionally assigned to the value ranges.These distributions ensure, for example, that rarely occurring extreme values from the value ranges are also less represented in the sensitivity analysis and thus do not distort the parameter influences.The value ranges and associated distributions are shown in Table 2.Besides value range boundaries, reference values are also given there.On one hand, these are the base values for the OVAT procedure and, on the other hand, they represent the values with the highest occurrence for the specified distribution functions.
The parameter d po is a special case.Definition of the range of values for d po is based on manufacturer specifications.There it can be seen that the wall thicknesses t p of commercially available pipes are usually uniquely determined by the respective pipe diameters.Therefore, the wall thickness was not considered as a separate parameter, but was defined as a dependent parameter of the pipe diameter.For pipe diameters of 32 mm, the wall thickness is usually 2.9 mm and for pipe diameters of 40 mm

Initial and boundary conditions
the thickness is 3.7 mm.To take into account the wall thickness increasing together with the pipe diameter, the following relationship was determined for the sensitivity analysis: Results for the sensitivity of the pipe diameter therefore refer to the simultaneous change of pipe diameter and wall thickness.
The parameters concerning the refrigerant vary in a range between pure water and the values of refrigerants commonly used in BHEs, as mentioned by Stober and Bucher (2014) following Rosinski and Zapp (2007).The most commonly used mixture has an amount of 25% of monoethylene glycol, which is used as reference.The values are constant during the simulation, so that temperature dependencies of the parameters can be neglected.
The parameters annual heat demand Q and variation of the undisturbed subsurface temperature T geo are set with a seasonal curve in the simulation model.For the annual heat demand, this is achieved by a monthly varying heat load, which is shown in Fig. 4.This is based on a curve determined using degree days from the Institut Wohnen und Umwelt GmbH (2021) (IWU).To get a representative value for Germany, the degree days of the cities of Hamburg, Berlin, Erfurt, Düsseldorf, Lahr and Munich were averaged.
T geo as a parameter for the variation of the undisturbed subsurface temperature repre- sents in this study both the initial condition of the temperature in the subsurface and the surface temperature.The surface temperature is incorporated into the model as a seasonal curve through a boundary condition at the ground surface.The temperature curve is based on data from the German Weather Service-Deutscher Wetterdienst (DWD) at the locations mentioned above for the years 2014 to 2018. Figure 5 shows the temperature curve.
On the other hand, the subsurface temperature is incorporated into the model as an initial condition in the form of a temperature profile.This profile corresponds to the curve shown schematically in Fig. 6.The reference profile, which is varied by the parameter T geo , has a temperature of about 10.73 °C in the neutral zone.The part of the curve representing the seasonal zone was determined in a long-term simulation of 100 years based on the influence of the surface temperature presented above. (1) Fig. 4 Annual curve of the heat load (evaporator power) in the sensitivity analysis As mentioned earlier, depending on the parameter constellation, there is a groundwater flow in the model which crosses the model in a horizontal direction.The influence of the flow was investigated by means of three parameters.Two of these parameters characterize the position of the layer with groundwater flow.These are shown schematically in Fig. 7.
The depth from the ground surface to the top of the layer through which the groundwater flows is described by the parameter z aqf .Furthermore, the parameter t aqf describes the vertical expansion (thickness) of the layer with groundwater flow.The Parameters z aqf and t aqf are varied independently.Due to technical limitations of the simulation, the groundwater flow can only start below the seasonal zone of the temperature profile in order not to affect this zone.Therefore, the value range for z aqf is defined such that the minimum depth is 15 m below the top of the subsurface and the maximum value corresponds to the current BHE length.For t aqf the value range starts at 0 m and reaches a maximum of the current BHE length reduced by 15 m.Since multiple parameters are varied simultaneously during parameter sampling, the BHE length also changes parallel to the position and dimension of the layer with groundwater flow.However, mainly groundwater flow in the area around the BHEs is of interest for the investigation.Therefore the parameters are referred to L BHE * , the BHE length reduced by 15 m.Since parameters are to be investigated independently of each other, parameter constellations in which parts of the layer with groundwater flow lie below the BHEs cannot be avoided due to the nature of the model.An example of this is shown in Fig. 7.The third parameter characterizing the groundwater flow is the Darcy velocity v f .

Results of parameter screening
The Parameter screening provides a first impression of the parameter sensitivities.The OVAT method is based on 41 numerical calculations: two calculations per parameter and one simulation with base values.In Fig. 8, results from these calculations are shown in a so-called tornado plot.
The bars in the diagram show the difference of the mean fluid temperature after 1 year of simulation T 365 as the quantity of interest resulting from the calculation with the maximum and minimum value of the respective parameter value range compared to the calculation with reference value (see also Table 2).A significant split between parameters with large and those with small influence cannot be observed.Rather, the decrease in sensitivity is more or less steady.As a criterion for distinguishing between significant and insignificant parameters, the total change in the quantity of interest between the calculation with the maximum value and the calculation with the minimum value was therefore considered.To select the parameter as significant, it was specified that the parameter must cause a minimum change of 0.5 K.According to this criterion, all parameters below the parameter T geo in Fig. 8 can be considered insignificant.
The results from the Morris method are based on 420 calculations in OGS.They are shown in Fig. 9 by the result quantities µ ⋆ and σ of the Morris method.
For further investigations, parameters are relevant that show high influence due to a large value for µ * compared to other parameters.Furthermore, parameters that have a comparatively high value for σ despite a low µ * value is also interesting, as they may have unrecognized influences on the quantity of interest T 365 through interactions with other parameters.The results show that parameters with a large value for µ * , also have large values for σ , except for the parameter T geo , where σ is small compared to the others.A smaller value for σ is an indication of a more linear influence of the parameter on the quantity of interest and of little to no interaction with other parameter influences.In the result plots, parameters with a small influence leading to a low value for µ * do not show any noticeable values for σ either.
Due to its approach, the Morris method offers a more nuanced view on parameter sensitivities than the OVAT method.However, apart from a few exceptions, the sensitivities of the parameters also decrease rather continuously.Thus, the ten least influential parameters were first considered insignificant.These are the same parameters that were also assessed insignificant in the OVAT procedure.In addition, the Morris procedure would assign the parameters µ rf and T geo to the insignificant category.Yet, since the two parameters were assessed as significant in the OVAT procedure, they are not excluded from further investigations.
It was already clear from the parameter definition that, due to some value constellations of z aqf and t aqf , it cannot be ruled out that parts of the layer with groundwater flow lie below the BHEs.In the results of the Morris method, the parameter for the depth z aqf , where the layer with groundwater flow begins, is more influential than the thick- ness t aqf .Therefore, it can be assumed that the influence of z aqf is mainly caused by a reduction of the part of t aqf that lies in the area of the BHEs.To gain a better under- standing regarding the mutual influence of z aqf and t aqf , a separate investigation was first carried out with OGS simulations.In doing so, different constellations of these two parameters were considered, while all other parameters were set to the reference values specified in Table 2. Results of this investigation are shown with respect to the mean fluid temperature T 365 in Fig. 10.
The dashed parts of the curves indicate the curve section where the sum of z aqf and t aqf is greater than L BHE * and thus parts of the layer with groundwater flow lie below the BHE (see also Fig. 7).In these cases, an increase in one of the two parameters z aqf or t aqf only results in an increase in the thickness of the layer with groundwater flow below the BHE, while the fraction of the BHE that is effectively surrounded by groundwater flow remains the same.The almost horizontal course of these curve sections shows that the influence on the BHE hardly changes.Comparing the curve for z aqf = 0 L BHE * with the curve for z aqf = L BHE * , it becomes obvious that the main influence of z aqf on T 365 results from the reduction of the effective fraction of t aqf in the area of the BHEs.Ultimately, with this modeling approach, both the influences of the effective fraction of t aqf and z aqf cannot be clearly separated from each other.Similarly, the influence of t aqf is only par- tially related to the BHE.Therefore, the parameter z aqf was also excluded for further investigations.So in the following investigations, only the thickness t aqf of the layer with groundwater flow is varied.In the model, this layer always starts at a depth of -15 m and extends maximally to the bottom tip of the BHE at t aqf = L BHE * .Consequently, this parameter now describes the share of the BHE that is surrounded by groundwater.Finally, the 11 parameters from Table 3 remain for further investigations.

Results of global sensitivity analysis
With this reduced number of parameters, the sensitivity analysis is carried out with Sobol' indices.As explained, the parameter sampling points in this study are computed with a proxy model trained with results from OGS calculations.For the training, 2000 variations of parameter value constellations were calculated using OGS.For this purpose, parameter constellations were created using a LHS.To check the quality of the approximation by the proxy model, 208 additional OGS models were computed with parameter sampling points from a Sobol' sequence, which were not part of the proxy training.Subsequently, these parameter sampling points were also evaluated with the proxy model to compare the proxy results with the OGS results.The coefficient of determination ( R 2 ) and the root mean square error (RMSE) were used for this purpose.R 2 ranges between 0 and 1.Values close to 0 represent a larger deviation of the proxy output from the actual simulated output temperatures and values close to 1 indicate more accurate proxy outputs.For the RMSE, this is the other way around: the higher the value, the bigger the deviation from the OGS results.By creating proxy models trained with a subset of the 2000 sampling points from OGS, the convergence of the proxy model quality as a function of the training sample size could also be checked in this step.The corresponding results are presented in Fig. 11.Especially the curves for the coefficient of determination demonstrate the convergence with increasing number of training sampling points.With a number of 2000 training sampling points, the coefficient of determination is 99.8%, while the RMSE value has stabilized at around 0.5 K.
In this study, the Sobol' indices are calculated based on parameter variations that were generated using a Sobol' sequence.For the computation of first-and secondorder Sobol' indices, using the Sobol' Sequence, the sample size n is determined by the following relation (Saltelli 2002): Here j is the number of parameters and N the number of sampling points, where N should be a power of two according to Owen (2021).Due to the low computational cost of the proxy model, Sobol' indices based on different sample sizes could be considered in a convergence analysis.This analysis was carried out for N = 2 6 to N = 2 14 .Figure 12 shows Sobol' indices of the 11 parameters considered in dependence on the sample size.
Based on the shown convergence of the Sobol' indices, a Sobol' sequence with N = 2 14 was used to obtain the sensitivity analysis results.For 11 parameters, a Sobol' sequence with N = 2 14 corresponds to 393,216 parameter constellations within the value ranges, which were evaluated using the proxy model.In addition to the value ranges, the distribution functions from Table 2 were also taken into account.The greatest influence is exerted by the temperature conditions T geo at the site.The second largest influence is the extracted heat energy Q.Third largest is the parameter t aqf for the share of the BHE surrounded by groundwater flow.The influence of the Darcy velocity v f is slightly less than half that of the thickness t aqf .In contrast, the specific volumetric heat capacity of the geological subsurface c vs , the geothermal gradient T geo and the viscosity of the BHE's refrigerant µ rf show little to hardly any influence compared to the other parameters considered here.Considering confidence intervals, the thermal conductivity of the geological subsurface s and the Darcy velocity v f show a comparable influence.
Interactions between the parameters and their influences on the quantity of interest can be observed by means of second-order Sobol' indices.These are shown in Fig. 14.
The greatest influence due to parameter interaction is revealed for s and t aqf .In addition, also the second-order Sobol' indices for the thermal conductivity s in conjunction with Darcy velocity v f and the annual heat demand Q in conjunction with parameter t aqf stand out.
For the two parameter combinations with s , the interactions were investigated in more detail.As for the results shown in Fig. 10, parameter variations were calculated with OGS models in these investigations.While the two respective parameters are varied, all other parameters were set to reference values specified in Table 2. Figure 15 displays results for the two parameter combinations, which have the greatest Sobol' indices in Fig. 14.
The results for the parameter combination s and t aqf are shown in Fig. 15a.In the graph it can be seen that for larger values of s the curves have a flatter shape and for larger values of t aqf the distance between the curves becomes smaller.These proper- ties indicate that the influence of both parameters t aqf and s decreases as the value of the other parameter increases.Figure 15b illustrates the interaction between s and v f .The results show that at v f = 0 m s the largest temperature changes due to variations in s are obtained.At larger values of v f the curves become flatter, similar to what can be seen in Fig. 15a in relationship with t aqf .Also, it can be seen that for larger values of s , the curves for different values of v f are closer together.In the models, the groundwater flow is characterized by the parameters t aqf and v f .From first-order Sobol' indices in Fig. 13, it was seen that in this study the Darcy velocity v f has only slightly less than half the influence compared to the share of the BHE surrounded by groundwater t aqf .Figure 14 shows a smaller second-order Sobol' index for the combination of these two parameters than for the previously examined parameter combinations.Nevertheless, the second-order Sobol' index for t aqf and v f stands out from the majority of the other parameters.The relationship between these two parameters was also examined in more detail and is shown in Fig. 16.
As expected from the modeling, there is only an influence from groundwater flow if both parameter values are greater than zero.To create the graph, the parameter values of

Uncertainty quantification
Uncertainties in the designing process of geothermal systems essentially manifest in the fact that parameter values cannot be determined exactly.Instead, they usually lie within a certain interval-the uncertainty range.Examples of this are geological parameters, which are usually only investigated at discrete points or have to be estimated based on existing data collected at some distance from the site.If one wants to take these uncertainties into account, for instance, in the context of a planning process, this can be done by an uncertainty quantification.Uncertainty quantifications serve to analyze the sum of uncertainties due to the input variables on a result quantity.The aim of the investigation presented in this section is to show by example how uncertainty quantification can be carried out in the designing process of a geothermal system.This also includes elaborating the resulting uncertainties exemplarily.The uncertainty quantification builds on the results of the sensitivity analysis presented earlier and focuses on the uncertainty of the parameters that have shown the greatest impact on BHE performance.
As described in the introduction, an uncertainty quantification for BHEs in shallow subsurface is primarily of interest in terms of improving the BHE design.The uncertainty is assessed based on the BHE fluid temperature, as in the sensitivity analysis.An interesting aspect of the evaluation of the uncertainty quantification is whether the required limits are still met when the uncertainties are taken into account.For the designing of BHEs, the required values to be complied with often refer to lower temperature limits.During operation of a BHE, the annual course of the fluid temperature as well as of the subsurface decreases the most during the first few years, then only slightly to finally reach a kind of static state.For the uncertainty quantification, the temperature at the end of the operating lifetime is therefore most interesting.According to VDI 4640 (2020), a German guideline for shallow GSHP systems, an operating lifetime of 50 years should normally be assumed.To reduce the high computational effort required for numerical simulations over a period of this magnitude, only the first 5 years were considered for this study.This way, the largest changes in temperature could be taken into account while limiting the computational effort.In the fifth year, at a certain time the temperature of the fluid entering the BHE reaches a minimum T min in .The uncertainty is determined based on this lowest fluid temperature T min in .Due to the annually repeating heating curve, this is also the lowest fluid temperature occurring in the entire simulation period.This is always reached at the same time in the OGS model, regardless of the variation in parameter values.

Modeling uncertainty ranges
The uncertainty quantification is based on the preliminary design of a BHE for a fictional site using the software Earth Energy Designer (EED) in accordance with common design practice.An annual heat demand of 11 MWh was assumed, which is to be covered by a double U-tube BHE.The site was modeled oriented to the conditions in the region of Leipzig.The design with EED, assuming a sandstone-dominated subsurface, resulted in a BHE length of 100 m.The setup of the FE model used for the uncertainty quantification is basically the same as for the sensitivity analysis, but with different parameter settings.Therefore, in contrast to the calculation in EED, this model has a layer with groundwater flow.Furthermore, the temperature profile corresponding to Fig. 6 was again used as the initial condition for the temperature in the FE model.In EED, the subsurface temperature was 10.6 • C , whereas the mean subsurface temperature in the model domain of the OGS model is 9.9 • C , because of the more complex initial condition used in the FE- model.The temperature at the surface is defined in the FE model according to the curve shown in Fig. 17, while in EED it is set to constant 8.8 • C.
Figure 18 shows a comparison, of temperature curves with the values of the original design for five years.Figure 19 additionally shows the entire curve of mean fluid temperature calculated with EED.
Based on this design, uncertainty ranges were defined for the 11 parameters, which were already considered in the global sensitivity analysis.In addition, the parameter z aqf for the depth at which the layer with groundwater flow is located was included.Due to smaller ranges of values (uncertainty ranges, respectively) for t aqf and z aqf , no conflicts arise in the modeling, as they did in the sensitivity analysis.An uncertainty range was defined for each parameter.This is the value range containing the true parameter value, but which cannot be determined more precisely, or only with greater effort.The ranges are shown in Table 4.The values applied in the designing with EED are marked as reference values.
Parameter names are the same as presented for the sensitivity analysis in Table 1.The uncertainty range for the BHE length L BHE is an assumption which can be regarded as a measurement deviation.Since the BHE pipes are usually supplied by the manufacturer in the correct length, only minor measurement inaccuracies can be assumed here.
For grout material, the use of thermally enhanced materials was assumed.The associated thermal conductivity is given by gr = 2 W m −1 K −1 to 2.2 W m −1 K −1 (Stober and Bucher 2014).The lower value for gr represents the case of faulty backfilling of the BHE.Errors in the grout zone can be caused on the one hand by an improper execution during installation.Stober and Bucher (2014) for example write of errors due to use of too high water concentration when mixing the suspension in favor of better ease of pumping.On the other hand, errors can result from conditions in the subsurface, for example, if part of the grout is washed away or also from freezing in the grout zone.These faults have the effect of poorer heat transfer in the grout zone.In the uncertainty modeling, it was therefore assumed that up to 15% of the grout zone of the BHE is not filled with grout material or that there is no direct contact between rock and grout material.It was assumed that 10% of the faulted area is filled with groundwater and 5% is filled with air or air-like gas mixtures.Using a simple interpolation of the thermal conductivity over the volume fractions, the error was estimated in terms of a lower thermal conductivity gr .The uncertainty range used for the viscosity µ rf is based on the manufacturer's data for a commercially available BHE refrigerant consisting of a water-glycol mixture.Here, the uncertainty refers to viscosity values that can arise due to the varying fluid temperature during BHE operation.As a lower limit, a temperature of − 5 °C was chosen here.According to VDI 4640 (2020), this temperature is the lower limit.The upper limit temperature was set at 20 °C, which should not be exceeded as the upper limit for shallow geothermal systems according to VDI 4640 (2020).
For this study, the flow rate Vrf was assigned an assumed uncertainty of ±5 % .Fluc- tuations are conceivable here, for example, due to the temperature-dependent viscosity change of the refrigerant mentioned above or also when using a modulating heat pump.
Regarding the annual demand for thermal energy Q, two scenarios were considered.In the first scenario, hereinafter referred to as scenario 1, an uncertainty of the demand due to weather fluctuations was considered.Here, heating days from the years 2011 to 2021 for Leipzig were taken into account, which are published by the

IWU (Institut Wohnen und Umwelt GmbH 2021
).In the long-term average (2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021), the number of heating days was 252.7 days.For uncertainty modeling, the maximum and minimum percentage deviation of heating days from the mean number of heating days was considered and applied to the annual heat demand.The variation here was − 12.1% and + 6.2%.In the second scenario, hereafter referred to as scenario 2, greater uncertainty was assumed, which could result from a change in the use of the facility, for example due to a change of users.For this, a possible change of ±33.33 % was assumed.This is based on the specification of 1200-2400 annual full load hours in the VDI 4640 (2020), which is given there as a reference for parameter design of geothermal systems up to 30 kW heat load.
Uncertainty ranges for thermal conductivity s and specific volumetric heat capacity c vs in the geological subsurface are also based on data in VDI 4640 (2020).The char- acteristic values are already given there in the form of value ranges.If in the planning phase of a geothermal system there is only knowledge about which rock types are present at the site, it is the responsibility of the planning engineer to select a value from the value ranges given in the literature.Therefore, the ranges of values specified there were assumed to be the uncertainty ranges of this analysis.
With respect to groundwater, the depth and the thickness of the aquifer layers through which groundwater flows can be determined from contour maps of water table and the geological layer structure.Groundwater levels are values that are interpolated between measuring points.The thickness of the groundwater body can be estimated from the location of the nearest subjacent aquitard, which forms the groundwater bottom.This location of the corresponding aquifer is modeled here as a constant depth, even though this is an idealized simplification of the practice (Hölting and Coldewey 2013).Therefore, the thickness cannot be precisely determined with reasonable effort for the entire area under consideration.In addition to uncertainty due to interpolation between groundwater monitoring wells, groundwater levels may also change over the years, altering the location and thickness of the layer with groundwater flow.These sources of uncertainty are considered here with an uncertainty range of ±5 m .On the other hand, the uncer- tainty for parameter z aqf , which defines the depth of the layer with groundwater flow in the model, was estimated based on data on the change in groundwater levels from 1991 to 2017 in the region of Leipzig (Stadt Leipzig Amt für Umweltschutz 2017).Besides the location and thickness of the layer with groundwater flow, the direction and velocity of the groundwater flow at the site are also determined based on groundwater levels.For the simulation in OGS, the groundwater velocity is defined using the Darcy velocity.It is calculated according to the following equation: where k f represents the coefficient of hydraulic conductivity and i denotes the ground- water gradient.In practice, the groundwater gradient can be determined using groundwater levels at groundwater measuring points in the respective area.For the site, the groundwater gradient was determined from a groundwater level map of Leipzig with i = 1.3 • 10 −3 m m .Ranges of values for the permeability coefficient k f are given in litera- ture.These value ranges extend over several powers of ten and are only roughly categorized into rock types.The uncertainty range for Darcy velocity was therefore defined (3) based on the value range for the permeability coefficient in Hölting and Coldewey (2013) at constant groundwater gradient.Here, the values for k f range between 10 −4 and 10 −8 .
In the model, the initial temperature distribution is described by the temperature T geo and the geothermal gradient T geo .Uncertainties for these two parameters were assumed to be ±0.5 K for T geo and ±0.5 K m −1 for T geo based on experience.The param- eter T geo basically defines the undisturbed subsurface temperature.For larger BHE sys- tems, this is often determined together with the geothermal gradient as part of a TRT on a pilot borehole.Inaccuracies can arise if, for example, the measuring is carried out to soon after the test drilling (Rumohr 2021a, b).

Results of uncertainty quantification
To evaluate as large a number of parameter variations as possible in a reasonable amount of time, a proxy model is again built as previously described.In accordance with the approach of the sensitivity analysis, results from simulations with OGS are used to train the proxy model.For this purpose, 1000 parameter constellations were calculated in OGS specified with an LHS.The proxy model was used to increase the sample size of parameter variations up to 100,000.The quality of the proxy model was checked in the same way as for the sensitivity analysis.The R 2 value and the RMSE value were checked based on 100 sample points not included in the training.The results for scenario 1 are shown in Fig. 20.
With a number of 1000 sample points for training, the coefficient of determination is 99.99%, while the RMSE value has stabilized at around 0.02 K. Compared to the results for the quality of the proxy model in the sensitivity analysis (Fig. 11), the RMSE values of this proxy model are smaller and the R 2 values are higher.
To evaluate the results of uncertainty quantification, the previously mentioned cumulative distribution function p UQ (T min in ) of the minimum fluid temperature at the inlet to the BHEs in the fifth year of the simulation T min in is used here.The graph shows the distribution of all temperature values, which may result from the uncertainty ranges of the input parameters.Based on the graph, it can therefore be seen what percentage of cases the mean fluid temperature is below or above a given temperature value.Figure 21 shows the results for both scenarios (see also Table 4).Also, temperature T EED ref is marked in the graphs.This is the minimum fluid temperature entering the BHE in the fifth year, which was calculated in the original designing with EED.In this software, only a temperature curve for the mean fluid temperature can be determined.Therefore, the required supply temperature was formed by assuming a temperature difference of 3 K between the in-and outflow in the BHE.Furthermore, the minimum fluid temperature occurs in EED somewhat earlier in time than in OGS (see Fig. 18).Therefore, unlike the simulations in OGS, this minimum temperature was evaluated at that earlier time.In the graphs, the temperature T OGS Nevertheless, in both scenarios there is a proportion that lies below the reference temperatures.Regarding the initial design with EED, the uncertainty for temperatures lower than T EED ref is 4.6% in scenario 1. Referring to the OGS simulation with the BHE design parameter values, the uncertainty for lower values than T OGS ref is 15.0 %.The uncertainty for temperatures T min in ≤ 0 • C , in comparison, is 0.3 %.The lowest tem- perature for T min in is − 1.39 • C and the largest is 5.43 • C .In scenario 2, a larger uncer- tainty range for the annual heat demand Q was considered.Therefore, in Fig. 21b with the results for this scenario, a larger temperature range can be seen for the distribution of T min in .The smallest value here is − 3.77 • C and the maximum temperature for T min in in this observation is 6.34 • C .For lower fluid temperatures than in EED, the Scenario 1. Scenario 2.
Fig. 21 Results of uncertainty quantification using cumulative distribution functions p UQ (T min in ) with respect to the scenarios for the uncertainty in the annual heat demand.The uncertainty p UQ indicates the probability that T min in is below the corresponding temperature value uncertainty in this scenario is 13.7 %.Compared to T OGS ref , the uncertainty for lower temperatures in this case is 24.3 %.
In Germany, the design of BHE parameters is usually based on the guideline VDI 4640 (2020).It is specified there that for geothermal systems, the fluid temperature at the inlet of a BHE may reach −5 • C for a short period of time and should not fall below 0 • C as a monthly average.Applying this standard to the results shown in Fig. 21, it can be concluded that the design of BHE parameters is sufficient with respect to peak load in the first five years of BHE operation.However, no statement can be made about the monthly average temperature based on the evaluation of T min in .In the annual pattern of thermal energy demand applied here, the greatest demand and thus the lowest monthly mean temperature are in February.Therefore, additional proxy models were trained using the existing 1000 OGS calculations with the monthly mean fluid temperature at BHE inlet for the month February.Subsequently, these proxy models were also evaluated with the previously used LHS 100,000 constellations of values for the parameters.Using the cumulative distribution of the result values for the average fluid temperature in February, the uncertainty p UQ for a monthly mean temperature falling below 0 • C was determined.Figure 22 shows these uncer- tainty values for the first 5 years in both scenarios considered.
For scenario 1, the uncertainty is 0 % in the first year and remains at 0.6 % in the fifth year.For scenario 2, on the other hand, uncertainty reaches 1.3 % in the first and 6.9 % in the fifth year.
In addition to evaluating the uncertainty by looking at statistical distributions, Sobol' indices were calculated to evaluate sources of uncertainty.This consideration was again based on the minimum fluid temperature T min in entering the BHE in the fifth year of the simulation.For this purpose, the generated proxy model was used to evaluate 106,496 parameter variations based on a Sobol' sequence, and first-order Sobol' indices were calculated.The results are shown in Fig. 23.
From the graph it can be seen that the uncertainty at the considered site for both scenarios is mainly determined by the parameters Q, s and T geo .For scenario 1, the largest source of uncertainty in the considered case is the thermal conductivity of the geological subsurface s , whereas in scenario 2 the annual demand for thermal energy Q dominates as a source of uncertainty.

Discussion
The following section critically reviews the results presented earlier and discusses possible reasons and conclusions for the results obtained.Since the sensitivity and uncertainty quantification serve different purposes, this is done one after the other for the two studies despite their similarity.First, the results of the sensitivity analysis are discussed, followed by those of the uncertainty quantification.

Sensitivity analysis
In the sensitivity analysis, the relevance of 20 parameters for the design of BHEs in shallow geothermal systems is investigated.As a result, the undisturbed subsurface temperature T geo was identified as the parameter with the greatest influence on the performance of a BHE.This parameter represents the temperature conditions in the geological model, which are basically due to the undisturbed subsurface temperature.As a consequence of this result, special attention must be paid to the determination of temperature conditions at the site.In addition to careful exploration of undisturbed subsurface temperatures, the seasonal influences at the site of the planned BHE should also be appropriately considered to improve the design process.Furthermore, the temperature conditions in the subsurface are described by the geothermal gradient T geo .Compared to the other parameters of the analysis, however, this one showed a subordinate influence on the mean fluid temperature in the BHE.This is probably mainly due to the small depth of the shallow systems considered in this study, where the geothermal gradient can hardly exert its effect.
In the study, the annual heat demand Q had the second largest influence on the operation of the BHEs.For this reason, the demand should also be determined as precisely as possible from the perspective of facility systems and, if necessary, appropriate reserves should be planned.Furthermore, peak loads should also be considered for the design of a specific system, since they are relevant especially with respect to the temperature limits.In this study, these were not considered because their magnitude and occurrence are often highly stochastic and difficult to generalize for a global sensitivity analysis.For future considerations, cooling demand is also an important point in this context.Particularly due to climate change, this is becoming increasingly relevant.In the context of this study, the use of geothermal systems for cooling was not taken into account.In principle, however, a positive influence of cooling can Fig.23 First-order Sobol' indices for the parameters of the uncertainty quantification be assumed to a certain extent, since it has a regenerating effect on the temperature resources in the subsurface surrounding the BHE.If there are energy requirements for cooling and heating, there is a compensating effect in the sensitivity analysis with regard to the parameter influence due to the opposing demands.Thus, a lower influence of the parameter Q can be expected.Taking cooling into account, a change in the importance of parameters characterizing heat transport and heat storage in the geological subsurface would also be conceivable.For such considerations, other quantities may be of importance for the evaluation of parameter influences, such as integral values representing the fluid temperature over the entire operating time considered.
The parameter t aqf for the thickness of the layer with groundwater flow or, more precisely, the share of the BHE which is flowed around by groundwater showed the third largest influence on the performance of the BHE.In comparison, the Darcy velocity of the groundwater v f only had an influence that was slightly less than half of t aqf .For the design of a BHE, this suggests that, with respect to the performance of the BHE, special attention should be paid to identifying the areas where groundwater flows around the BHE.The Darcy velocity of the groundwater flow was less crucial.
In an isolated analysis it could be shown that Darcy velocities in the lower part of the considered value range are already sufficient to obtain a positive influence on the BHE.Increasing the Darcy velocity above a certain value did not show any significant improvement of the positive influence.With respect to these results, Darcy velocities in the lower range of values ( v f = 1 • 10 −8 m s to 3 • 10 −6 m s ) should be determined more precisely.In particular for larger geothermal systems, such as for quarters, it is often necessary to consider the propagation of temperature plumes during the planning process.This serves the purpose of being able to estimate the influence of neighboring BHEs and to be able to comply with limits for temperature changes around the BHE due to applicable regulations.In the case of such larger geothermal systems, the determination of flow direction and Darcy velocity has great importance in the design process.
When considering the performance, the influence of the thermal conductivity s is comparable to v f .An isolated comparison of these two parameters showed that in the case of a large value for one parameter, the influence of the other parameter gets smaller.An analogous behavior was also observed for the parameters s and t aqf .This behavior is due to the fact that all three parameters characterize a form of heat transport in the subsurface and, therefore, result in similar effects with respect to the performance of the BHE.
Comparing the 11 parameters from Sobol' sensitivity analysis, especially specific volumetric heat capacity c vs , dynamic viscosity of refrigerant µ rf and geothermal gradient T geo showed little influence.However, when considering the results pre- sented, it should be noted that the 11 parameters in the analysis with Sobol' indices are already the parameters with the greatest influence.Ultimately, all of these parameters are relevant to the design of a BHE.The results are not meant to be an encouragement to neglect all parameters presented as having little influence in the context of the sensitivity analysis.For example, parameters such as the thermal conductivity of the grout material gr , which in comparison is rather in the lower middle range of the Sobol' indices, should not be disregarded.The choice of thermally improved grout materials with higher thermal conductivities is definitely recommended as well as careful execution of the backfilling under no circumstances should be underestimated.
It should be noted that the results of a sensitivity analysis are only valid in the context of the model assumptions made and only for the defined value ranges.The parameter value ranges used in this study are intended to represent the region of Germany.The conditions for the studies were selected in such a way that the highest possible general validity could be achieved.
With the selected simulation duration, it was also possible to capture the change in parameter influences over an annual cycle.However, it is quite conceivable that the sensitivity of some parameters may still change if more than 1 year is considered.Interesting results in further research could therefore be the consideration of the change of parameter sensitivities in the life cycle of a BHE.

Uncertainty quantification
In contrast to the more general consideration in the sensitivity analysis, the uncertainty quantification in this paper is a site-specific investigation.Therefore, other sites may have different uncertainties.The site considered here is a fictitious site for an exemplary consideration.However, the result of this uncertainty quantification provides interesting insights for planning practice and may provide important hints for handling uncertainty considerations for shallow geothermal systems.
For practical application, the consideration of uncertainties regarding the regulations to be complied with is of particular interest.In this work, temperature limits given in the German guideline VDI 4640 (2020) were compared to the results of the uncertainty quantification.This showed that the limits for the fluid temperature at peak load are complied with in both scenarios considered, at least within the first 5 years.However, further examination of the mean monthly fluid temperature at inlet of the BHE showed that it reaches temperatures below 0 °C regardless.In the fifth year of the simulation, in February, the monthly mean fluid temperature at the inlet of the BHE showed an uncertainty of 0.6% in scenario 1 and 6.9% in scenario 2 for values below 0 °C.These values for the uncertainties can be directly interpreted as the respective probability of undersizing the BHE in the particular scenario.With regard to compliance with the temperature limits, however, it should be noted here again that only the first 5 years of BHE operation were considered or, more precisely, only 4 full years, as shown in Fig. 18.VDI 4640 (2020), on the other hand, specifies an operating lifetime of 50 years, which is generally assumed.When considering the total operating life of a geothermal system of 50 years, lower fluid temperatures and therefore larger uncertainties regarding undersizing of the BHE can be expected.However, the largest changes in annual fluid temperature are expected to occur within the first few years, as seen in Fig. 19 of this paper.
Besides the uncertainty regarding too low temperatures, the results of the uncertainty quantification in the example show a high probability for higher fluid temperatures than expected from the original design.If there is still a large probability of higher fluid temperatures when the entire operating time is taken into account, this can be interpreted as a large probability of oversizing the system at the location under consideration.
In addition to a detailed investigation of the most influential parameters from the sensitivity analysis, oversizing must be addressed primarily by reducing the uncertainty ranges of the input parameters.For this purpose, it is important to understand which parameters contribute most to the uncertainty of the result.For the example site, this was determined by means of Sobol' indices.The results showed that the uncertainty for both scenarios is mainly due to four parameters.These are the annual demand for thermal energy Q, the thermal conductivity of the geological subsurface s , the Darcy velocity v f of the groundwater flow and parameter T geo , which essentially describe the undisturbed subsurface temperature.
In scenario 1, the largest uncertainty in the performance of the BHE is caused by the thermal conductivity of the geological subsurface s .The defined uncertainty range for this parameter is based on the assumption that no measured value is available at the site and therefore a value from literature (here VDI 4640 ( 2020)) had to be selected with respect to the rock type.TRTs offer a possibility to determine the thermal conductivity more precisely.Geothermal systems often require multiple BHEs, so that such a test can be carried out on a pilot borehole.In the example site, however, only one BHE is planned, so a TRT can only check the value for s after the installation of the BHE.In the case of only one borehole, it is worthwhile to investigate carefully in the planning process whether, for example, a more accurate estimate of the value can be obtained using geoinformation systems based on existing data from boreholes in the surrounding area.On the other hand, TRTs themselves can also be a source of uncertainty (Richter et al. 2023a, b).This is mainly due to the fact that the starting point for the evaluation procedure must be chosen by the user and that an averaged value for the thermal conductivity is determined.
For scenario 2, the annual demand for thermal energy Q constitutes the largest source of uncertainty.In principle, the uncertainty of the energy demand can be reduced with a demand analysis that is as detailed as possible, including for example a simulation of the facility systems.However, in the scenario presented here, we consider an uncertainty that may arise due to a change in use.This is difficult to predict during the design of a system and is not reasonable from an economic point of view in most cases.Therefore, it might be more appropriate to address such uncertainties by reviewing the original design of the system in the event of a change of user behavior and, if necessary, to compensate for additional needs by expanding the system.
As mentioned before, the uncertainty quantification showed a high probability for higher fluid temperatures than expected from the original design with EED.Nevertheless, the probability of falling below reference temperatures is significant in both scenarios, which means that uncertainty is still not negligible.For scenario 1, the uncertainty for lower temperatures was about 4.6% and for scenario 2 about 13.7%.In the simulation with OGS on which the uncertainty quantification results are based, the modeling approach for subsurface temperature resulted in a lower mean temperature of the subsurface in the model compared to the EED design.Therefore, the reported uncertainties for lower temperatures are a rather conservative estimation in this case.On the other hand, the conservative nature of EED designing is nevertheless evident from the results.For example, the simulation with the parameter values of the EED design in OGS resulted in a higher fluid temperature than in EED despite the lower subsurface temperature.This is mainly due to the fact that groundwater flow at the site cannot be taken into account in the calculation with the EED software.
As mentioned above, the presented uncertainty quantification is intended to give an exemplary insight into its possibilities and methodology.No claim is made to have included all sources of uncertainty.Due to software-related problems, for example, it was not possible to take into account uncertainties resulting from the position of the individual pipes within the BHE.The position of the pipes is usually assumed to be ideal in planning processes but is difficult to control or predict in reality when installing a BHE.There are also uncertainties due to the lack of knowledge of the exact layer structure in the heterogeneous geological subsurface.Further investigations could therefore consider more detailed information on the lithology and also the permeability of the layers in which groundwater flows.Oftentimes, information from other boreholes is only available at greater distance from the proposed facility.The exact path of the BHE in the subsurface is also subject to uncertainties, which Steinbach et al. (2021), for example, considers in an uncertainty quantification.Furthermore, it is likely that uncertainties in the parameters also affect other areas that were not considered.For instance, the evaluation by means of Sobol' indices shows that the viscosity of the BHE fluid has only a small influence on the performance.However, the actual value of viscosity Vrf lies in a large uncertainty range (see also Table 4).Sto- ber and Bucher (2014) notes about the large change in viscosity during BHE operation that it has an impact on the required electric pump current, thus additionally affecting the performance of the geothermal system.In the case of poor design due to uncertainties, impacts on groundwater are also conceivable, as considered in other studies (Griebler et al. 2016;Zhu et al. 2017;Vienken et al. 2019).However, these were not the subject of the investigation here, as the focus is on the performance of the BHE.
The pursued approach for uncertainty quantification represents a good method to obtain a differentiated overview of the uncertainties.It can be noted here that compared to the sensitivity analysis, in the uncertainty quantification the proxy model was able to achieve a better fit to the OGS results even with a smaller sample size for training.The reason for this is most likely that the uncertainty ranges have a smaller span than the value ranges used in the sensitivity analysis.However, this uncertainty quantification is not yet suitable for use in planning processes in practice due to the high time requirements of numerical simulations despite the use of machine learning.Further research should therefore address how to achieve comparable results in a reasonable time for planning practice.In addition, the computation time should be short enough to allow simulations of 50 years to be considered in order to make statements about the uncertainty for the entire operating life of the geothermal system.One obvious way to address this is to reduce the size of the OGS model.The basic model used here was originally designed for sensitivity analysis, where larger temperature plumes may occur due to the larger parameter value ranges and thus larger model size was required.Therefore a smaller model size would have been possible.Other ideas include the use of Monte Carlo simulations, such as those used by Watanabe et al. (2010) for example, or the reduction of the considered parameters.As mentioned before, the result of this analysis already showed that the uncertainty is mainly due to a few parameters.

Conclusion
In this paper, both a global sensitivity analysis for 20 design parameters of shallow geothermal systems and an uncertainty quantification for a fictitious example site, considering twelve parameters, were performed.A main motivation of this study was to investigate parameter sensitivities for design practice and how uncertainty analyses as a tool could be beneficial for practitioners.The aim of designing a BHE is to achieve the best possible performance.Therefore, both analyses consider the influence on performance measured by changes in fluid temperature, in this study.The analyses are mainly based on numerical simulations with OGS in combination with machine learning using Gaussian process regression.Using this machine learning method, so-called proxy models were trained, which in turn were used for the evaluation of a large number of parameter constellations.Although small losses in precision had to be conceded in this way, the computational effort could be reduced tremendously.To our knowledge, this is the first time that sensitivity analyses respectively uncertainty quantifications for the shallow geothermal field have been performed on such a large scale and based on such detailed numerical simulations.
The global sensitivity analysis is based on a careful preparation of value ranges, which represent the country of Germany as well as possible.By choosing a simulation period of 1 year, seasonal influences could also be fully considered and represented.To be able to examine the effects of the design parameters in isolation, a rather simple model was constructed.For the analysis the performance is measured by the mean fluid temperature at the end of the year.In a first investigation, sensitivities were evaluated for a parameter screening with two methods (OVAT and Morris method).The screening allowed to exclude nine parameters that showed less influence on the simulation results compared to all other considered parameters.This reduced the calculation effort for the actual global sensitivity analysis.For the evaluation of parameter constellations in the global sensitivity analysis with Sobol' indices, a proxy model was used, which was previously trained with 2000 OGS calculations.As a result, parameter sensitivities regarding the performance of a BHE could be identified, which provide important information for the design of BHEs in shallow geothermal systems in practice.The most important findings are summarized below: • The initial temperature conditions at the site showed the greatest influence on the BHE performance, which suggests that in practice special attention should be paid to the determination of the undisturbed subsurface temperature and, if possible, to detailed modeling in the designing process.• The second most important parameter is the heat demand.Thus this should be determined as precisely as possible in planning practice.
• In third rank is the share of the BHE that is surrounded by groundwater flow, so in practice greater attention should be paid to the layers in which groundwater flow occurs.• The size of the share of the BHE that is surrounded by groundwater flow has a greater influence than the Darcy velocity of the groundwater.• Darcy velocity and thermal conductivity showed a comparable influence.
In the presented uncertainty quantification, the uncertainty of the design of a BHE with respect to its performance for a fictitious site was examined.On the one hand, the aim was a methodical consideration of uncertainty quantification in the context of planning processes for shallow geothermal systems.On the other hand, the analysis should give an exemplary impression of how existing uncertainties regarding the site conditions and the system itself can affect the performance of a geothermal system.Two scenarios were considered, which differ in their uncertainty range for the annual heat demand.In scenario 1, an uncertainty of the heat demand based on weather fluctuations was incorporated.Scenario 2 takes into account a larger range of uncertainty for the heat demand, which might occur due to a change of user behavior.Effects of the uncertainty in the input parameters on the operation of the BHE were considered analyzing the lowest inlet fluid temperature in the fifth year of the simulation.The result of the uncertainty quantification is a statistical distribution that provides information about the probability that the fluid temperature is above or below certain values.The starting point for the uncertainty quantification was a design of the BHE for the site based on a calculation with the software EED, which is a commonly used tool for the design of such systems.Even though the uncertainty quantification is based on simplified model assumptions, it provides a good insight into the potential of such an analysis.On the one hand, this shows a profound view with regard to regulations to be complied with and, on the other hand, it provides information that can be used to reduce oversizing.The key findings from the uncertainty quantification are summarized below: • In both example scenarios, uncertainties could be found for falling below the fluid temperature expected according to the design (depending on the scenario and the reference temperature, the uncertainty for lower temperatures than expected lies between 4.6 % and 24.3 %).• After 5 years, the probability that the fluid temperature falls below 0 °C in the monthly mean of the month with the highest heat demand was 0.6 % in scenario 1 and 6.9 % in scenario 2. • The probability for oversizing the BHE is much larger than the uncertainty for undersizing, suggesting that uncertainty quantification could also be useful for minimizing drilling meters.• The comparison of the two scenarios showed that in this example the change in user behavior causes greater uncertainties than weather changes.• Sobol' indices can be used to identify the parameters that contribute most to uncertainty and can therefore help reduce oversizing.
• In the example, the Sobol' indices have shown that the uncertainty is mainly caused by four parameters (the thermal conductivity of the subsurface, the heat demand, the temperature conditions at the site and the Darcy velocity of the groundwater flow).
Overall, it can be concluded that the application of uncertainty quantification in planning practice can be very beneficial.To make this possible, future studies have to address the reduction of the computational effort of such analyses.

Fig. 1
Fig. 1 Workflow for the sensitivity analysis

Fig. 2
Fig. 2 Workflow for the uncertainty quantification

Fig. 3
Fig. 3 Schematic structure of the FE model according to the dual-continuum approach

Fig. 5 Fig. 6 Fig. 7
Fig. 5 Annual curves and value range of the surface temperature for different values of T geo

Fig. 8
Fig. 8 Change of mean fluid temperature after 1 year with respect to the calculation with reference values

Fig. 9
Fig. 9 Results of parameter screening with the Morris method

Fig. 10
Fig. 10 Influence of parameter t aqf on mean fluid temperature T 365 after 1 year of BHE operation for different values of z aqf shows these Sobol' indices (of first order), which were calculated from the (2) n = N (2j + 2).

Fig. 12
Fig. 12 Convergence of Sobol' indices depending on the sample size.Each curve represents the Sobol' index of 1 of the 11 parameters

Fig. 14
Fig. 14 Second-order Sobol' indices of parameters used in the global sensitivity analysis Fig. 15 Detailed examination of relationships between the parameter combinations with the greatest second-order Sobol' indices for the result quantity T 365

Fig. 16
Fig. 16 Detailed view of the relationships between parameters t aqf and v f for the influence on the result quantity T 365 .Solid lines represent 11 values of v f , which are evenly distributed over the value range.In addition, dashed lines represent calculated values of v f , as specified in the legend

Fig. 17
Fig. 17 Reference curve of surface temperature in the OGS model for one year

Fig. 19
Fig. 19 Mean fluid temperature of the BHE from the calculation with EED

Fig. 20
Fig. 20Quality of the proxy model for scenario 1 of the uncertainty quantification compared to OGS results using the coefficient of determination R 2 and the root mean square error RMSE for the quantity of interest T min in depending on the training sample size ref denotes the minimum fluid temperature T min in obtained from the simulation with OGS using the parameter values derived from designing with EED.Temperature T EED ref is equal to 0.72 • C and temperature T OGS ref is 1.55 • C .In the graph, the values on the ordi- nate show how the result values are distributed.The uncertainty p UQ indicates the probability that T min in is below the corresponding temperature value.At first glance, it is obvious that the probability of falling below the temperatures T EED ref and T OGS ref is smaller in both scenarios than the probability of exceeding them.

Fig. 22
Fig. 22 Uncertainty p UQ that the monthly average fluid temperature is below 0 • C in the month of February within the first 5 years

Table 1
Preselection of input parameters for sensitivity analysis

Table 2
Value ranges and statistical distributions of the input parameters for the sensitivity analysis *Relative value related to the length of the BHE reduced by 15 m L BHE * **Parameter refers to the change of the curve shown in Fig. 17 as well as the temperature profile of the initial condition for subsurface temperature

Table 3
Parameter selection for the sensitivity analysis with Sobol' indices Quality of the proxy model for sensitivity analysis compared to OGS results using the coefficient of determination R 2 and the root mean square error RMSE for the quantity of interest T 365 depending on the training sample size

Table 4
Uncertainty ranges and reference values for uncertainty quantification *Relative value related to the length of the BHE reduced by 15 m L BHE * **Parameter refers to the change of the curve shown in Fig. 17 as well as the temperature profile of the initial condition for subsurface temperature