International Journal of Heat and Mass Transfer 91 (2015) 1119–1127

Contents lists available at ScienceDirect

International Journal of Heat and Mass Transfer

journal homepage: www.elsevier.com/locate/ijhmt

The effects of borehole thermal resistances and fluid flow rate

on the g-functions of geothermal bore fields
Massimo Cimmino ⇑
Department of Mechanical Engineering, McGill University, 817 Sherbrooke St. W., Montreal, Quebec H3A 0C3, Canada

a r t i c l e i n f o a b s t r a c t

Article history: A model to calculate g-functions accounting for the axial variation of borehole wall temperatures and
Received 31 March 2015 heat extraction rates due to thermal interaction between the fluid and the borehole wall is presented.
Received in revised form 12 August 2015 Boreholes are divided into segments, with each borehole segment modeled as a finite line source.
Accepted 14 August 2015
Borehole wall temperatures are calculated from the spatial and temporal superposition of the finite line
Available online 3 September 2015
source solution. A delta-circuit of thermal resistances is used to consider the heat transfer between the
fluid and the borehole wall. A system of equations is built in the Laplace domain and solved for an equal
Keywords:
inlet fluid temperature for all boreholes in the bore field. The values of the g-function tend toward the
Boreholes
Geothermal
g-function obtained using a uniform heat extraction rate boundary condition as the borehole thermal
g-function resistance increases and tend toward the g-function obtained using a uniform temperature boundary
Thermal response factor condition as the borehole thermal resistance decreases. It is shown that modeling the thermal interaction
Ground heat exchangers between the fluid and the borehole wall is required to accurately predict the g-function.
Ó 2015 Elsevier Ltd. All rights reserved.

1. Introduction of the boreholes and equal for all boreholes. This assumption often
leads to the overestimation of Eskilson’s g-functions for bore fields
The design and simulation of geothermal heat pump systems with a large amount of closely spaced boreholes [5]. Methods have
coupled to vertical boreholes relies on the prediction of the fluid been proposed to consider the variation of the heat extraction rates
and ground temperatures during operation. The simulation of the among boreholes [6,7] and along the length of the boreholes to
ground temperatures in a geothermal bore field can be achieved obtain a uniform borehole wall temperature [8].
using thermal response factors, also known as g-functions. Ther- Although Eskilson’s g-functions are often used as reference for
mal response factors give the time variation of the borehole wall judging the accuracy of thermal response factors for geothermal
temperatures in a bore field for a constant rate of heat extracted bore fields, the validity of the uniform temperature boundary
from the ground. condition – or that of the uniform heat extraction rate boundary
The concept of g-functions was initially introduced by Eskilson condition – has yet to be validated. In practice, for a field of bore-
[1]. Eskilson’s g-functions were calculated numerically, under the holes connected in parallel, boreholes receive a heat carrier fluid at
assumption that the temperature along the wall of the boreholes a common inlet temperature and both the borehole wall tempera-
is uniform and equal for all boreholes in a bore field. The basis of tures and the extraction rates vary along the length of the bore-
this assumption is that the difference between the inlet and outlet holes. Fig. 1 shows a representation of the borehole wall
fluid temperatures in the boreholes is sufficiently small so that the temperature and heat extraction rate profiles along the length of
temperature variation along the length of the boreholes is also one borehole for the uniform heat extraction rate boundary condi-
small. Only the thermal properties of the ground surrounding the tion (UQ), the uniform temperature boundary condition (UT), and
boreholes and the dimensions of the boreholes are taken into the case were both the heat extraction rate and the borehole wall
account for the calculation of the g-functions. temperature vary along the length of the borehole due to the ther-
The finite line source analytical solution has been proposed to mal interaction between the fluid and the borehole wall (present
approximate the g-functions [2–4]. The finite line source solution model).
assumes that the heat extraction rate is uniform along the length This paper provides an extension to the model of Cimmino and
Bernier [8] to calculate the g-functions of geothermal bore fields,
using an equal inlet fluid temperature for all boreholes. The
⇑ Tel.: +1 438 275 5887. method accounts for the thermal interaction between the fluid
E-mail address: massimo.cimmino@polymtl.ca and the borehole wall and incorporates the borehole thermal

http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.08.041
0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.
1120 M. Cimmino / International Journal of Heat and Mass Transfer 91 (2015) 1119–1127

numerical simulations. Each borehole was modeled in a 2D radial-

axial grid. Spatial superposition of the temperature fields around
the boreholes was used to calculate the total temperature variation
at the borehole walls. A constant temperature was maintained at
the ground surface by superimposing the temperature fields of
image boreholes above the ground surface. At each time step of
the simulations, the heat extraction rates required to obtain a uni-
form temperature at the borehole walls, equal for all boreholes,
were evaluated. Although the interaction between the fluid loop
and the borehole wall was included in Eskilson’s Superposition
Borehole Model [20], its effect on the g-functions was not assessed.
The concept of g-functions was reviewed by Cimmino et al. [6] and
Bernier [21]. The definition of the g-function is given by [1]:

Q
Tb ¼ Tg
 gðt=t s ; r b =H; B=H; D=HÞ ð1Þ
2pks

where Tb is the borehole wall temperature, Tg is the undisturbed

ground temperature, Q
is the average heat extraction rate per unit
borehole length, ks is the ground thermal conductivity and g is the
g-function. The g-function varies according to the non-
dimensional time t/ts, where ts = H2/9as is the characteristic time
of the bore field, H is the borehole length and as is the ground ther-
mal conductivity; the non-dimensional borehole radius rb/H; the
non-dimensional borehole spacing B/H; and the non-dimensional
buried depth D/H.
The g-functions of a field of 4  4 boreholes is shown on Fig. 2.
Fig. 1. Schematic representation of the temperature and heat extraction rate
profiles. g-Function curves increase with time, indicating that the borehole
wall temperature drops with the extraction of heat from the
ground. A constant value is reached at a time ln(t/ts)  2.3 when
resistance, the short-circuit thermal resistance between the down- heat conduction between the boreholes, the ground and the
ward and upward U-tube pipes and the fluid flow rate into the ground surface reaches steady-state conditions. The value of the
evaluation of the g-functions, using the analytical solution of Eskil- g-function increases when the spacing between boreholes
son [1] and Hellström [9] for the fluid temperature variation in the decreases. The g-function of a field of infinitely spaced boreholes
boreholes. The method is aimed at the calculation of long-term (i.e. B/H ? 1) is equal to the g-function of a single borehole.
g-functions and neglects the thermal capacity of the borehole as g-Function curves are usually presented for single values of the
well as the travel time of the fluid through the boreholes. non-dimensional borehole radius and buried depth.
g-Functions obtained using this new methodology are compared The finite line source solution, obtained by the spatial integra-
to the g-functions obtained using the uniform temperature and tion of the point source solution, was proposed by Eskilson to
uniform heat extraction rate boundary conditions. It is shown that approximate the g-functions. Zeng et al. [2] calculated the average
the borehole thermal resistance affects the value of the g-functions temperature variation along the borehole lengths using the finite
and that the use of the uniform temperature and uniform heat line source solution. Simplified solutions for the evaluation of the
extraction rate boundary conditions can both lead to errors when average temperature using the finite line source solution were pro-
compared to the g-functions obtained with the presented posed by Lamarche and Beauchamp [3] for the case D = 0 and
methodology. Claesson and Javed [4] for the case D P 0. However, due to the
difference in the boundary condition at the borehole walls used
to calculate the g-functions, the finite line source solution was
2. Literature review shown to overestimate Eskilson’s g-functions when thermal
interaction between boreholes becomes significant [5].
Early contributions to the calculation of thermal response fac-
tors of geothermal heat exchangers were made by Ingersoll and
Plass [10] and Ingersoll et al. [11,12]. The infinite line source ana-
lytical solution and Carslaw and Jaeger’s [13] cylindrical heat
source solution were used to predict the ground temperature vari-
ations. Series solutions and correlations for the evaluation of the
cylindrical heat source solution were proposed by Beaudoin [14],
Cooper [15], Bernier [16,17] and Bernier and Salim-Shirazi [18].
One-dimensional models such as the infinite line source and the
cylindrical heat source do not account for axial conduction in
the ground caused by the influence of the ground surface and the
ground below the boreholes. The upper limit of validity of these
models is t = ts/10 (ln(t/ts) =
2.3) according to Eskilson [1]. At this
time, the error attributed to the use of one-dimensional models is
about 3% [19].
The pioneering work of Eskilson [1] led to the development of
g-functions. Eskilson’s g-functions were obtained by finite difference Fig. 2. g-Functions of a field of 4  4 boreholes.
 M. Cimmino / International Journal of Heat and Mass Transfer 91 (2015) 1119–1127 1121

A semi-analytical method to calculate the g-functions was The g-function of the bore field is obtained by calculating the
developed by Cimmino and Bernier [8]. Each borehole was divided time variation of the average borehole wall temperature T b caused
into a series of finite line source segments. Spatial superposition by the simultaneous extraction of heat by all boreholes at an aver-
and temporal superposition was used to build a system of
(Eq. (1)).
age heat extraction rate per unit length Q
equations in the Laplace domain. The system is solved to obtain
a uniform temperature at the borehole walls, equal to the 3.1. Finite line source
g-function. The g-functions obtained with the method were
compared to g-functions obtained with a uniform heat extraction Each borehole is divided into nq segments and each borehole
rate over the length of the boreholes. The differences between segment is modeled by a line source of finite length. The tempera-
the g-functions obtained using the two methods increased with ture variation at the wall of a borehole segment caused by the
the number of boreholes in the bore fields. The authors developed extraction of heat from another segment is given by the finite line
a preprocessor that calculates the g-functions of bore fields based source solution [8]. For a constant heat extraction rate at the uth
on user inputs of the borehole individual positions, lengths and segment of the ith borehole:
buried depths [7]. The g-function of a single borehole was
Q i;u
validated by experiments on a small-scale geothermal borehole DT i!j;u!v ðtk Þ ¼
 hi!j;u!v ðtk Þ ð2Þ
in laboratory conditions [22]. 2pks
Zanchini and Lazzari [23] developed a 2D axisymmetric finite Z
2 2 
1
1 1
element model to calculate g-functions for geothermal bore fields. hi!j;u!v ðtÞ ¼ pffiffiffiffiffiffiffi s2 expð
dij s Þ ierf ððDj;v
Di;u þ Hj;v ÞsÞ
g-Functions are calculated for a uniform heat extraction rate along 2Hj;v 1= 4as t

the length of the boreholes. Monzó et al. [24] developed a 3D finite
ierf ððDj;v
Di;u ÞsÞ þ ierf ððDj;v
Di;u
Hi;u ÞsÞ
element model to calculate the g-functions of geothermal bore
ierf ððDj;v
Di;u þ Hj;v
Hi;u ÞsÞ
fields. Boreholes are modeled as cylinders of highly conductive
þ ierf ððDj;v þ Di;u þ Hj;v ÞsÞ
ierf ððDj;v þ Di;u ÞsÞ
material connected above the ground surface. A constant heat is
injected into the ground through the highly conductive material. þ ierf ððDj;v þ Di;u þ Hi;u ÞsÞ

Due to the high thermal conductivity, the temperature tends
ierf ððDj;v þ Di;u þ Hj;v þ Hi;u ÞsÞ ds; ð3Þ
towards uniformity at the borehole walls and the model replicates
Z X
Eskilson’s g-functions accurately. 0 1
ierf ðXÞ ¼ erf ðx0 Þdx ¼ Xerf ðXÞ
pffiffiffiffi ð1
expð
X 2 ÞÞ; ð4Þ
Other solutions have been proposed for inclined boreholes 0 p
[25–27], boreholes connected in series [28], boreholes in multiple
(
layers of ground with unequal thermal properties [29], geothermal rb for i ¼ j
piles [30–36] and boreholes under the influence of groundwater dij ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ð5Þ
2 2
flow [37–40].
ðxi
xj Þ þ ðyi
yj Þ for i – j

where DT i!j;u!v is the temperature variation at the wall of the vth

3. Model segment of the jth borehole caused by the extraction of heat at
the wall of the uth segment of the ith borehole at a constant rate
A field of nb = 3 boreholes is shown on Fig. 3. All the boreholes Qi,u, dij is the distance separating borehole i from borehole j located
have the same length H, the same radius rb and are buried at a dis- at coordinates (xi,yi) and (xj,yj), hi!j;u!v is the segment-to-segment
tance D from the ground surface. Each borehole i is positioned at response factor of the uth segment of the ith borehole on the vth
coordinates (xi,yi). Boreholes are connected in parallel and are fed segment of the jth borehole given by the finite line source solution,
the heat carrier fluid at the same inlet temperature Tf,in and mass Hi;u ¼ nHq and Di;u ¼ D þ ðu
1ÞH
nq
are the length and buried depth of the
flow rate m._ The heat carrier fluid exits each borehole i at a temper- uth segment of the ith borehole.
ature Tf,out,i. The ground is homogenous and isotropic, has a thermal The total temperature variation at the wall of a borehole segment
conductivity ks, a thermal diffusivity as and is initially at a uniform is obtained by the spatial and temporal superposition of the finite
temperature Tg. line source solution for all borehole segments in the bore field:
X k X nb Xnq
Q i;u ðtp Þ
DT b;j;v ðtk Þ ¼
Dhi!j;u!v ðt k
t p
1 Þ ð6Þ
p¼1 i¼1 u¼1
2pks

where DT b;j;v ¼ T b;j;v
T g is the total temperature variation at the

wall of the vth segment of the jth borehole, Q i;u ðtp Þ is the heat
extraction rate per unit length of the uth segment of the ith
borehole from t ¼ t p
1 to t ¼ tp , Dhi!j;u!v ðtk Þ ¼ hi!j;u!v ðt k Þ

hi!j;u!v ðt k
1 Þ is the segment-to-segment response factor increment,
nb is the number of boreholes in the bore field and nq is the number
of segments per borehole.
As per the definition of the g-function (Eq. (1)), the total heat
extraction rate in the bore field is constant:
Pnb Pnq
u¼1 Q i;u ðt k Þ
:
i¼1
¼Q ð7Þ
nb  nq
Eq. (6) can be evaluated for all segments of all boreholes. Eqs.
(6) and (7) form a system of nb  nq þ 1 equations with 2  nb  nq
unknown variables (Q i;u and DT b;j;v ). The system of equations is
completed by considering the interactions between the fluid loop
Fig. 3. Field of 3 boreholes. and the borehole walls.
1122 M. Cimmino / International Journal of Heat and Mass Transfer 91 (2015) 1119–1127

qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi  
3.2. Internal resistances 2 b0 þb0
ðb01 þ b02 Þ =4 þ b012 ðb01 þ b02 Þ, d ¼ c10 b012 þ 1 2 2 and z is the axial

The borehole cross-section is modeled as a delta-circuit of ther- distance from the top of the borehole.
mal resistances between the downward and upward moving heat The outlet fluid temperature is given by:
carrier fluid and the borehole wall as shown on Fig. 4. An energy Z
f 1 ð1Þ þ f 2 ð1Þ 1 H

balance on the downward and upward moving heat carrier fluid T f 2;i ð0; tÞ ¼ T f 1;i ð0; tÞ  þ T b;i ðz0 ; tÞ
f 3 ð1Þ
f 2 ð1Þ H 0
0   
yields: 0 0 0
f 1
zH þ f 5 1
zH 0
@T f 1;i ðT b
T f 1;i Þ ðT f 2;i
T f 1;i Þ  4 dz : ð17Þ
_ p
mc ¼ þ ; ð8Þ f 3 ð1Þ
f 2 ð1Þ
@z RD1 RD12
Assuming a uniform wall temperature along the length of each
borehole segment, Eqs (10), (11) and (17) become:
@T f 2;i ðT b
T f 2;i Þ ðT f 1;i
T f 2;i Þ    
_ p

mc ¼ þ ; ð9Þ u u
@z RD2 RD12 T f 1;i;u ðt k Þ ¼ T f 1;i;0 ðt k Þ  f 1 þ T f 2;i;0 ðt k Þ  f 2
nq nq
Z ðu
v þ1ÞnH  00 
where T f 1;i and T f 2;i are the temperatures of the downward and 1X u
q 0 z 00
upward moving fluid in the ith borehole, RD1 and RD2 are the delta- þ T b;i;v ðt k Þ f4 dz ; ð18Þ
H v ¼1 ðu
v ÞnH H
q
circuit thermal resistances between the downward and upward
moving heat carrier fluid and the borehole wall, RD12 is the delta-    
u u
circuit short-circuit thermal resistance between the downward T f 2;i;u ðt k Þ ¼
T f 1;i;0 ðtk Þ  f 2 þ T f 2;i;0 ðt k Þ  f 3
nq nq
and upward moving heat carrier fluid, m _ is the fluid mass flow rate Z ðu
v þ1ÞnH  00 
and cp is the specific heat of the heat carrier fluid. 1X u
q 0 z 00

T b;i;v ðt k Þ f5 dz ð19Þ
The solution to Eqs. (8) and (9) for an arbitrary borehole wall H v ¼1 ðu
v Þn
H H
q
temperature profile is obtained using Laplace transforms and is
given by Eskilson [1] and Hellström [9]: f 1 ð1Þ þ f 2 ð1Þ
z z T f 2;i;0 ðt k Þ ¼ T f 1;i;0 ðt k Þ 
f 3 ð1Þ
f 2 ð1Þ
T f 1;i ðz; tÞ ¼ T f 1;i ð0; tÞ  f 1 þ T f 2;i ð0; tÞ  f 2 Z nq
v þ1  0  0
Z z
H   H 1 Xnq
T b;i;v ðt k Þ nq H
0 z 0 z 00
1 0 z z0 0
þ f4 þ f5 dz ; ð20Þ
þ 0
T b;i ðz ; tÞ  f 4
dz ð10Þ H v ¼1 f 3 ð1Þ
f 2 ð1Þ nq
v
nq H
H H
H 0 H H
where T f 1;i;u and T f 2;i;u are the temperatures of the downward and
z z
T f 2;i ðz; tÞ ¼
T f 1;i ð0; tÞ  f 2 þ T f 2;i ð0; tÞ  f 3 upward moving heat carrier fluid at the bottom of the uth segment
H   H of the ith borehole, T f 1;i;0 and are the inlet and outlet fluid temper-
Z z
1 0 z z0 0 atures of borehole i and T b;i;v is the uniform wall temperature along

T b;i ðz0 ; tÞ  f 5
dz ð11Þ
H 0 H H the length of the vth segment of the ith borehole. A schematic rep-
resentation of a borehole i along with the temperatures of the fluid
z  z h  z  z i
f1 ¼ exp b0 cosh c0
d sinh c0 ; ð12Þ and the borehole wall and the heat extraction rate at each borehole
H H H H segment is presented on Fig. 5.
The heat extraction rate per unit length of the uth segment of
z  z  b0  z
f2 ¼ exp b0 12
sinh c0 ; ð13Þ the ith borehole is given by an energy balance on the heat carrier
H H c 0 H fluid:
z  z h  z  z i _ p
mc
Q i;u ðt k Þ ¼ ðT f 1;i;u ðtk Þ
T f 1;i;u
1 ðt k Þ þ T f 2;i;u
1 ðt k Þ
f3 ¼ exp b0 cosh c0 þ d sinh c0 ; ð14Þ H=nq
H H H H

T f 2;i;u ðtk ÞÞ: ð21Þ
z  z   z  b0 b0
  z 

0
f4 ¼ exp b0 b01 cosh c0
db01 þ 2 0 12 sinh c0 ; From Eqs. (18)–(21):
H H H c H     
Q i;u ðt k Þ H u u
1
ð15Þ ¼ T f 1;i;0 ðt k Þ X 1
X1
mc_ p nq nq nq
    

z  z   z  b0 b0
  z 
X 1 ð1Þ u u
1
0
f5 ¼ exp b0 b02 cosh c0 þ db02 þ 1 0 12 sinh c0 ;
X2
X2
H H H c H X 2 ð1Þ nq nq
   

ð16Þ
X u
u
v þ1 u
v
þ T b;i;v ðt k Þ X 3
X3
v ¼1 nq nq
where b01 ¼ mcH
, b02 ¼ mcH
and b012 ¼ mc H
are the non-dimensional    

_ RD p 1 _ RD p 2 _ RD p 12 X
u
1
u
v u
v
1
b02
b01
T b;i;v ðt k Þ X 3
X3
thermal conductances of the borehole, b0 ¼ , c0 ¼ nq nq
2 v ¼1
   
X 2 nuq
X 2 u
1 nq Xnq

T b;i;v ðtk Þ
X 2 ð1Þ v ¼1
   

nq
v þ 1 nq
v
 X3
X3 ; ð22Þ
nq nq
z z z
X1 ¼ f1 þ f2
H H H
 z   z  b0 þ b0  z 

¼ exp b0 cosh c0
1 0 2 sinh c0 ; ð23Þ

Fig. 4. Internal resistances of a borehole.
H H 2c H
 M. Cimmino / International Journal of Heat and Mass Transfer 91 (2015) 1119–1127 1123

tk ¼ t k
Dt is needed prior to the evaluation of the functions

expð
rtÞ and expðrtÞ during the direct and inverse Laplace trans-
forms. The damping coefficient r is calculated according to Wede-
pohl’s criterion [43]:

lnðNt Þ
r¼2 ; ð30Þ
t max
where Nt is the total number of time steps and tmax is the maximum
value of the time variable t  .
Eqs. (6), (7) and (22) are rewritten in non-dimensional form in
the Laplace domain:
Fig. 5. Fluid temperatures, borehole wall temperatures and heat extraction rates of
X
nb X
nq
a borehole with nq = 3 segments.
LðHb;j;v Þ ¼ ~ i;u Þ  LðDhi!j;u!v Þ;
LðQ ð31Þ
i¼1 u¼1
z z
z
X2 ¼ f3
f2 Pnb Pnq ~
H H H u¼1 LðQ i;u Þ
 z   z 
¼ Lð1Þ; ð32Þ
i¼1
 z  b0 þ b0 nb  nq
0
¼ exp b cosh c0 þ 1 0 2 sinh c0 ; ð24Þ
H H 2c H     
~ i;u Þ 2pks H u u
1
z 1 Z h z  z i LðQ
_ p
¼ LðHf 1;0 Þ X 1
X1
mc nq nq
X3 ¼ f4 þ f5 dz     

H H H H X 1 ð1Þ u u
1
 z  b0 þ b0  z
X2
X2
X 2 ð1Þ nq nq
¼ exp b0 1 2
sinh c0 ð25Þ    

H c0 H X u
u
v þ1 u
v
þ LðHb;i;v Þ X 3
X3
Eq. (22) can be evaluated for all segments of all boreholes. Eqs. v ¼1 nq nq
(6), (7) and (22) form a system of 2  nb  nq þ 1 equations with    

X
u
1
u
v u
v
1
2  nb  nq þ nb unknown variables (Q i;u , DT b;j;v and T f 1;i;0 ). The sys-
LðHb;i;v Þ X 3
X3
v ¼1 nq nq
tem of equations is completed by imposing an equal inlet fluid    
temperature T f 1;i;0 ¼ T f ;in for all boreholes. X 2 nq
X 2 nq X
u u
1 nq

LðHb;i;v Þ
X 2 ð1Þ v ¼1
3.3. Non-dimensional equations in the Laplace domain    

nq
v þ 1 nq
v
 X3
X3 ; ð33Þ
The system of equations can be difficult to solve due to the nq nq
summation involving the time variable tp in Eq. (6). Marcotte and v T
T
where Hb;j;v ¼
b;j; g

=2pks is the non-dimensional temperature at the
Pasquier [41] suggested the use of the Fourier transform to solve Q

the convolution product between the heat extraction rate and ~ i;u ¼ Q =Q
wall of the vth segment of the jth borehole, Q
is the nor-
i;u
the response factor. Cimmino et al. [6] and Cimmino and Bernier malized heat extraction rate increment of the uth segment of the ith
[8,22] used the numerical Laplace transform to simplify the T
T
borehole and Hf 1;0 ¼
Qf
;in=2pkgs is the non-dimensional inlet fluid tem-
convolution product in Eq. (6). The Laplace transform pairs are
perature, equal for all boreholes.
expressed as:
The solution of the system of equations (Eqs. (31)–(33)) gives
Z 1
the non-dimensional wall temperatures, the normalized heat
FðsÞ ¼ Lðf ðtÞÞ ¼ f ðtÞ expð
stÞdt; ð26Þ
0 extraction rates of all borehole segments and the non-
dimensional inlet fluid temperature. The g-function of the bore
Z rþj1 field is given by the average of the non-dimensional temperatures
1
f ðtÞ ¼ L
1 ðFðsÞÞ ¼ FðsÞ expðstÞds; ð27Þ of all borehole segments:
2pj r
j1
Pnb Pnq
Hb;i;u ðt k Þ
where L and L
1 are the direct and inverse Laplace transforms, f is gðt k Þ ¼ i¼1 u¼1
: ð34Þ
an arbitrary function in the time domain and F the corresponding nb  nq
function in the Laplace domain, s is the complex frequency in the
pffiffiffiffiffiffiffi
Laplace domain and j ¼
1 is the imaginary number. 3.4. Other boundary conditions
The Laplace transform is obtained from the Fourier transform
using a variable change s = r + jx [42]: Two additional boundary conditions at the borehole walls are
Z 1 considered for comparison with the proposed model: the boundary
FðsÞ ¼ ½f ðtÞ expð
rtÞ expð
jxtÞdt ¼ F ðf ðtÞ  expð
rtÞÞ; ð28Þ condition of uniform and equal heat extraction rates (UQ) and the
0
boundary condition of uniform and equal temperature at the bore-
Z þ1 hole walls (UT).
expðrtÞ
f ðtÞ ¼ FðsÞ expðjxtÞdx ¼ expðrtÞ  F
1 ðFðsÞÞ; ð29Þ Boundary condition UQ is obtained by having the heat extrac-
2p
1
tion rates of all borehole segments equal to the average heat
where x is the angular frequency in the Fourier domain, F and F
1 extraction rate in the bore field:
are the direct and inverse Fourier transforms and r is a real positive ~ i;u Þ ¼ Lð1Þ:
LðQ ð35Þ
constant.
The Laplace transform is calculated numerically using a fast The solution to Eqs. (31) and (35) gives the non-dimensional
Fourier transform (FFT) algorithm. Since the heat extraction rates temperature variations at the wall of the borehole segments. The
Qi,u(tk) are only defined for t P t1 (t 1 ¼ Dt), a variable change g-function is then obtained from Eq. (34).
1124 M. Cimmino / International Journal of Heat and Mass Transfer 91 (2015) 1119–1127

The boundary condition UT is obtained by having the non-

dimensional temperature variation at the wall of borehole seg-
ments be equal for all segments of all boreholes:

Hb ðt k Þ ¼ Hb;j;v ðt k Þ: ð36Þ
The solution to Eqs. (31), (32) and (36) gives the non-
dimensional temperature variation at the wall of the borehole seg-
ments, equal to the g-function.

4. Results

The g-functions of a field of one (1  1) borehole and a field of

4  4 boreholes are calculated using the proposed methodology.
The effect of the variation of the thermal resistances of the bore-
holes and the fluid flow rate per borehole on the g-functions is first Fig. 7. Steady-state values of the g-functions of the 1  1 bore field.

studied. All results are calculated using nq = 24 segments.

The boreholes have a length H = 150 m, a radius rb = 0.075 m
The steady-state values of the g-function are shown to decrease
(rb/H = 0.0005) and are buried at a depth D = 4 m (D/H = 0.0267).
slightly with increasing values of the non-dimensional internal
The boreholes are equally spaced at a distance B = 7.5 m
conductance and decrease with decreasing values of the non-
(B/H = 0.05) on a rectangular grid. The ground thermal conductivity
dimensional borehole thermal resistance. The values are mostly
is ks = 2 W/m-K. The fluid flow rate is varied from m_ = 0.0625 kg/s
  confined within the values obtained with the UQ and the UT
per borehole 2pks H
_ p ¼ 7:54
mc
to m _ = 5.0 kg/s per borehole boundary conditions. The maximum difference between all values,
  including those obtained with the UQ and UT boundary conditions,
2pks H
_ p ¼ 0:094 . The fluid specific heat is cp = 4000 /kg-K. The
mc
is 1.1%. In the case of a bore field of a single borehole, the use of any
RD
borehole thermal resistance Rb ¼ 21 is varied from Rb = 0.001 m-K/W of the studied boundary conditions provide an accurate estimation
to Rb = 1000 m-K/W. The borehole internal resistance is varied of the g-function values.
from Ra = 0.10 m-K/W to Ra = 0.75 m-K/W. The short-circuit Fig. 8 shows the steady-state values of the g-functions of the
thermal resistance is given by [9]: 4  4 bore field. In this case, the non-dimensional internal thermal
conductance has negligible effect on the steady-state value of the
2Ra RD1 g-function. Once again, the steady-state values of the g-functions
RD12 ¼ : ð37Þ
2RD1
Ra decrease with decreasing values of the non-dimensional borehole
thermal resistance.
The range of variation of the g-functions of the 1  1 and 4  4
bore fields for the range of parameters indicated above are pre-
sented on Fig. 6. The g-function curves tend to follow the same 4.2. Borehole thermal resistance
shape and the differences are better studied by the steady-state
values of the g-functions. The non-dimensional borehole thermal The dependence of the g-function on the non-dimensional bore-
resistance Xb ¼ 2pks Rb , the non-dimensional internal thermal con- hole thermal resistance is investigated on Fig. 9. The steady-state
values of the g-function of the 4  4 bore field are calculated for
ductance b0a ¼ mc H 2pks H
_ p Ra and the non-dimensional fluid flow rate mc
_ p
a single value of the non-dimensional internal thermal conduc-
proved to be the better suited parameters for the analysis of the
tance b0a ¼ 0:300 and different values of the non-dimensional fluid
steady-state values. pks H
flow rate 2mc
_ p . It is shown that the steady-state value of the

g-function converges to that obtained with the UT boundary condi-

4.1. Internal short-circuit conductance
tion when the non-dimensional borehole thermal resistance
approaches zero and converges to that obtained with the UQ
Fig. 7 shows the steady-state values of the g-functions of the
boundary condition when the non-dimensional borehole thermal
1  1 bore field. The steady-state values are compared to
resistance approaches infinity. The effect of the non-dimensional
the steady-state values obtained with the UQ boundary
fluid flow rate is negligible.
condition (Eq. (35)) and the UT boundary condition (Eq. (36)).

Fig. 6. g-Functions of the 1  1 and 4  4 bore fields. Fig. 8. Steady-state values of the g-functions of the 4  4 bore field.
 M. Cimmino / International Journal of Heat and Mass Transfer 91 (2015) 1119–1127 1125

Fig. 9. Dependence of the g-function on the borehole thermal resistance and the
fluid flow rate.

The steady-state value of the g-function of the 4  4 bore field

pks H
_ p = 0.9425 is 25.82 when ln(Ob) =
4.38 and 29.69 when
with 2mc
ln(Ob) = 9.44, while the steady-state values are 25.78 and 29.69
for the UT and UQ boundary conditions. For a ground thermal con-
ductivity ks = 2 W/m-K and a borehole thermal resistance
Rb = 1.0 m-K/W (ln(Ob) = 2.53), the steady-state value of the
g-function is 28.02. Using the UQ boundary condition overesti-
mates this value by 6.0% and using the UT boundary condition
underestimates this value by 8.0%. For a lower borehole thermal
resistance Rb = 0.1 m-K/W (ln(Ob) = 0.23), the steady-state value Fig. 10. (a) Dependence of the g-function on the borehole thermal resistance and
of the g-function is 26.31. Using the UQ boundary condition over- the size of the bore field and (b) range of variation of the g-functions.
estimates this value by 12.8% and using the UT boundary condition
underestimates this value by 2.0%. This shows that there is an
important variation of the g-function for different values of the boreholes in the bore field increased. However, as shown by
non-dimensional borehole thermal resistance. The thermal inter- Cimmino et al. (2013) and Cimmino and Bernier (2014) and
action between the fluid loop and the borehole wall, b0a taken into illustrated on Fig. 10b, the difference between the g-function
account in the present model, should be modeled in the evaluation values evaluated with the UQ and UT conditions increases when
of the g-function. the number of boreholes increases.
For very high borehole thermal resistances, the difference For the 10  10 bore field, the steady-state values of the
between the fluid temperatures and the borehole wall tempera- g-function obtained with the UQ and UT boundary conditions are
tures is a lot larger than the axial variation of the fluid and bore- 92.80 and 59.99, respectively. For a ground thermal conductivity
hole wall temperatures and the heat extraction rates thus tend to ks = 2 W/m-K and a borehole thermal resistance Rb = 1.0 m-K/W
be close to uniform and equal for all boreholes. For very low bore- (ln(Ob) = 2.53), the steady-state value of the g-function is 73.74.
hole thermal resistances, the borehole wall temperatures tend to Using the UQ boundary condition overestimates this value by
be close to uniform and equal for all boreholes. The non- 25.8% and using the UT boundary condition underestimates this
dimensional fluid flow rate has a very small impact on the value by 18.6%. For a lower borehole thermal resistance
steady-state value of the g-function. Rb = 0.1 m-K/W (ln(Ob) = 0.23), the steady-state value of the
g-function is 62.71. Using the UQ boundary condition overestimates
4.3. Number of boreholes this value by 48.0% and using the UT boundary condition underes-
timates this value by 4.3%.
The normalized steady-state values of the g-functions of fields
of 4  4, 7  7, and 10  10 boreholes are shown on Fig. 10. The 4.4. Borehole spacing
steady-state values are presented for a single value of the non-
dimensional internal thermal conductance b0a ¼ 0:300 and of the The normalized steady-state values of the g-functions of a 4  4
pks H
non-dimensional fluid flow rate 2mc
_ p = 0.9425 as a function of the
bore field are shown on Fig. 11 for different borehole spacing to
non-dimensional borehole thermal resistance. length ratios B/H. The steady-state values are presented for a single
value of the non-dimensional internal thermal conductance
The normalized steady-state values of the g-functions are given
by: b0a = 0.300 and of the non-dimensional fluid flow rate
2pks H
_ p = 0.9425 as a function of the non-dimensional borehole ther-
gðt ! 1Þ
g UT ðt ! 1Þ mc
g~ðt ! 1Þ ¼ ; ð38Þ mal resistance. In this case, the steady-state value of the g-
g UQ ðt ! 1Þ
g UT ðt ! 1Þ
function gets closer to the steady-state value obtained with the
where g~ is the normalized g-function, g UT is the g-function calcu- UQ boundary condition as B/H increases.
lated with the uniform temperature boundary condition and g UQ
is the g-function calculated with the uniform heat extraction rate 4.5. Temperature and heat extraction rate profiles
boundary condition.
The steady-state values of the g-function are shown to be closer The temperature and heat extraction rate profiles along the
to that obtained with the UT boundary condition as the number of length of a corner borehole in a field of 4  4 boreholes for the 3
1126 M. Cimmino / International Journal of Heat and Mass Transfer 91 (2015) 1119–1127

Fig. 13. Variation of the heat extraction rates among boreholes in a 4  4 bore field.

Table 1
Time required for the calculation of g-functions.

Boundary condition Calculation time

11 44 77 10  10
borehole boreholes boreholes boreholes
Proposed model 1.33 s 19.3 s 2 min 42 s 13 min 10 s
Fig. 11. (a) Dependence of the g-function on the borehole thermal resistance and (nq = 24)
the borehole spacing and (b) range of variation of the g-functions. Uniform temperature 1.31 s 18.7 s 2 min 30 s 12 min 19 s
(nq = 24)
Uniform heat extraction 0.11 s 0.39 s 0.95 s 1.95 s
rate (nq = 1)
boundary conditions are shown on Fig. 12. The profiles are calcu-
lated for a non-dimensional internal thermal conductance
b0a = 0.300, a non-dimensional fluid flow rate 2mc
pks H
_ p = 0.9425, and
~ = 1). At steady-state, the normalized
to be closer to the average (Q
for non-dimensional borehole thermal resistances Ob = 1.26 and
Ob = 12.6. The temperature and heat extraction rate profiles are heat extraction rate of the corner boreholes is Q ~ = 1.31 and
closer to those obtained using the UT condition for a low non- ~ = 1.13 for Ob = 1.26 and Ob = 12.6, respectively, while the UT
Q
dimensional borehole thermal resistance and get closer to those ~ = 1.41.
boundary condition predicts Q
obtained using the UQ boundary condition when the non-
dimensional borehole thermal resistance increases. 4.6. Calculation time
Similar results are obtained for the time variation of the nor-
malized heat extraction rates of the boreholes. Fig. 13 shows the The time required for the calculation of the g-functions of the
normalized heat extraction rates of 2 boreholes in a field of 4  4 fields of 1  1, 4  4, 7  7 and 10  10 boreholes, using nq = 24
boreholes. It is shown that corner boreholes (i.e. borehole 1) have segments and 27 time steps, is presented on Table 1. Calculations
a larger heat extraction rate than center boreholes (i.e. borehole were done on a desktop PC equipped with an Intel i5-2500K
2). As the non-dimensional borehole thermal resistance increases, (3.30 GHz) processor and 24 GB of RAM. The calculation time for
the normalized heat extraction rates of individual boreholes tend the proposed model and the uniform temperature boundary condi-
tion are similar: the calculation time for the field of 10  10 bore-
holes was 13 min 10 s and 12 min 19 s for the proposed model and
the uniform temperature boundary condition, respectively.

5. Conclusion

An analytical model for the calculation of g-functions is pre-

sented. Boreholes are divided into segments and the temperatures
at the wall of the borehole segments are obtained from the finite
line source solution. The spatial and temporal superpositions are
used to build a system of linear equations with the temperature
and heat extraction rates at the wall of the borehole segments as
the unknowns. The system of equations is completed using an ana-
lytical model describing the heat transfer between the fluid and
the borehole wall. The system of equations is solved for an equal
inlet fluid temperature for all boreholes in the bore field. The solu-
Fig. 12. Temperature and heat extraction rate profiles along the length of a corner tion to the system of equations gives the g-functions of the bore
borehole in a 4  4 bore field. field.
 M. Cimmino / International Journal of Heat and Mass Transfer 91 (2015) 1119–1127 1127

The effects of the internal thermal resistances of the boreholes, [12] L.R. Ingersoll, O.J. Zobel, A.C. Ingersoll, Theory of earth heat exchanger for the
heat pump, in: Heat Conduction: With Engineering, Geological and Other
the number of boreholes and the borehole spacing on the
Applications, McGraw-Hill, 1954, pp. 240–271.
g-function are reported. The steady-state values of the g-functions [13] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, vol. 2, Clarendon Press,
obtained using the proposed model are compared to the steady- Oxford, 1959.
state values of the g-functions obtained using the uniform temper- [14] A. Baudoin, Stockage intersaisonnier de chaleur dans le sol batterie
d’échangeurs baionnette verticaux; modèle de prédimensionnement (Ph.D.
ature and uniform heat extraction rate boundary conditions. It is Thesis), Université de Reims Champagne-Ardenne, France, 1988.
shown that the differences between the proposed model and the [15] L.Y. Cooper, Heating of a cylindrical cavity, Int. J. Heat Mass Transfer 19 (1976)
other two boundary conditions are mostly influenced by the non- 575–577.
[16] M.A. Bernier, A review of the cylindrical heat source method for the design and
dimensional borehole thermal resistance. The uniform tempera- analysis of vertical ground-coupled heat pump systems, in: presented at the
ture and uniform heat extraction rate boundary conditions are only 4th International Conference on Heat Pumps in Cold Climates, 2000.
valid for very low and very high non-dimensional borehole thermal [17] M.A. Bernier, Ground-coupled heat pump system simulation, ASHRAE Trans.
107 (2001) 605–616.
resistances, respectively. For a field of 10  10 boreholes with [18] M.A. Bernier, A. Salim-Shirazi, Solar Heat Injection into Boreholes: a
ground conductivity ks = 2 W/m-K and borehole thermal resistance Preliminary Analysis, in: presented at the 2nd Canadian Solar Buildings
Rb = 1.0 m-K/W (Ob = 12.6), differences of 25.8% and 18.6% were Conference, Calgary, 2007.
[19] M. Philippe, M. Bernier, D. Marchio, Validity ranges of three analytical
observed between the steady-state value obtained using the pro- solutions to heat transfer in the vicinity of single boreholes, Geothermics 38
posed model and each of the uniform heat extraction rate and uni- (2009) 407–413.
form temperature boundary conditions, respectively. The observed [20] P. Eskilson, Superposition Borehole Model: Manual for Computer Code,
University of Lund, Lund, Sweden, 1986.
differences demonstrate that modeling the thermal interaction
[21] M. Bernier, Sizing and simulating bore fields using thermal response factors
between the fluid and the borehole wall is required to properly (Keynote presentation), in: presented at the 11th IEA 2014 Heat Pump
predict the long-term temperature variations in bore fields. Conference, Montreal, Canada, 2014.
The non-dimensional borehole thermal resistance was also [22] M. Cimmino, M. Bernier, Experimental determination of the g-functions of a
small-scale geothermal borehole, Geothermics 56 (2015) 60–71.
shown to influence the distribution of the heat extraction rates [23] E. Zanchini, S. Lazzari, Temperature distribution in a field of long borehole heat
among boreholes and the temperature and heat extraction rate exchangers (BHEs) subjected to a monthly averaged heat flux, Energy 59
profiles along the length of the boreholes. The time required for (2013) 570–580.
[24] P. Monzó, P. Mogensen, J. Acuña, F. Ruiz-Calvo, C. Montagud, A novel numerical
the calculation of the g-function using the proposed model is com- approach for imposing a temperature boundary condition at the borehole wall
parable to the calculation time using the uniform temperature in borehole fields, Geothermics 56 (2015) 35–44.
boundary condition. [25] P. Cui, H. Yang, Z. Fang, Heat transfer analysis of ground heat exchangers with
inclined boreholes, Appl. Therm. Eng. 26 (2006) 1169–1175.
[26] D. Marcotte, P. Pasquier, The effect of borehole inclination on fluid and ground
Conflict of interest temperature for GLHE systems, Geothermics 38 (2009) 392–398.
[27] L. Lamarche, Analytical g-function for inclined boreholes in ground-source
heat pump systems, Geothermics 40 (2011) 241–249.
None declared. [28] D. Marcotte, P. Pasquier, Unit-response function for ground heat exchanger
with parallel, series or mixed borehole arrangement, Renewable Energy 68
Acknowledgements (2014) 14–24.
[29] S.L. Abdelaziz, T.Y. Ozudogru, C.G. Olgun, J.R. Martin Ii, Multilayer finite line
source model for vertical heat exchangers, Geothermics 51 (2014) 406–416.
The author received a postdoctoral research scholarship from [30] Y. Man, H. Yang, N. Diao, J. Liu, Z. Fang, A new model and analytical solutions
the Fonds de recherche du Québec – Nature et Technologie (FRQNT). for borehole and pile ground heat exchangers, Int. J. Heat Mass Transfer 53
(2010) 2593–2601.
Part of this work was done while the author received scholarships [31] P. Cui, X. Li, Y. Man, Z. Fang, Heat transfer analysis of pile geothermal heat
from the Canadian GeoExchange Coalition and the Natural Sciences exchangers with spiral coils, Appl. Energy 88 (2011) 4113–4119.
and Engineering Research Council of Canada (NSERC) and a Grant- [32] M. Li, A.C.K. Lai, New temperature response functions (G functions) for pile and
borehole ground heat exchangers based on composite-medium line-source
in-aid award from the American Society of Heating, Refrigerating theory, Energy 38 (2012) 255–263.
and Air-Conditioning Engineers (ASHRAE). [33] M. Li, A.C.K. Lai, Heat-source solutions to heat conduction in anisotropic media
with application to pile and borehole ground heat exchangers, Appl. Energy 96
(2012) 451–458.
References
[34] F. Loveridge, W. Powrie, Temperature response functions (G-functions) for
single pile heat exchangers, Energy 57 (2013) 554–564.
[1] P. Eskilson, Thermal Analysis of Heat Extraction Boreholes (Ph.D. Thesis), Dep. [35] F. Loveridge, W. Powrie, G-Functions for multiple interacting pile heat
of Mathematical Physics, University of Lund, Lund, Sweden, 1987. exchangers, Energy 64 (2014) 747–757.
[2] H.Y. Zeng, N.R. Diao, Z.H. Fang, A finite line-source model for boreholes in [36] A. Leroy, M. Bernier, Development of a novel spiral coil ground heat
geothermal heat exchangers, Heat Transfer – Asian Res. 31 (2002) 558–567. exchanger model considering axial effects, Appl. Therm. Eng. 84 (2015) 409–
[3] L. Lamarche, B. Beauchamp, A new contribution to the finite line-source model 419.
for geothermal boreholes, Energy Build. 39 (2007) 188–198. [37] M.G. Sutton, D.W. Nutter, R.J. Couvillion, A ground resistance for vertical bore
[4] J. Claesson, S. Javed, An analytical method to calculate borehole fluid heat exchangers with groundwater flow, J. Energy Res. Technol. 125 (2003)
temperatures for time-scales from minutes to decades, ASHRAE Trans. 117 183–189.
(2011) 279–288. [38] N. Diao, Q. Li, Z. Fang, Heat transfer in ground heat exchangers with
[5] M. Fossa, The temperature penalty approach to the design of borehole heat groundwater advection, Int. J. Therm. Sci. 43 (2004) 1203–1211.
exchangers for heat pump applications, Energy Build. 43 (2011) 1473–1479. [39] N. Molina-Giraldo, P. Blum, K. Zhu, P. Bayer, Z. Fang, A moving finite line source
[6] M. Cimmino, M. Bernier, F. Adams, A contribution towards the determination model to simulate borehole heat exchangers with groundwater advection, Int.
of g-functions using the finite line source, Appl. Therm. Eng. 51 (2013) 401– J. Therm. Sci. 50 (2011) 2506–2513.
412. [40] M. Tye-Gingras, L. Gosselin, Generic ground response functions for ground
[7] M. Cimmino, M. Bernier, Preprocessor for the generation of g-functions used in exchangers in the presence of groundwater flow, Renewable Energy 72 (2014)
the simulation of geothermal systems, in: Presented at the Proceedings of 354–366.
BS2013, Chambéry, France, 2013. [41] D. Marcotte, P. Pasquier, Fast fluid and ground temperature computation for
[8] M. Cimmino, M. Bernier, A semi-analytical method to generate g-functions for geothermal ground-loop heat exchanger systems, Geothermics 37 (2008) 651–
geothermal bore fields, Int. J. Heat Mass Transfer 70 (2014) 641–650. 665.
[9] G. Hellstrom, Ground Heat Storage: Thermal Analysis of Duct Storage Systems [42] P. Moreno, A. Ramirez, Implementation of the numerical laplace transform: a
PhD, Department of Mathematical Physics, University of Lund, Lund, Sweden, review, IEEE Trans. Power Delivery 23 (2008) 2599–2609.
1991. [43] L.M. Wedepohl, Power system transients: errors incurred in the numerical
[10] L.R. Ingersoll, H.J. Plass, Theory of the ground pipe heat source for the heat inversion of the Laplace transforms, in: 26th Midwest Symposium on Circuits
pump, Heating Piping Air Conditioning 20 (1948) 119–122. and Systems, Mexico, 1983, pp. 174–178.
[11] L.R. Ingersoll, F.T. Adler, H.J. Plass, A.C. Ingersoll, Theory of earth heat
exchangers for the heat pump, Heating Piping Air Conditioning 22 (1950)
113–122.