Uncertainty quantification of borehole heat exchanger design length based on a global sensitivity analysis

. Ground heat exchangers are essential components of ground-source heat pump systems. One of the most-used ground heat exchanger types is a borehole heat exchanger (BHE). Many parameters, such as the building thermal loads, ground thermal properties, and BHE components, influence the required BHE length. Therefore, understanding how individual parameter uncertainties are propagated to the output is counterintuitive. In this study, a global sensitivity analysis (GSA) was performed using the framework of the American Society of Heating, Refrigerating, and Air Conditioning Engineers (ASHRAE) BHE design method. The GSA was conducted to identify the subparameters behind the input parameters of ASHRAE design method, which indirectly determine the BHE lengths by impacting the direct input parameters in the ASHRAE formula. The results reveal that the zone infiltration, heat pump efficiency, ground temperature, effective borehole thermal resistance, and electric equipment heat gain are the most influential parameters on the BHE design length.


Introduction
The ground heat exchanger is the most crucial component in ground-source heat pump systems. The vertical closed-loop ground heat exchanger, which is called a borehole heat exchanger (BHE), is the most commonly used type.
Various factors influence the BHE length, such as the building load, ground-source heat pump efficiency, building geometry, internal gains, microclimate, BHE geometry, and thermal properties of BHE components and ground. Additionally, many factors cannot be fully controlled in the BHE construction process. Moreover, the effective thermal properties of the ground, a porous media, have significant uncertainties. These uncertainties can cause wrong decisions regarding the BHE sizing and lead to an increase in initial cost or operation failure. Thus, the probabilistic approach is reasonable for estimating the required BHE length because the effect of the uncertain factors propagating to the output can be quantified.
Therefore, we performed a global sensitivity analysis (GSA) in the framework of designing the BHE with the American Society of Heating, Refrigerating, and Air Conditioning Engineers (ASHRAE) method [1]. The individual parameters in the ASHRAE design formula and their associated latent variables determining the individual input value were included in the GSA for a holistic understanding of the uncertainty propagation. The results of the GSA provide insight into which parameters are more influential in the BHE design and which uncertainties should be managed to improve design accuracy. * Corresponding author: wonjun.choi@jnu.ac.kr

Methodology
The ASHRAE design method needs the following parameters: thermal loads assigned to the ground, effective ground thermal resistances, mean temperature of the circulating fluid in BHE, undisturbed ground temperature, effective borehole thermal resistance, and penalty temperature. First, the GSAs were carried out to screen the subparameters that underlie the input parameters of building loads and effective ground thermal resistances to select the influential parameters on the ASHRAE design method. The final GSA was then conducted to identify the uncertainty of the BHE design length, incorporating the screened subparameters and all the parameters required in the ASHRAE design method.

ASHRAE design method
Kavanaugh et al. [1] provided the formula for the BHE design length for cooling and heating operations as expressed in Eq. (1). The BHE design length is primarily influenced by the peak load [2]. In this study, only the BHE design for heating season is considered because the target building has a much higher peak heating load than the peak cooling load. Information for the target building is described in Section 3.
There are four types of ground loads ( , , 4ℎ , and ), which are differentiated by the time scale. The building load is not directly assigned to the BHEs as the coefficient of performance (COP) of the heat pump should be considered to obtain the ground loads. Eq. (2) expresses the relationship between the ground load and the building load. Hence, the COP is considered an uncertain factor in the GSA.
The term represents the annual average load. The terms , 4ℎ , and denote the peak loads for 720, 4, and 1 h, respectively. The effective ground thermal resistances ,10 , , , and ,4ℎ correspond to loads of 10 years, 1 month, and 4 hours, respectively. In addition, * represents the effective thermal resistance of the borehole. Further, indicates the undisturbed ground temperature, represents the mean fluid temperature of the heat pump, and denotes a temperature penalty that considers the heat interference of adjacent boreholes in the borefield. In this paper, the penalty temperature is considered a negligible parameter because we assumed that each borehole is positioned sufficiently far from the others.
The g-function is a nondimensional thermal response function influenced by the ground thermal conductivity, response time, borehole configuration, and arrangement [3]. Any response models can be used for obtaining g-function. We used the finite line source (FLS) model [4,5]. The FLS model accounts for the axial heat transfer at the end of the BHE, which is crucial for a long-term ground response. There are various FLS model equations regarding boundary conditions and computational speed. This study used the FLS model implemented by Claesson

Global sensitivity analysis
The sensitivity analysis is a technique that determines the influence of changes in input parameters on the model results. Sensitivity analyses enable a designer to make reliable decisions by quantifying the effects of uncertain parameters.
The sensitivity analysis is divided into local sensitivity analysis (LSA) and global sensitivity analysis (GSA). LSA involves varying one input parameter at a time to explore its effect on the output, providing insights into the effect of the parameter on a specific location of the parameter space. However, it may be challenging to examine the overall impact of a parameter on the output, especially when the parameter space is large and the causality between the input and output is nonlinear. On the other hand, GSA is a method that explores the entire input space and simultaneously varies all input parameters to examine input parameters' influence on the output. Therefore, we used GSA in this study.
The GSA distinguishes the effects of input parameter variations on the model output. This study applied Sobol's method based on the variance decomposition technique to consider the interaction between input parameters [8].
in Eq. (6) is the first-order sensitivity index of the input parameter . The conditional expectation of output , given all input parameters except , is denoted by ~( | ). The term is the variance of . The term ( ) is the variance of the model output. A high indicates a significant effect of the variations in the output parameters.
In Eq. (7), the variable represents the total order index of input parameter . The conditional variance of output , given all input parameters except , is denoted by ( |~). In other words, ( |~) is the variance caused by the interaction with other input parameters when changes within its parameter space. The outer expectation value ~, which is the expectations of all possible input parameters except , was calculated for ( |~).
In summary, a first-order index means a single effect of each input parameter, and the total order index is the summation of the first-, second-, and higher-order indices. The total order index represents the effect of the entire input parameter space influenced by each input parameter.

Target building
For the load calculation, a small office building in Chicago, one of the commercial reference buildings provided by the US Department of Energy, was used [9]. The building loads were calculated using EnergyPlus and employed in the ASHRAE method. The building has one story and is divided into six zones: east, west, south, north, core, and attic. The five zones, excluding the attic, are air-conditioned. The geometric configuration of the building is depicted in Fig. 1. The floor area of the building is 511 m 2 . The energy simulation used weather data from the Chicago O'Hare International Airport TMY3 [10]. The heating, ventilation, and air conditioning system in the original model file was replaced with an ideal load air system with 100% efficiency to calculate the building loads. Although the annual cooling load (51.9 GJ) is higher than the annual heating load (45.9 GJ), the peak heating load (122 MJ) is about twice as high as the peak cooling load (61 MJ). As the size of the BHE is significantly influenced by the peak load, only the heating load is considered for the BHE design [2].

Input variable screening
Before conducting Sobol's method on the borehole design length, the input parameters of the building loads and ground thermal resistance were screened through a GSA. The purpose of variable screening is to conduct a sensitivity analysis using only parameters significantly affecting the BHE design length.

Input variables of the building load
This study considered 17 parameters related to the thermal load, such as the building envelope performance (thermal conductivity, material thickness, solar heat gain coefficient (SHGC), and window U-value), internal heat gains (occupancy, lighting, and electric equipment), and zone infiltration rate. All input parameters considered for the GSA of thermal loads are listed in Table 1. The uncertainty range and probability distribution were based on prior research, literature, and standards [11][12][13][14][15][16][17][18][19].
In case of an office, heat gain from lighting and electric equipment depends on the number of occupants, but it does not fully reflect the variability of occupancy. Especially, we can guess that heat gain from lighting occurs within a range where reduction factor applied to the range of occupancy, because lights are mostly turned on to use occupied areas. In addition, it is reasonable to suppose that the dependency of electric equipment heat gain on the number of occupants is less than the lighting. Hence, we made an assumption that the uncertainty range of lighting and eletric equipment heat gain would be set at 50%, 60% of uncertainty range of occupancy heat gain, respectively, although this assumption may be somewhat arbitrarily.
For the GSA of thermal loads, Monte Carlo simulations were conducted using the input parameters listed in Table 1. EnergyPlus was executed with 5,000 sample sets generated using the Saltelli sampling method [20]. The influence of each variable was verified using Sobol's method. The average percentages of firstorder indices for four types of loads were 61% for the zone infiltration (with the strongest influence), 5.5% for the thermal conductivity of the attic floor insulation, 5% for the occupancy heat gain, 4.9% for the window Uvalue, 4.9% for the thermal conductivity of concrete, 4.8% for the electric equipment heat gain, and 3.6% for the lighting heat gain. The screened influential parameters are illustrated in Fig. 2 and highlighted in bold in Table 1.

Input variables of the effectivsse ground thermal resistances
The effective ground thermal resistances of various time scales ( ,10 , , , ,4ℎ ) were calculated using Eqs. (3)-(5) [1]. The ground thermal conductivity and the time scales 1 , 2 , and 3 are 3,650, 3,680, and 3,680.167 days, respectively. Before quantifying the uncertainty range of effective ground thermal resistances, the range of uncertain parameters affecting the effective ground thermal resistance was set based on prior research [21] and listed in Table 2.
The 15,000 sample sets to quantify the uncertainty range of the effective ground thermal resistances were generated using the Saltelli sampling method. Then, Monte Carlo simulations were performed using Eqs. (3)-(5) and the FLS model. Sobol's method was applied to determine the contribution of each parameter to the effective ground thermal resistances.
The percentages of first-order indices, which indicate the direct contribution of each parameter to the output, are presented in Table 3. The influence of the ground volumetric heat capacity (59%) on ,4ℎ has the highest effect, followed by the borehole radius (34.7%) and ground thermal conductivity (6.3%). In the case of , , the influence of the ground thermal conductivity (99.6%) is greater than that of the other parameters. For ,10 , the ground thermal conductivity is 97.3%, followed by the ground volumetric heat capacity at 2.7%. The uncertainty of the active borehole length has almost no effect. The screened input parameters are highlighted in bold in Table 3. Table 2. Parameters for effective ground thermal resistances and their uncertainty ranges. Table 3. First-order sensitivity index of screened variables for ,4ℎ , , , and ,10 .

Results
The necessary input parameters required for calculating the BHE design length are presented in Table 4. The 16,000 samples were generated using the Saltelli sampling method for input parameters listed in Table 4. Sobol's method quantified the contribution of the subparameters of the ASHRAE method, including previously screened input parameters. The calculated 95% confidence interval of BHE design length ranges from 377.5 to 547.9 m, as depicted in Fig. 3. The single effect of input parameters on the BHE design length was also assessed. The zone infiltration (29.4%), COP of the heat pump (17.8%), ground temperature (15.7%), effective borehole thermal resistance (13.5%), and electric equipment heat gain (7%) were the most influential factors. Figure 4 illustrates the first-order effect. The higher-order effect is negligible in the ASHRAE method. The results suggest that the designer must consider the building load subparameters, such as the zone infiltration and electric equipment heat gain, which are latent variables of the ASHRAE method but have a significant impact on the uncertainty of the BHE design length, in addition to the ground thermal properties and effective borehole thermal resistance.

Conclusion
This study quantified the contribution of each parameter affecting the BHE design length using the Sobol' method in the framework of the ASHRAE design method. Furthermore, the GSA was conducted to quantify the contribution of subparameters, which are the determining inputs for the ASHRAE method. As a result, the following subparameters: zone infiltration, COP of the heat pump, ground temperature, effective borehole thermal resistance, and electric equipment heat gain, which are not explicitly expressed in the ASHRAE deisgn formula, have a dominant influence on the BHE design length. This result demonstrates that the parameters determining the building load should be considered in addition to the thermal properties of the ground and borehole when estimating the BHE design length. The uncertainty of the weather data was not considered in estimating the building loads. Therefore, a more comprehensive understanding of the uncertainty in the BHE design can be achieved through the GSA considering the uncertainty of climate factors in the future.