Applied Energy 208 (2017) 1113–1127

Contents lists available at ScienceDirect

Applied Energy
journal homepage: www.elsevier.com/locate/apenergy

A computationally eﬃcient pseudo-3D model for the numerical analysis ofborehole heat exchangers

MARK

⁎

Giuseppe Brunettia, , Hirotaka Saitob, Takeshi Saitoc, Jiří Šimůnekd
aDepartment of Civil Engineering, University of Calabria, Rende, CS 87036, Italy
Graduate School of Agriculture, Tokyo University of Agriculture and Technology, 3-5-8 Saiwaicho, Fuchu, Tokyo, Japanc
Graduate School of Science and Engineering, Saitama University, 255 Shimo-Okubo, Sakura-ku, Saitama, Japand
Department of Environmental Sciences, University of California, Riverside, CA 92521, USAb

H I G H L I G H T S
eﬃcient modeling framework is proposed.
• AThecomputationallyapproach includes a widely used hydrological model.
• The modelingis validated against experimental data from two Thermal Response Tests.
• The model
model is used to perform a statistical sensitivity analysis.
• The validated• inﬂuence of groundwater and lithologic heterogeneities is examined.

A R T I C L E I N F O

A B S T R A C T

Keywords:
Geothermic
Heat transfer
Thermal response test
Modeling
Sensitivity analysis

Ground-Source Heat Pump (GSHP) systems represent one of the most eﬃcient renewable energy technologies.
Their eﬃciency is highly inﬂuenced by the thermal properties of the ground, which are often measured in-situ
using the Thermal Response Tests (TRTs). While three-dimensional mechanistic models oﬀer signiﬁcant advantages over analytical solutions for the numerical interpretation of TRTs, their computational cost represents alimiting factor. Moreover, most of the existing models do not include a comprehensive description of hydrological processes, which have proven to strongly inﬂuence the behavior of GSHP. Thus, in this study, we propose
a computationally eﬃcient pseudo-3D model for the numerical analysis and interpretation of TRTs. The numerical approach combines a one-dimensional description of the heat transport in the buried tubes of the exchanger with a two-dimensional description of the heat transfer and water ﬂow in the surrounding subsurface
soil, thus reducing the dimensionality of the problem and the computational cost. The modeling framework
includes the widely used hydrological model, HYDRUS, which can simulate the movement of water, heat, andmultiple solutes in variably-saturated porous media. First, the proposed model is validated against experimental
data collected at two diﬀerent experimental sites in Japan, with satisfactory results. Then, it is combined with
the Morris method to carry out a sensitivity analysis of thermal properties. Finally, the model is exploited toinvestigate the inﬂuence of groundwater and lithologic heterogeneities on the thermal behavior of the GSHP.

1. Introduction

to one of the most recent reviews by Lund and Boyd [1], the total installed capacity of ground-source heat utilization systems increased
46.2% worldwide from 2010 to 2014 and more than 600% in the past
20 years. Among diﬀerent ground-source heat utilization systems,
GSHP systems account for more than 70% of the installed capacity
worldwide. The heat exchange capacity of the heat exchanger, which
determines the performance of the GSHP systems, is mainly characterized by ground thermal properties, such as the thermal conductivity and the thermal capacity. As the optimum design and

Ground-source heat pump (GSHP) systems use ground or groundwater as a heating/cooling source by circulating a non-freezing solution
inside of a vertically or horizontally installed closed-loop tube totransfer heat to/from the ground. Because GSHP systems are based on
renewable energy and can achieve a much higher coeﬃcient of performance (COP) than conventional air-source heat pump systems, theuse of GSHP systems has been rapidly increasing worldwide. According⁎

Corresponding author.
E-mail address: Giusep.bru@gmail.com (G. Brunetti).

http://dx.doi.org/10.1016/j.apenergy.2017.09.042
Received 24 April 2017; Received in revised form 31 July 2017; Accepted 10 September 2017
Available online 18 September 20170306-2619/ © 2017 Elsevier Ltd. All rights reserved.

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

estimated from the slope m of the plot. This estimate depends largely on
the selection of the temporal data section. By selecting the time period,
during which there is an obvious linear relationship, early time data arediscarded from the analysis. Also, since the volumetric heat capacity C(or the speciﬁc heat c) that characterizes the heat storage in a substanceis not estimated [7], C must be obtained by another method.
In addition to ignoring the eﬀects of heterogeneity of the geologicalsystem and its three-dimensionality, the analytical solution based on
the line source solution of Kelvin cannot easily account for several
practical issues, such as variations in ambient temperatures during thetest [8] or variations or instability in the heating source at the inlet [9].
It should be emphasized that the heat transfer is not solely due to
conduction, but also a result of convection due to regional groundwater
ﬂow. In a region where groundwater ﬂow is relatively fast, one can
expect a higher heat exchange rate. Multiple studies have investigatedthe eﬀect of groundwater ﬂow on borehole heat exchangers (BHE). For
example, Liuzzo-Scorpo et al. [10] demonstrated that even very low
groundwater ﬂow rates can signiﬁcantly reduce the Inﬂuence Length ofthe exchanger. In another study, Wang et al. [11] conducted a thermal
performance experiment of a BHE under groundwater ﬂow in China
and showed that the presence of groundwater ﬂow enhanced thethermal performance of the borehole heat exchanger and inﬂuenced thetemperature proﬁle in the aquifer. Although recent studies have proposed complex analytical solutions [12], such eﬀects are not taken into
account in most of them considering the high nonlinearity and mutualinteractions of hydrological and thermal processes involved. These interactions were properly emphasized in a recent study of Liuzzo-Scorpo
et al. [13], where the authors proposed a fast graphical method to estimate the hydraulic conductivity of the aquifer using the TRT analysis.
It is thus necessary to develop more reliable and more robust approaches that can account for the above-discussed factors aﬀecting
ground heat exchange so that ground thermal properties may be moreaccurately estimated. For example, numerical solutions of the equations
governing the phenomena related to ground heat exchange processes
have been coupled with parameter estimation procedures to estimate
thermal properties from TRT data [14]. Numerical approaches provide
much more ﬂexibility in accounting not only for more physical processes, detailed geometry, and diﬀerent materials, but also for various
initial, boundary, and operational conditions that are usually oversimpliﬁed when analytical solutions are used. Numerical models,therefore, allow more detailed representations of the system and phenomena related to TRT. More importantly, by using numerical solutions, we do not need to omit early time data of TRT, which are usually
discarded, when estimating ground thermal properties using analytical
approaches [15]. By being able to use early data, which reﬂect thermal
properties of the tube and the grout, it is possible to estimate thermal
properties of diﬀerent materials composing the system. There havebeen numerous studies that have used numerical models. These studies
usually employed the Finite Diﬀerence (FD), Finite Elements (FE), or
Finite Volume (FV) methods. Florides et al. [16] applied a model thatcombined three-dimensional (3D) conduction with one-dimensional(1D) mass ﬂow and 1D convective heat transfer within the carrier ﬂuid.
The resulting 3D model was implemented in the FlexPDE® environment
and used to investigate the thermal performance of single and double Utube borehole heat exchangers in multiple-layer substrates. However,
the eﬀect of groundwater was neglected by Florides et al. [16]. A similar approach was used by Ozudogru et al.[17], who coupled thecommercial software COMSOL® for the description of the heat transfer
in the borehole with a 1D description of heat and water transport in thepipe. In another study, Han et al. [18] developed a 3D coupled Finite
Element Model, which was then utilized to simulate steady-state andtransient behaviors of a geothermal heat exchanger. Han et al. [18]validated the model and then used it to investigate the inﬂuence ofseveral factors on the borehole behavior by carrying out a sensitivity
analysis. However, one of the main disadvantages of 3D mechanistic
models is their computational cost, which limits their use inoperation of GSHP systems depend on the heat exchange capacity, it iscritical to know the ground thermal conductivity before the installationof the GSHP systems.
The apparent thermal conductivity of the ground can be obtained
from one of many in-situ thermal conductivity tests. In many in-situ
tests, the ground is considered to be homogeneous although it often
displays clear heterogeneities, such as layering. It is common to neglect
the eﬀects of layering or heterogeneity when estimating the groundthermal conductivity. This assumption was recently investigated by Luo
et al. [2], who analyzed two modeling approaches in which groundproperties were considered to be either homogeneous or stratiﬁed. The
comparison of ﬂuid outlet temperatures between these two cases
showed relatively minor diﬀerences, indicating that the assumption of ahomogeneous ground introduces a relatively low bias and that the
Thermal Response Tests or TRTs [3] are suitable to determine the effective thermal conductivity of the ground. TRTs are the in-situ tests
most frequently used to estimate the apparent ground thermal conductivity and the borehole thermal resistance. In TRTs, a non-freezing
solution or water is circulated inside of a closed-loop pipe installedvertically or horizontally in the subsurface while it is heated by aconstant heat source at the inlet. The inlet and outlet ﬂuid temperatures
are then continuously recorded. These temperatures are later used toestimate thermal properties of the ground and the borehole. A numberof in-situ mobile or immobile devices for TRT tests have been developed[4]. Furthermore, a number of data analysis methods, based on either
simple analytical solutions or more complex numerical solutions, havebeen used to estimate thermal properties of the ground.
When a vertically installed heat-exchange well (which is in manycases 50–100 m long) is used, the closed-loop pipe is much longer in thevertical direction than its diameter, which is at most several tens of mm.
For this reason, the data analysis typically uses an analytical solutionthat is either based on the line source function of Kelvin [5] or thecylindrical heat source theory [6]. The Kelvin's source function assumesa system in cylindrical coordinates with an inﬁnitely long line heatsource in the center. Under the assumptions that the heat ﬂux from theheat source is constant regardless of time or depth, the medium has
uniform thermal characteristics, and the heat transfer in the verticaldirection can be neglected, the following analytical solution is obtainedby solving the heat conduction equation that describes the change in
the ground temperature ΔT at a given time t and distance r. It is also
assumed that temperature at inﬁnity is constant:∞

ΔT = T (r ,t )−Ti = −

s=

Q′
Q′ ⎛
(−1) k − 1s k ⎞Ei (−s ) = −γ + lns− ∑
⎜
4πλ
4πλ
k · k! ⎟
k=1
⎝
⎠

r 2C
4λt

(1)
(2)

where Ei(−s) is the integral exponential function, Q′ is the heat per unit
length of the line heat source [M L T−3], λ is the thermal conductivity[M L T−3K−1], C is the volumetric heat capacity [M L−1 T−2 K−1], andγ is the Euler constant, which has a value of 0.5772 [5,6]. When t islarge and s is suﬃciently small (s < 0.05), Eq. (1) can be approximatedusing the thermal diﬀusivity parameter α (=λ /C) [L2 T−1] by the following equation:

(

)

r2

ΔT = T (r ,t )−Ti = − 4πλ′ ln 4αt + 0.5772
Q

Q
Q= 4πλ′ ln(t )− 4πλ′

(

)

r2
ln 4α + 0.5772= mln(t ) + b

(

r2

)

m = 4πλ′ , b = − 4πλ′ ln 4α + 0.5772
Q

Q

(3)

In practice, the average ﬂuid temperature between an inlet and an
outlet of the U-tube during TRT is taken as T(0, t) in Eq. (3) to estimate
λ . The most straightforward approach is to plot T(0, t) as a function of
the natural logarithm of time. The thermal conductivity λ can then be1114

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

conjunction with parameter estimation algorithms. Thus, more researchis needed to develop accurate and computationally eﬃcient models.
When the numerical approach is adopted to analyze the TRT data,
special attention needs to be paid to a dramatically diﬀerent dimensionality of the tube and the ground. While the thickness of the
tube is on the order of several millimeters, the tube length and thegeological system are usually on the order of a hundred or even several
hundred meters. Signorelli et al. [19] used the ﬁnite-element-based
numerical code FRACTure [20] to evaluate the eﬀect of factors notconsidered by analytical solutions since his code allowed one to use acombination of lower and higher dimension elements. In their analysis,
a one-dimensional model was used for heat transport in tubes, while athree-dimensional model was used to simulate heat transfer in thesurrounding ground. Bozzoli et al. [21] took a similar modeling approach where both one- and three-dimensional equations were solvedfor diﬀerent domains to inversely estimate thermal properties of the soil
and the grout. In their analysis, the heat transfer between the ﬂuid
inside the tube and the solid phase (e.g., tube or grout) was modeledusing the Robin boundary condition with the heat transfer coeﬃcient,which is a function of ﬂuid parameters. Instead of using the Robinboundary condition, Christodoulides et al. [22] used a source/sink termin the governing equation to account for the heat transfer between theﬂuid and the tube, assuming the heat transfer within the tube body was
instantaneous because the thickness of the tube was extremely thin. An
interesting modeling approach was proposed by Kim et al. [23], whoapplied a Model Order Reduction (MOR) technique to develop a computationally eﬃcient model for the numerical analysis of a verticalborehole heat exchanger. Kim et al. [23] combined a vertical slice
model with a Finite Element discretization of horizontal cross sectionsand validated it against analytical solutions. The authors reported acomputational time reduction of 95% compared to a fully 3D model.
However, the proposed model lacked a proper description of hydrological processes in the vadose zone, which proved to signiﬁcantly affect the thermal performance of the borehole heat exchanger.
In this view, another factor that is usually excluded from the analysis is the eﬀect of the mass and energy transport in the unsaturated
zone, which usually surrounds a signiﬁcant portion of the verticallyinstalled well. In the case of horizontally installed U-tubes, surrounding
soils are often unsaturated. The computational cost of evaluating ﬂow
and transport processes in the unsaturated zone is signiﬁcantly higher
than in the saturated zone. This may become an issue when numericalapproaches are used, especially for three-dimensional models.
Developing a numerically eﬃcient method that can account for detailed
saturated-unsaturated processes is therefore essential in promoting a
numerical data analysis of borehole heat exchangers. Since GSHP arestrongly inﬂuenced by hydrological processes occurring in the vadosezone, it is necessary to bridge the gap and develop models that are able
to describe the thermal and hydrological behavior of the entire groundborehole system.
Thus, the ﬁrst objective of this study was to develop a computationally eﬃcient pseudo-3D model for the numerical analysis and interpretation of TRTs. The numerical approach combines a one-dimensional description of the heat transport in the buried tubes of theexchanger with a two-dimensional description of the heat transfer inthe surrounding subsurface soil, therefore reducing the dimensionalityof the problem and the computational cost. Furthermore, the modelingframework includes the widely used hydrological model, HYDRUS [24],which can simulate the movement of water, heat, and multiple solutes
in variably-saturated porous media. The combination of a reducedorder numerical approach with a hydrological model guarantees a highmodeling ﬂexibility and represents a new contribution to this ﬁeld. The
second objective of this study is to validate the newly developed model
against experimental data collected at two diﬀerent experimental sites
in Japan and to use the validated model to carry out a sensitivity
analysis, based on the Morris method [25], that would investigate theinﬂuence of diﬀerent materials on the TRT response. Finally, the third

Fig. 1. A schematic of the experimental setup used to perform the thermal response test.

objective is to demonstrate the capability of the model to handle various hydrological processes in the vadose zone by investigating theinﬂuence of groundwater ﬂow and aquifer thickness on the boreholeheat exchanger. A further scenario has been used to demonstrate theability of the proposed modeling framework to simulate the eﬀect oflayered soil proﬁles on the heat transfer process.
2. Materials and methods
2.1. Site description
The Thermal Response Tests (TRT) were performed at the campusesof Tokyo University of Agriculture and Technology (TUAT) (Fuchu city,
Tokyo; 35°40′59″N, 139°28′58″E) and Saitama University (SU)
(Saitama city, Saitama; 35°51′44″N, 139°36′34″E) in Japan. A schematic of the experimental setup is shown in Fig. 1. Two 50-m long
vertical boreholes with an inner diameter of 180 mm were drilled in2011 at TUAT to be used as borehole heat exchangers for a groundsource heat pump system installed for a 25 m2 constant temperatureroom [26]. Without ﬁnishing with casing, two U-shaped 3-mm thick 20mm diameter pipes made of high-density polyethylene (HDPE) wereinserted into the borehole along with space-ﬁlling grout, which consisted mainly of silicate sand. Two U-shaped pipes were inserted together into the borehole since the bottom of both tubes were connected.
The groundwater table at TUAT was 10–12 m below the ground surface.
At SU, a 50-m long vertical borehole with an inner diameter of 180 mm
was also drilled in 2011 to be used for a long-term heating and coolingexperiment [27]. Similar to TUAT, two U-shaped 3-mm thick HDPE
pipes were installed in the borehole. Silica sand was used as a spaceﬁlling grout material. The groundwater table at SU was 1–4 m below
the ground surface.
Fig. 2 shows the distribution of particle sizes, water contents, andbulk densities measured at boring core samples collected at both BHE[28]. The particles greater than 2 mm (i.e., gravel) were excluded fromthe particle size distribution analysis. While the elevation of the TUAT
site, which is located at the Tachikawa terrace, is about 60 m above themean sea level (MSL), that of the SU site, which is in the Arakawa
Lowland near the Tokyo metropolis, is about 10 m above MSL. At the
TUAT site, there are three layers of gravel located at depths of 5–11 m,
32–34 m, and 43–45 m, where particle size distribution data were notavailable (or where water content values were very small due to low
water retention of gravel). The average gravel content in each layer was
0.785 kg kg−1, 0.848 kg kg−1, and 0.808 kg kg−1, respectively. At theSU site, a gravel layer was observed only around the depth of 40 m. On
top of the gravel layer, silt layers and sand layers were deposited alternately.
During TRT, a non-freezing solution of 20% propylene glycol washeated to keep the temperature diﬀerence between the inlet and outlet
at 4.5 °C with a constant power of 2.5 kW for 48 h while the ﬂuid was1115

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

Fig. 2. Particle size distributions, water contents,
and dry bulk densities from core samples collected at TUAT (the upper plot) and SU (the lowerplot), respectively.

circulated at a constant rate of 9.0 L min−1 in one of the U-tubes in one
HEB at TUAT. The temperatures of inlet and outlet ﬂuid were continuously measured using a platinum temperature sensor whose measurement accuracy was less than 0.1 °C. As a result, the average heatexchange rate at the TUAT site was estimated to be 28.2 W m−1. At SU,
tap water was circulated in both U-tubes for 48 h with a constantheating power of 2.5 kW at a constant ﬂow rate of 13.5 L min−1. As aresult, the average heat exchange rate at SU was estimated to be 56.4 W
m−1. The inlet–outlet temperature diﬀerence was kept at 3 °C. In thefollowing analysis of the SU data, instead of using the double U-tubes,
we assumed a single U-shape tube and used an eﬀective radius of
3.7 mm which resulted in the surface area equal to the original double
U-shape tubes.

1116

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

Fig. 3. A schematic of the modeled transport domain consisting of the subsurface soil and two branches of the heatexchanger.

borehole wall, while it is less pronounced for cement-grouted boreholeheat exchangers, as reported by Choi and Ooka [33]. Choi and Ooka[33] also emphasized how the eﬀect of the natural convection in theborehole was mainly dominated by the hydrothermal properties of thesurrounding soil, rather than the ﬁlling material. In particular, thepermeability of the soil is directly related to the magnitude of thebuoyancy ﬂow [34].
The modeling approach proposed in this study should be usedcarefully in scenarios when the borehole is ﬁlled with groundwater andwhen relatively high heat injection rates are used as these could lead tobiased results. Under such circumstances, more complex 3D models
should be adopted.

2.2. Modeling theory
2.2.1. Modeling approach
A computationally eﬃcient pseudo 3D-model for the numericalanalysis of heat transfer in borehole heat exchangers has been developed in this study. The proposed approach combines a one-dimensionaldescription of the heat transport in the buried tubes of the exchanger
with a two-dimensional description of the heat transfer in the surrounding subsurface soil. The main assumption is that the heat transfer
in the soil domain occurs only in the horizontal plane. Similar approaches have been proposed by Yavuzturk et al. [29] and Austin et al.
[30]. However, the main novelty of the present study is the use of thewell-established hydrological model, HYDRUS, in the modeling framework.
The vertical three-dimensional soil proﬁle is discretized in N twodimensional horizontal cross-sections. Each cross-section contains two
nodes of the one-dimensional ﬂuid domain (Fig. 3), which represent thetwo branches of the tube. The heat transfer between the two domains,i.e., between two-dimensional cross sections and a one-dimensional
tube, is represented by a Robin boundary condition. The reduced dimensionality of the proposed approach guarantees a signiﬁcant computational gain compared to a fully three-dimensional model while
maintaining a comparable accuracy. However, it must be emphasizedthat the proposed approach neglects the eﬀects of natural convection(and other vertical, parallel with the borehole, heat ﬂuxes) on theborehole that may occur under saturated conditions and that may enhance the heat transfer and reduce the thermal resistance of the borehole. The magnitude of the natural convection depends on the temperature diﬀerence between the ﬂuid in the U-tube and the surrounding
ground and it is thus directly related to the heat injection rate. Gustafsson and Gehlin [31] reported a decrease of more than 10% in theborehole thermal resistance when the heat injection rate increased from
40 W/m to 80 W/m. As pointed out by Spitler et al. [32], this process israther signiﬁcant for groundwater-ﬁlled borehole heat exchangers inwhich groundwater ﬁlls the annular space between the U-tube and the2.2.2. One-dimensional heat transport in the vertical domain
The heat transport in both tubes of the heat exchanger is assumed to
be one-dimensional and purely advective [14]:

ρf cf

∂Tf
∂t

= ρf cf u

∂Tf
(4)

∂z

where ρf is the density of the ﬂuid [M L−3], cf is the speciﬁc heat capacity of the ﬂuid [L2 T−2 K−1], Tf is the temperature of the ﬂuid [K], u
is the ﬂuid velocity [L T−1] (positive for the descending tube and negative for the ascending tube), and z is the vertical coordinate [L] assumed positive downward. The convective heat transfer between the
tube and the ﬂuid was described by a source term added to Eq. (4),which becomes:

ρf cf

∂Tf
∂t

= ρf cf u

∂Tf
∂z

−

h
Af

∫P [Tf (z,t )−Tp (p,z,t )] dp
−3

−1

(5)

where h is the heat transfer coeﬃcient [M T K ], P is the pipe innerperimeter [L], Tp is the temperature on the inner surface of the pipe [K],and Af is the area of the tube [L2]. The integral in the last term of (5)averages the temperature between the ﬂuid and the HDPE pipe along its
inner perimeter. Eq. (5) can be written also as:1117

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

ρf cf

∂Tf
∂t

= ρf cf u

∂Tf
∂z

+

h
Af

∫P Tp (p,z,t ) dp− 4dhp Tf (z,t )where δij is the Kronecker delta function, and λL and λT are the longitudinal and transverse thermal dispersivities [L]. λ0(θ) is describedwith the equation [36]:(6)

where dp is the inner diameter of the pipe [L] (Fig. 2).
The convective heat transfer coeﬃcient can be estimated using therelation:

h=

λf
dp

Nu

λ 0 (θ) = b1 + b2 θ + b3 θ 0.5

where b1, b2, and b3 are empirical parameters [M L T−3 K−1].
The convective heat transfer between the ﬂuid and the solid domain
has been described using a Robin boundary condition (Eq. (13)):

(7)
−3

−1

where λf is the thermal conductivity of the ﬂuid [M L T K ], and Nu
is the Nusselt number [–], which is expressed using the Dittus-Boeltercorrelation as:

λij

where Re is the Reynolds number [–], Pr is the Prandtl number [–], andn is a coeﬃcient assumed to be 0.33 for cooling (the wall hotter thanﬂuid) and 0.4 for heating (the wall colder than ﬂuid).
Eq. (6) has been solved using the Finite Volume Method (FVM). A
second-order accurate QUICK (Quadratic Upstream Interpolation forConvective Kinematics) scheme [35] has been adopted due to its capability to reduce the artiﬁcial numerical dispersion and achieve a good
accuracy of the solution. A semi-implicit time stepping scheme has beenused. The ﬁrst (i.e., advective component) and second terms on theright-hand side of Eq. (6) have been solved explicitly in time, while thethird term has been solved implicitly. The linearization of the second
term, which is a function of the temperature of the tube wall, allows fora signiﬁcant reduction of the computational cost. An implicit treatmentof this term would have required running a two-dimensional model forthe soil domain several times at each time step, signiﬁcantly increasing
the simulation time.
The problem deﬁnition was completed by the speciﬁcation of theinitial and boundary conditions. The condition of a constant power
supplied to the ﬂuid has been implemented using a periodic Dirichletcondition for the descending tube:

T fdescending (0,t ) = T fascending (0,t ) + ΔT

−λij

λij

2.2.3. Two-dimensional heat transport in the horizontal domain
The HYDRUS (2D/3D) software [24] has been used to simulate theheat transport in the horizontal domains. The software contains a twodimensional ﬁnite element model HYDRUS-2D for simulating the
movement of water, heat, and multiple solutes in variably-saturated
porous media. Neglecting the eﬀects of water vapor diﬀusion, two-dimensional heat transport can be described as:

qj qi
|q|

(10)

+ λ 0 (θ) δij

(14)

∂Ts
ni = 0∂x j(15)

1. A common time step is set for both models;
2. The simulation begins with the numerical discretization of the vertical domain for both legs of the U-tube. Python reads the temperatures on the border of the cross-sectional domains and solves
Eq. (6), updating the values of the ﬂuid temperature in diﬀerent
nodes;
3. A for loop is used to iterate through diﬀerent horizontal cross-sections. At each iteration, the initial condition is ﬁrst updated and thecalculated ﬂuid temperatures are passed to the HYDRUS solver forthe computation of the Robin boundary condition. HYDRUS-2D isthen executed and the resulting calculated solid temperatures are
stored in a 2D matrix for the next time step. This matrix contains
temperature distributions for all horizontal cross sections and is
updated at each time step.

where C(θ) and Cf are the volumetric heat capacities of the solid andliquid phases, respectively [M L−1 T−2 K−1], θ is the volumetric watercontent [–], λij(θ) is the apparent thermal conductivity of the soil
[M L T−3 K−1], and qi represents the water ﬂux in the porous media[L T−1]. The ﬁrst term on the right-hand side of Eq. (10) represents theheat transfer by conduction in the soil, while the second term accounts
for heat being transported by ﬂowing water. The apparent thermal
conductivity λij(θ) combines the thermal conductivity λ0(θ) of theporous medium in the absence of ﬂow and the macrodispersivity, whichis assumed to be a linear function of the velocity. It can be expressed as:

λij (θ) = λT Cw |q|δij + (λT −λL) Cw∂Tsni + TCf qi ni = T0 Cf qi ni∂x j

2.2.4. Models coupling strategy
The one-dimensional vertical domain and the two-dimensionalcross-sectional domains have been coupled using a Python script, which
acted as a “glue” between the two models. The script simultaneously
solves the 1D advective heat transport in the tube, interacts with
HYDRUS-2D, and exchanges data between the two models. As mentioned previously, a horizontal HYDRUS-2D cross section was located ateach node of the 1D ﬂuid domain. The coupling strategy is summarized
below:

(9)

∂Ts
∂ ⎛
∂Ts ⎞
∂Ts
=
⎜λ ij (θ )⎟−Cf qi∂t∂x i ⎝∂x j ⎠∂x i(13)

in which T0 is the temperature of the incoming ﬂuid. When theboundary is impermeable or when water ﬂow is directed out of theregion, Eq. (14) reduces to a Neumann type BC of the form:

where ΔT is a constant temperature increment [K]. The U-tube connection was modeled by imposing the same temperature on the two legsof the tube at z = L, which was numerically accomplished by using a
Dirichlet boundary condition for the ascending pipe. The one-dimensional heat transport in the vertical domain has been solved numerically in the Python programming language.

C (θ)

∂Ts
ni = h·[Tf −Tpi]∂x jwhere Tpi is the temperature on the inner surface of the tube [K].
Since the Robin BC is not included in the standard version of
HYDRUS, an ad hoc version has been developed. Once the ﬂuid temperature and the convective heat transfer coeﬃcient are speciﬁed,
HYDRUS automatically updates the value of the heat ﬂux at theboundary. HYDRUS-2D uses an implicit adaptive time stepping to solve
Eq. (10), thus minimizing the computational cost and guaranteeing theaccuracy of the solution.
A Cauchy-type boundary condition was speciﬁed on the remaining
boundaries of the domain (green1 line in Fig. 2). The Cauchy BC prescribes the heat ﬂux along the boundary, and can be expressed as:(8)

Nu = 0.023Re 0.8 Pr n

(12)

An important feature of this approach is that each horizontal crosssection can have its own hydraulic and thermal properties reﬂecting
diﬀerent geological layers and their saturations, which guarantees asigniﬁcant modeling ﬂexibility.
1
For interpretation of color in Fig. 2, the reader is referred to the web version of thisarticle.

(11)
1118

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

2.3. Model validation

Table 1
Geometric characteristics, thermal properties of diﬀerent materials, and input parameters
used in the numerical simulation.

2.3.1. Objective function
To assess its accuracy, the model was validated against experimental data collected at two experimental facilities described above. In
particular, the measured and simulated outlet temperatures werecompared. The use of the ﬂuid outlet temperature as a reference datafor validating the model is widespread in conventional TRTs. However,it must be emphasized that limited information content of the outlet
temperatures does not allow for a precise estimation of thermal properties of diﬀerent soil layers. Beier et al. [37] investigated the usefulness of diﬀerent experimental data sets for the TRT analysis and model
validation. Beier et al. [37] also considered, in addition to usual measurements of outﬂow temperatures, temperatures measured on theborehole wall and in the surrounding soil. The analysis revealed thatboth data sets can be used to better constrain the parameter estimationand provide independent validation measurements. For example,
thermal probe data can provide an independent measurement of thesoil thermal conductivity. This information can be used to validate
optimized parameters or as prior information in the Bayesian inference
framework.
In the present study, the Root Mean Square Error (RMSE) was used
to quantify the overall quality of the ﬁt:

Internal diameter of the U-tube, dp (m)U-tube thickness, tp (m)U-tube length, L (m)Fluid velocity, uf (m s−1)
Heat transfer coeﬃcient, h (W m−2 K−1)
Initial temperature, T0 (K)Inlet-outlet temperature diﬀerence, ΔT (K)
Fluid density, ρf (kg m−3)
Fluid speciﬁc heat, cf (J kg−1 K−1)
HDPE density, ρp (kg m−3)
HDPE speciﬁc heat, cp (J kg−1 K−1)
HDPE thermal conductivity, λp (W m−1 K−1)
Diameter of the heat exchange well, dw (m)∑i = 1 (Tobs,i−Tmod,i )2(16)nwhere Tobs,i is the ith measured temperature [K], Tmod,i is the ith modeled temperature [K], and n is the number of observations. The RMSE isa common choice in borehole model validation and the TRT analysis[38].
2.3.2. Numerical domain and boundary conditions
The thermal properties were taken from a previous study of Saito
et al. [39], which was carried out at the same experimental facilities.
Only the thermal properties of the grout for the TUAT experimental sitewere manually adjusted to increase the overall accuracy of the ﬁt. The
thermal properties that were used for model validation in both casestudies are listed in Table 2. Since FVM schemes can be aﬀected by false
numerical diﬀusion when the velocity of the simulated process is relatively high and when the mesh size is large, a preliminary sensitivity
analysis has been carried out to investigate the inﬂuence of the vertical
discretization on the numerical solution. Usually, a combination ofhigher order schemes and ﬁne mesh sizes is needed to reduce the false
numerical diﬀusion. However, the use of a very ﬁne vertical discretization would have increased the computational cost of the analysis.
Thus, the main aim of the preliminary mesh sensitivity analysis was toestablish a balanced vertical discretization, guaranteeing both a highaccuracy and a reasonable computational cost of the numerical solution.
The domain was discretized into two-dimensional triangular elements using the MESHGEN tool of HYDRUS-2D. The mesh was reﬁned
near the pipes to accommodate the signiﬁcant temperature gradientsParameter

TUAT

SU

HDPE

Thermal conductivity λp‡ (W/m2 K)
Volumetric heat capacity Cp (MJ/m3 K)

0.51
1.71

0.51
1.71

Grout

Thermal conductivity λg (W/m2 K)
Volumetric heat capacity Cg (MJ/m3K)

1.70
1.20

2.80
3.10

Soil

Thermal conductivity λs (W/m2 K)
Volumetric heat capacity Cs (MJ/m3 K)

2.28
3.10

1.70
3.10

‡

0.020
0.003
50
0.15
834.7
290.15
4.5
1000
4182
950
1800
0.51
0.178

0.037
0.003
50
0.21
1038.1
289.35
3
1000
4182
950
1800
0.51
0.178

2.4. Sensitivity analysis
2.4.1. Morris method
The main aim of developing an accurate and reliable model is topotentially exploit it to analyze diﬀerent components of the simulated
process. A classic example is to use the model to investigate the eﬀects
of diﬀerent parameters on some variables of interest, namely a sensitivity analysis. For example, Florides et al. [40] used a previously validated three-dimensional ﬁnite element model to examine the inﬂuence of the ground properties on the thermal behavior of a boreholeheat exchanger. In particular, authors highlighted how the soil thermal
conductivity plays a major role toward dissipating the heat into theground as well as when this sensitivity starts increasing with time. In
another study, Wagner et al. [41] used a three-dimensional ﬁnite element model to investigate the eﬀect of a shank spacing, a non-uniform
initial thermal distribution, and the thermal dispersivity on a thermal
response test. However, the main limitation of such types of sensitivity
analyses is their lack of a rigorous statistical base.
The Morris method [25] belongs to the class of Screening methods(SM). SMs aim to provide qualitative sensitivity measures for diﬀerent
factors using a relatively small number of model evaluations. In general, the Morris method is a one-factor-at-a-time (OAT) local method,since it computes the elementary eﬀect by changing only one factor at atime. However, it can be viewed as a Global method, since it averages
several elementary eﬀects computed at diﬀerent points in the parameter space, providing a qualitative measure of the importance andinteractions of diﬀerent factors.
In this study, the modiﬁed version of the Morris method proposed

Table 2
Thermal properties (thermal conductivity λ and thermal capacity C) of various materials(HDPE is a pipe material, high-density polyethylene) used in the model validation.
Material

SUgenerated by the convective heat transfer, therefore increasing theoverall accuracy of the solution. The generated FE mesh had 1114 and1234 nodes, and 2117 and 2346 two-dimensional elements for the
TUAT and SU experimental site, respectively. The Robin BC was used tosimulate the heat transfer between the carrier ﬂuid and the pipes, whilea Cauchy BC was used for the remaining boundaries. The initial temperatures and input parameters are summarized in Table 1.
Considering that only the inlet and outlet ﬂuid temperatures were
measured, the soil proﬁle was assumed to be homogeneous, thus neglecting the eﬀects of diﬀerent soil textures. Moreover, the heat transfer
in the porous medium was assumed to be independent of the volumetricwater content, i.e., the contribution of the terms with the parameters b2
and b3 to the thermal conductivity (see Eq. (12)) was neglected. The
inclusion of these parameters in the inverse analysis would increase thedimensionality of the inverse problem and pose a parameter identiﬁcation problem, considering that only one type of measured data (i.e.,
outlet temperature) was available for the analysis. The overall analysis
was conducted with a laptop equipped with a CPU Intel® Core i7-4700
MQ 2.40 GHz processor and 8 GB of RAM.

nRMSE =

TUAT

Subscripts p, g, and s refer to HDPE, grout, and soil, respectively.

1119

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

by Campolongo et al. [42] has been used to investigate the inﬂuence ofthermal properties of ground, grout, and pipe on the TRT response. σ
and μ∗ are the two sensitivity measures calculated in the Morris method.
While the former summarizes the interaction eﬀect, the latter reﬂects
the overall importance of a particular parameter. For a detailed description of the method, refer to Morris [25] and Campolongo et al.
[42]. To interpret the results by simultaneously taking into accountboth sensitivity measures, Morris suggested their graphical representation in the (μ∗–σ) plane.
The sensitivity analysis has been performed for two diﬀerent simulation durations of 10,800 and 30,000 s to examine the dynamics ofthe sensitivity measures. Although, for the sake of simplicity, the analysis has been applied only to the TUAT experimental site, its results canbe generalized. The RMSE has been used as the objective function in thesensitivity analysis. Considering the computational cost of the modeland the intent of the present analysis, which was targeted to a qualitative evaluation of the parameters inﬂuence, the sample size was set to5, with a total of 35 numerical simulations for each scenario.

2.4.2. Groundwater eﬀects
As reported by Signorelli et al. [19], groundwater ﬂow can signiﬁcantly aﬀect borehole heat exchangers when groundwater ﬂuxes areequal to or above 10−1 m/s, which can occur in permeable aquifers or
regions with high hydraulic gradients. In particular, groundwater canenhance the performance of the system by dissipating heat aroundboreholes faster than unsaturated soil. Moreover, groundwater can introduce signiﬁcant nonlinearities in the measured temperatures whichcannot be interpreted accurately with a simpliﬁed model based on an
inﬁnite line source. It must be emphasized that while previous modeling studies have already investigated the eﬀect of groundwater[19,43,44], none of these studies could fully exploit the features of acomplex hydrological model, such as HYDRUS, which can give acomprehensive description of a variety of subsurface physical processes.
The proposed modeling framework has been used to investigate theinﬂuence of groundwater on the borehole heat exchanger. In particular,
the eﬀects of the aquifer thickness and groundwater ﬂow have been
examined. A 50-m long U-tube borehole heat exchanger, similar to theone used at the TUAT experimental site, has been simulated. The soil
proﬁle was assumed to be homogeneous and composed of sand. The
thermal properties of a typical sandy soil are reported in Table 3.
Three diﬀerent aquifer thicknesses were simulated: 10, 20, and
30 m. The Python code automatically overwrites the initial pressure
head at each time step depending on the vertical coordinate of the
Finite Volume calculated. A pressure head of −1 m, which corresponds
to θ ≈ 0.05 for sandy soils, is imposed above the water table. A positive
pressure head (full saturation) is set below the water table.
Similarly, the eﬀect of groundwater ﬂow in a 30-m deep unconﬁned
aquifer on the TRT response has been simulated. Three diﬀerent waterﬂow values were considered: 0.0, 2.5 · 10−6, and 5 · 10−6 m/s. Waterﬂow was simulated in HYDRUS-2D using the constant ﬂux BC, with awater temperature of 290.15 K. The input parameters (i.e., the U-tube
thickness, ﬂuid velocity, the inlet-outlet temperature diﬀerence, etc.)
and other thermal properties (i.e., HDPE, grout) were the same as thoseused for the TUAT model validation. A schematic of the analyzed scenario is reported in Fig. 4.

Fig. 4. A schematic of the modeling scenario used to investigate the inﬂuence ofgroundwater ﬂow. The water ﬂow is perpendicular to the borehole heat exchanger.

2.4.3. Layered soil proﬁle
The assumption of a homogeneous and isotropic medium may not
be valid for particular geological settings. Heterogeneities in subsurface
thermal properties can signiﬁcantly inﬂuence the heat transfer between
the borehole heat exchanger and the surrounding soil [45]. However, inorder to estimate thermal properties of each layer, it is necessary toinclude diﬀerent types of measurements in the optimization process.
When applying the conventional TRT, only the inlet and outlet ﬂuid
temperatures are usually available. In such circumstances, an inverse
parameter estimation of thermal properties of multiple soil layers
would cause parameter identiﬁcation problems and increase the parameter uncertainty. If diﬀerent types of measurements were available(e.g., ﬂuid temperatures in the U-tube), it would be possible to better
constrain the optimization problem, considering increased information
content of experimental data. In such case, a proper modeling framework would include the use of a multi-objective optimization algorithm
to obtain a trade-oﬀ between diﬀerent likelihood functions assigned toeach measurement. In this view, optical ﬁbers represent a promising
tool to retrieve important experimental data from borehole heat exchangers.
Optical ﬁber thermometers in geothermal wells started to be used in
the late 1980s. Optical ﬁbers can be used to measure the ﬂuid temperature along the buried pipe of the exchanger. For example, Fujii
et al. [46] conducted a series of TRTs using optical ﬁbers in a U-tube
borehole heat exchanger in Fukuoka, Japan. The measured ﬂuid temperatures were then used to estimate the vertical distribution of thermal
conductivities. In Fujii et al. [46], a cylindrical source function coupledwith a nonlinear regression technique was applied. However, the ﬂuid
temperature was assumed to be constant during the circulation period,
thus simplifying the problem. In another study, Acuña and Palm [47]proposed the use of optical ﬁbers to measure the borehole wall temperature. The technique consisted of a ﬁber optic cable pressed against
the borehole wall with a ﬂexible plastic pipe that is ﬁlled with water
once it is installed in the ground.
To demonstrate the capability of the newly-developed model toprovide a comprehensive description of the heat transfer in a layered
proﬁle, a speciﬁc scenario has been simulated. In particular, a lithologic

Table 3
Thermal properties of sand [36].
λL (m)λT (m)

b1 (W/mK)b2 (W/m K)

b3 (W/mK)

Cs (MJ/m3
K)

Cf (MJ/m3
K)

0.05

0.01

0.228

−2.406

4.909

3.100

4.182

1120

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

numerical solution, indicating that the model is relatively insensitive of the vertical mesh discretization;– Considering the small numerical error and relatively low computational cost, the model with a coarse vertical mesh can be used when
the aim of the analysis is to test diﬀerent thermal properties ofdiﬀerent materials (i.e., pipe, grout, and soil) and to simulate theoutlet temperature;– When the numerical analysis considers the eﬀects of diﬀerent layers,
the groundwater table depth, and water ﬂow on the thermal behavior of the borehole heat exchanger, a ﬁner mesh should be used. Insuch cases, an additional preliminary mesh sensitivity analysis isrequired.
– Our analysis revealed that a value of N = 5 provided a good tradeoﬀ between computational cost and numerical accuracy. As a result,the model with N = 5 was used in the following modeling scenarios.

Table 4
Simulated lithologic composition and soil distribution.
Depth (m)Soil type

0–10
10–20
20–30
30–40
40–50

Loam
Sand
Clay
Sand
Clay

column composed of a series of sandy, clayey, and loamy soil layers has
been considered. The exact soil distribution along the vertical depth isreported in Table 4.
A 50-m long U-tube borehole heat exchanger, similar to the oneused at the TUAT experimental site, has been simulated. The pressure
head was assumed to be -1 m and constant in the entire lithologic
proﬁle. Values reported in Chung and Horton [36] for loam, sand, andclay were used to characterize thermal properties. The input parameters(i.e., the U-tube thickness, ﬂuid velocity, the inlet-outlet temperature
diﬀerence, etc.) and other thermal properties were again the same as in
the TUAT model validation. In order to highlight diﬀerences between
homogeneous and heterogeneous systems, the results for the layered
proﬁle were compared with those for a homogeneous lithologic proﬁlecomposed of sand.
For the sake of clarity, diﬀerent modeling scenarios with their assumptions are summarized in Table 5.

The comparison between the measured and modeled outlet temperatures for both experimental sites is shown in Fig. 6. In particular,
the simulated temperatures of the ﬂuid at the outlet calculated using
the values listed in Table 2 (red lines in Fig. 6) and the value of thermal
conductivity obtained through the Line Source approach (Eq. (3))(black lines in Fig. 6) are reported, respectively. The calculated RMSEs
were 0.16 and 0.25 °C for the TUAT and SU experimental sites, respectively, when the best ﬁt values were used. The low RMSE values
conﬁrmed the satisfactory accuracy of the proposed model in reproducing the behavior of the system. In both scenarios, the quality of
the ﬁt was suﬃcient, with a slightly better performance for the TUAT
site. In particular, the model for the SU site slightly overestimatedtemperatures during the ﬁrst 60,000 s of the simulation. This overestimation could be related to thermal properties of the grout material,
whose inﬂuence is signiﬁcant during the ﬁrst part of the test. It must beemphasized that an improved ﬁt could have been obtained by using an
automatic optimization algorithm (e.g., Particle Swarm or Genetic Algorithm). However, this was not the aim of this study, which mainly
focused on the development of a computationally eﬃcient model. Future applications of the model can include the global optimization ofthermal properties of diﬀerent materials, as well as the sensitivity anduncertainty analysis.
The behavior was diﬀerent when the values of thermal conductivity
calculated using the Line Source approach (Eq. (3)) were used in the
numerical simulations. In particular, the estimated RMSEs were 0.19
and 1.93 °C for the TUAT and SU experimental sites, respectively. In
both cases, there was a degradation in terms of ﬁtting quality, which
was more evident for the SU site. For the TUAT, the model reproduced
the thermal behavior of the BHE with suﬃcient accuracy. The assumption of homogeneous thermal properties led to a slight underestimation of the ﬂuid temperature in the ﬁrst part of the simulated
period, which was mainly related to an overestimated value of thethermal conductivity of the grout. This bias tended to vanish for longer
simulation periods. Conversely, the simulated thermal behavior of theborehole at the SU experimental site was quite diﬀerent. In this case,the assumption of homogeneous thermal properties led to a signiﬁcant
overestimation of the ﬂuid temperatures at the outlet, indicating a
potential underestimation of the thermal conductivity of the grout. The
deviation between the two models at t = 120,000 s is almost 2 °C. The
analysis suggests an important role of the grout in the thermal behaviorof the BHE. This aspect has been further investigated in the present
study.
As shown in Fig. 7, the proposed model can calculate the temperature distribution in both the ﬂuid and soil domains and at diﬀerent
time steps, providing an overall picture of the physics of the problem.
At the beginning of the test (t = 30,000 s), it is evident that the heattransfer mainly aﬀects the pipe and the grout, having only a limited
inﬂuence on the surrounding soil. On the contrary, the soil is signiﬁcantly heated after t = 90,000 s, especially near the surface where
the eﬀect of the heating system is important. The temperature

3. Results and discussion
3.1. Model validation
Before proceeding with the model validation and further modelingscenarios, a mesh sensitivity analysis has been carried out. Results are
summarized in Fig. 5, in which the outlet temperatures (Fig. 5, left) are
reported for diﬀerent values of N (the number of vertical elements) andcorresponding computational costs (duration of simulations) (Fig. 5,
right). Generally, ﬁner discretizations are associated with more accurate numerical solutions since false numerical diﬀusion tends to vanish.
At the ﬁrst inspection, it is evident how the outlet temperature wasrelatively insensitive to the vertical mesh, underlining a general robustness of the proposed model. Simulated temperatures for N = 5 andN = 10 were almost the same (overlapped), indicating that a further
mesh reﬁnement would lead to negligible accuracy gains. Conversely,
the results for N = 2 exhibited a slight overestimation of outﬂow
temperatures, which can be related to the false numerical diﬀusion.
However, the temperature diﬀerence at t = 60,000 s for N = 2 andN = 5 was only about 0.2 °C, which is still a relatively small error. On
the other hand, the computational cost (Fig. 5, right) for N = 2 wasonly about half of that for N = 10, and also signiﬁcantly lower than for
N = 5. As a results of this analysis, a number of conclusions can be
drawn:– A coarse vertical mesh generally only introduces a low bias in the
Table 5
Summary of modeling scenarios and their assumptions.
Scenario

Lithologic proﬁleHorizontal
domainInﬂuence of watercontent

Model validation
Morris sensitivity
analysis
Groundwater ﬂow
Groundwater depth
Layered soil proﬁleHomogeneous
Homogeneous

Circular
Circular

No
No

Homogeneous
Homogeneous
LayeredSquaredCircular
Circular

Yes
Yes
Yes

1121

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

Fig. 5. Simulated outlet temperatures (left) for diﬀerent values of N (the number of vertical elements) with their computational cost (right).

Fig. 6. Measured and modeled ﬂuid temperatures at the tube/pipe outlet. Circles represent the measured ﬂuid temperatures, black lines the modeled ﬂuid temperatures calculated using
the thermal conductivity obtained with the Kelvin line source equation, and red lines the modeled ﬂuid temperatures calculated using the values reported in Table 2. (For interpretation of
the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)

λp, which exhibited the highest σ value, indicating its important role ininteractions with other factors. The volumetric heat capacity of the
HDPE had a negligible eﬀect on the model’s response as highlighted byits low values of σ and μ∗, an expected behavior considering the limited
diameter of the pipe. The thermal properties of the ground, λs and Cs,
only had a limited inﬂuence on the output. Such behavior is intuitive,
considering that in the ﬁrst part of the test the ground is only partiallyinvolved in the heat transfer process, which is mainly driven by thegrout and the pipe. This is emphasized in Fig. 8, in which two distinctgroups can be identiﬁed: while parameters λg, Cg, and λp signiﬁcantly
inﬂuence the model’s response, parameters λs, Cs, and Cp have only
limited eﬀects on the output.
The sensitivity measures signiﬁcantly changed at t = 30,000 s (rightplot in Fig. 8). As expected, the thermal properties of the soil increased
their inﬂuence on the model’s output, as indicated by their relatively
high values of σ and μ∗. However, the grout and the pipe still had signiﬁcant eﬀects on the model’s response. Contrary to what was highlighted at t = 10,800 s, it was very diﬃcult to identify groups ofparameters with distinct eﬀects on the model behavior. The only

distribution is asymmetric in the upper part of the domain, as it is expected for a U-tube borehole heat exchanger. It should be emphasizedthat the ﬂuid temperature in both legs of the pipe appears smooth anddevoid of nonlinearities. This behavior is directly related to the assumption of homogeneous soil. In this case, the convective heat transfer
regulated by the Robin BC is not aﬀected by inhomogeneities typical oflayered lithologic proﬁles.

3.2. Sensitivity analysis
3.2.1. Morris method
The results of the sensitivity analysis, carried out using the Morris
method, are reported in Fig. 8. In particular, the scatter plots μ∗–σ fortwo selected simulation times are shown.
At t = 10,800 s (left plot in Fig. 8), the thermal properties of thegrout exhibited the highest μ∗, indicating a signiﬁcant inﬂuence on themodel’s response. More speciﬁcally, the thermal conductivity of thegrout λg was the most inﬂuential parameter, followed by the grout
volumetric heat capacity Cg. The third most inﬂuential parameter was
1122

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

Fig. 7. The temperature distribution in both the ﬂuid and horizontal two-dimensional domains at two diﬀerent simulation times.

visible at t = 30,000 s.
The sensitivity analysis provided important information regarding
the roles of diﬀerent materials at diﬀerent simulation times. In particular:

irrelevant factor was Cs, while the other parameters exhibited a relatively large impact.
The inﬂuence of λg and λs on the model’s output has been further
investigated in Fig. 9, which reports the outlet temperatures for different values of the grout (the left plot in Fig. 9) and ground (the rightplot in Fig. 9) thermal conductivities. Blue arrows and dotted lines indicate increasing values of the thermal conductivity and the simulation
times chosen for the Morris sensitivity analysis, respectively. The results
were obtained by changing the value of the thermal conductivity while
maintaining all other factors ﬁxed. Fig. 8 conﬁrmed the ﬁndings of theMorris sensitivity analysis. The inﬂuence of the grout was signiﬁcantduring the simulation, with a preeminent role in the ﬁrst part of theexperiment. This behavior is particularly evident when observing
temperature values at two dotted lines. At t = 10,800 s, the model’s
output was mainly sensitive to λg, which generated a relatively highvariance of temperatures. After 6000–8,000 s, the thermal conductivityof the ground started to have more important eﬀects, which were more

– The volumetric heat capacity of the pipe generally had a negligible
eﬀect on the model’s output. In the optimization framework, thisfactor could be ﬁxed to any feasible value in the parameter space
without signiﬁcantly aﬀecting results.
– The thermal properties of the grout had a signiﬁcant inﬂuence on
the model’s response during the ﬁrst part of the simulation. Thissuggests that the thermal properties of the grout could be determined with reasonable accuracy by limiting the analysis to theﬁrst few hours of the simulation, in which the eﬀect of the ground israther limited. Considering only a short simulation period could
speed up the numerical optimization, even if the model is computationally expensive. Moreover, the independent estimation of the1123

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

Fig. 8. Scatter plots of the Morris sensitivity measures for various thermal parameters (λ is the thermal conductivity, C is the thermal capacity, and subscripts p, g, and s refer to HDPE,
grout, and soil, respectively) at two diﬀerent simulation times (left t = 10,800 s and right t = 30,000 s).

correlated with the outlet temperatures of the ﬂuid, indicating an increased heat dissipation in soil. This eﬀect is more pronounced for increasing aquifer thicknesses, which indicates that large portions of thesubsurface domain become saturated. In such circumstances, the soil
becomes more conductive (Eq. (12)) and the heat transfer between thepipe and the soil is enhanced. Since saturation signiﬁcantly aﬀects
thermal properties of soil, its eﬀect becomes evident as soon as the soil
is involved in the heat transfer process (t ≈ 10,000 s), similarly to what
was indicated in the previous Morris sensitivity analysis. At
t = 70,000 s, the largest deviation between simulated outlet temperatures is almost 1 °C and is kept almost constant for the remainder of thesimulation.
This behavior is shown in detail in Fig. 11, in which temperature
diﬀerences between the envelope curves (red and black lines in Fig. 10)
are shown. The temperature diﬀerence increases signiﬁcantly at the
beginning of the simulation, indicating a signiﬁcant inﬂuence of theaquifer thickness on the heat transfer process, and stabilizes aroundthermal properties of the grout could signiﬁcantly reduce the dimensionality of the subsequent optimization problem. This could befocused on the determination of the thermal properties of theground, which is the core of the thermal response test.
– The eﬀects of thermal properties of the ground on the model’s
output increased with time. It is important to temporally extend thenumerical simulation to identify the thermal properties of theground and limit the inﬂuence of the grout and the pipe. A shortsimulation time could lead to a biased estimation of the thermal
properties.

3.2.2. Groundwater eﬀects
The inﬂuence of groundwater on the TRT is shown in Fig. 10. In
particular, the upper and lower ﬁgures show the eﬀects of groundwater
ﬂuxes (0.0, 2.5 · 10−6, and 5 · 10−6 m/s) and aquifer thicknesses (10,20, and 30 m), respectively. At the ﬁrst inspection, it is evident thatboth aquifer thicknesses and groundwater ﬂuxes are negatively

Fig. 9. The ﬂuid temperature at the outlet for diﬀerent values of the grout (left plot) and ground (right plot) thermal conductivities. Blue arrows indicate increasing values of the thermal
conductivity. Dotted lines indicate the simulation times examined in the sensitivity analysis. (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the
web version of this article.)1124

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

Fig. 11. Temperature diﬀerences between the two envelope curves (red and black lines inFig. 6) for the aquifer thickness (solid line) and groundwater ﬂow (dashed line). (Forinterpretation of the references to color in this ﬁgure legend, the reader is referred to the
web version of this article.)exchanger. It must be noted that the newly-developed model can provide a comprehensive description of the heat transfer in the entire
system at a reduced computational cost compared to a fully 3D model.
3.2.3. Layered soil proﬁle
The results of simulated scenario evaluating the eﬀects of soil heterogeneity are presented in Fig. 13, which displays ﬂuid temperaturesin the U-tube for a homogeneous (lower plot in Fig. 13) and a layered(upper plot in Fig. 13) soil proﬁle. At the ﬁrst inspection, it is evident
that the homogeneous proﬁle, which is composed of sand, dissipates
more heat compared to the layered scenario. The temperature diﬀerence is higher than 1 °C almost everywhere along the pipe. While theﬂuid temperature proﬁle for the homogeneous soil is smooth, it exhibits
some nonlinearities for the layered scenario. In particular, at the interface between clay and sand, the drift from the theoretical temperature proﬁle for the homogeneous soil is signiﬁcant. At a vertical depthof 15 m, the sandy layer causes a small increase in the concavity of thetemperature curve, indicating an increase in the heat dissipation. On
the contrary, at a vertical depth of 25 m, the temperature curve exhibits
an inﬂection point due to the presence of the low-conductive clay layer,
which in turn lowers the dissipation capacity of the surrounding soil.
It must be emphasized that, for the sake of simplicity and computational eﬃciency, only 5 Finite Volumes were used to discretize thelithologic proﬁle. This example was meant to demonstrate the capabilities of the proposed modeling framework to handle layered proﬁles. Future applications can include the application of the model toreal geological settings, with ad hoc mesh discretizations.

Fig. 10. Fluid temperatures at the outlet for diﬀerent values of groundwater ﬂuxes (theupper plot) and aquifer thicknesses (the lower plot). Blue arrows indicate increasing
values of ﬂux and thickness, respectively. (For interpretation of the references to color inthis ﬁgure legend, the reader is referred to the web version of this article.)t = 70,000 when it is almost 1 °C. This behavior is signiﬁcantly different when examining the eﬀect of groundwater ﬂow. As shown in
Fig. 10, the eﬀect of groundwater ﬂow becomes visually appreciable
only after t = 20,000 s. The temperature diﬀerence increases almost
linearly with time. At t = 40,000 s, its value is below 0.2 °C, indicating
a limited inﬂuence on the system during the ﬁrst part of the test. At
t = 90,000 s, the temperature diﬀerence is about 0.4 °C. However, it
must be emphasized that contrary to what was observed for the aquifer
thickness, the temperature diﬀerence does not stabilize at this time,rather still exhibits an increasing trend. This indicates that the dissipation eﬀect of groundwater ﬂow is enhanced when the system becomes warmer.
The eﬀect of diﬀerent groundwater ﬂow velocities on the heattransfer process is shown in Fig. 12, which displays three HYDRUS-2D
cross sections for three diﬀerent water ﬂuxes. For q = 0 m/s, the distribution of temperatures is symmetric with a maximum value of
25.9 °C on the pipe surface. As the ﬂux increases, the temperature ﬁeld
becomes asymmetric with a considerable amount of heat following
saturated ﬂow. The convective heat transport dissipates a signiﬁcantamount of heat, thus lowering the temperature of the borehole heat

4. Conclusions and summary
The objective of the present study was to develop a computationallyeﬃcient pseudo-3D model for the analysis of Thermal Response Tests.
The proposed modeling framework combines a one-dimensional description of the heat transport in the buried tubes of the exchanger witha two-dimensional description of the heat transfer in the surrounding
subsurface soil, thus reducing the dimensionality of the problem andthe computational cost. One of the main novelties of the study is the1125

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

Fig. 12. HYDRUS-2D horizontal cross sections indicating the temperature distribution for three diﬀerent values of groundwater ﬂow at a depth of 45 m.

thermal properties of the ground. Multi-objective optimization algorithms and Monte Carlo uncertainty procedures represent valuable tools
for the inverse estimation of soil thermal properties, also in conjunction
with surrogate models, which have proven to be rather useful in otherengineering applications (e.g., Brunetti et al. [48]). Furthermore, acomputational benchmark between existing models for borehole heatexchangers is suggested in order to compare the eﬃciency of diﬀerent
modeling frameworks.

inclusion in the model of a widely used hydrological model, HYDRUS,
which can describe the simultaneous movement of water, heat, andsolutes in porous media. The proposed modeling approach is computationally eﬃcient and the use of HYDRUS adds several importantfeatures that can be exploited not only in TRT but also in other engineering applications. The analysis demonstrated that the proposedmodel can reproduce the thermal behavior of borehole heat exchangers
with good accuracy and that it can be exploited to investigate a variety
of processes occurring in the vadose zone that aﬀect the heat exchange(i.e., groundwater, ﬂow, variably-saturated conditions, inhomogeneous
lithologic proﬁles, etc.). Furthermore, the combination of the modelwith speciﬁc statistical techniques can clarify the inﬂuence of diﬀerent
factors on the thermal behavior of the borehole heat exchanger. In thisview, the application of the Morris method represents a step forward
compared to the traditional One-factor-at-a-time (OAT) techniques typically applied in this ﬁeld since it guarantees a global exploration of
the parameter space. Future work should investigate the use of diﬀerent
types of measurements for a more accurate determination of the

Acknowledgements
The experimental part of this research was supported by Core
Research for Evolutionary Science and Technology (CREST) from Japan
Science and Technology Agency (JST). All TRT data analyzed in thisstudy were collected with technical support from Mr. Eiji Yamashita of
Agricluster Co Ltd, Japan.

Fig. 13. Temperature proﬁles of the ﬂuid in the U-tube for the layered (the upper plot) and homogeneous (the lower plot) lithologic proﬁles. The red and black lines indicate thedescending and ascending legs of the U-tube, respectively.

1126

Applied Energy 208 (2017) 1113–1127

G. Brunetti et al.

of trace elements in shallow groundwater in Fuchu in western Tokyo Metropolis.
Japan. Environ Earth Sci 2016;75:559.
[27] Saito T, Hamamoto S, Ueki T, Ohkubo S, Moldrup P, Kawamoto K, et al.
Temperature change aﬀected groundwater quality in a conﬁned marine aquifer
during long-term heating and cooling. Water Res 2016;94:120–7.
[28] Saito T, Hamamoto S, Ei Mon E, Takemura T, Saito H, Komatsu T, et al. Thermal
properties of boring core samples from the Kanto area, Japan: development ofpredictive models for thermal conductivity and diﬀusivity. Soils Found
2014;54:116–25.
[29] Yavuzturk C, Spitler JD, Simon PE, Rees J. A transient two-dimensional ﬁnite volume model for the simulation of vertical U-Tube ground heat exchangers. ASHRAE
Trans 1999;105(2):465–74.
[30] Austin W, Yavuzturk C, Spitler JD. Development of an in-situ system for measuring
ground thermal properties. ASHRAE Trans 2000;106(1):365–79.
[31] Gustafsson AM, Gehlin S. Inﬂuence of natural convection in water-ﬁlled boreholes
for GCHP. ASHRAE Trans 2008;114(PART 1):416–23.
[32] Spitler JD, Javed S, Ramstad RK. Natural convection in groundwater-ﬁlled boreholes used as ground heat exchangers. Appl Energy 2016;164:352–65.
[33] Choi W, Ooka R. Eﬀect of natural convection on thermal response test conducted insaturated porous formation: comparison of gravel-backﬁlled and cement-grouted
borehole heat exchangers. Renew Energy 2016;96:891–903.
[34] Hellström G, Tsang C-F, Claesson J. Buoyancy ﬂow at a two-ﬂuid interface in a
porous medium: analytical studies. Water Resour Res 1988;24:493–506.
[35] Leonard BP. A stable and accurate convective modelling procedure based on
quadratic upstream interpolation. Comput Methods Appl Mech Eng 1979;19:59–98.
[36] Chung S-O, Horton R. Soil heat and water ﬂow with a partial surface mulch. Water
Resour Res 1987;23:2175–86.
[37] Beier RA, Smith MD, Spitler JD. Reference data sets for vertical borehole groundheat exchanger models and thermal response test analysis. Geothermics
2011;40:79–85.
[38] Zhang C, Guo Z, Liu Y, Cong X, Peng D. A review on thermal response test of
ground-coupled heat pump systems. Renew Sustain Energy Rev 2014;40:851–67.
[39] Saito H, Saito T, Kohgo Y, Hamamoto S, Moldrup P, Komatsu T. Eﬀect of reducing
the test period on estimating eﬀective thermal conductivity from Kelvin’s line
source function in thermal response test: a numerical study. J Jpn Soc Soil Phys
2014;128:11–20.
[40] Florides GA, Christodoulides P, Pouloupatis P. An analysis of heat ﬂow through aborehole heat exchanger validated model. Appl Energy 2012;92:523–33.
[41] Wagner V, Bayer P, Kübert M, Blum P. Numerical sensitivity study of thermal response tests. Renew Energy 2012;41:245–53.
[42] Campolongo F, Cariboni J, Saltelli A. An eﬀective screening design for sensitivity
analysis of large models. Environ Model Softw 2007;22:1509–18.
[43] Gehlin SEA, Hellström G. Inﬂuence on thermal response test by groundwater ﬂow invertical fractures in hard rock. Renew Energy 2003;28:2221–38.
[44] Raymond J, Therrien R, Gosselin L, Lefebvre R. Numerical analysis of thermal response tests with a groundwater ﬂow and heat transfer model. Renew Energy
2011;36:315–24.
[45] Fujii H, Okubo H, Nishi K, Itoi R, Ohyama K, Shibata K. An improved thermal
response test for U-tube ground heat exchanger based on optical ﬁber thermometers. Geothermics 2009;38:399–406.
[46] Fujii H, Okubo H, Itoi R. Thermal response tests using optical ﬁber thermometers.
Geotherm Resour Counc Trans 2006;30:545–51.
[47] Acuña J, Palm B. Distributed thermal response tests on pipe-in-pipe borehole heatexchangers. Appl Energy 2013;109:312–20.
[48] Brunetti G, Šimůnek J, Turco M, Piro P. On the use of surrogate-based modeling forthe numerical analysis of low impact development techniques. J Hydrol
2017;548:263–77.

References[1] Lund JW, Boyd TL. Direct utilization of geothermal energy 2015 worldwide review.
Geothermics 2016;60:66–93.
[2] Luo J, Rohn J, Bayer M, Priess A, Xiang W. Analysis on performance of boreholeheat exchanger in a layered subsurface. Appl Energy 2014;123:55–65.
[3] Gehlin S. Thermal Response Test. Lulea University of Technology; 2002.
[4] Spitler JD, Gehlin SEA. Thermal response testing for ground source heat pumpsystems—an historical review. Renew Sustain Energy Rev 2015;50:1125–37.
[5] Carslaw HS, Jaeger JC. Conduction of heat in solids. 2nd ed. Oxford Clarendon
Press; 1959.
[6] Ingersoll LR, Adler FT, Plass HJ, Ingersoll AC. Heat conduction with engineering,
geological, and other applications. New York: McGraw-Hill; 1954.
[7] Wagner R, Clauser C. Evaluating thermal response tests using parameter estimation
for thermal conductivity and thermal capacity. J Geophys Eng 2005;2:349–56.
[8] Bandos TV, Montero Á, Fernández de Córdoba P, Urchueguía JF. Improving parameter estimates obtained from thermal response tests: eﬀect of ambient air temperature variations. Geothermics 2011;40:136–43.
[9] Hu P, Meng Q, Sun Q, Zhu N, Guan C. A method and case study of thermal response
test with unstable heat rate. Energy Build 2012;48:199–205.
[10] Liuzzo-Scorpo A, Nordell B, Gehlin S. Inﬂuence of regional groundwater ﬂow on
ground temperature around heat extraction boreholes. Geothermics
2015;56:119–27.
[11] Wang H, Qi C, Du H, Gu J. Thermal performance of borehole heat exchanger undergroundwater ﬂow: a case study from Baoding. Energy Build 2009;41:1368–73.
[12] Hu J. An improved analytical model for vertical borehole ground heat exchanger
with multiple-layer substrates and groundwater ﬂow. Appl Energy
2017;202:537–49.
[13] Liuzzo Scorpo A, Nordell B, Gehlin S. A method to estimate the hydraulic conductivity of the ground by TRT analysis. Groundwater 2017;55:110–8.
[14] Rainieri S, Bozzoli F, Pagliarini G. Modeling approaches applied to the thermal
response test: a critical review of the literature. HVAC & R Res 2011;17:977–90.
[15] Shonder JA, Beck JV. A new method to determine the thermal properties of soil
formations from in situ ﬁeld tests; 2000.
[16] Florides GA, Christodoulides P, Pouloupatis P. Single and double U-tube groundheat exchangers in multiple-layer substrates. Appl Energy 2013;102:364–73.
[17] Ozudogru TY, Olgun CG, Senol A. 3D numerical modeling of vertical geothermalheat exchangers. Geothermics 2014;51:312–24.
[18] Han C, Yu X. (Bill). Sensitivity analysis of a vertical geothermal heat pump system.
Appl Energy 2016;170:148–60.
[19] Signorelli S, Bassetti S, Pahud D, Kohl T. Numerical evaluation of thermal response
tests. Geothermics 2007;36:141–66.
[20] Kohl T, Hopkirk RJ. “FRACTure” – a simulation code for forced ﬂuid ﬂow andtransport in fractured, porous rock. Geothermics 1995;24:333–43.
[21] Bozzoli F, Pagliarini G, Rainieri S, Schiavi L. Estimation of soil and grout thermal
properties through a TSPEP (two-step parameter estimation procedure) applied toTRT (thermal response test) data. Energy 2011;36:839–46.
[22] Christodoulides P, Florides GA, Pouloupatis P, Messaritis V, Lazari L. Ground heatexchanger modeling developed for energy ﬂows of an incompressible ﬂuid. Int J
Mech Aerospace, Inductrial, Mechatron Manuf Eng 2012;6:640–4.
[23] Kim E-JJ, Roux J-JJ, Rusaouen G, Kuznik F. Numerical modelling of geothermalvertical heat exchangers for the short time analysis using the state model size reduction technique. Appl Therm Eng 2010;30:706–14.
[24] Šimůnek J, van Genuchten MT, Šejna M. Recent developments and applications of
the HYDRUS computer software packages. Vadose Zo J 2016;15:25.
[25] Morris MD. Factorial sampling plans for preliminary computational experiments.
Technometrics 1991;33:161–74.
[26] Thuyet DQ, Saito H, Saito T, Moritani S, Kohgo Y, Komatsu T. Multivariate analysis

1127