Accepted Manuscript

Evaluation of a transient borehole heat exchanger model in dynamic simulation of a

ground source heat pump system

Hassan Biglarian, Madjid Abbaspour, Mohammad Hassan Saidi

PII: S0360-5442(18)30037-9
DOI: 10.1016/[Link].2018.01.031
Reference: EGY 12143

To appear in: Energy

Received Date: 13 April 2017

Revised Date: 6 December 2017
Accepted Date: 5 January 2018

Please cite this article as: Biglarian H, Abbaspour M, Saidi MH, Evaluation of a transient borehole heat
exchanger model in dynamic simulation of a ground source heat pump system, Energy (2018), doi:
10.1016/[Link].2018.01.031.

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to
our customers we are providing this early version of the manuscript. The manuscript will undergo
copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please
note that during the production process errors may be discovered which could affect the content, and all
legal disclaimers that apply to the journal pertain.
 ACCEPTED MANUSCRIPT
Evaluation of a transient borehole heat exchanger model in
dynamic simulation of a ground source heat pump system

Hassan Biglarian, Madjid Abbaspour, Mohammad Hassan Saidi*

PT
Center of Excellence in Energy Conversion (CEEC), School of Mechanical Engineering,
Sharif University of Technology, Tehran, Iran

RI
Abstract

SC
The performance of a vertical ground source heat pump system (GSHPS) largely depends on
the fluid temperature leaving the borehole heat exchanger (BHE) that may be affected by the
short-term behavior of the BHE. Although considerable research has been carried out to

U
analyze the short-term transient response of the BHEs, few studies have investigated its
AN
impact on dynamic simulation of GSHPS. Therefore, this paper presents a numerical
approach based on a transient BHE model to evaluate the performance of a residential
M

GSHPS over short and long timescales. The numerical results are compared with the results
of EnergyPlus software. It is shown that the proposed model can appropriately predict the
D

dynamic behavior of the system. Moreover, effect of borehole thermal capacity on the
performance of the GSHPS is investigated in comparison with a quasi-steady state model. It
TE

is found that including the borehole thermal capacity substantially affects the design borehole
length. Using the transient model instead of the quasi-steady state model leads to a 16%
EP

reduction in the required borehole length.

Keywords: Ground source heat pump; Borehole heat exchanger; Short-term response; Part-
C

load fraction; Transient model; Thermal capacity.

AC

* Corresponding author. Tel.: +98 21 66165522; fax: +98 21 66000021.

E-mail addresses: saman@[Link] (M.H. Saidi); biglarianh@[Link] (H. Biglarian).

1
 ACCEPTED MANUSCRIPT
Nomenclature
fluid (or liquid)
 coefficients of Eq. (21)  fan
 area (m2)  grout
 ,  and  coefficients of Eqs. (13)-(15) 
ground
 specific heat capacity (J kg-1 K-1)  middle grout
 thermal capacity (J m K ) -1 -1
 indoor
  coefficient of performance  coil inlet flow
half shank spacing (m)  index
 and  coefficients of Eqs. (16) and (17) ℎ heating

PT
 liquid loop heat balance error (W) ℎ heating coil
ℎ specific enthalpy (J kg-1)  hydraulic
 thermal conductivity (W m-1 K-1)  heat pump

RI
 mass flow rate (kg s ) -1
 motor
 maximum heat pump cycling rate (s-1)  outdoor
!

SC
fraction of on-cycle power use pipe
 power consumption (W)  pipe inner surface
" part-load fraction  pipe outer surface
" $ adjusted part-load fraction  pump

U
"% part-load ratio !&
reference
' heat transfer rate (W) (& sensible
AN
! radius (m) ()!& source side coil
% thermal resistance (mK W-1) (*( system
+ seasonal performance factor + shaft
, ,,
M

time (s) total

- temperature (ºC) ./ wet bulb
0 volumetric flow rate (m3 s-1) Superscripts
1 , /
D

depth (m) inner and outer resistances

Greek symbols
downward
3 thermal diffusivity (m2 s-1) 
TE

iteration index
∆ pressure drop (Pa) ) upward
∆, time step (s) 5 downward or upward
6 accuracy (ºC) Abbreviations
EP

7 efficiency BHE borehole heat exchanger

8 heat pump time constant (s) COP coefficient of performance
9 run time fraction EER energy efficiency ratio
C

Subscripts GSHPS ground source heat pump system

 air ITMY Iranian typical meteorological year
AC

/ borehole PLF part-load fraction

 building PLR part-load ratio
 cooling RC resistance-capacity
 cooling coil SPF seasonal performance factor

/ dry bulb TRCM thermal resistance-capacity model

2
 ACCEPTED MANUSCRIPT
1. Introduction
Buildings are responsible for a considerable proportion of the global total energy use and CO2
emissions [1]. Therefore, there are significant opportunities for the use of renewable energy
in building sector [2-6]. In recent years, ground source heat pump system (GSHPS) has
received increasing attention for space heating and cooling. Since the ground temperature is
almost constant throughout the year, it can provide stable conditions for the heat pump

PT
operation. The most common type of ground heat exchanger is borehole heat exchanger
(BHE) that consists of one or more vertical boreholes. Each borehole contains single or

RI
double U-tube through which a fluid circulates and exchanges heat with the surrounding
ground. Due to high initial cost of GSHPS which is mainly attributed to the borehole

SC
construction, accurate design of BHE is crucial. The main goal of the BHE design is to retain
the borehole leaving fluid temperature within an appropriate range to provide optimal
operating conditions for the heat pump. Therefore, in addition to evaluation of the system

U
performance, the BHE models are used for designing GSHPS.
AN
A considerable number of models ignore thermal capacity of the borehole components versus
the surrounding ground to simplify the modeling. The main assumption of these models is
M

that the heat transfer inside the borehole occurs in a steady state and hence, they are only
suitable for gradually varying heat loads. However, to meet the variable load demand, heat
D

pumps are usually controlled through on-off cycling that may result in rapid variations in the
TE

BHE load. In this situation, transient models should be employed to account for the short-
term response of the BHE. Several studies have been carried out based on the analytical or
numerical approaches. Using a 2D numerical model, Yavuzturk and Spitler [7] extended
EP

Eskilson’s long time-step g-functions [8] to include the short-term response of BHE. This
model was then extended by Fisher et al. [9] for variable heat loads and implemented in
C

EnergyPlus software [10]. Several models have been proposed [11-14] in which the radial
AC

heat transfer in the borehole was analyzed using equivalent diameter approximation. Bauer et
al. [15] and Pasquier and Marcotte [16] suggested thermal resistance-capacity models
(TRCMs) to investigate the short-term behavior of BHE. These models are well adaptable to
various types of BHE e.g. single and double U-tubes. In addition, Minaei and Maerefat [17]
presented an analytical solution based on a TRCM. An analytical composite-medium line-
source model [18] was also suggested that includes the heat capacity of the grouting material.
Kim et al. [19] proposed a hybrid reduced model, combining a modified g-function with a 2D
numerical model to account for the transient heat transfer inside the borehole.

3
 ACCEPTED MANUSCRIPT
The models mentioned earlier are based on 1D or 2D heat transfer analysis in the horizontal
cross-section of the borehole and therefore do not consider the fluid transport along the U-
tube and axial heat conduction in the pipes and grout. However, in dynamic simulation of
GSHPS with on-off cycling, these effects may be important to the system performance. In
addition, these models assume an equally distributed loading versus depth. However, in
actual situation, the borehole loading and flow temperature are not uniform along the

PT
borehole axis especially when the heat pump operates at low capacity. Consequently, these
assumptions may affect the borehole leaving fluid temperature as a main design parameter of

RI
the BHE. In this circumstance, the models which consider axial heat transfer in the BHE, can
provide better estimate. Zarrella et al. [20] suggested a capacity-resistance model for a double

SC
U-tube BHE. Also, the TRCMs were extended to include the fluid advection and axial
conductive effects [21, 22]. Rees [23] developed an extended 2D numerical model to address
the effects of fluid circulation along the U-tube. Maestre et al. [24] proposed a hybrid

U
approach, combining g-functions with a resistance-capacity model. Ma et al. [25] used the
AN
composite-medium line-source model [18] in a quasi-3D approach to account for the fluid
temperature variations along the U-tube.
M

On the other hand, several 3D models have been proposed based on the finite volume
approach. Li and Zheng [26] used an unstructured grid for heat conduction in the grout and
D

ground. Moreover, heat convection in the fluid domain was described by a 1D equation.
TE

Similar numerical model was proposed by Rees and He [27] except that multi-block
structured grid was applied to all domains including the fluid, pipe, grout and ground. Lyu et
al. [28] investigated some parameters influencing the performance of a single U-tube within a
EP

borehole using Fluent software. However, the main drawback of these models is their
computational costs associated with the use of fine grids which might limit their application
C

in the long-term simulation. In addition, several research works have been devoted to the
AC

experimental study of GSHPS. Montagud et al. [29] and Sebarchievici and Sarbu [30]
compared experimental data with results of TRNSYS software. Pärisch et al. [31] compared
three BHE models of TRNSYS with experimental measurements. They also suggested a
modification to improve the short-term response of the steady-state models by an adiabatic
pipe, representing the borehole thermal capacity. Villarino et al. [32] evaluated the
performance of a GSHPS for an office building using experimental data and compared it with
other technologies.

4
 ACCEPTED MANUSCRIPT
Recently, Biglarian et al. [33] developed a numerical model based on the finite volume
approach to predict the thermal response of BHE over both short and long time scales. The
model takes into account the transient 2D heat conduction in the ground and transient 3D heat
transfer inside the borehole using a resistance-capacity (RC) approach. The numerical model
was validated by experimental data of Beier et al. [34] and a fully discretized 3D model with
a multi-block fine grid simulated by ANSYS CFX. The full Navier-Stokes equations along
with the standard  − 6 turbulence model were applied to the fluid domain in the ANSYS

PT
CFX simulation. Even though having appropriate accuracy, the proposed model requires very

RI
low computational cost in comparison with the fully discretized model.

Despite considerable studies have been devoted to development of the transient models as

SC
discussed earlier, less attention has been paid to the impact of the short-term response of
BHEs on the performance of GSHPS. Kummert and Bernier [35] simulated a residential

U
GSHPS with double U-tube BHE using TRNSYS and compared the system performances by
transient and steady-state models. They employed EWS model of Wetter and Huber [36] for
AN
transient simulation of the BHE. This model uses only one lumped capacity in each layer for
the filling material and neglects the pipe thermal capacity. He et al. [37] simulated a water-to-
M

water heat pump coupled to a single BHE using a 3D numerical model. The results were
compared with those of a 2D numerical model and the short time-step g-function model [7]
D

of EnergyPlus [10]. However, the transient behavior of the heat pump under cyclic conditions
TE

was not considered.

In this paper, the transient BHE model [33] is used for dynamic simulation of a GSHPS under
EP

cyclic conditions. This model is coupled to a water-to-air heat pump model and the governing
equations are iteratively solved at each time step. The numerical results are compared with
the results of EnergyPlus over short and long time scales. Furthermore, the transient model is
C

compared with a quasi-steady state model and the effects of borehole thermal capacity on the
AC

system performance and borehole sizing are studied. In addition, the impact of time step size
on the results is investigated.

2. Mathematical description
Fig. 1 depicts the schematic diagram of a unitary ground source water-to-air heat pump
system investigated in this paper. The air loop consists of a mixing box, fan, heating coil,
cooling coil, building and a supplemental heating coil. A blow through configuration is
considered where the fan is located before the heat pump coils. In the mixing box, the outlet

5
 ACCEPTED MANUSCRIPT
air from return duct is mixed with the outdoor fresh air and the leaving air moves through the
fan to the coils. The fan intermittently turns on and off along with the compressor to supply
the building’s heating and cooling loads. The liquid loop, on the other hand, consists of a
circulating pump and a BHE in which a fluid, namely water or aqueous antifreeze solution
circulates through the U-tube to exchange heat with the ground. In dynamic simulation, the
governing equations for all components of the air and liquid loops are simultaneously solved

PT
and the system performance is evaluated accordingly. The mathematical formulations for
each component are presented in the following sections.

RI
2.1 Borehole heat exchanger

SC
As stated earlier, the transient model [33] is employed for thermal analysis of the BHE shown
in Fig. 2. In this model, transient heat conduction in the ground surrounding the borehole is

U
governed by a 2D radial-axial equation [33]:
AN
; ?@=> ; ? ?@=> ? G @=>
= D! E+ (1)
<=> ?A C ?C ?C ?H G

where ! is the radial distance from the borehole axis and 1 is the depth below the ground
M

surface. In addition, transient heat transfer inside the borehole is described by a 3D model
D

which accounts for thermal capacity of the borehole components including the fluid, pipes
and grout. Also, the fluid circulation through the U-tube and heat conduction in the fluid,
TE

pipes, and grout are included. The 3D model was developed based on a detailed thermal RC
circuit with arbitrary number of grout thermal capacities [33]. Fig. 3a illustrates the borehole
cross-section in which the grout area is divided into three sub-sections, namely middle (II ),
EP

upward (IK ) and downward (IL ). Also, Fig. 3b depicts the RC circuit for the case of one
C

capacity for the middle sub-section and three capacities for the upward and downward sub-
sections. This is the optimal discretization scheme in terms of the accuracy and computational
AC

cost as remarked by Biglarian et al. [33]. In this case, heat transfer inside the borehole is
governed by a set of equations as follows [33]:

?@NO ?@NO @QO R@NO ? G @NO

M ± M M =S + M M (2)
?A ?H N TS U
Q ?H G

?@QO @NO R@QO O R@ O

@QW ? G @QO
V =S + + V V
Q
U X (3)
?A N TSQ SQ ?H G

6
 ACCEPTED MANUSCRIPT
O R@ O
@QW O R@ O
@QW O R@
@QW
+ + =0
Q =,Y ==
X U U (4)
SQ S=,Y S==

O
?@=,Y O R@ O
@QW O R@ O
@=,G O
? G @=,Y
I,; = + + I I,;
=,Y =,Y
U X TS U (5)
?A S=,Y S=,Y =,G ?H G

O
?@=,G O R@ O
@=,Y O R@ O
@=,\ O
? G @=,G
I,[ = + + I I,[
=,G =,G
X TSU X TSU (6)
?A S=,Y =,G S=,G =,\ ?H G

O
?@=,\ O R@ O
@=,G O
@X R@=,\ O
? G @=,\
I,] = SX + + I I,]
=,\
(7)

PT
?A U X ?H G
=,G TS=,\ S=,\

?@== > R@
@QW ^ R@
@QW ? G @==
II = + + I II
== ==
U X (8)
?A S== S== ?H G

RI
SC
where superscript 5 could be
or ) referred to the domains associated with the downward-
and upward-flowing fluid, respectively, as shown in Fig. 3b. The sign ± in Eq. (2) is positive
for 5 =
and negative for 5 = ). Equations (1)-(8) are coupled and solved using the finite

U
volume approach. The fully implicit scheme is employed for discretization of the equations in
AN
order to guarantee the stability of the numerical model. The convection term in Eq. (2) is
discretized using the first-order upwind scheme. In the present work, the transient model is
M

also referred to as “RC model”.

As mentioned before, the quasi-steady state models neglect the thermal capacity of the fluid,
D

U-tube and grout. In this situation, more oscillations in the fluid temperature would be
TE

expected under variable heat loads since the thermal capacity has damping effects.
Fortunately, the RC model can readily cover this type of analysis by setting the values of
thermal capacities (M , V , I, and II ) in Eqs. (2)-(8) to zero. This case is referred to as
EP

“resistance model” in the present study. More details about the RC and resistance models can
be found in [33].
C
AC

2.2 Pump
A single speed pump is considered for circulating the heat carrier fluid through the liquid
loop as shown in Fig. 1. The pump turns off when there is no demand by the building. The
hydraulic, shaft, and motor powers are calculated as follows:

_, = 0M ∙ ∆ (9)

7
 ACCEPTED MANUSCRIPT
_,a = _, ⁄7_ (10)

_,c = _,a ⁄7c (11)

where ∆ is the pressure drop caused by the friction losses. The power transferred to the fluid

PT
will eventually be converted to heat due to friction in the liquid loop. Assuming no losses
through the pump casing, the heat gained by the fluid is equal to the shaft power:

RI
SC
'_ = _,a (12)

2.3 Heat Pump

U
The water-to-air heat pump’s capacity and power consumption at full-load conditions can be
AN
evaluated by the equation-fit model of Tang [38]. Using five non-dimensional equations, this
model predicts the steady-state performance of the heat pump in cooling and heating mode
M

[38]:
D

• Cooling mode
TE

deWe,ff @U,iX,jk @ 
mU mN
= ; + [ + ] @N,jk + l m + n m (13)
deWe,ff,ghN @ghN ghN U,ghN N,ghN
C EP

dohk,ff @U,>X,jk @U,iX,jk @ 

mU mN
dohk,ff,ghN
= ; + [ @ghN
+ ] @ghN
+ l @N,jk + n m + p m (14)
ghN U,ghN N,ghN
AC

_ff @U,iX,jk @ 
mU mN
_ff,ghN
= ; + [ @ghN
+ ] @N,jk + l m + n m (15)
ghN U,ghN N,ghN

• Heating mode

8
 ACCEPTED MANUSCRIPT
dqf @U,>X,jk @N,jk 
mU mN
= ; + [ + ] + l + n (16)
dqf,ghN @ghN @ghN mU,ghN mN,ghN

_qf @U,>X,jk @N,jk 

mU mN
= ; + [ + ] + l 
+ n 
(17)
_qf,ghN @ghN @ghN mU,ghN mN,ghN

PT
where the reference temperature, -CrM , is 283 K [38]. Also, other reference parameters are the
conditions at which the heat pump delivers the maximum capacity [38]. However, it should

RI
be noted that the reference conditions for the heating mode may be different from the
indicated conditions for the cooling mode [38]. Equations (13)-(17) imply that the heat pump

SC
performance is dependent upon the temperature and volumetric flow rate of the liquid and air
entering to the heat pump coils. The coefficients  ,  ,  ,  and  in Eqs. (13)-(17) are

U
obtained from the manufacturer’s performance data by applying the generalized least square
AN
method [38]. The data should include the capacity and power consumption of the heat pump
versus a range of temperatures and volumetric flow rates of the air and liquid.
M

Assuming that heat loss through the compressor casing is negligible, the source side heat
transfer rate for the cooling and heating mode is given by [38]:
D
TE

'stKCur,u = 'AtA,uu + uu (18)

EP

'stKCur,v = 'vu − vu (19)

C
AC

Equations (13)-(19) represent the heat pump performance at steady-state conditions.

However, heat pumps are cycling on and off to meet the variable load demand. This transient
behavior has negative effect on the heat pump performance. In the present study, part-load
fraction (PLF) model [39] is employed to account for the performance deterioration under
cyclic conditions. The PLF is defined as the energy efficiency ratio (EER) at part load
divided by the full-load EER [38] that can be calculated using the following iterative formula
[39]:

9
 ACCEPTED MANUSCRIPT

_S R;
" wT; = 1 − 48 D1 − E |1 − &} ~ …† (20)
_z{ l€U‚ ƒ;R_S⁄_z{ „

where PLR is the part-load ratio defined as the average capacity of the heat pump at part load

PT
divided by the capacity at full-load conditions [38]. Fig. 4 shows the variations of PLF versus
PLR resulted from Eq. (20) for a time constant (8) of 60 s and a maximum cycling rate
( ) of 2.5 cycles per hour. These values were recommended by Henderson et al. [40] for

RI
a “typical” condition. As shown in Fig. 4, a polynomial trend-line of order 3 perfectly
matches the results obtained from Eq. (20). Therefore, instead of using Eq. (20) which

SC
requires time-consuming iterations, a regression equation can be employed:

" = ‡ + ; "% + [ "% [ + ] "% ]

U
AN
(21)
M

Bonne et al. [41] showed that considering the off-cycle power significantly affects the PLF at
low PLRs. They suggested a correction to adjust the PLF as follows [41]:
D
TE

_S
" $ = ˆ‰Š ˆ‰Š (22)
TD;R EVC
ˆ‰‹ ˆ‰‹
EP

2.4 Liquid loop

C

The heat balance for the liquid loop should be preserved at each time step. In heating mode,
AC

heat extracted from the ground together with heat gain from the pump is transferred to the
source side coil and eventually absorbed by the refrigerant circulating through the coil.
Conversely, in cooling mode, heat rejected to the fluid in the source side coil and heat
resulting from the pump are transferred to the BHE and dissipated in the ground. Assuming
no losses, the liquid loop heat balance for the cooling and heating mode is given by:

'Œ,u = 'stKCur,u + '_ (23)

10
 ACCEPTED MANUSCRIPT
'Œ,v = 'stKCur,v − '_ (24)

2.5 Air loop

The basis for the air loop shown in Fig. 1 is to solve the mass and energy balance equations
for each component. The building may be considered as a black box with specified heating

PT
and cooling loads which can be calculated using building energy simulation programs, e.g.
EnergyPlus [10]. The building loads, indoor and outdoor conditions are the boundary

RI
conditions needed for the solution of the mass and energy balance equations at each time
step. As mentioned before, the fan operates through on-off cycling along with the

SC
compressor. Hence, the air flow rate, fan power and heat pump capacity are adjusted by the
PLR in the course of simulation [38]. The PLR can be calculated by the following closed-
form expressions derived from the air loop energy balance for the cooling and heating mode:

U
AN
d,eWe,f
"%u = d (25)
ohk,ff R_‹ RU,W ƒvU,W RvU,j „
M

d,q
"%v = d
D

(26)
qf T_‹ TU,W ƒvU,W RvU,j „
TE

where ℎ is the specific enthalpy of air evaluated by the psychrometric correlations and ,t
EP

is the full-load fresh air mass flow rate entering the mixing box. Also, 'srŽ,uu and 'vu are the
full-load sensible cooling and heating capacities at operating conditions calculated by the
C

equation-fit model [38].

AC

The compressor power is adjusted by the run-time fraction defined as [38]:

_S
9 = _z (27)

The heat pump coefficient of performance (COP) is also defined for the cooling and heating
mode as follows:

11
 ACCEPTED MANUSCRIPT

deWe,ff
 _,u = _ff
(28)

dqf
 _,v = (29)
_qf

PT
where the heat pump’s capacity and power consumption should be adjusted for part-load

RI
conditions. Also, the whole system performance can be determined by the net COP,
accounting for the pump and fan powers:

U SC
deWe,ff R_‹
 ss,u = (30)
_ff T_‹ T_ˆ
AN
dqf T_‹
 ss,v = _ (31)
M

qf T_‹ T_ˆ

The COP measures the instantaneous performance of the heat pump or system that does not
D

reflect the performance over a time interval. Therefore, a seasonal performance factor (SPF)
TE

may be defined to evaluate the performance under varying operational conditions over a
certain time span. The general expression for the SPF is given by
EP

e
‘e G dLA
Y
+ = e (32)
‘e G _LA
Y
C

which can be defined for the heat pump and system in cooling or heating mode.
AC

Fig. 5 displays the logic and sequences of all calculation steps given in the section 2 in the
form of a flowchart. All governing equations are solved using an in-house computer code
developed in C# programming environment and the results are printed in output files. The
input data include the building loads, indoor and outdoor conditions i.e. dry bulb temperature,
relative humidity and barometric pressure for each time step. As shown in Fig. 5, the
calculations of both the liquid and air loops at each time step are performed by an iterative
procedure to maintain the accuracy of the heat pump entering liquid temperature about
6 = 0.01 ºC.

12
 ACCEPTED MANUSCRIPT

3. System specifications
In order to study the performance of the proposed model, a residential GSHPS is considered.
The system specifications are presented in the following sections.

3.1 Building characteristics

PT
A typical single-family house with a 130 m2 of floor area is considered that is assumed to be
located in the city of Tehran, Iran. The building has been modeled as a single thermal zone

RI
with internal loads for a typical family of four. Table 1 details the specifications of the
building. The building load calculation has been carried out by EnergyPlus 8.5 for a 1-year

SC
period based on a 10-min time step. In this regard, the weather data for the city of Tehran
have been taken from Iranian Typical Meteorological Year (ITMY) data set [42]. Fig. 6

U
shows the building load profile over a 1-year period resulted from the simulation. Note that
the negative values represent the total cooling load, while the positive values imply the
AN
heating load. As it is shown in Fig. 6, the peak heating and cooling loads are 5.2 kW and 3.0
kW, respectively. Also, the annual heating and total cooling energy demands are 3242 kWh
M

and 1812 kWh, respectively, calculated by integrating the building load profile numerically
over the 1-year period.
D
TE

3.2 Heat pump specifications

A 1.5-ton water-to-air heat pump has been selected based on the peak heating load of the
EP

building. The rated heating and cooling capacities of the heat pump are 5.3 kW and 5.4 kW,
respectively. Also, the rated and reference volumetric flow rates of the working liquid are
C

0.265 dm3 s-1 and 0.315 dm3 s-1, respectively. Moreover, the rated and reference volumetric
flow rates of the air are 269 dm3 s-1 and 323 dm3 s-1, respectively. As noted in section 2.3, the
AC

heat pump steady-state performance is calculated by Eqs. (13)-(17). Table 2 summarizes the
coefficients of these equations for the specified heat pump. Furthermore, the parameters of
the PLF model including the time constant, maximum cycling rate and fraction of on-cycle
power use ( !) are respectively taken as 60 s, 2.5 cycles per hour and 0.01 corresponding to a
typical condition [40]. Furthermore, the volumetric flow rate of the fresh air entering the
mixing box is considered as 9 dm3 s-1 per person at full-load conditions. The power of the
supplemental heater is also taken zero.

13
 ACCEPTED MANUSCRIPT

3.3 BHE specifications

A single borehole with diameter of 12.5 cm and length of 100 m is considered. The
undisturbed ground temperature is assumed to be 17 °C which is approximately equal to the
mean annual air temperature of the city of Tehran. The volumetric flow rate of the working
liquid (water) is considered as 0.265 dm3 s-1 which is equal to that of the heat pump at rated

PT
conditions. The motor’s power and efficiency and also the circulating pump efficiency are
taken as 50 W, 0.88 and 0.6, respectively. Other parameters of the BHE are given in Table 3.

RI
SC
4. Results and discussion
The dynamic simulation of the GSHPS described in the section 3 has been performed using

U
the proposed model. Also, the integrated system has been separately simulated by EnergyPlus
AN
8.5 and the obtained results are compared with those from the proposed model. Moreover, the
effects of borehole thermal capacity and time step size on the system performance are
investigated. In addition, design of the BHE is carried out using the specified BHE models.
M

4.1 Comparison of the proposed model with EnergyPlus

D

The GSHPS simulation has been performed by EnergyPlus and the proposed model for a 1-
TE

year period. The RC model is used for the BHE simulation in the proposed approach. On the
other hand, the short time-step g-function data for the single borehole used in EnergyPlus
model have been taken from [7]. The time step size is taken as 10 minutes for both two
EP

models. A constant on-demand flow control is considered for the circulating pump in both
two models through which full flow is provided whenever there is a load. Fig. 7 shows the
C

heat pump’s load and source side heat transfer rates ('vu and 'stKCur,v ) and the compressor’s
AC

power consumptions for a typical heating day (January 17) by the proposed model and
EnergyPlus. It is observed that there is a good agreement between the two models. Also, Fig.
7 shows that the heat pump turns off when there is no heating demand (22:00-22:40). Fig. 8
compares the borehole heat transfer rates by the proposed model and EnergyPlus over the
same period. It can be seen that there are notable discrepancies between the two models.
Since the heat gain resulting from the circulating pump is constant and small (44 W) during
its operation, on the basis of Eq. (24), it is expected that the borehole heat transfer rates
should follow the same trend as the source side heat transfer rates shown in Fig. 7. Although

14
 ACCEPTED MANUSCRIPT
the proposed model shows such trend, EnergyPlus predicts a different pattern especially
when there is a sudden rise in the heating demand (6:00-6:40 and 22:40-23:10). In this
regard, the liquid loop heat balance error for the two models can be calculated based on Eq.
(24) as:

 = “'Œ,v + '_ − 'stKCur,v | (33)

PT
As shown in Fig. 8, the heat balance error for the proposed model is negligible over the

RI
whole period. Indeed, the iterative procedure used in the RC model maintains the error below
0.1 W at each time step throughout the simulation. On the other hand, the maximum error for

SC
EnergyPlus is 1675 W at 6:05 a.m. This error could influence the borehole entering and
leaving fluid temperatures and borehole wall temperature as shown in Fig. 9. In addition, the

U
short-time step g-function model does not account for the fluid transport and fluid thermal
capacity [37]. These factors may also cause the discrepancy observed in Fig. 9. It should be
AN
noted that the minimum borehole leaving fluid temperature over the year occurs on January
17 at 6:45 a.m. (Fig. 9). The heat extracted from the ground over this day by the proposed
M

model and EnergyPlus is 38.2 kWh and 39.7 kWh, respectively. Also, Fig. 9 indicates that
when the circulating pump turns off (22:00-22:40), the RC model accounts for the thermal
D

recovery process of the BHE. Hence, the fluid and borehole experience a temperature rise
during this period. Fig. 10 shows COPs and PLRs of the heat pump by the proposed model
TE

and EnergyPlus over the same period. The PLRs calculated by the two models, perfectly
match each other. This implies that Eq. (26) yields the same results as the iterative procedure
EP

of EnergyPlus. The proposed model, on the other hand, predicts slightly lower COP than
EnergyPlus. It is caused by the lower borehole leaving fluid temperature by the proposed
C

model as shown in Fig. 9. The maximum COP over this period by the proposed model and
AC

EnergyPlus is 4.27 and 4.39, respectively, occurring at 6:05 a.m. In addition, Fig. 10 reveals
that the COP improves with increase in the PLR and vice versa. Despite the borehole leaving
fluid temperature rises during 6:45 to 22:00 (Fig. 9), the heat pump COP declines over the
same period due to reduction in the PLR as shown in Fig. 10. Therefore, it can be concluded
that the cycling losses significantly affect the heat pump performance and hence, should be
taken into consideration in the GSHPS simulation.

15
 ACCEPTED MANUSCRIPT
Fig. 11 shows the monthly average fluid and borehole wall temperatures by the proposed
model and EnergyPlus. Herein, the average fluid temperature represents the mean
temperature of the fluid entering and leaving the BHE. It is apparent that there is a close
agreement between the two models for the heating months. However, for the cooling months,
the proposed model predicts slightly higher fluid and borehole wall temperatures than
EnergyPlus. The annual heat rejected to the ground in the cooling mode is 2850 kWh and

PT
2664 kWh for the proposed model and EnergyPlus, respectively. The more heat rejection
results in the higher fluid and borehole wall temperatures which lead to the reduction in the

RI
heat pump cooling performance. Fig. 12 compares the monthly cooling SPFs of the heat
pump by the proposed model and EnergyPlus. As mentioned, due to the higher fluid

SC
temperature over cooling months (Fig. 11), the proposed model predicts slightly lower SPF
than EnergyPlus. However, the discrepancy between the two models is negligible as shown in
Fig. 12.

U
Table 4 summarizes the annual energy performances for the cooling and heating mode
AN
resulted from the annual simulations with the proposed model and EnergyPlus. The
maximum difference between the two models is 6.5%. In addition, Table 4 shows the
M

cumulative loop heat balance errors over the cooling, heating, and whole periods. It is
observed that the proposed model appropriately preserves the liquid loop heat balance over
D

the long-term periods.

TE

4.2 Effects of time step and borehole thermal capacity

EP

The time step is one of the main parameters in the numerical simulation of GSHPS. The
larger time steps result in the larger truncation errors which may affect the model accuracy.
C

The small time steps, on the other hand, lead to the vast computational cost which is not
AC

suitable for the long-term simulation. Therefore, a trade-off between the accuracy and
computational cost should be considered. For the system under study, the time step may
directly influence the building loads and the BHE thermal response due to the transient
numerical approaches used for their calculations. Due to slow variations of building’s indoor
and outdoor conditions, the hourly-analysis is common practice. Therefore, the building loads
calculated by EnergyPlus based on a 10-min time step (Fig. 6) are fairly accurate. Since the
objective is to assess the impact of time step on the BHE model and the consequence effects
on the system performance, the building load profile is assumed to be unchanged under

16
 ACCEPTED MANUSCRIPT
different time steps. Nevertheless, to account for smaller time steps, the building loads are
estimated by the interpolation method from the specified load profile. Also, for larger time
steps, the loads are determined by averaging the load profile over each interval. In this
respect, several time step values are considered. Fig. 13 shows the effect of time step on the
borehole leaving fluid temperature over the same day as the previous section (January 17). It
is observed that increasing the time step from 1 to 30 minutes, does not lead to significant

PT
difference in the results of the RC model. The discrepancies are less than 0.3 °C between the
corresponding points. This is an outstanding feature of the RC model which makes it suitable

RI
for the long-term simulation of GSHPS. Fig.13 also reveals that with a 60-min time step, the
discrepancy increases when there is a sudden temperature variation. However, the deviations

SC
do not exceed 0.5 °C compared with the corresponding values associated with the 1-min time
step. In addition, Fig. 13 shows the borehole leaving fluid temperature by the resistance
model. As noted before, this model ignores the thermal capacity of the fluid, pipes and grout.

U
In view of this assumption, the fluid temperature shows more fluctuations. The discrepancy
AN
between the RC and resistance models is more evident when the temperature varies sharply
as shown in Fig. 13. The maximum temperature difference between the two models with the
M

same time step is 2.44 °C, occurring at 6:05 a.m.

D

Table 5 compares the annual performances by the resistance and RC models with several
TE

time steps. It can be found that the discrepancies in the annual energy performances by the
two models regarding different time steps are less than 1% which are negligible. It was
EP

expected that the annual performance predicted based on the RC model would not differ
significantly versus time step since the borehole leaving fluid temperature is only slightly
influenced by the time step as discussed earlier. Surprisingly, the resistance model predicts
C

approximately the same results except for the maximum and minimum borehole leaving fluid
AC

temperatures. Despite the discrepancy in the borehole leaving fluid temperatures (Fig. 13)
could influence the system performance, its effect is only notable at short timescales. The
oscillation behavior of the fluid temperature may cause these errors to cancel over the longer
periods. Therefore, the use of the resistance model may not remarkably affect the annual
ground heat exchange, electricity consumption, and SPF when compared with the RC model.
However, as shown in Table 5, there is a significant difference in the minimum (and
maximum) borehole leaving fluid temperatures by the two models which should be carefully
considered in the design of the BHE. It can be found that the resistance model tends to

17
 ACCEPTED MANUSCRIPT
underestimate the minimum borehole leaving fluid temperature by 1.52 °C. Table 5 also
presents the computation time and speedup for all simulations which have been carried out on
a PC with quad-core 2.40 GHz CPU and 8 GB of memory. Herein, the speedup is the ratio
between the program run time using a 1-min time step based on the RC model and the
computation time required for each case. It is observed that the use of a 1-min time step leads
to a large computation time. However, the obtained results do not considerably differ from

PT
those predicted by the RC model with larger time steps. Also, it is worth noting that the
resistance model requires more computation time than the RC model with the same time step.

RI
Indeed, the resistance model needs more iterations at each time step to converge.

SC
The analysis presented in this section indicates that the RC model can appropriately predict
the system performance over both short and long time scales especially when the time step is

U
taken below 30 minutes. However, it should be noted that the results were obtained under the
specified building loads with a 10-min time step. If the effect of time step on the building
AN
load being considered, the situation might be different. Therefore, since the RC model is less
sensitive to the time step, it is suitable to select the time step for the whole system simulation
M

equal to that required for the load calculations. If the building conditions vary suddenly,
smaller time steps are needed and for slow variations, larger time steps can be used. A time
D

step between 10 to 30 minutes may be reasonable in terms of accuracy and computational

TE

cost.
EP

4.3 Design of the BHE

The maximum and minimum entering liquid temperatures to the heat pump are the main
C

parameters for the design of BHEs. In the previous sections, a borehole with length of 100 m
AC

was chosen. It was shown that there are some discrepancies in the heat pump entering liquid
temperatures predicted by the RC, resistance and the short-time step g-function models. Now,
the objective is to determine the appropriate borehole length for the specified GSHPS using
these models. Since the building under study is heating dominated, therefore the design is
carried out based on the heating mode. According to Kavanaugh and Rafferty [43], for an
optimal design in heating mode, the entering liquid temperature to the heat pump should be 5
to 8 °C lower than the undisturbed ground temperature. Therefore, herein a value of 10 °C is
chosen as the minimum heat pump entering liquid temperature which is 7 °C less than the

18
 ACCEPTED MANUSCRIPT
undisturbed ground temperature (17 °C). Table 6 compares the design borehole lengths by
the RC, resistance and g-function models obtained from trial and error processes. All
simulations have been carried out for a 10-year design life using a 10-min time step. It can be
found that the RC model shows a 16% reduction in the required borehole length in
comparison with the resistance model. Therefore, neglecting the borehole thermal capacity
leads to overestimation of the borehole length. Table 6 also shows that the short time-step g-

PT
function model of EnergyPlus predicts a smaller borehole than the RC model. The difference
could be attributed to the liquid loop heat balance issue as discussed earlier.

RI
SC
5. Conclusions
In this study, a numerical approach is proposed for dynamic simulation of a GSHPS. The RC
model is used to predict the transient response of the BHE. This model accounts for the

U
borehole thermal capacity, fluid circulation through the U-tube and heat conduction in the
AN
fluid, pipes, grout and ground. In addition, the equation-fit model along with the PLF model
is used to calculate the transient performance of the heat pump. The proposed numerical
M

model is employed to simulate a residential ground source water-to-air heat pump system in
the climatic conditions of Tehran, Iran. The predicted results are generally in good agreement
D

with those from EnergyPlus software. However, there are some discrepancies in the borehole
heat transfer rates and fluid temperatures that could be attributed to the liquid loop heat
TE

balance of the two models. Also, the results obtained based on the RC model are only slightly
affected by increasing the time step size from 1 to 30 minutes. Hence, larger time steps can
EP

be used to reduce the computational cost. In addition, the results show that ignoring the
borehole thermal capacity has a remarkable impact on the borehole leaving fluid temperature
C

at short timescales, while its effect on the annual energy performance of the GSHPS is less
impressive. The minimum borehole leaving fluid temperature by the resistance model is 1.52
AC

°C less than that predicted by the RC model. This temperature disparity would result in
different design values for the borehole length. It is found that including the borehole thermal
capacity through the RC model leads to a 16% reduction in the design borehole length.
Consequently, much lower capital costs would be expected when transient models instead of
quasi-steady state models are used for the design of BHEs.

Acknowledgments

19
 ACCEPTED MANUSCRIPT
Financial support from INSF, through project No. 940017 is greatly acknowledged.

References
[1] Li DHW, Yang L, Lam JC. Zero energy buildings and sustainable development implications – A
review. Energy. 2013;54:1-10.
[2] Li G. Energy and exergy performance assessments for latent heat thermal energy storage
systems. Renewable and Sustainable Energy Reviews. 2015;51:926-54.
[3] Li G. Sensible heat thermal storage energy and exergy performance evaluations. Renewable and

PT
Sustainable Energy Reviews. 2016;53:897-923.
[4] Li G, Zheng X. Thermal energy storage system integration forms for a sustainable future.
Renewable and Sustainable Energy Reviews. 2016;62:736-57.

RI
[5] Li G, Hwang Y, Radermacher R. Cold thermal energy storage materials and applications toward
sustainability. In: Zhang X, Dincer I, editors. Energy Solutions to Combat Global Warming: Springer
International Publishing; 2017, p. 67-117.

SC
[6] Li G, Hwang Y. Energy Storage Systems for Buildings. In: Gonzalez JE, Krarti M, editors. Handbook
of Integrated and Sustainable Buildings Equipment and Systems, Volume I: Energy Systems. New
York: ASME; 2017, 74 p.
[7] Yavuzturk C, Spitler JD. A short time step response factor model for vertical ground loop heat

U
exchangers. ASHRAE Transactions. 1999;105(2):475-85.
[8] Eskilson P. Thermal Analysis of Heat Extraction Boreholes. Doctoral Thesis. Sweden: University of
AN
Lund; 1987.
[9] Fisher DE, Rees SJ, Padhmanabhan SK, Murugappan A. Implementation and validation of ground-
source heat pump system models in an integrated building and system simulation environment.
M

HVAC&R Research. 2006;12(3a):693-710.

[10] Crawley DB, Lawrie LK, Winkelmann FC, Buhl WF, Huang YJ, Pedersen CO, et al. EnergyPlus:
creating a new-generation building energy simulation program. Energy and Buildings.
D

2001;33(4):319-31.
[11] Lamarche L, Beauchamp B. New solutions for the short-time analysis of geothermal vertical
boreholes. International Journal of Heat and Mass Transfer. 2007;50:1408-19.
TE

[12] Javed S, Claesson J. New analytical and numerical solutions for the short-term analysis of
vertical ground heat exchangers. ASHRAE Transactions. 2011;117(1):3-12.
[13] Shirazi AS, Bernier M. Thermal capacity effects in borehole ground heat exchangers. Energy and
EP

Buildings. 2013;67:352-64.
[14] Gordon D, Bolisetti T, Ting DSK, Reitsma S. Short-term fluid temperature variations in either a
coaxial or U-tube borehole heat exchanger. Geothermics. 2017;67:29-39.
[15] Bauer D, Heidemann W, [Link]-Steinhagen, Diersch H-JG. Thermal resistance and capacity
C

models for borehole heat exchangers. International Journal of Energy Research. 2011;35(4):312-20.
[16] Pasquier P, Marcotte D. Short-term simulation of ground heat exchanger with an improved
AC

TRCM. Renewable Energy. 2012;46:92-9.

[17] Minaei A, Maerefat M. A new analytical model for short-term borehole heat exchanger based
on thermal resistance capacity model. Energy and Buildings. 2017;146:233-42.
[18] Li M, Lai ACK. Analytical model for short-time responses of ground heat exchangers with U-
shaped tubes: Model development and validation. Applied Energy. 2013;104:510-6.
[19] Kim E-J, Bernier M, Cauret O, Roux J-J. A hybrid reduced model for borehole heat exchangers
over different time-scales and regions. Energy. 2014;77:318-26.
[20] Zarrella A, Scarpa M, Carli MD. Short time step analysis of vertical ground-coupled heat
exchangers: The approach of CaRM. Renewable Energy. 2011;36:2357-67.
[21] Pasquier P, Marcotte D. Joint use of quasi-3D response model and spectral method to simulate
borehole heat exchanger. Geothermics. 2014;51:281-99.

20
 ACCEPTED MANUSCRIPT
[22] Bauer D, Heidemann W, Diersch H-JG. Transient 3D analysis of borehole heat exchanger
modeling. Geothermics. 2011;40:250-60.
[23] Rees SJ. An extended two-dimensional borehole heat exchanger model for simulation of short
and medium timescale thermal response. Renewable Energy. 2015;83:518-26.
[24] Maestre IR, Gallero FJG, Gomez PA, Perez-Lombard L. A new RC and g-function hybrid model to
simulate vertical ground heat exchangers. Renewable Energy. 2015;78:631-42.
[25] Ma W, Li M, Li P, Lai ACK. New quasi-3D model for heat transfer in U-shaped GHEs (ground heat
exchangers): Effective overall thermal resistance. Energy. 2015;90, Part 1:578-87.
[26] Li Z, Zheng M. Development of a numerical model for the simulation of vertical U-tube ground

PT
heat exchangers. Applied Thermal Engineering. 2009;29(5–6):920-4.
[27] Rees SJ, He M. A three-dimensional numerical model of borehole heat exchanger heat transfer
and fluid flow. Geothermics. 2013;46:1-13.
[28] Lyu Z, Song X, Li G, Hu X, Shi Y, Xu Z. Numerical analysis of characteristics of a single U-tube

RI
downhole heat exchanger in the borehole for geothermal wells. Energy. 2017;125:186-96.
[29] Montagud C, Corberán JM, Ruiz-Calvo F. Experimental and modeling analysis of a ground source
heat pump system. Applied Energy. 2013;109:328-36.

SC
[30] Sebarchievici C, Sarbu I. Performance of an experimental ground-coupled heat pump system for
heating, cooling and domestic hot-water operation. Renewable Energy. 2015;76:148-59.
[31] Pärisch P, Mercker O, Oberdorfer P, Bertram E, Tepe R, Rockendorf G. Short-term experiments

U
with borehole heat exchangers and model validation in TRNSYS. Renewable Energy. 2015;74:471-7.
[32] Villarino JI, Villarino A, Fernández FÁ. Experimental and modelling analysis of an office building
AN
HVAC system based in a ground-coupled heat pump and radiant floor. Applied Energy.
2017;190:1020-8.
[33] Biglarian H, Abbaspour M, Saidi MH. A numerical model for transient simulation of borehole
heat exchangers. Renewable Energy. 2017;104:224-37.
M

[34] Beier RA, Smith MD, Spitler JD. Reference data sets for vertical borehole ground heat exchanger
models and thermal response test analysis. Geothermics. 2011;40(1):79-85.
[35] Kummert M, Bernier M. Sub-hourly simulation of residential ground coupled heat pump
D

systems. Building Service Engineering Research and Technology. 2008;29(1):27-44.

[36] Wetter M, Huber A. TRNSYS Type 451: Vertical Borehole Heat Exchanger EWS Model, Version
TE

3.1-Model description and implementing into TRNSYS. 1997.

[37] He M, Rees S, Shao L. Simulation of a domestic ground source heat pump system using a three-
dimensional numerical borehole heat exchanger model. Journal of Building Performance Simulation.
2011;4(2):141-55.
EP

[38] TANG CC. Modeling Packaged Heat Pumps in a Quasi-steady State Energy Simulation Program.
Master's Thesis. Stillwater, OK: Oklahoma State University; 2005.
[39] Henderson HI, Rengarajan K. Model to predict the latent capacity of air conditioners and heat
C

pumps at part-load conditions with constant fan operation. ASHRAE Transactions. 1996;102(1):266-
74.
AC

[40] Henderson H, Huang YJ, Parker D. Residential Equipment Part Load Curves for Use in DOE-2.
Berkeley, CA: Lawrence Berkeley National Laboratory; 1999.
[41] Bonne U, Patani A, Jacobson RD, Mueller DA. Electric-driven heat pump systems: Simulations
and controls II. ASHRAE Transactions. 1980;86(1):687-705.
[42] Iranian Typical Meteorological Year (ITMY) data set, [Link]
location/asia_wmo_region_2/IRN//IRN_Tehran-Mehrabad.407540_ITMY; 2017 [accessed March
2017].
[43] Kavanaugh S, Rafferty K. Geothermal Heating and Cooling: Design of Ground-Source Heat Pump
Systems. Atlanta, GA: ASHRAE; 2014.

21
 ACCEPTED MANUSCRIPT
Table captions:

Table 1 Summarized parameters of the building.

Table 2 Coefficients of the equation-fit model for the selected heat pump.

Table 3 Design parameters of the BHE.

Table 4 Comparison of the annual performances by the proposed model and EnergyPlus.

PT
Table 5 Comparison of the annual performances by the resistance and RC models with several time
steps.

RI
Table 6 Design of the BHE based on different models.

U SC
AN
M
D
TE
C EP
AC

22
 ACCEPTED MANUSCRIPT
Table 1 Summarized parameters of the building.

House characteristics
Floor area (m2) 130
Conditioned volume (m3) 448.5
Window-wall ratio (%) 12.93
Thermal characteristics
External walls U-value (W m-2 K-1) 0.510
External roof U-value (W m-2 K-1) 0.223

PT
Floor U-value (W m-2 K-1) 2.945
Window U-value (W m-2 K-1) 1.984
Set points

RI
Heating (6:00-22:00) (°C) 22
Heating (22:00-6:00) (°C) 18
Cooling (7:00-22:00) (°C) 23
Cooling (22:00-7:00) (°C) 26

U SC
AN
M
D
TE
EP
C
AC

23
 ACCEPTED MANUSCRIPT
Table 2 Coefficients of the equation-fit model for the selected heat pump.
; [ ] l n -
-1.354 3.340 -1.240 0.144 0.013 -
; [ ] l n p
0.380 17.688 -17.564 -0.561 0.472 0.005
; [ ] l n -
-5.561 -2.371 9.243 -0.011 -0.203 -
; [ ] l n -
-3.948 -0.789 5.549 0.049 0.098 -

PT
; [ ] l n -
-5.248 5.067 1.332 -0.369 0.019 -

RI
U SC
AN
M
D
TE
C EP
AC

24
 ACCEPTED MANUSCRIPT
Table 3 Design parameters of the BHE.

Parameter Value
Ground thermal conductivity (W m-1 K-1) 2.1
Ground volumetric heat capacity (J m-3 K-1) 1.9E+6
Grout thermal conductivity (W m-1 K-1) 1.0
Grout volumetric heat capacity (J m-3 K-1) 3.3E+6
Pipe thermal conductivity (W m-1 K-1) 0.4
Pipe volumetric heat capacity (J m-3 K-1) 2.28E+6

PT
Pipe outer diameter (m) 0.032
Pipe inner diameter (m) 0.026
Half shank spacing (m) 0.03
Buried depth (m) 5.0

RI
U SC
AN
M
D
TE
EP
C
AC

25
 ACCEPTED MANUSCRIPT
Table 4 Comparison of the annual performances by the proposed model and EnergyPlus.

Proposed model EnergyPlus Difference (%)

Extracted heat from ground (kWh) 2494 2600 4.3
Rejected heat to ground (kWh) 2850 2664 6.5
Pump electricity consumption in heating mode (kWh) 154.9 154.2 0.5
Pump electricity consumption in cooling mode (kWh) 106.2 105.9 0.3
Fan electricity consumption in heating mode (kWh) 68.6 67.9 1.1
Fan electricity consumption in cooling mode (kWh) 47.5 46.3 2.5

PT
Heating coil electricity consumption (kWh) 976.4 961.5 1.5
Heating coil heating energy (kWh) 3604 3579 0.7
Cooling coil electricity consumption (kWh) 638.2 610.8 4.3
Cooling coil total cooling energy (kWh) 2116 2067 2.3

RI
Cooling coil sensible cooling energy (kWh) 1992 1944 2.4
Heat pump heating SPF 3.691 3.722 0.8
System heating SPF 3.061 3.081 0.7

SC
Heat pump cooling SPF 3.316 3.384 2.1
System cooling SPF 2.612 2.648 1.4
Liquid loop heat balance error in heating mode (kWh) 2.8 -15.1 -
Liquid loop heat balance error in cooling mode (kWh) -2.8 16.9 -

U
Annual liquid loop heat balance error (kWh) 0.0 1.8 -
AN
M
D
TE
C EP
AC

26
 ACCEPTED MANUSCRIPT
Table 5 Comparison of the annual performances by the resistance and RC models with several time
steps.

RC model Resistance model

Time step (min) 1 10 30 60 10
Extracted heat from ground (kWh) 2492 2494 2492 2488 2489
Rejected heat to ground (kWh) 2855 2850 2855 2859 2855
Total electricity consumption in heating mode (kWh) 1203 1200 1204 1211 1203
Total electricity consumption in cooling mode (kWh) 794.8 791.9 794.7 796.7 798.6

PT
Heat pump heating SPF 3.687 3.691 3.686 3.677 3.678
System heating SPF 3.051 3.061 3.050 3.033 3.051
Heat pump cooling SPF 3.312 3.316 3.311 3.309 3.279
System cooling SPF 2.602 2.612 2.603 2.597 2.588

RI
Maximum borehole leaving fluid temperature (°C) 25.49 25.47 25.43 25.43 26.06
Minimum borehole leaving fluid temperature (°C) 9.12 9.23 9.32 9.42 7.71
Computation time (s) 2547 299 117 72 370

SC
Speedup 1 8.5 21.8 35.4 6.9

U
AN
M
D
TE
C EP
AC

27
 ACCEPTED MANUSCRIPT
Table 6 Design of the BHE based on different models.

RC model Resistance model g-function model

Minimum heat pump entering liquid temperature (°C) 9.998 10.005 10.007
Required borehole length (m) 115 137 102

PT
RI
U SC
AN
M
D
TE
C EP
AC

28
 ACCEPTED MANUSCRIPT
Figure captions:

Fig. 1 Schematic diagram of a ground source water-to-air heat pump system with a blow through configuration.

Fig. 2 Schematic of a borehole heat exchanger.

Fig. 3 (a) Borehole cross-section (b) corresponding thermal RC circuit.

Fig. 4 Part-load fraction versus part-load ratio.

PT
Fig. 5 Flowchart of calculations.

Fig. 6 Load profile of the building.

RI
Fig. 7 Heat transfer rates and power consumption of the heat pump by the proposed model and EnergyPlus for
the specified heating day (Jan. 17).

SC
Fig. 8 Comparison of borehole heat transfer rates by the proposed model and EnergyPlus for the specified
heating day (Jan. 17).

U
Fig. 9 Fluid and borehole wall temperatures by the proposed model and EnergyPlus for the specified heating
AN
day (Jan. 17).

Fig. 10 COPs and PLRs by the proposed model and EnergyPlus for the specified heating day (Jan. 17).
M

Fig. 11 Monthly average fluid and borehole wall temperatures by the proposed model and EnergyPlus.

Fig. 12 Comparison of monthly heat pump cooling SPFs by the proposed model and EnergyPlus.
D

Fig. 13 Effects of time step and borehole thermal capacity on the borehole leaving fluid temperature for the
specified heating day (Jan. 17).
TE
C EP
AC

29
 ACCEPTED MANUSCRIPT

Fan
Source side Source side
cooling coil heating coil

PT
RI
SC
Supplemental

U
AN
Cooling coil Heating coil heating coil

M
D
TE
EP
Mixing box Building

C
AC
BHE Pump

Outside air Relief air Air loop (Load side) Liquid loop (Source side)
 ACCEPTED MANUSCRIPT
4.5 1.0
Proposed model, Heat pump COP
EnergyPlus, Heat pump COP
4.3 0.8

PT
Proposed model, Heat pump PLR

RI
EnergyPlus, Heat pump PLR

SC
4.1 0.6

PLR
COP

AN
M
3.9 0.4

D
TE
EP
3.7 0.2

C
AC
3.5 0.0
00:00 04:00 08:00 12:00 16:00 20:00 24:00
Time (hour)
 ACCEPTED MANUSCRIPT
22

20

PT
Temperature (°C)

RI
SC
18

U
AN
M
16
Proposed model, Average fluid

D
TE
EnergyPlus, Average fluid
14
EP
C Proposed model, Average borehole wall
AC

EnergyPlus, Average borehole wall

12
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
Time (month)
 ACCEPTED MANUSCRIPT

Proposed model EnergyPlus

3.5

PT
Heat pump cooling SPF

3.2

RI
SC
2.9

U
AN
M
2.6

D
TE
2.3 EP
C
AC

2.0
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
Time (month)
Borehole leaving fluid temperature (°C) 15 ACCEPTED MANUSCRIPT

PT
13

RI
U SC
AN
11
RC model, Δt=1 min

M
D
RC model, Δt=10 min

TE
9 RC model, Δt=30 min

EP
RC model, Δt=60 min

C
AC
Resistance model, Δt=10 min
7
00:00 04:00 08:00 12:00 16:00 20:00 24:00
Time (hour)
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
M
D
TE
EP
C
AC
 ACCEPTED MANUSCRIPT Tb
Rbg,3 Rbg,3
Cg,3 Cg,3
d u
Tg,3Tg,3
Rag,3 Rag,3

PT
Adg Aug b
Rbg,2
Cg,2 Rg,2

RI
Cg,2
Tdg,2 Tug,2

SC
Ap 2D Ap Rag,2 Rag,2

U
AN
Rbg,1 Rbg,1
Cg,1 Cg,1
Agg

M
po
Tug,1
2rpi

Af Af Tdg,1

D2r
Rag,1 Rag,1
Tgg

TE
d
Tpo Tupo
Rbp Ragg Rbgg
EP Rbp
Cp Tpd Cgg Tup Cp
C

Rpa Rpa
AC

rb

Rf Rf
Cf Tdf Tuf Cf

(a) (b)
 ACCEPTED MANUSCRIPT

PT
RI
U SC
AN
M
D
TE
EP
C
AC
 Start

Read input data including building

loads, indoor and outdoor conditions

Set time step, initial conditions and design parameters

BHE module for grid generation and setting coefficient matrices

ACCEPTED MANUSCRIPT

Update old temperature matrices

Calculate conditions of air entering and leaving the fan using

mass and energy balances for the mixing box and fan

PT
RI
Set an initial guess for Tf,in
U SC

cooling heating
mode
AN
M

Equation-fit model, Eqs. (13) Equation-fit model, Eqs. (16)

to (15) off and (17)
D
TE

PLF model, Eqs. (21), (22), PLF model, Eqs. (21), (22),
(25) and (27) (26) and (27)
C EP

Calculate QP , QSource,c & Qb,c , Calculate QP , QSource,h & Qb,h ,

Eqs. (12), (18) & (23) Eqs. (12), (19) & (24)
AC

Qb = 0

BHE model, solve Eqs. (1) to (8)

Updatee liquid loop temperatures

no yes
k+1 k k k-1
Tf,in = Tf,in |Tf,in - Tf,in |< t=t+ t

yes no
Print t > t max
End
results
 ACCEPTED MANUSCRIPT
6

PT
Building load (kW)

RI
SC
2

U
AN
M
0

D
TE
EP
-2 C
AC

-4
Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec Jan
Time(month)
 ACCEPTED MANUSCRIPT
6 6
Proposed model, Load side heat transfer rate
EnergyPlus, Load side heat transfer rate
5 5
Proposed model, Source side heat transfer rate

PT
Heat transfer rate (kW)

EnergyPlus, Source side heat transfer rate

RI
4 4
Proposed model, Compressor power

SC

Power (kW)
EnergyPlus, Compressor power

U
AN
3 3

M
D
2 2

TE
EP
1 1

C
0
00:00 04:00 08:00
AC 12:00 16:00 20:00
0
24:00
Time (hour)
 ACCEPTED MANUSCRIPT
6
Proposed model, Borehole heat transfer rate
5 EnergyPlus, Borehole heat transfer rate

PT
Heat transfer rate (kW)

Proposed model, Loop heat balance error

RI
4 EnergyPlus, Loop heat balance error

U SC
AN
3

M
D
2

TE
EP
1

C
AC
0
00:00 04:00 08:00 12:00 16:00 20:00 24:00
Time (hour)
 ACCEPTED MANUSCRIPT
15

13

PT
Temperature (°C)

RI
SC
11

U
AN
Proposed model, Borehole wall

M
9
EnergyPlus, Borehole wall

D
TE
Proposed model, Outlet fluid

EP
7 EnergyPlus, Outlet fluid

C
Proposed model, Inlet fluid

AC
EnergyPlus, Inlet fliud
5
00:00 04:00 08:00 12:00 16:00 20:00 24:00
Time (hour)
 ACCEPTED MANUSCRIPT

Highlights (Revised)

• A transient BHE model is used for the long-term simulation of a GSHPS.

• The transient model shows that the design borehole length can be reduced by 16%.

PT
• Steady-state model underestimates minimum fluid temperature leaving BHE by 1.52 °C.

• Neglecting the borehole thermal capacity affects annual energy performance by <1%.

RI
U SC
AN
M
D
TE
C EP
AC