ADVANCES IN MODELING OF
GROUNDSOURCE REAT
PUMP SYSTEMS
By
ANDREW D. CHIASSON
Bachelor of Applied Science
University of Windsor
Windsor, Ontario, Canada
1989
Master of Applied Science
University of Windsor
Windsor, Ontario, Canada
1992
Submitted to the Faculty of the
Graduate College of the
Oklahoma State University
in partial fulfillment of
the requirement for
the degree of
MASTER of SCIENCE
December 1999
ADVANCES IN MODELING OF
GROUNDSOURCE HEAT
PUMP SYSTEMS
Thesis Approved:
___UJ.......;;' ~~Ofl~ ~l~~e
ij
ACKNOWLEDGEMENTS
This thesis would not have been possible without the contribution of everal
people. To all of these people, I wish to express my gratitude.
First and foremost, I wish to thank my thesis adviser, Dr. Jeffrey Spitler. Hi
expertise and respected reputation in the thermal sciences field allowed for the funding of
many projects and the opportunity for me to gain knowledge and experience in everal
geothermal topics.
I also wish to thank to Dr. Simon Rees who served as an unofficial coadviser
along the course of this thesis. His expertise in computational methods was quite
valuable in helping me to achieve the project goals.
Thanks also to Dr. Marvin Smith who provided input and exp rimental data
regarding the supplemental heat rejecter studies
The review and contribution to this thesis by the committee members, Dr. R.
Delahoussaye and Dr. A. Ghajar are appreciated.
Last, but not least, I thank my colleagues in BlO, both here and gone, namely
Cenk Yavuzturk and Nagendra Jain, for their ideas and suggestions along the way.
111
Chapter
TABLE OF CONTENTS
Page
1. Introduction 1
1.1. Overview of GroundSource Heat Pump Systems 1
1. 1.1. GroundWater Heat Pump Systems 4
1.1.2. GroundCoupled Heat Pump Systems 5
1.1.2.1. Vertical GroundCoupled Heat Pump Systems 6
1.1.2.2. Horizontal GroundCoupled Heat Pump Systems 10
1.1.3. SurfaceWater Heat Pump Systems 12
1.104. Standing Column Well Systems 13
1.2. Thesis Objectives and Scope 15
1.3. The Overall Modeling Approach 17
2. A Preliminary Assessment of the Effects of GroundWater Flow
on ClosedLoop GroundSource Heat Pump Systems 20
2.1. Introduction 20
2.2. Coupled GroundWater Flow and Heat Transport 24
2.2.1. GroundWater Flow 24
2.2.2. Heat Transport in Ground Water. 26
2.2.3. Typical Hydraulic and Thermal Properties of Soils and Rock 27
2.2.4. Conduction versus Advection in Geologic Material 30
2.2.5. Numerical GroundWater Flow and Heat Transport Models 33
2.3. The Numerical Model 37
2.3.1. Model Description 37
2.3.2. The Finite Element Mesh 38
2.3.3. Boundary Conditions 39
2.3.4. Validation of the Numerical ModeL 41
2.4. Results and Discussion 43
2.4.1. SingleBorehole Simulations 43
2.4.2. Simulated InSitu Thermal Conductivity Tests 45
2.4.3. Borehole Field Simulations 49
2.5. Concluding Remarks and Recommendations for Future Work 57
3. A Model for Simulating the Performance of a Shallow Pond
as a Supplemental Heat Rejecter with ClosedLoop
GroundSource Heat Pump Systems 60
IV
Chapter Pag
3.1. Introduction 60
3.2 Heat Tran fer In Ponds ,61
3.2.1. General Overview 6 t
3.2.2. Existing Pond and Lake Models 62
3.3. Experimental Methods 66
3.3.1. Pond Description and Data Collection 66
3.3.2. Weather Instrumentation and Data Collection 67
3.4. Model Development 68
3.4.1. Governing Equations 68
3.4.1.1. Solar Radiation Heat Gain 69
3.4.1.2. Thennal Radiation Heat Transfer 71
3.4.1.3. Convection Heat Transfer at the Pond Surface 72
3.4.1.4. Heat Transfer to the Ground 75
3.4.1.5. Heat Transfer Due to Ground Water Seepage 76
3.4.1.6. Heat Transfer Due to Evaporation 77
3.4.1.7. Heat Transfer Due to the Heat Exchange Fluid 79
3.4.1.8. Solving the Overall Energy Balance Equation 83
3.4.2. Computer Implementation 84
3.5. Results and Discussion 87
3.5.1. Model Comparison to Experimental Results With No Heat Rejection 87
3.5.2. Model Comparison to Experimental Results With Heat Rejection 89
3.5.3. Model Application 92
3.6. Concluding Remarks and Recommendations for Future Work 97
4. A Model for Simulating the Performance of a Pavement
Heating System as a Supplemental Heat Rejecter with
ClosedLoop GroundSource Heat Pump Systems 99
4.1. Introduction 99
4.2. Heat Transfer In Pavement Slabs 100
4.3. Experimental Methods 101
4.3.1. Test Slab Description and Data Collection 10 I
4.3.1.1. Bridge Deck Test Section 10 I
4.3.1.2. Parking Lot Test Section 104
4.3.2. Weather Instrumentation and Data Collection 105
4.4. Model Development 106
4.4.1. Governing Equations 106
4.4.2. The Finite Difference Grid 108
4.4.3. Boundary Conditions 108
4.4.3.1. Solar Radiation Heat. Flux 110
4.4.3.2. Thermal Radiation Heat Flux 111
4.4.3.3. Convection Heat Flux at the Pavement Surfaces 111
4.4.3.4. Heat Flux Due to Rain and Snow I J3
4.4.3.5. Heat Transfer Due to the Heat Exchange Fluid 116
v
Chapter Page
4.4.4. Computer Implementation J 18
4.5. Results and Discussion 121
4.5.1. Model Comparison to Experimental Re ult With No Heat Rej ction 121
4.5.2. Model Comparison to Experimental Results With Heat Rejection 121
4.5.3. Model Application 126
4.6. Concluding Remarks and Recommendation for Future Work 131
5. Summary and Conclusions 133
References 140
APPENDIX A  Model and Experimental Uncerta.inty Analysis 145
vi
Table
LIST OF TABLES
Page
21. Typical Values of Hydraulic and Thermal Properties of Soil and Rock 29
22. Peelet Numbers Corresponding to Typical Value of
Hydraulic and Thermal Properties of Soils and Rocks 32
23. Case Summary of Simulated InSitu Ground Thermal Conductivity Te t 47
24. Summary of Borehole Field Design Parameters for Each Te t Case 50
31. Summary of Pond Performance for Example GSHP System Simulation 96
41. Summary of Pavement Heating System Performance for
Example GSHP System Simulation 131
VII
Figure
LIST OF FIGURES
Page
11. Schematic of cycles in a GSHP system in cooling and heating mode 2
12. A schematic of a groundwater heat pump ystem 5
13. A schematic of a vertical borehole and a horizontal
groundcoupled heat pump system 7
] 4. A schematic of a surfacewater heat pump system 13
15. A schematic of a standing column well system 14
21. Schematic of typical apparatus for insitu
ground thermal conductivity testing 2 I
22. Finite element mesh representing a single borehole 38
23. Finite element mesh representing a 16 borehole field 39
24. Average borehole fluid temperature vs. time for
coarse sand, fine sand, and shale 44
25. Average borehole fluid temperature for the 12 simulated
insitu ground thermal conductivity test cases 48
26. Predicted effective ground thermal conductivity values ver u
groundwater flow velocity for the simulated insitu thermal
conductivity tests 49
27. Hourly ground loads for the test building 5l
28. Design borehole depths versus groundwater flow velocity for
the simulated insitu thermal conductivity tests 52
VIIl
Figure Page
29. Annual maximum and minimum average borehole fluid
temperatures for the t 6 borehole field simulation S4
210. Temperature contours for Case 8 for the end of September
of years 1, 2, 5, and 10 S6
31. Heat transfer mechanisms in shallow ponds 63
32. Pond model component configuration 85
33. Pond model computer algorithm 86
34. Comparison of observed and simulated average pond temp~ratures
with no heat rejection for the (a) 2feet deep pond and
(b) 3.5 feet deep pond 88
35. Comparison of observed and simulated heat exchange fluid return
temperatures for the (a) 2feet deep pond and (b) 3.5 feet deep pond 90
36. Comparison of observed and simulated heat rejected to the
(a) 2feet deep pond and (b) 3.5 feet deep pond 91
37. System schematic for the example model of a GSHP system with a
shallow pond supplemental heat rejecter 92
38. Building thermal loads for the example building in Tulsa, OK. 94
39. Entering heat pump water temperatures for the example GSHP
system simulation with no pond and with a 2feet deep, 6000 ft2 pond 9S
41. Serpentine pipe configuration in a hydronicallyheated pavement
slab in (a) plan view and (b) crosssectional view 102
42. Slinky pipe configuration in a hydronicallyheated pavement
slab in (a) plan view and (b) crosssectional view 102
43. Heat transfer mechanisms in hydronicallyheated pavement slabs with
(a) no bottom exposure to the atmosphere and (b) bottom expo me
to the atmosphere 103
44. Finite difference cell geometry and notation 107
45. Model domain showing the finitedifference grid and boundary
conditions 109
IX
Figure Page
46. Pavement heating model component configuration ] 19
47. Pavement heating model computer algorithm 120
48. Comparison of observed and simulated slab surface temperature
for the parking lot test section with no heat rejection to the slab ]22
49. Comparison of observed and simulated slab surface temperature and
heat exchange fluid return temperatures for the bridge deck test ection
for the nights of (a) March 78,1996 and (b) March 2423, ]996 124
410. Comparison of observed and simulated (a) heat exchange fluid return
temperatures and (b) heat rejected to the parking lot te t section 125
411. System schematic for the examp}e model of a GSHP sy tern with a
pavement heating system supplemental heat rejecter. l27
412. Entering heat pump water temperatures for the example GSHP
system simulation with no pavement heating and with a 24,000 ft2
(2230 m2
) parking lot with pavement heating 129
x
Symbols
a
f3
cr
f.i
f.i'
v
B
P
p'
(j
A
C,Cp
d
D
D*
ff
g
h
H
I
K
k
L
f
m
• II m
NOMENCLATURE
= thermal diffusivity  ft2lhr (m2/s);
= solar absorptivity ()
= coefficient of thermal expansion  R1 (K1
)
= emissivity coefficient ()
=Euler's constant (0.577215664901)
= dynamic viscosity  lbrr/fthr (Pas)
= extinction coefficient  ft I (m I
)
= kinematic viscosity ~ fr2fhr (ro2/s)
= angle of incidence of sun's rays (radians)
= density  Ib/ft3 (kg/m3
)
= reflectance of pond surface ()
= StephanBoltzmann constant
= 0.1714 x 108 Btu/hrft2
OR
4 (5.67 x 1O!S W/m2_K4
)
= transmittance of solar radiation ()
= area  ft2 (m2)
= specific heat capacity  Btu/lbm<>P (J/kg°C)
= depth  ft (m)
= diffusion coefficient  ft2/s (m2/s)
= effective thermal diffusivity  ft2/s (m2/s)
= fouling factor  fr2oFlhrBtu (m2OCIW)
= acceleration due to gravity  ftls 2 (m/s2
)
= heat or mass transfer coefficient  Btu/hrft2OF (W/m2_oC);
= hydraulic head  ft (m)
= borehole depth  ft (m)
= hydraulic gradient  fUft (m/m)
= solar radiation flux on horizontal  Btu/hrft2 (W/m2
)
= hydraulic conductivity  ftls (mls)
= thermal conductivity  Btu/hrftOF (WIm°C)
= characteristic length  ft (m)
= length term described by Eskilson (1987)  ft (m);
= distance between nodes or thickness  ft (m)
= mass flow rate lbrn/hr (kgls)
= mass flux  Ibrnlhrft2 (kg/sm2
)
xi
Le
n
N
Nu
P
Pe
Pr
q
q*,
q
Q
Q*
r
R
R*
Ra
Re
Ss
T
t
Tom
U
V
v
V*
w
X,Z
,,1t
Subscripts
i
II
a
AB
b
c
circuit
d
= Lewis Number ()
= porosity ()
= quantity ()
= Nusselt Number ()
= perimeter  ft (m);
= pressure  atmospheres
= Peelet number ()
= Prandtl Number ()
= heat transfer rate  Btuthr (W);
= specific discharge  ftls (m/s)
= ground thermal load  Btuthr (W]
= heat flux  Btu/hrft2 (W/m2
)
= volumetric flow rate  ft3/s (m3/s)
= heat source/sink term  op/s (oC/s)
= radius  ft (m);
= unpolarized solar radiation component ()
= retardation coefficient ();
= thennal resistance  ft2
 OPlhrBtu (m2
°CfW)
= groundwater recharge SI
= Rayleigh Number ()
..... Reynolds Number ()
= specific storage coefficient  ft 1 (m I
)
= temperature  OF (OC)
= time  s
= undisturbed mean temperature of the ground  oP (oC)
=overall heat transfer coefficient  Btu/hrft2
OF (W/m2_oC)
= volume  fr3 (m3
)
= average linear groundwater velocity  ftls (m/s)
= velocity  ftls (m/s)
= volumetric flow rate  fe/s (m3/s)
= humidity ratio  Ibm water/Ibm dry air (kg water/kg dry air)
= rectangular coordinates
= time step  s
= perpendicular component
= parallel component
= absorbed component of solar radiation
= transfer from material A (water) to material B (air)
= beam radiation;
= borehole
= convection
= now circuit or spool
= diffuse radiation;
XII
= diffusion
db = dry bulb
dp = dew point
eff = effective
fg = latent heat of vaporization
fluid = heat exchange fluid
= nodal index;
= pipe inside
i,j coordinate indices
if latent heat of fusion
in = inlet
R = liquid phase
m,n = coordinate locations
0 = pipe inside;
= background or initial
out = outlet
r = thermal radiation;
= refraction
s = solid phase
surf = surface
sW,b = steadystate, borehole wall
t = current time step:
= total
tL1t = previous time step
W = water;
= injected/extracted water
wb = wet bulb
X,Z = coordinate indices
(x,) ) =coordinate index for a surface node
XIII

1. Introduction
1.1. Overview of GroundSource Heat Pump Systems
Groundsource heat pump (GSHP) systems (also referred to as geothermal heat
pump systems, earth energy systems, and GeoExchange system) have received
considerable attention in the recent decades as an alternative energy source for re idential
and commercial space heating and cooling applications. GSHP application are one of
three categories of geothermal energy resources as defined by ASHRAE (1995). The e
categories are: (l) hightemperature (>302 of (> 150 °C) ) electric power production, (2)
intermediate and lowtemperature «302 OF «150 °C» directuse applications, and (3)
GSHP applications (generally <90 OF «32 °C». The GSHP application are
distinguished from the others by the fact that they operate at relatively low temperature.
The term "groundsource heat pump" has become an allin~lu~iv.eJerm to
describe a heat pump system that u es the earth, ground water, or surface water as a heat
source and/or sink. GSHP systems consist of three loops or cycles as shown in Figure II.
The first loop is on the load side and is either an air/water loop or a water/water loop,
depending on the application. The second loop is the refrigerant loop inside a watersource
heat pump. Thermodynamically, there is no difference between the wellknown
vaporcompression refrigeration cycle and the heat pump cycle; both sy. terns absorb heat
at a low temperature level and reject it to a higher temperature level. The difference
between the two systems is that a refrigeration application i only concerned with the low
temperature effect produced at the evaporator, while a heat pump may be concerned with
both, the cooling effect produced at the evaporator as well as the heatin effec produced
,at the condenser. In these dualmode GSHP systems, a reversing valve is used to switch
between heating and cooling modes by reversing the refrigerant flow direction. The third
loop in the system is the ground loop in which water or an antifreeze solution exchanges
heat with the refrigerant and the earth. t J ( I'
r
To
H.... LOId
~.
Ii_ Elu:IlouaCYcle
t. l
'Hal
P1IIIIp i
To
CooiDt! LOId

1 , .
<a) (b)
Figure 11. Schematic ofcycles in a GSHP system in (a) cooling mode
and (b) heating mode.
Efficiencies ofGSHP systems are much greater than conventional airsource heat
pump systems. A higher COP (coefficient ofperformance) can be achieved by a GSHP
because the source/sink earth temperature is relatively constant compared to air
temperatures. Additionally, heat is ab orbed and rejected through water which i a more
desirable heat transfer medium because of its relatively high heat capacity.
GSHP systems rely on the fact that, under normal geothermal gradi _nts of about
0.5 °FIlOO ft (30 °CIkm) (Grant et aI., 1982), the earth temperature i roughJy con tant in
a zone extending from about 20 ft (6.1 m) deep to about 150 ft (45.7 m). deep (Hart and
Couvillion, 1986). This constant temperatura interval within the earth i the re ult of a
complex interaction of heat fluxes from above {the un and the atmo phere) and from
below (the earth interior). As a result, the temperature of this interval withjn th earth i
approximately.:qu_al !~~ av ra eannuaJ air tern e tgJ~JHart and Couvillion, 1986).
Above this zone (less than about 20 feet (6.1 m) deep), the earth temperature is a damp~_
version of the air temperature at the earth's surface. Below this zone (greater than about
150 ft (45.7 m) deep), the earth temperature begins to rise according to the natural
geothermal gradient
ASHRAE (1995) groups GSHP systems into three categorie ba ed on th heat
source/sink used. A fourth category is added here for the sake of com l~~
.  . These

categories are: (1) groundwater heat pump (GWHP) systems, (2) groundcoupled heat
pump (GCHP) systems, (3) surface water heat pump (SWHP) systems, and (4) standing
column well (SCW) systems. Each of these is discussed in the foHowing sub ections.
4
1.1.1. GroundWater Heat Pump Systems
GWHP systems, also referred to a openloop y tern are th original type of
GSHP system. The first GWHP system were installed in the late 1940 K vanaugh and
Rafferty, 1997). GWHP systems are oot the~f~__ofthis the i I 0 they will only be
brieflydescribed here.
A schematic of a GWHP system is shown in Figure 12. In GWHP sy tern ,
c.?}!v otional water wells and ~_~.u._PU.ID.I2S are used to su I round water to a heat pump
or directly to some application. Corrosion protection of the heat pump may be nece sary
._~~~......;:ality i~oor. The "used" ground water i typically di charged
to a suitable receptor, such as back to an aquifer, to the unsaturated. zone (a shown in
Figure 12), to a surfacewater bod , or to a sewer. Desi consid ratiQI) for GWHP
~' 
systems arefairly..w.e:ll:.established; welldrilling technologie and welltesting m th
have been wellknown for decades. Design con iderations include: graundwat r
availability, groundwater chemical quality, and groundwater di asal method.
The main advantage of GWHP systems is their low cost, simplicity, and small
amount of ground area required relative to other GSHP and conventional ysterns.
Disadvantages include limited availability and poor chemical quality of ground water in
some regions. With growing environmental concerns over recent decade, many legal
.~ ... 
issues have arisen over ground water withdrawal and reinjection in some locahtie .
s
I l
Figure 12. A schematic ofa groundwater heat pump system.
1.1.2. GroundCoupled Heat Pump Systems
GORP systems, also referred to as closedloop groundsource heat pump y tern ,
were ieered in the 1970s. Their main advantage over their waterwell Jm~eGeS~;s.ls
that they eliminate the problems ~as::.:;,.c::.;i~...... w.;it,.... ground water quality and availability and
they generally require much less pumping energy than water well systems because there
is less elevation head to overcome. GCHP systems can be installed at any location where
drilling or earth trenching is feasible.
In GCHP systems, heat rejection/extraction is accom l' . circulating a heat
exchange puid through a piping system buried in the earth. This fluid is either pure wat~
or an antifree solution and is typically circulated through highdensity polyeth lene ...
(HOPE) pipe installed in vertical boreholes r horizontalGh .' shown in Figur 13.
Thus, these systems are further subdivided into vertical GCHP y tern and horizontal
GCHP systems.
1.1.2.1. Vertical GroundCoupled Heat Pump Systems
Vertical borehole GCHP systems are the primary focus ofthi entir the i .
Therefore, they are described in some detail here and their d..:.es::;;;i..........;;;.;;.;.=~!..:..::explained,
la in the foundation for the motivation of this study.
In vertical borehole GCHP systems, ground heat exchanger configurations
typically consist of one to tens of boreholes each containing a Vshaped pipe through
which the heat exchange fluid is circulated. Some Swedish applications use borehole
inclined from the vertical. A number of borehole to borehole plu bin m t are
possible. Typical Vtubes have a diameter in the range of 3A in. (19 mm) t 1 Yz in. (38
mm) and each borehole is typically 100 ft (30.5 m) to 300 ft (91.4 m) d ep with
diameter ranging from 3 in. (76 mm) to 5 in. (127 mm). The borehole annulu i
generally backfilled with a material that prevents contamination of ground water.
7
(a> I t
(b)
Figure 13. A schematic of (a) a vertical borehole groundcoupled heat
pump system and (b) horizontal groundcoupled heat pump
system.
8
The design of vertical ground heat exchanger i complicat d by th variety of
geological formations and ro rti~ that affect their thermal performance (ASHRAE,
1995). Proper subsurface cha cteri~is not economically fe ibJe for every project.
One of the fundamental tasks in the design of a reliable GCHP y tern i prop rly izing
the groundcoupled heat exchanger length (i.e. depth of boreholes). Particularly for large
systems, an extensive effort is made to design the ground loop heat xchanger 0 that
they are not too large (resulting in too high of a first cost) or too smaJl (re ulting in th
building's thennalload not being met).
In the early days of GCHP technology, the task of sizing the groundloop heat
exchanger was accomplished using rules of thumb (i.e. 250 feet of bore length per ton of
heating or cooling capacity). These rules were slightly modified on a_casebyca e basi
using some estimates of thermal conductivity of the formation or using previou design
experience, but additional costs of more detailed testing or calculation wa judg d to
outwei h the costs of a conservative design. This relatively imple approach prov d to
be successful in most residential and other small applications, but in larger cal
commercial and institutional applications, some groundloop heat exchanger fai led to
meet their design loads after the first few years of operation. Further, the ractice of
greatly overdesi nin large GCHP systems was found to be unacceptable becau e the
~
first costs were sim Iy not competitive with the first co ts of conventional y tern .
Consequently, intensive research regarding methods to optimize groundloop heat
exchanger design has been ongoing for the last decade.
Simple approaches to sizing the groundloop heat e cb nger in larg  cale
applications are inadequate mainly becau e the beat tran fe proce: '"'
in th gr und are
com licated by thermally interacting borehole and hourly periodic heat
extraction/injection pulses. Annual heat rejection and. beat extraction are ually not
equal and peak temperatures rise or fall over a number of year . A a r~t, groundloop
heat exchanger designers need to consider hourly heating and cooling load of tb
building and need to perlorm some simulation of the groundloop temperature over the
lifecycle of the building. Recent research efforts have produced several methods and
computer software'programs for this purpose. However, none of the program con ider
the effects of ground water flow on groundloop beat exchanger performance~ tbe
effects have not been well understood, perhaps becau,e of the lack of relevant
investigations. This is the topic of Chapter 2 of tbi thesis.
Another challenge in the design of GCHP system ari es from the fact that mo t
commercial and institutional buildings, even in moderate climates, are generally coolingdominated
and therefore reject more heat to the ground than they extract over the annual
cycle. This load imbalance may require the heat exchanger length to be significantly
greater than the length required if the annual loads were balanced. A a result, the GSHP
system may be eliminated from consideration early in the design phase of the project due
to excessive first cost. This has given ri _to the conce f "supplemental heat rejecter"
~
or socalled "hybrid GSHP systems".
10
Supplemental heat rejecter have been 'Dte rated into buiLdi
~balance the ground loads and therefore reduce the nece ar J~n ~ of the
groundloop heat exchanger. In some applications, the excess ~t that would otherwi e
build up in the ground may be used for domestic hot water heaters, car washes, and
pavement heating systems. In cases where the exce heat cannot be u ed b n ficiaJly,
~
conventional cooling towers or shallow ponds can provide a co teffective mean to
reduce heat exchanger length.
Design of these supplemental components adds to the challenge of de igning the
overall hybrid GCHP system because of their highly transient nature. Heat rejection
systems are likely to operate more often during the nighttime hours or when the building
is not in use. Therefore, it is essential that the hourly (or less) behavior of these systems
be examined during their design phase. These are the topics of Chapters 3 (shallow
ponds) and Chapter 4 (pavement heating systems) of this the i .
1.1.2.2. Horizontal GroundCoupled Heat Pump Systems
In horizontal GCHP systems, ground heat exchanger configuration typically
consist of a series of parallel pipe arrangements laid out in dug trenches or horizontal
boreholes about 3 ft (0.91 m) to 6 ft (l.83 m) deep. A number of piping arrangements are
possible. "Slinky" configurations (as shown in Figure 13 (b» are popular and imple to
install in trenches and shallow excavations. In horizontal boreholes, straight pipe
configurations are installed. Typical pipes have a diameter in the range of % in. (19 mm)
11
to 1 Y2 in. (38 mm) aDd about400ft(12L.9 m) to 600 ft (182.9 m) of pipe i in tall d P r
ton of heating and cooling capacity.
The thennal characteristics of horizontal GCHP sy terns are imjlar to tho of
vertical ODes. The maiD difference is that horizontal groundLoop heat exchange ar
more affected by weather and air temperature fluctuations because of their roximity to
the earth's surface. This may result in Larger loop temperature fluctuation and therefore
lower heat pump efficiencies. Recent researc activit~s have focus ed on u iog th e
systems as supplemental heat rejecters with vertical borehole GCHP system . A P cific
application (i.e. the use of a shallow pavement heating system) is the focu of Chapter 4
of this thesis.
Aside from the invention of the Slinky coil itself and the use of these system as
suppLemental heat rejecters, horizontal systems have received much Ie s att ntion than
vertical systems with res ect to recent research efforts. _This.may be due to the fact that  ~ ....... ~ ...
vertical systems tend to be preferred in larger applications since much les ground area i
required. Also, since horizontal systems are installed at shallow depth, geologic site
characterization is relativeLy easier because soils can be readily seen and sampled.
Further, overconservative designs are not as cost prohibitive as with vertical borehole
designs because of the relatively low installation costs of the heat exchanger pipe.
1
1.1.3. SurfaceWater Heat Pump Systems
The third category of GSHP systems is the surfacewater heat pump (SWHP)
system. A specific application of SWHP systems (i.e. the use of a shallow pond as a
supplemental heat rejecter in vertical GCHP systems) i~focu of Chapter 3 of thi
thesis.
A schematic of a SWHP system is shown in Figure 14. The surfacewater heat
exchanger can be a closedloop or openloop type. Typical clo edloop configuration
are the Slinky coil type (as shown in Figure 14) or the loose bundle coil type. In the
closedloop systems, heat rejection/extraction is accomplished by circulating a heat
exchange fluid through HDPE pipe positioned at an adequate depth within a lake, pond,
reservoir, 0~uita~9 en cha.D1!~l. Typical pipe diameters range from % in. (19
mm) to I Y2 in. (38 mm) and a length of 100 feet (30.48 m) to 300 feet (91.44 m) per ton
of heating or cooling capacity is recommended by ASHRAE (1995b), depending on the
climate. In openloop systems, water i extracted from the urfacewater body through a
screened intake area at an adequate depth and is discharged to a suitable recepLor.
Heat transfer mechanisms and the thermal characteristics of surfacewater bodies
are quite different than those of soils and rocks. This subject will be further discus ed in
Chapter 3 of this thesis. At the present time, design tools for surfacewater heat pump
systems arein their developmental infancy (Kavanaugh and Rafferty, 1997). However, ._~.
many successful instaJlations are currently in operation and some guideline do exist. In
short, the loop design im
specifying adequate dian
and locating the coil at a
Figure 14
1.1.4. Standing Columa
The fourth categor
system. These systems an
recently receiving much a1
only briefly discussed hen:
A schematic ofm sew system is shown in Figure 15. This type ofGSHP draws
water to a heat pump from a standing column ofwater in a deep well bore and returns the
water to the same well These systems, primarily installed in hard rock areas, use
uncased bo eholes with typical diameters of about 6 in. (15.24 cm) and depths up to l500
feet (4~7.2 m). The uncased. borehole allows the. heat exchange fluid to be in direct
contact with the earth (unlike closedloop heat exchangers) and allows groupd water
~filtrati9n over the entire length of the borehole. Properly sited and designed, sew
systems have been shown to have significant installation cost savings over closedloop
GCHP systems. Design guidelines for sew systems are currently under development.
I.
Figure 15. A schematic ofa standing column well system.
1.2. Thesis Objectives and Scope
This s . deals with the modeling of vertical do edloop and
source heat pump systems. The challenges associated ith the de ign
were discussed in the previous section. A considerable amount of rese
decade has bee
....c.._~
this study is part of those efforts.
There are thTee primary objectives of this study. These are to:
(1) examine the effects of groundwater flow on closedloop G~
(2) develop a design and simulation tool for modeling the perfo
shallow pond as a supplemental heat rejecter with closedIoc
and
(3) develop a design and simulation tool for modeling the p rfo:
hydronic pavement heating system as a supplemental heat re
closedloop GSHP systems.
Chapter 2 of this thesis ad'd.r....~_r the first objective. Given the t
constraints of finding a suitable site with significant ground waterflow,
that site, and collecting and analyzing data at the site over a time period
computer modeling study was conducted as a preliminary assessment of
groundwater flow on closedloop GSHP systems. Hydraulic and therm
soils and rocks were compiled and AQUA3D, a commerciallyavailable
groundwater flow and heat transport model developed by Vatn . kil Con ultin
Engineers, Inc., Reykjavik, Iceland, was used to simulate the impact of ground\!
flow on the average heat exchange fluid temperature in single and multiple bor
systems. The impact of groundwater flow on the estimation of oil/rock th rma
conductivity from insitu test data w=as::....=a1=s;:.o_=='~
Chapter 3 of this thesis addresses the second objective. The development
validation of a design and simulation tool for modeling the perfonnance of a hal
pond as a supplemental heat rejecter with closedloop GSHP systems is presentee
.:...
model has been developed in the TRNSYS modeling environment and can theref(
coupled to other GSHP system component models for shorttime step (hourly or
minutely) system analyses. TRNSYS, developed at the University of Wi consin
Madison, is a widelyused modular transient thermal systems simulation program
each system component is described mathematically by a FORTRAN ubroutin .
simple language, the components are connected together in a manner analogou t<  ducting, and wiring in a physical system (Duffie and Beckman, 1991). An examp
application of the pond model is also presented.
Chapter 4 of this thesis addresses the third objective. The development ani
validation of a design and simulation tool for modeling the performance of a hydn
pavement heating system as a supplemental heat rejecter with closedloop GSHP "
is presented. This model has also been developed in the TRNSYS modeling envir,
17
for use in shorttime step system analyse . An example application of the pavem nt
heating model is al 0 presented.
Finally, Chapter 5 of this tbesi~e ~~mclu ion f the~
studies.
1.3. The Overall Modeling Approach
Each of the objectives of this thesis deals with the development and/or application
of a model. Since the term "model" can be used loosely, some definitions and  approaches as applicable to this study are described here.
~
A ''model'' is a physical or a mathematical representation of an actual system.
This study deals only with mathematical representations of system . American ociety
for Testing and Materials (ASTM) defines a mathematical model a "mathematical
equations expressing the physical system behavior and including simplifying
assumptions". Mathematical models are solved analytically or numerically u ing manual
or computer methods.
The overall modeling approach consists of six stages: (1) define the purpose of
the model, (2) develop a conceptual model of the system, (3) develop or define the
mathematical model of the system, (4) implement the solution method, (5) validate the
model, and (6) apply the model. Each of these is described in the following paragraphs.
1
The first stage in the modeling approach i to clearly d fin the purpo and
objectives of model. This helps to detennine the level of detail and accuracy desired by
the model and helps in making decisions regarding the resource needed.
The second stage in the modeling approach is to develop a conceptual model. of
the system. ASTM defines a conceptual model as "an interpretation or working
description of the characteristics and dynamics of the physical system". The purpose of
the conceptual model is to describe the system by a set of assumptions and concepts that
can be evaluated mathematically.
The third stage in the modeling approach is to develop a mathematical model. In
this stage, the conceptual model is translated to mathematical equations that can be
solved for the desired unknowns. The solution method and limiting a sumption or
simplifications are also identified.
The fourth stage in the modeling approach is to implement the solution method to
solve the mathematical equations. With respect to this thesis, this stage involved using
computer methods with a combination of commerciallyavailable software and
FORTRAN code developed specifically for this study.
The fLfth stage in the modeling approach is to validat th mod j. With r p t to
this thesis, this stage involved comparing model re uLts, where applicable, to an
analytical solution or to experimental data.
The sixth stage in the modeling approach is to finally use the model to analyze the
performance and behavior of actuaJ thermal systems.
20
2. A Preliminary Assessment of the Effects of GroundWater Flow on ..
ClosedLoop GroundSource Heat Pump Systems
2.1. Introduction
One of the fundamental tasks in the design of a reliable groundcoupled heat
pump system is properly sizing the groundcoupled heat exchanger length (i.e. depth of
boreholes). Recent research efforts have produced several method and commerciallyavailable
design software tools for this purpose (Ingersoll, 1954; Kavanaugh, 1984:
Eskilson, 1987; IGSHPA, 1991; Spitleret aI., 1996; and Kavanaugh and Rafferty, 1997).
All of these design tools are based on principles of heat conduction and rely on some
estimate of the ground thermal conductivity and volumetric specific heat. These
parameters are perhaps the most critical to the system design, yet adequately determining
them is often the most difficult task in the design phase.
Methods of determining the thermal properties of the ground have al 0 be n the
subject of considerable recent research (Eklof and Gehlin, 1996; and Au tin et aI., 2000).
Current methods range from estimating values from published data to conducting
laboratory experiments on soil/rock samples to conducting singleborehole, in. itu field
tests. In general, thermal property values derived from insitu field tests are mo t
representative because the values are sitespecific and a larger volume of material i
evaluated under more realistic conditions than is possible in the laboratory. The typical
field procedure in insitu tests is to measure the temperature response of a fluid flowing
through a ground heat exchanger in a single borehole. A schematic of the typical field
21
apparatus is howD in Figure 21. Tbe fluid temperature at eguJar time iot rv and beat
added to the fluid stream are recorded over the course of the test.
Pipe
Insulation
Temperature
Sensol'S
Portable
Unit
Flow
Meier
Waler
Circulating
Pumps
UTube Heat Exchanger
Figure 21. Schematic of typical apparatus for insitu ground thermal
conductivity testing.
Determination of thermal conductivity from temperaturetime data i an inverse
problem. Several analytical and numerical methods exist for interpreting the data et.
These methods will not be disoussed in detail here, but include the "cylinder ource"
analytical solution (Carslaw and Jaeger, 1946), the "linesource" analytical olution
(Kelvin, 1882; Ingersoll, 1954), numerical solutions (Mei and Emerson, 1985; Muraya et
aI., 1996; and Rottmayer et aI., 1997), and numerical solutions with parameter e timation
22
(Shonder and Beck, 1999, Au tin et al., 2000). Each of th e method i ba ed on
Fourier's Law of heat conduction. .1 ' 't
A further complication in the design of groundcoupled heat pump yst ms i the
presence of ground water. Where ground water is present, flow will occur in re ponse to
hydraulic gradients and the physical process affecting heat transfer i.n the ground i
inherently a coupled one of heat diffusion (conduction) and heat advection by moving
ground water. In general, groundwater flow can be expected to be beneficial to the
thermal perfonnance of closedloop ground heat exchangers ince it will have a
moderating effect on borehole temperatures in both heating and cooling mod
Complications in the borehole field design process due to the pre ence of flowing
ground water arise from the fact that both current insitu conductivity te t data
interpretation methods and groundloop heat exchanger design method. are based on
models that only consider heat conduction. Therefore, groundwater flow may impact the
design process in two ways: (I) thermal conductivities derived from iin itu test may
appear artificially high and (2) borehole fields designed from artificially high thermal
conductivity values may be over or underdesigned. An unusually high thermal
conductivity value was determined from insitu test data at a site in Minnesota where
significant groundwater flow was believed to occur (Remund, 1998).
The objectives of the work presented in this chapter have been to make a
preliminary examination of the effects of ground water flow on both insitu ground

thennal conductivity measurements and longterm borehole field performance. Thi h
been attempted firstly by examining the range of hydrogeological condition that might
be expected and estimating the order of magnitude of the corresponding groundwater
flows. A simple method of examining the importance of heat advection from groundwater
flow is then presented.
A finiteelement numerical groundwater flow and heat transport model ha been
used to simulate and observe the effects of groundwater flow on the average fluid
temperature in a single Vtube borehole in various geologic materials. The model wa
used to simulate several insitu ground thermal conductivity tests, and thermal
conductivities were derived from these data using a standard approach. For each te t
case, the derived thermal conductivities, along with the thennal loads from an actual
building, were used to design a hypothetical multiborehole field by employing
conventional design tools and procedures. For different sets of hydrogeological
conditions, a numerical model of the whole borehole field was u ed to simulate its longterm
performance. Conclusions are presented on the ability of conventional design
procedures to correctly predict the longterm performance of closedloop ground heat
exchangers under different groundwater flow conditions.
24
2.2. Coupled Ground Water Flow and Heat Transport
2.2.1. GroundWater Flow
Underground water occurs in two zones: the unsaturated zone and the aturated
zone. The tenn "ground water" refers to water in the saturated zone. The urface
separating the saturated zone from the unsaturated zone is known as the "water table". At
the water table, water in soil or rock pore spaces is at atmospheric pre sure. In the
saturated zone (below the water table), pores are fully saturated and water exist at
pressures greater than atmospheric. In the unsaturated zone, pore are only partially
saturated and the water exists under tension at pressures less than atmospheric. In this
paper, we deal only with water in the saturated zone.
Ground water is present nearly everywhere, but it is only available in u able
quantities in aquifers. An "aquifer" is defined by Driscoll (l986) as a formation, group
of formations, or part of a fonnation that contains sufficient saturated permeable material
to yield economical quantities of water to wells and springs. Aquifer are de cribed a
being either confined or unconfined. Unconfined aquifers are bounded at their upper
surface by the water table. Confined aquifers are bounded between two layers of lower
penneablity materials. In practice, the boreholes of groundloop heat exchanger may
partially penetrate several geologic layers.
The governing equation describing flow through porous media is Darcy's Law:
dh
q=Kdx
25
(21)
• I
where q is the specific discharge (volume flow rate per unit of cross ectional area), K i
the hydraulic conductivity, and h is the hydraulic head. The specific di charge i related
to average linear ground water velocity, v, by:
v =!l..
n
(22)
where n is the porosity and is introduced to account for the difference between the unit
crosssectional area and the area of the pore spaces through which the ground water flows
(Freeze and Cherry, 1979; Fetter, 1988).
By applying the law of conservation of mass to a control volume and by making
use of Darcy's Law (Equation 21), an equation defining the hydraulic head di tribution
can be derived. Transient groundwater flow with constant den ity can then be expre ed
in Cartesian tensor notation as:
s ah~(K.'~J=R·
or at ax. IJ ax I /
(23)
Since ground water at 1lOoF (43.3°C) (an extreme temperature limit expected in GSHP
applications) has a specific gravity of approximately 0.99 J , the assumption of constant
density flow for lowtemperature geothermal applications may be considered valid.
26
2.2.2. Heat Transport in Ground Water
Heat can be transported through a saturated porou medium by the following thre
processes:
temperature of rock/water matrix. It is the second term in Equation 24 that repre ents
applying the law of conservation of energy to a control volume, an equation for heat
differential equation of the advectiondispersion type (Freeze and Cherry, 1979). By
nR (24)
aT
+v aT ~(D aTJ=Q* at I ax ax. IJ ax I I I
transport in ground water can be found and can be expressed as:
(3) heat transfer through the liquid phase by advection.
where the velocity, Vi is determined from the solution of Equation 23 and T is the
(2) heat transfer through the liquid phase by conduction, and
(1) heat transferthrough the solid phase by conduction,
advection of heat by the ground water and couples Equations 23 and 24 together. If the
ground water velocity is zero, Equation 24 reduces to a form of Fourier's Law of heat
The governing equation describing mass or heat transport in groundwater i a partial
conduction.
The diffusion coefficient tensor Dij is modeled here as an effecti ve thermal
diffusivity given by:

DO = kerr
PlCl
• t
(25)
27
The effective thermal conductivity k~ff is a volumeweighted average th rmal conductivity
of the saturated rock matrix and can be expressed using the porosity as:
(26)
It is necessary to distinguish between the conductivity and thermal capacity of the water
and soil/rock in this way to account for the fact that heat is stored and conducted through
both the water and soil/rock, but heat is only advected by the water. Similarly, it i
necessary to define a retardation coefficient R accounting for retardation of the thermal
plume which results from differences in the liquid and solid volumetric heat capacities:
R = 1+(1 n)csps
n Clpe
2.2.3. Typical Hydraulic and Thermal Property Values for Soils and Rocks
(27)
In assessing the significance of groundwater flow to closed loop heat exchanger
performance, the question arises as to what locations have significant groundwater flow.
Darcy's Law indicates that flow is dependent on both the local hydraul ic gradient and the
hydraulic conductivity of the geologic material. Heat transfer is dependent on the flow
velocity and the thermal properties of the material. It is therefore useful, in making a
preliminary assessment of the significance of groundwater flow, to consider the range of
naturallyoccurring soil and rock properties and possible values of hydraulic gradient.
2
Naturallyoccurring ranges of value of hydraulic and thermal prop rti of oil
and rocks are summarized in Table 21. Values of hydraulic gradient are omewhat more
sitespecific; the United States Environmental Protection Agency (1996) report a typical
range of hydraulic gradient values of 0.000 1 to 0.05.
Some specific examples of natural ground water velocities include: 1796 ft/yr
(547.5 m/yr) to 7185 ft/yr (2190 rnlyr) under a hydraulic gradient of 0.002 to 0.012 in the
Snake River Group basalt, Idaho, USA (Lindholm and Vaccaro, 1988); 361 ft/yr (110
rnIyr) in the High Plains sand and gravel aquifer, western central USA (Week and
Gutentag, 1988); and 1.3 x 103 ftlyr (4.0 x 104 rnIyr) to 1.50 x 102 ft/yr (4.6 x 103
rnIyr) in glacial clay soils in Southern Ontario, Canada (Stephenson et aI., 1988). Local
pumping activities may further increase groundwater flow rates in aquifers.
The thermal properties of soils and rocks are function of mineral content,
porosity, and degree of saturation. Of these, porosity may be considered the mo t
important property simply because of the origin and nature of soils and rock. Rocks
originate under higher heat and pressure environments than soils and consequently
generally possess lower porosities. Lower porosities in rocks result in higher contact area
between grains and therefore higher thermal conductivities than soils, regardless of
mineral content. For saturated materials, increased porosity results in increased heat
capacities and therefore lower thermal diffusivities.
TABLE 21.
Typical Values of Hydraulic and Thermal Properties of Soils and Rocks
Porous MedIum Hydraulic Properties Thermal ProDertles
Hydraulic Conductivity' Porosity' Velocity" Thermal Conductivity'" Volumetric Heat Capacity'"
(K) (n) (v) (k) (Psc.)
ft/s H ft/yr Btulhr.ftof BtulffOF
Im/s) (mivr) rN/mOC) (J/m3OC)
Range Geometric Range Arithmetic Range Arithmetic Range Arithmetic
Averaae Averaae Averaae Averaae
So/ls
Gravel 9.84E04  9.84E02 9.84EQ3 0.24  0.38 0.31 1.00E+04 0.40  0.52 0.46    2.09E+01
3.00E04  3.00E02 3.00E03 3.05E+03 (0.70)  (0.90) (0.80)    (1.40E+06)
Sand (coarse) 3.0E06  2.0E02 2.4E04 0.31  0.46 0.385
.
1.98E+02 0.40  0.52 0.46    2.09E+Ol
I (9.0E071  (6.0E03) (7.3E051 1(6.01 E+Oll (0.701  (090) (0.80)    (1.40E+06)
Sand (fine) 6.6E07  6.6E04 2.1E05 0.26  0.53 0.40 1.66E+Ol 0.40  0.52 0.46    2.09E+01
I (2.0E07)  (2.0E04) (6.3E06) I(5.05E+00\ (0.70\  (0.90\ (0.801  .  (1.40E+06)
Silt 3.3E09  6.6EOS 4.6E07 0.34  0.61 0.475 3.D8E01 0.69  1.39 1.04 3.58E+01  4.92E+Ol 4.25E+Ol
I(1.0E09)  (2.0E05) (1.4E07) (9.40E02) (1.20)  (2.40) (1.80) (2.40E+061  (3.30E+061 (2.85E+06)
Clay 3.3E11  1.SE08 7.1El0 0.34  0.60 0.47 4.78E04 0.49  0.64 0.56 4.47E+Ol  5.37E+Ol 4.92E+Ol
I(1.0Elll  (4.7E·091 (2.2El0) (1.46E041 (0.85)  (1.10) (0.98) (3.00E+06)  (3.60E+06) (3.3E+061
Rocks
Limestone, Dolomite 3.3E09  2.0E05 2.5E07 0  0.20 0.10 8.02EOl 0.87  1.91 1.39 3.17E+02  8.20E+Ol 1.99E+02
1(1.0E09)  (6.0E·OS) (7.7E08) (2.44EOl) (1.50)  (3.30) (2.40) (2.13E+071  (S.50E+06) (1.34E+07)
Karst Limestone 3.3EQS  3.3E02 3.3E04 0.05  0.50 0.275 3.76E+02 1.44  2.48 1.96 3.17E+02  8.20E+Ol 1.99E+02
1(1.0E06)  (1.0E02) (1.0E041 I(1.15E+02) (2.501  (4.30) (3.40) (2.13E+07)  (5.50E+06) (1.34E+07)
Sandstone 98El0  2.0E05 1.4E07 0.05  0.30 0.18 2.51 EOl 1.33  3.76 2.54 3.17E+Ol  7.46E+Ol 5.31E+Ol
I(3.0El01  (6.0E061 (4.2E08) (7.6SE021 (2.30\  (6.50\ (4.401 (2.13E+06)  (5.00E+061 (3.56E+061
Shale 3.3E13  6.6E09 4.6E11 0  0.10 0.0525 2.79E04 0.87  2.02 1.44 3.54E+Ol  8.20E+01 5.87E+Ol
I(1.OE13)  (2.0E09) (1.4E11) (8.50ED5) (1.50)  (3.50) (2.50) (2.38E+06)  (5.50E+06) (3.94E+06)
Fractured Igneous 2.6E08  9.8EQ4 5.1 EOS 0  0.10 0.05 3.21E+Ol 1.47  3.83 2.65   3.28E+Ol
and Metamorphic I(8.0E09)  (3.0E04) (1.5E06) I(9.78E+00) (2.50)  (6.60) (4.581    (2.20E+06)
Unfractured Igneous 9.8E14  6.SEl0 8.0E12 0  0.05 0.025 1.01E04 1.47  3.83 2.65    3.28E+Ol
and Metamorphic I(3.0E13) • (2.0E10\ (2.4E121 (3.09EOSl (2.501  (6.60) (4.581    (2.20E+06)
Notes: Thermal conductivity values are taken to represent that of materials in the dry condition .
• hydraulic conductivity and porosity data from Domenico and Schwartz (1990).
" v is the average linear ground water velocity based on an assumed gradient of 0.01 fVft (mlm) .
••• thermal property data from Hellstrom (1991). For sedimentary rocks, Hellstrom lists only Cs. In these cases, a density of 2500 kglm3 is assumed.
The porosity of soils and rocks can al 0 be an importan controlling'nflu n on
hydraulic conductivity (Freeze and Cherry, 1979). Materials with higher poro ity
generally also have higher hydraulic conductivity. However, thi correlation doe not
hold for finegrained soils (see Table 21.). Porosity and hydraulic oonductivity of oil
and rocks can be increased by socalled "secondary porosity" which is attributJ d to
solution channels (i.e. in karst limestone) or to fracturing (i.e. in rock and cohe iv
soils).
2.2.4. Conduction Versus Convection in Geologic Materials
o
It has already been noted that it is the presence of advection that distingui hes the
heat transfer regime under groundwater flow conditions from that of heat conduction
alone. Some assessment of the significance of the flow can be made by considering the
order of magnitude of the advection of heat compared to conduction (diffusion).
A dimensionless parameter describing conduction versus convection is th Peelet
number (Pe). In this application, the Peelet number expresses the transport of heat by
bulk fluid motion to the heat transported by conduction. Domenico and Schwartz (1990)
define Pe for heat transport in ground water as:
Pe = ptctqL
k eff
The term L is defined as some characteristic length dependent on the situation.
(28)
According to Bear (1972), L can be chosen as any length dimension, so long as it is
consistent with other comparisons. In principle, advection becomes significant when Pe
3
is of order one. The exact value of Pe at which advection becom significant i lightly
dependent on the choice of L.
The Peelet number has been used to quantify the relative importance of
mechanical (or advective) dispersion versus molecular diffusion in mas tran port in
ground water. Many studies have been conducted and the data have been ummarized by
Bear (1972). In short, when the characteristic length was chosen as mean grain ize,
diffusion is the process controlling mass transport at Peelet number Ie s than about 0.4.
At Peelet numbers in the range of 0.4 to 5, a transition occurs where mechanical
dispersion (or advective dispersion) and diffusion are of the same order of magnitude.
Above a Peelet number of about 5, mechanical dispersion (or advective dispersion)
dominates. No similar studies conducted for heat transport have been found.
An analysis of the Peelet number using the typical hydraulic and thennal values
of soils and rocks presented in Table 21 may be used to assess the role of ground water
flow in the design of elosedloop ground heat exchangers. The characteristic length
could conceivably be chosen as (I) a typical borehole spacing or (2) the length of the
borehole field in the direction of flow. The calculated Peelet numbers are listed in Table
22 using a typical borehole spacing of 14.8 ft (4.5 m) and assuming the fluid property
values of Pt, ct , and k( as 62.4lb/ft3 (1000 kglm\ 1.0 Btu/lboF (4180 J/kgOC), and
0.347 Btu/hrftoF (0.60 W/moC).
2
A review of the data presented in Table 22 reveals that heat advection by ground
water flow is significant process contributing to heat transfer in coar egrained oil
(sands and gravels) and in rocks exhibiting secondary porositie (fracturing and olution
channels). When the characteristic length is defined as the borehole pacing, Peelet
numbers exceeding 1 exist only for sands, gravels, and karst limestones. It i pos ible
however, that even when the Peelet number is of order one or higher, the effect of the
groundwater flow on the temperature response may not be seen within the normal time
scale of an insitu thennal conductivity test. This is one of the reasons for conducting
numerical borehole field simulations for the duration of several year .
TABLE 22.
Peclet Numbers Corresponding to Typical Values of Hydraulic and Thermal
Pr.operties of Soils and Rocks
Porous Medium Peclet Number
where L =a typical
borehole spacing of
14.8 ft (4.5 m)
[]
Soils \
Gravel 5.72E+02
Sand (coarse) 1.34E+01
Sand (fine) 1.15E+OO
Silt 1.28E02
Clay 3.24E05
Rocks
Limestone, Dolomite 5.92E03
Karst Limestone 5.28E+OO
Sandstone I 1.77E03
Shale 1.05E06
Fractured Igneous 6.32E02
and Metamorphic
Unfractured Igneous 1.00E07
and Metamorphic
3
2.2.5. Numerical GroundWater Flow and Heat Transport Models
In assessing the effects of groundwater flow on Vtube heat exchanger
performance, one is mainly interested in the temperature of the heat exchange fluid.
Therefore, modeling of the Vtubes in some detail is important in this problem. The heat
exchange fluid temperature is affected by the transient building thermal loads .in addition
to the heat transfer in the porous medium around the horehole. Consequently, this
problem is characterized by an irregular model domain with timevarying boundary
conditions and is best handled by a numerical model.
Numerous commerciallyavailable and public domain numerical software codes
exist for modeling mass and/or heat transport in ground water. Of these, the following
8 were selected for a more detailed review for potential application to thi project:
• 3DFEMFAT (3Dimensional Finite Element Method Flow and Transport) by G.
Yeh, Pennsylvania State University. This code was developed to imulate ma s
transport in variablysaturated porous media. Densitydependent flow can also be
simulated.
• AQUA3D by Vatnaskil Consulting Engineers, Reykjavik, Iceland. This code is aL 0 a
threedimensional, finiteelement code. It was developed mainly for simulation of
masstransport problems but allows easy adaptation of boundary conditions to model
heat transport without densitydependent groundwater flow.
, I
• FEFLOW (Finite Element FLOW) by WASY Institute for Wate
Planning and Systems Research, Ltd., Berlin, Germany. Thi code i al 0 a thr 
dimensional, finiteelement code. It is capable of imulating both m and h at
transport in variabledensity groundwater flow systems.
• Flowpath II by Waterloo Hydrologic, Inc. (WHD, Waterloo, Ontario. Tbi code i a
twodimensional finite difference code. It was developed originally for imuLation of
groundwater flow problems only; contaminanttransport simulation capabiLitie ,
mainly in the horizontal plane, have been recently added.
• HST3D (Heat and Solute Transport in 3 Dimensions) by USGS, Denver, Colorado.
This code is a threedimensional finitedifference code. It is capable of simulating
mass and heat transport in variabledensity groundwater flow ystem . It was
developed mainly for simulating problems involving wa te injection into aqui~ r .
• MT3d6 (Mass Transport in 3 Dimensions) by S.S. Papadopulos & As ociates, Inc.,
Bethesda, Maryland. This code is also a threedimensional finitedifference code. It
solves the mass transport equation only and requires a solution to the groundwater
flow equation from another code. It was developed for simulating contaminanttransport
problems in ground water.
5
• SUTRA (SaturatedUnsaturated TRAnsport) by United Sta e Geological Sur y
(USGS), Denver, Colorado. This code is a twodimen ioncU finjt~ I m n cod . It i
capable of simulating mass or energy transport in variablysaturated variabl den ity
groundwater flow systems. It was mainly developed as a cros  ectional model for
simulating saltwater intrusion into freshwater aquifers.
• SWIFT (Sandia Waste Isolation Flow and Transport) by HSI GeoTran ,Sterling,
Virginia. This code is a threedimensional finitedifference code. It is capable of
simulating mass and heat transport in variabledensity groundwater flow sy tern in
porous or fractured media. It was developed mainly for simulating problem
involving deepwell radioactive waste injection into geologic repositories.
In the code selection process, particular attention was paid to the following item : (I)
the type of boundary conditions handled by the code, (2) the solution cherne mployed
by the code, (3) verification of the code, and (4) cost. Each of the e point i described
in more detail below.
The type of boundary conditions handled by the code was perhap the most important
consideration in the code selection process. For simulation of periodic heat extraction or
heat rejection to the ground, the selected code needed to be capable of handling timevarying,
heat flux boundary conditions. Further, since a relatively large number of timevarying
data were to be used as input, the selected code needed to be capable of reading
timevarying conditions from an external file.
6
The software codes that were reviewed for this project employ either finitediff rence
methods (FDM) of finiteelement methods (FEM) to solve the partial differential
equations describing heat/mass transport in ground water (Equation 23 and 24). The
solution scheme was considered important for two main reasons. First, FEM offer an
advantage over FDM in the ability to represent complex or irregular geometries (i.e.
circular Vtubes in a rectangular domain). Second, there has been controver y in the
literature over advantages of FEM over FDM in solving the advectiondispersion
equation (Equation 24). In general, experience has shown the FEM to be generally
superior to FDM in solution stability (Wang and Anderson (1982) and Mercer and Faust
(1986». Consequently, codes employing an FEM solution scheme were preferred.
Documented verification of the code was an important consideration since it often
requires years to find and fix bugs in these types of software programs. All of the 8
codes listed above have been originally developed in the 1980s and many validation
examples exist.
The cost of the code was also a main consideration in the selection proce since the
project had an allocated budget for software.
7
2.3. The Numerical Model
2.3.1. Model Description
The computer code selected for this study was AQUA3D. It i a commerciallyavailable
software package originally developed in 1983 by Vatnaskil Consulting
Engineers of Reykjavik, Iceland. The partial differential Equations 23 and 24 are
discretized spatially by a Galerkin finite clement method using triangular element with
linear weighting functions (Vatnaskil, 1998). The temporal term of the equation i dealt
with by first order backward differencing in time. AQUA3D does not allow for the
explicit representation of the heat transport equation, but provides a general form of the
mass transport equation (EquatioD 24). Temperatures were in fact calculated by suitable
choice of the coefficients of the mass transport equation and corresponding adaptation of
the boundary conditions.
The finite element groundwater flow and mass/heat transport model wa used in
this study as the primary means of assessing the effects of groundwater flow on c1osedloop
heat exchangers. Use of a numerical model allows a wide range of conditions to be
examined and is the only practical means of modeling a whole borehole field. In each
test case, a unidirectional flow field was imposed over the whole numerical domain. As
the flow was assumed to be fuJlysaturated and within homogeneous geologic material, it
was only necessary to use a twodimensional model.
3
2.3.2. The Finite Element Mesh
Finite element meshes for a single borehole geometry, and for a complete
borehole field geometry have been constructed using triangular element . Nodal pacing
was kept relatively fine around the pipe walls where the steepe t temperature gradient
were expected. The mesh for the single borehole geometry was constructed within a
square domain and consisted of 465 nodes, as shown in Figure 22.
I I
I I
I I
~ 3.46 in. (8.8 em) +!
,<III 14.4 It (4.4 m) ~,
Figure 22. Finite element mesh representing a single borehole.
A mesh for a fourbyfour configuration borehole field w con tructed by u iug
the single borehole mesh (Figure 22) as the basis for the me h at each borehole and by
expanding the mesh in the direction of groundwaterflow a shown in Figure 23. This
mesh consisted of 4532 nodes.
1.;.4 1640 ft(SOOm) .t
.~..  57.71t117.8m) .t~!
Figure 23. Finite element mesh representing a 16 borehole field.
2.3.3. Boundary Conditions
Two sets of boundary conditions are required: one set for the flow model and one
set for the transport model.
o
In the flow model, firsttype or fixedvalue (Dirichlet) boundary condit,ion were
set on the lefthand and righthand boundaries in order to impo e a fiJ(l d hydraulic
gradient across the domain. Secondtype or fixedgradient (Neumann) boundary
conditions were set on the upper and lower boundaries of the flow domain and were
specified as zero flux. This work assumes that no groundwater recharge take place
across the water table within the model domain. In the transport model, Dirichlet
boundary conditions were set on all four sides of the model domain. The e conditions
represent fixed background or farfield temperatures.
In order to impose the ground thermal loads as boundary condi.tion at the Vtube
pipe walls, some adaptation of the usual boundary condition was required. Thi ari es
from the use of the mass transport equation to model heat transport. First, a zero flux
condition for the mass (heat) transport equation was applied at each of the sixteen nodes
forming each pipe wall. The required heat flux is imposed using a source term in the
groundwater flow equation at these nodes (representing injection of the heating/cooling
water). The flow injected, V·, was set negligibly small (3.53 x 10 19 fe/s [1.0 x 10.20
m3/s]) so as not to disturb the groundwater flow field. The temperature of this injection
flow, Tw was set to achieve the required heat input (the ground thermal load), 0 that,
T q* :=::;=
W PlCl v* (29)
The values of ptand ctare taken as constants of 62.4 Jb/ft3 (1000 kg/m3
) and 1.0 Btu/lb
OP (4180 J/kg°C). The average temperature of the heat exchange fluid in each borehole
is taken as the average of the nodal temperatures of the 32 nodes defining the Vtube pipe
41
in each borehole. Where single borehole cases were simulated, the heat input per pipe
node, q., was set at a fixed value representative of insitu thermal conductivity te t
conditions. Where the whole borehole field was modeled, q. was determined from the
timevarying building loads.
2.3.4. Validation of the Numerical Model
In order to check the accuracy of the AQUA3D model and the implementation of
the boundary conditions, an appropriate analytical solution was sought. Numerous
analytical solutions have been developed for the advectiondispersion equation (Equation
24). However, these are mostly specific to pollutanttransport problems (T i replaced
by solute concentration in Equation 24) involving point or line sources with uniform
concentration in time. Fetter (1988) and Bear (1972) summarize solutions for boundary
and initial conditions describing situations that are commonly found in nature. The
literature gives little to no attention to analytical solutions describing the explicit
transport of heat in ground water.
The most appropriate analytical solution found was that described by E kil on
(1987) for steadystate heat extraction from a borehole in a groundwater flow field.
However, Eskilson's solution contains some approximations.
Eskilson (1987) describes the steadystate temperature at the borehole wall (Tsw,b)
as:
T T +q' {1[ In(H Jp (H )~} (210)
.!W.b  om H 2n:k
s
2r
b
IV £
The logarithmic term gives the steadystate resistance and Pw ( ~ }s a correction term
accounting for groundwater flow which is approximated as:
42
error between the two solutions is 4.9%. This error was considered to be acceptable.
Comparisons of the steadystate temperature at the wall of a pipe were made
(211)
Q~m
0
(212)
tll
where £»rb
using both the AQUA3D and the Eskilson (1987) solution. The following data were u ed
P (.H) =InH + r11n(3)
IV £ U 2
for f. defined as:
depth of 250 ft (76.2 m), soil thermal property values of sand from Table 21, water
properties as described above, a farfield temperature of 63°F (17.2°C), and a heat flux of
as inputs: a pipe diameter of 0.787 in. (2 cm) (i.e. a single leg of a typical Vtube), a pipe
Eskilson's (1987) solution were 101.4°F (38.54°C) and 99.59~ (37.55°C), respectively.
8530 Btulhr (2500 W). The temperature predicted by AQVA3D at steady tate time and
In tenns of temperature increase above the farfield temperature, the percent difference in

4'
2.4. Results And Discussion
2.4.1. SingleBorehole Simulations
The numerical model was initially used to determine average borehole
temperatures for a range of soil and rock types over a twoyear simulation time. The
objective was to examine trends in heat exchanger performance with increa ing Peelet
number. The hydraulic and thermal property inputs are those listed in Table 21 and a
hydraulic gradient of 0.01 was assumed. A constant heat flux of 8530 Btu/hr (2500 W)
was applied on a Vtube in a 250 ft (76.2 m) deep borehole. The initial temperature and
firsttype boundary conditions were set at 63°F (17.2°C). The model domain is that
shown in Figure 22. A simulation time of two years with a time step of 5 days was used
for these simulations. For comparison purposes, simulations were made for each case but
with zero groundwater flow.
Plots of average borehole fluid temperature versus time for three example
geologic materials are shown in Figure 24. A review of these plots reveals that a
"typical" ground water flow rate in a coarse sand dramatically lowers the average
borehole fluid temperature when compared to the zeroflow case. After a oneyear
period, the average fluid temperature in the borehole is approximately 15°F (8.3°C) lower
than the average fluid temperature in the borehole where no groundwater flow was
simulated, and appears to have reached a steady state. A small reduction in peak
temperature is shown for the case of fine sand. However, "typical" ground waterflow
44
(a) coarse sand
120 rrr,rr,.,......r 48.9
110 ttC7'"4F'+tttIIt 43.3
15,6
0.50 0,75 1.00 1,25 1.50 1.75 2.00
Time (yeers)
0,25
~ ~ 3 ~. ~
~ ~>oo f'jIIf++fIrf 37,8 ~ E
o f ~ ! 'i a • 3
~! 90 32,2 ~ =
&~ ~~
!!. 80 .111 ~::~.."N~o=Grou'nd.,.,W.,cat"er·cF=loW 26.7 !! ~
~ ~ ~ ~
c( __Ground Water Row Ra1e 0<
70 tII =196.8 ftJyr (60 m/yr) 111 21.1
60 I+IT+r+4_T~'"____l
0.00
(b) fine sand
120 48,9
15.6
1,75 2.00
~
'5
37.8 ~ E
o .,
J:. ~ e :>
32.2 0 e III ., &e I+f 26. 7 ~ ~
~
0,50 0,75 1.00 1,25 1.50
Time (yeara)
Ground Water Row Rate
70·t+j = 16,6 ftlyr (5.05 m/yr) 111 21.1
1 r 1
601lI+I~++l
0.00 0.25
110 ++:7""""""fIf++ft 43,3
~ .tIi:" 100 .hF/II+++ft+I o • J:. ~
~ ~ 90 ItI+++tf+i
III •
e& .E 80 1t _+f':":·:···:··:·N:o~~Grrooou:n1dd\WNa8tie~r=Filmo.w
~ ~
0<
(c) shale
:2_
37,8 ~ E
.! !
o a
32,2 'i e
~ II
Ilil a.
. E
26.7 ~ ~
0<
48.9
43,3
21.1
15.6
0.75 1.00 1.25 1.50 1,75 2.00
Time (yeara)
 .....• ' No Ground Waler Row I
Ground Water Flow Rate = I
I _rE
4 ~(8'~1'5 miT I
60
0.00 0.25 0.50
120
110
Figure 24. Average borehole fluid temperature vs. time for (a) coarse sand,
(b) fine sand, and (c) shale.
45
rates in rocks such as granites, limestones, dolomite , and shales were found to have a
negligible effect on the average borehole fluid temperature.
The trends shown in these results are in agreement with the previou Peelet
number analysis. At Peelet numbers of order one or higher, advection of heat by flowing
ground water is a significant process contributing to heat transfer in the ground. At
Peelet numbers of order less than one, conduction is the dominant heat transfer proce s
and enhancement to the heat exchanger performance is negligible.
2.4.2. Simulated InSitu Thermal Conductivity Tests
The second objective of the singleborehole simulations was to determine the
effects of groundwater flow (in a material where groundwater flow is expected to be
significant) on the interpretation of data from insitu ground thermal conductivity t t.
The previous results showed the effects of groundwater flow to be most significant in the
cases of gravel and eoarse sand. Accordingly, the simulated insitu thermal conductivity
test calculations have been based on coarse sand properties.
As previously discussed, insitu thermal conductivity tests involve the application
of a steady heat flux to a test borehole along with the measurement of the temperature
response of the circulating water. These data are used either with an analytical model or
with a numerical model and parameter estimation technique to arrive at a value of
thermal conductivity of the soil! rock formation. Here, the borehole temperature response
46
is calculated for a range of groundwater flows using the AQUA3D model. The da
have been analyzed in exactly the same way as if the data had been measured in itu.
Hence, 'effective' thermal conductivities have been estimated for different flow
conditions.
Insitu test conditions were modeled by applying a constant heat flux of 8530
Btulhr (2500 W) on a Utube in a 250 ft (76.2 m) deep borehole. The simulation time
periods were 50 hours and one week, corresponding to typical duration of in itu ground
thermal conductivity tests. The model time step was 2.5 minutes. The initial temperature
and firsttype boundary conditions were set at 63°F (17.2°C) and the model domain i that
shown in Figure 22. Model input hydraulic and thermal property values are those listed
in Table 21 for a coarse sand, except the ground water flow velocity was varied from a
"typical" value of 196.8 fttyr (60 m/yr) to a more extreme value of J968.5 fttyr (600
m/yr) by adjusting the hydraulic conductivity value. Twelve case wer simulat d a
listed in Table 23.
Resulting temperature responses for the 12 cases are plotted in Figure 25. A
review of Figure 25 shows that ground water flow in a coarse sand significantly impacts
the average borehole fluid temperature over the time scales of an insitu ground thermal
conductivity test. Two noteworthy conclusions can be drawn from these simulation : (1)
as ground water velocity increases, the time to reach steadystate conditions decrease
and (2) as ground water velocity increases, the steadystate temperature decrease . Also,
the deviation from the zeroflow condition can be seen to be dependent on the duration of
47
the test; temperatures are further reduced with increasing duration. Hence, the duration
of the test could be expected to have an influence on the estimated thermal conductivity
derived from an insitu test.
TABLE 23
Case Summary of Simulated InSitu Ground Thermal Conductivity Tests
Case Simulation Time Period Ground Water Flow Velocity
1 50 hours No Ground Water Flow
2 50 hours 196.8 ft/yr (60 m1yr)
3 50 hours 393.7 ft/yr (120 m1yr)
4 50 hours 787.4 ft/yr (240 m1yr)
5 50 hours 1574.8 ftlyr (480 rn/yr)
6 50 hours 1968.5 ft/yr (600 rn/yr)
7 1 week No Ground Water Flow
8 1 week 196.8 ft/yr (60 rn/yr)
9 1 week 393.7 ft/yr (120 rn/yr)
10 I week 787.4 ft/yr (240 m1yr)
] 1 ] week ]574.8 ft/yr (480 rn/yr)
]2 ] week 1968.5 ft/yr (600 rn/yr)
The effective thermal conductivity values predicted by the Austin et al. (2000)
model are plotted against the corresponding groundwater flow velocity for each of the
two insitu test simulation times (50 hours and 1 week) in Figure 26. The actual values
are listed by case number in Table 24. A review of these results shows that as ground
water flow velocity increases, the predicted effective thermal conductivity value from a
conductionbased model are significantly different, depending on the duration of the
simulated test. These values are "effective" values since they include the effects of
ground water advection. However, at this stage of the design process, it is not clear if the
I
50hour data set or the Iweek data set produces more representative value , or if either
data set produces representative values at all.
(a)
Case 2
Case 1_
~
...,0.; ~'" "
~F Case 4 ./ Case 5
, ~ "
I Case 6 '"

Case '
~ase 8 r  ./ ~ Case 9
~~
 
,R' Case 0/ Cse11,
r ./   1 Ca~ 12 "
f
 I
60 80 100 120 140 160
Time (hrs)
35.0
32.2
15.6
32.2
G' QI 0
29.4 "0 ;
.r:. ..
QI ::s
26.7 0 ~
~ 8
23.9 C'l E
l!! III
III ~
21.1 ~ ~
u:: 18.3
29.4 G'
tP't..,.
'0 GI 26.7 .c :s
! 1i
23.9 ~ ~
GI E
C'l GI
21.1 ~ ~
> 'tI
c( 'S
18.3 u:
15.6
40 50
(b)
20 30
Time (hrs)
40
10
20
90
60
o
95
ii:'
~ ~ 85 .oc .C.P.
GI ::s o iii 80
lEI i
& E 75
to GI
~ ~
~ ~ 70
it
65
60
o
90
ii:' 85
GI~
:g ~ 80
of l!!
lD GI 75 CD Q.
C'l E
~ ~ 70
> 'tI
c( oS
u: 65
Figure 25. Average borehole fluid temperatures for the 12 simulated insitu
ground thermal conductivity test cases in a coarse sand with ground
water velocities ranging from 0 to 1968 ft/yr (600 m/yr) for (a) 50
hours and (b) one week.
/
4
In order to investigate the effects of groundwater flow on borehole field
performance and system design procedures further, the predicted ground thermal
Figure 26. Predicted effective ground thermal conductivity values versus
groundwater flow velocity for 50hour and Iweek imulated insitu
thermal conductivity tests.
24.2 0
0
20.8 E
~
17.3 CI) ::s
13.8 a;
>
..Ill:
10.4 'C
CI) 6.9 u
'C
3.5
.C.I.) D.
600
27.7
0.0
656.2 984.2 1312.3 1640.4 1968.5
Ground Water Velocity (ftlyr)
Ground Water Velocity (mlyr)
100 200 300 400 500
328.1
/
/ f .50 hour simulated insitu test / f ___ 1 week simulated insitu test
 /
./
/'
/ /
.,/' v
~.......
8
6
4
2
o
0.0
CI) ::s
~
~
"CI) .~
".C.I.) c.
o
_ 16
LL.
0, 14
=..,=.
~ 12 m 10
2.4.3. Borehole Field Simulations
conductivities were used to design a borehole field for a test building. The test building
was an actual building located in northcentral Oklahoma. This building i a ingle tory
office building with 8 thermal zones and has a predominant demand for cooling. The
hourly building loads were determined for one year using building energy simulation
software (BLAST, 1986). The building loads were then converted to ground loads under
the assumption that all heat pumps in the system have a constant coefficient of
performance of 4.0. The ground loads for this building are shown in Figure 27.
TABLE 24
Summary of Borehole Field Design Parameters for Each Test Case
Case Simulation GroundWater Ground Thermal Design Borehole
Number Duration Flow Rate Conductivity Depth
Predicted by Predicted by
Numerical Model of Design Software of
Austin et al. (2000) Spitler et al. (1996)
(hours) ftlyr BtulhrftOF ft
(m/yr) (W/mOC) (m)
1 50 0 0.643 239.98
(1.11 ) (73.15)
2 50 196.85 0.650 238.56
(60.00) (1.12) (72.71 )
3 50 393.70 0.731 224.10
(120.00) (1.26) (68.31)
4 50 787.40 1.146 171.56
(240.00) (1.98) (52.29)
5 50 1574.80 3.657 87.24
(480.00) (6.33) (26.59)
6 50 1968.50 6.074 61.58
(600.00) (10.51) (18.77)
7 168 0 0.625 243.86
(1.08) (74.33)
8 168 196.85 0.691 230.86
(60.00) (1.20) (70.37)
9 168 393.70 0.962 191.58
(120.00) (1.66) (58.39)
10 168 787.40 2.250 115.91
(240.00) (3.89) I (35.33)
11 168 1574.80 8.229 , 4B.02
(480.00) (14.24) (14.64)
12 168 1968.50 15.107 26.90
(600.00) (26.14) (8.20)
5
I
51
250 59
200 44
~.s::. 150 29
":i ~
 100  ! 15 III
III 50 'D
'D III
III 0 .0.J 0 0 ..J 'D
'D  15 §
I:: 50 e :I
0 29 Cl (; 100
150 ~4 44
200 f...,,' 59
o 2190 4380 6570 8760
Time (hours)
Figure 27. Hourly ground loads for the test building. Heating load is shown
negative, representing heat extracted from the ground; cooling load
is shown positive, representing heat rejected to the ground.
Borehole field designs were produced for each of the twelve test cases. This was
done with GLHEPRO, a commerciallyavailable groundloop heat exchanger design
software tool developed by Spitler et al. (1996). A 16 borehole field (fourby four
boreholes in a square pattern) was deemed adequate for the test building ground load
(Figure 27). The monthly loads and peak hourly loads are generally input in the design
software. For this study, no peak hourly loads were specified for the sake of the
computational time required for the subsequent borehole field simulations (see discu sion
below). Peak design entering fluid temperatures to the heat pump were specified at 90°F
(32.2°C) maximum and 3SoF (1.7°C) minimum. The borehole depths were sized for 20
years of operation.
5
For each test case, the corresponding effective thermal conductivity hown in
Figure 26 and Table 24 was input into the groundloop heat exchanger de ign oftware
(GLlIEPRO). The borehole depths predicted by GLHEPRO are plotted again t the
corresponding ground water flow velocity for each of the two insitu te t imulation times
(50 hours and I week) in Figure 28. The values are listed by case number in Table 24.
Figure 28. Design borehole depths versus groundwater flow velocity for 50hour
and Iweek simulated insitu thermal conductivity test.
o 0.0
0.0 328.1 656.2 984.3 1312.3 1640.4 1968.5
Ground Water Velocity (ftlyr)
E
30.5
61.0 .c:: Co
CIl
C
.s!
o
s:
~
o
ED
c
Cl
'iii
CIl c
15.2
Ground Water Velocity (mlyr)
o 100 200 300 400 500 600
250 t+ttttj 76.2
g
;: 200 +l'c''''<:t1 _ 1 week simulated insitu test
Co
CIl
C
.s! 150 ++~t~..,____ttj45.7
o
.c::
~ 100 ++t""'d:+~,.__+l
ED
c
Cl "iii 50 ++I++~a:+1
CIl
C
AQUA3D was further used to simulate the longterm performance of each
borehole field designed from the simulated insitu ground thermal conductivity te teases.
The model domain was that previously described for the multiborehole field simulations
shown in Figure 23. The simulation time for all cases was 10 year and the model time
step was 5 days. The simulated heat flux at the internal boundary node defining the Utube
pipes was a timevarying source corresponding to the monthly ground loads for the

s
test building. Hydraulic and thermal property input for each borehole i ld c number
is the same as the corresponding singleborehole case number, except for tbe borebol
depths, which are those listed in Table 24. Each lOyear simulation required
approximately 60 hours of computation time on a personal computer with a 233
megahertz pentium II processor.
Annual maximum and minimum peak temperatures are plotted for each case in
Figure 29. Examination of the cases with no groundwater flow (cases I and 7) how
annual rises in peak temperature typical of coolingdominated buildings. After the
second year, all of the cases with groundwater flow show no changes in minimum and
maximum temperatures from year to year.
Some notable differences can be seen between the borehole field designs based on
50hour test data compared to one week test data. Thi is shown by ca 'e 5 and 6 which
represent thermal conductivity values determined from a 50hour te t at ground water
flow velocities of 1574.8 ft/yr (480 m/yr) and 1968.5 ft/yr (600 m1yr), respectively, and
by cases 11 and 12 which are for the same flow rates but based on thermal conductivity
values determined from oneweek test data. The thermal conductivity values determined
in cases 11 and 12 are unrealistically high and consequently the design borehole depth
are too shallow; the result is that the maximum peak temperature of the simulated
borehole field in both cases exceeds the maximum design temperature during the fir t
year. This implies that for insitu test cases where the average borehole fluid temperature
reaches steadystate in a relatively short time (as demonstrated by case 4/10, case 5/11,
4
and case 6/12 in Figure 25), increasing the duration of the in itu te t produce d d
confidence in the accuracy of the effective thermal conductivity value determined from
the test.
(a) (b)
TIm. (y.....) Tim. (yosr.)
25 >+ I _.J ·3.9
o 9 W
'~
''..
''..
'...... :•''..
18.3 . 12.8 :;
lii ;;
D.
E
7.2 ~
E
"J§
__ Mni'rumOesIgnT'''''''''alure 1.7 i
43.3 55
~ ~ · · 37.8 Z :; 55 • lii
;; ;;
D. D.
E E
32.2 :!. l· 45
E E
E" J"§
26.7 :I £ 35
:II :Ii
II I,I
! . t
I !. I i I I ! i
00 3J;1;f1mt
80 I; !~I I ,; I;' I
70 i1 j! I 21.'
o 9 10
110
•~ 100
•D.
E
~
E
"E
:=II
Icasal +casa::........ Cas. 3 ....cas.4 __eas. 5__Cas. aI Feasel +case2:"cas.3 ':"'C;se4 __Cas. 5 eas;;8j
''..
>fWUUNUuUU)(
twBxlrrum Oeslgn T.",,",elur. _
~~::
• • • a A a A • a •
110
t
•~ tOO
.;;.E
:!. 90
E
"E
i 80
:II
o
(e)
.. .. 43.3 55
~ f
• •
37.8 ~ :; 55
10
OJ ;;
D. D.
E E
322 ~ l· 45
E E
E" "E
·26.7 i "2 35·
:II i
(d)
___________ 18.3
~~12.8;
.,. !
l:.
/
1 ,. II Ie Ie I( II Ie l( E
7.2 :!.
E JI. · • • · · · .. .§ ...,Irrum OIlSIgn Terrpe,alur.            1.7 i
..
:'
'f
70 >___~_.l21.1
o 3 10
Time (y••n)
25 l___ .J ·3.9
o 3 6 8 9 10
TIm. (y..,.)
Icase 7+Cas.8 ~eas.9 ....~"'0 __ease 11 __Cas. 12.] t:Cas. 7+ea;;; 8....... eas.9 ....cas.l0__ease 11 __eas.,2 ,
Figure 29. Annual maximum (a) and (c) and minimum (b) and (d) average
borehole fluid temperatures for the 16 borehole field simulations.
Except for cases 11 and 12, the annual maximum and minimum temperatur fell
within the design conditions. Having followed conventional de ign procedur
interesting to note from Figure 29 (a) that it is the cases where the groundwater flow i
moderate (cases 2, 3, and 4) that are most overdesigned. These cases have maximum
peak temperatures of about 74°F (23.3°C), some 16°F (8.9°C) below the maximum design
temperature. Considerable drilling cost savings could be seen in cases like this where
shallower borehole depths could have been adequate. It is at higher flow rates (cases 5
and 6) that the peak temperature was closest to the original design condition after ten
years. This illustrates the nonlinearHy introduced into the design problem by the
presence of advection. It also illustrates the difficulty in adapting conventional design
methods to accurately size closedloop ground heat exchangers in ca es with ignificant
groundwater flow.
The temperature field predicted by the numerical model for case 8 is hown in
Figure 210 in the fonn of a series of contour plots over the to year imulation period.
The data are plotted for the end of September, when the average annual ground
temperatures are the greatest (i.e. at the end of the cooling season).
, ,
ground water flow rate = 60 m/yr (197 ft/yr)
•
end of Sept. Year 1
17,500 _0
18,(0)
18,500 c:J
19,(0)
19,500
20,000
20,500
21,000
21,500
22,(0)  22,500
'Iemperature ey
\C)
t
100m
t
< 500m
end of Sept. Year 5
end of Sept. Year 10
I:_~
Figure 210. Temperature contours for Case 8 for the end of September ofyears 1, 2, 5,
and 10.
57
A further feature that is shown in the predicted temperature field (Figure 2l 0) i
the development of a peak in the temperatures immediately down tr am of the bor hole
field after year 2. This arises from the advection of the heat rejected to the ground at the
boreholes during the previous year. In the contours plotted for year 10, thermal plumes
from previous years can be distinguished.
2.5. Concluding Remarks and Recommendations for Future Work
Using a compilation of "typical" hydraulic and thermal properties of oil and
rocks, a preliminary analysis of the effects of groundwater flow on the design and
performance of closedloop groundcoupled heat pump ystem has been made. A
simple but useful method of assessing the relative importance of heat conduction in the
ground versus heat advection by moving ground water is demonstrated through the use of
the dimensionless Peclet number.
A finiteelement numerical groundwater flow and heat tran port model wa u d
to simulate and observe the effects of groundwater flow on the heat transfer from a
single Utube closedloop ground heat exchanger in various geologic material . From
these simulations and from a Peelet number analysis, it appears tbat it i only in geologic
material with high hydraulic conductivities, such as coarsegrained soil (sand and
gravels) and in rocks exhibiting secondary porositie such fractures and olution
channels, that groundwater flow could be expected to have a significant effect on c1osedloop
heat exchanger performance.
:'J ',.
........ .'. '.
......
':~
5
The effect of groundwater flow on in itu thermal conductivity te t re ult h
been examined by numerically simulating test conditions around a ingle bOI bol under
different flow conditions. These data were analyzed as if it came from real in itu
sources to arrive at effective thennal conductivity values. As expected in all ca of
groundwater flow, these values were artificially high. Results from one week te t data
have been shown to be no more reliable than data from 50hour test .
The finiteelement numerical ground water flow and heat transport model was
also used to simulate the 10year performance of borehole fields de igned from
application of conventional design procedures using the derived thennal conductivity
data. For coarsegrained sands, the presence of moderate groundwater flows had the
effect of removing the yearbyyear increase in ground temperature found in systems
where there is no groundwater flow. The borehole fields de igned using conv ntional
methods and the derived effective thennal conductivities were generally overde igned.
However, in some cases at very high groundwater flow rate , temperatures were found
to rise above the design criteria.
From this preliminary assessment of the effects of groundwater flow, it appears
difficult to adapt results from current design and insitu measurement method to fully
account for groundwater flow conditions. Over the last decade, considerable progress
has been made in developing both insitu test methods and design procedure for
borehole field design for situations where there is no groundwater flow. Research would
3. A Model for Simulating the Performance of a Shallow Pond as a
Supplemental Heat Rejecter with ClosedLoop GroundSource Heat
Pump Systems
3.1. Introduction
Commercial buildings and institutions are generally coolingdominated, and
therefore reject more heat than they extract over the annual cycle. In order to adequately
dissi ate the imbalanced annual loads, the required groundloop heat exchanger length
are significantly greater than the required length if the annual loads were balanced.
Consequently, under these c' cumstances, groundsource heat pump sy terns may be
eliminate~fr~~onside..!,!tion during the feasibility study phase of the HVAC de ign
process because of excessive first cost.
To effectively balance the ground loads and reduce the nece 'sary size of the
ground loop heat exchanger, supplemental components can be integrated into the groundloop
heat exchanger design. GSHP systems that incorporate some type of upplemental
heat rejecter are commonly referred to as hybrid GSHP system. In ome application,
the excess heat that would otherwise build up in the ground may be used for dome tic hot
   ~. 
water heaters, car washes, and pavement heating systems. In cases where the exce s heat
cannot be used ~!lJili~iall ,shallow ponds can provide a costeffective_mea~ to balance
the thermal loading to the ground and reduce heat exchanger length.
The objective of the work presented in this chapter has been to develop a design
and simulation tool for modeling the performance of a shallow pond that can be useful1y
6,1
and costeffectively integrated into a groundsource heat pump y tern a a upplemental
heat rejecter. The pond model has been developed in the TRNSYS modeling
environment (SEL, 1997) and can be cou led to other GSHP sy tern component model
for shorttime step (hourly or minutely) system analy es. The model ha b n validated
by c.o..mparing simulation results to experimen_~a! data. As an example of the model'  ._ ~.  .~ " ..
applicability, GSHP system simulation results are presented for a commercial building
located in Tulsa, Oklahoma with a h othetical closedloop GSHP y tern with and
without a shallow pond supplemental heat rejecter.
3.2. Heat Transfer In Ponds
3.2.1. General Overview
Pertinent concepts of heat transfer in ponds and lakes have been summarized by 
many sources. Dake and Harleman (1969) conducted studie of thermal tratification in
lakes and addressed the overall thermal energy distribution in lakes. ASHRAE ( 1995),
ASHRAE (1995b), and Kavanaugh and Rafferty (1997) describe heat tran fer in lakes in
relation to their use as heat sources and sinks.
Solar energy is identified as the ~~~ing~for ponds and lakes.
The ~n f.9.~lin~mechan!§!!l is evaporation. Thermal radiation can also account for a
significant amount of cooling during night hours. Convective heating or cooling to the
atmosphere is less significant. Natural convection of water due to buoy~c:y effects. i the
primary mechanism for heat transfer within a surface water body. Conduction heat r'
transfer to the ground is generally a relatively insignificant proce ,except in c
the water surface is frozen.
Shallow ponds are generall thermally unst~at~i~d. Natural stratification of
deeper ponds and lakes is due to buoyancy forces and to the fact that water i at it
62
where
greatest density at 39.2°F (4°C). Therefore, over the annual cycle, water in deeper ponds
will completely overtum. Thermal stratification in ponds is also dictated by inflow and
outflow rates or ground water seepage rates. If inflow and outflow rates are high enough,
the pond will not stratify. Consequently, thennal stratification occur only in pond and
lakes that are relatively deep, generally greater than 20  30 ft (6.1  9.1 m), with low
inflow rates. The relevant heat transfer mechanisms occurring within shallow ponds are
illustrated in Figure 31.
3.2.2. Existing Pond and Lake Models
Several mathematical and computer models have been developed for simulation
of lakes used as heat sinks/sources and for solar ponds.
Raphael (1962) develo ed a numerical model for determining the temperature of
surface water bodies as heat sinks for power plants. Thermal stratification of the water
body was not considered. Input data to the model included weather data and inflow and
outflow data for the water body. Raphael reported that the model successfully predicted  
the temperature changes in a ri ver used as a heat sink for a power plant.
6
Figure 31. Heat transfer mechanisms in shallow pond.
Thermal RadiatiCl1
Evaporation
~ ~ ~
~~miilltl~~
Convection
Wind
Jobson (1973) develo ed a mathematical model for water bodie u ed a heat 
sinks for power plants. Thermal stratification of the water body Wa not con id r d. The
results of that work showed that the heat transfer at the water/air interface i highly
.. 
dependent on the natural water temperature and the wind speed..
Cantrell and Wepfer (1984) devel ed a numerical model for evaluating the
potential of shallow ponds for dissipating heat from buildings. The model takes weather
data and building cooling load data as inputs and computes the steady tate pond
temperature using an energy balance method. Thermal stratification of the pond was not
considered. The model showed that a 3acre (12,141 m2
), IOfeet (3.048 m) deep pond in
64
Cleveland, Ohio could reject 1000 ton (3516 kW) of thermal energy wjth a maximum
increase in pond temperature of about SOp (2.78°C) over a daily cycle.
Rubin et al. (1984) develo ed a model for solar ponds. The purpo e of a olar
pond is to concentrate heat energy from the sun at the pond bottom. Thi i accomplished
..... . 
by suppressing natural convection within the pond induced by bottom heating, u ually by
adding a brine layer at the pond bottom. As a res~lt, solar ponds have three distinct zones
as described by Newell (1984): (1) a top layer which is stagnated by some method and
acts as a transparent layer of insulation, (2) a middle layer which is usually allowed to be
mixed by natural convection, and (3) a lower layer where solar energy is collected. The
model of Rubin et al. (1984) applied an implicit finite difference scheme to solve a onedimensional
heat balance equation on a solar pond. Largescale convecti.ve currents in
the pond were assumed to be negligible while smallscale convective currents were
handled by allowing the coefficient of heat diffusion to vary through the pond depth.
Solar radiation was modeled as an exponentially decaying function through the pond
depth. The model successfully predicted seasonal variation in solar pond temperatures.
Srinivasan and Guha (1987) ~p.e.d a model similar to the model of Rubin et
al. (1984) for solar ponds. The Srini vasan and Guha (I 987) model con isted of three
coupled differential equations, each describing a thermal zone within the solar pond.
Solar radiation in each zone is computed as a function of depth. The model ~I~o_
~ucc~~sful.1¥PFediGted seasonal variations in solar pond temperatures with various heat
extraction rates.
65
Pezent and Kavanaugh (1990)rdeveloped a model for lake u ed a heat our
or sinks with watersource heat pumps. The model es entially combined the model of
Srinivasan and Guha (1987) to handle stratified cases and of Raphael (1962) to handle
unstratified cases. As such, thermal stratification of a lake could be handled in the 
summer months, when lakes are generally most stratified, and neglected in the winter
months, when lakes are generally unstratified. The model i driven by monthly av rage
bin weather data and handles both heat extraction and heat rejection. With no heat
extraction or rejection, the model favorabl predicted a lake temperature profile in
Alabama. The temperatures within the upper zone of the lake (the epiUmnion) and the
lower zone of the lake (the hypolimnion) were predicted to within 4<>P (2.22°C) and
approximately lOp (0.55°C), respectively. However, the model had orne difficulty in
matching the intermediate zone (the thermocline), perhaps due to the fact that thi zone
possesses moving boundaries (unlike the boundaries of a solar pond which are more
distinct). As concluded by Pezent and Kavanaugh (1990), a numerical method i
necessary to more accurately predict the thermocline rof .
The model presented in this paper is based on the as urnption that thermal
gradients in shallow ponds are negligible, especially during times of heat rejection. Thi
model is developed in the TRNSYS modeling environment and can be coupled to other
component models for larger system simulations. Purther, this model allow the pond
performance to be simulated on hourly or minutely time scales.
3.3. Experimental Methods : .. _.  
3.3.1. Pond Description and Data Collection
~on~ construction and data collection for this study have ~~~. conduc~~~
researc~ers affiliated with _~~ Division of En ineering Technology at the Oklahoma State
University. 
Two ponds were constructed on a test site at the north end of the cam u . The   _.
ponds are rectangular with a Ian area of 40 ft (12.19 m) by 3 ft (0.91 m). Each pond was
constructed with vertical sidewalls with one of the ponds being 2 feet (0.61 m) deep and
the other being 3.5 ft (1.07 m) deep. The walls and the bottom of each pond were
constructed of reinforced concrete, approximately 8 in. (20.3 cm) thick.
Heat was rejected to each pond by circulating heated water through a "slinky"
heat exchanger (a pipe coiled in a circular fas ion such that each loop overlaps the
adjacent loop) installed in each pond. Each slinky pipe was made of highdensity
polyethylene plastic and is 500 feet (152.40 m) long with a nominal diameter of % in.
(0.019 m). The pipe was coiled such that the resulting slinky heat exchanger wa 40 ft
(12.19 m) long with a diameter of 3 ft (0.91 m) and a 1Oin. (0.254 m) pitch (the
separation distance between the apex of each successive loop). In the 2ft (0.61m) deep
pond, the slinky heat exchanger was positioned horizontally within the pond at a depth of
approximately 10 in. (0.254 m). In the 3.5ft (1.07m) deep pond, the slinky heat
....
exchanger was positioned vertically within the pond along the centerline of the long axi
of the pond.
The temperature of the pond water was measured by tbennistor po itioned at
four locations witbin the pond: (1) near the pond surface at the center of tbe linky, (2)
below the slinky at its center, (3) near the pond surface at the end 0 0 ite from the
supply end, and (4) below the slinky at the end of the pond opposite from the upplyend.
Slinky supply and return water temperatures were measured by thermistor embeddeg, in
the slinky header. Each system also included a flow meter, a water heating element, and
a watt transducer. All sensor information was recorded by the data acquisition y te
time intervals of 6 minutes.
The tests were controlled to maintain a set supply water temperature by heating
the supply water if the temperature fell below a set point. Two et point temperatures
were used in this study: 90<>P (32.2°C) in the summer season and 75°F (23.9°C) in the
winter season.
3.3.2. Weather Instrumentation and Data Collection
Weather data for this study were obtained from the Oklahoma Mesonet
(mesoscale network), which is a weather station network consisting of weather
monitoring sites scattered throughout Oklahoma. Depending on the weather parameter,
)
I)
.)
".
data are recorded at time' ,=~=ingJrom3 seconds to 30 seconds and averaged
over 5minute observation intervals.
6
Weather data at ISminute intervals for the Stillwater monitoring tation w re
acquired for the time periods of interest for this study. The Stillwater station is located
approximately I mile from the test pond site. Data for the following parameter w re
obtained: wind speed, wind direction, air temperature, relative humidity, and oJar
radiation. Further details of the weather station network may be found in Elliott et aL.
(1994).
3.4. Model Development
3.4.1. Governing Equations
The governing equation of the model is an overall energy balance on the pond
using the lum ed ca acitance (or lum ed ar~ler) approach:
)
2.)
:>
(31) ~,
.:.I.
where qill is the heat transfer to the pond, quul is the heat transfer from the pond, V is the
pond volume, p is the density of the pond water, cp is the specific heat capacity of the
pond water, and dT is the rate of change of temperature of the pond water. This
dt
approach assumes that tern erature radients within the ~~~~.2.9_c!Y..~~~_~~Eligi~I~:.....
Considering the heat transfer mechanisms shown in Figure 31, Equation 31 can be
expressed to describe the rate of change in average pond temperature as:
6
dT = q solar + q IhermLl( +q convection + q ground + q groundwater +q YOpOTUIWn + q fluid
dt VPCp
(32)
where, qsolar is the solar radiation heat gain to the pond, qtlJemUlI is the thermal radiation
heat transfer at the pond surface, qconvection is the convection heat tran fer at the pond
surface, qgroulw is the heat transfer to/from the ground in contact with the pond, qgroutwlVaJer
is the heat transfer due to groundwater inflow or outflow, qevaporalion i the heat/rna
transfer due to evaporation at the pond surface, and (jjluid is the total heat tran fer to/from
the heat exchange fluid flowing in all spools or coils in the pond. Each of the heat
transfer terms used in the above equation is defined below.
3.4.1.1. Solar Radiation Heat Gain
Solar radiation heat gain (qsolar) is the net solar radiation absorbed by the pond. It
is assumed that all solar radiation incident on the pond urface becom s heat gain xcept
for the portion reflected at the surface.
To determine the reflected component of solar radiation, the angl of incidence
(8) of the sun's rays is first computed at each time step from correlation giv n by
Spencer (1971), Duffie and Beckman (1991), and ASHRAE (1997). The angle of
refraction (Or) of the un's rays at the pond surface i given by Snell' Law. The
reflectance (p') is computed from:
b
(33)
70
where ris the transmittance of solar radiation by the pond orface and the ub cript 'a
refers to the absorbed component. These are computed after Duffie and Beckman (1991)
as:
and
(34)
(35)
where 11' is the extinction coefficient for water, d is the pond depth, nl represents the
parallel component of unpolarized radiation, and r.1 represents the perpendicular
component of unpolarized radiation which are computed after Duffie and Beckman
(1991) as:
)
)I)
... )... )
J:.
tan 2 (8r 8)
r. =_,:'_C..
II tan 2 (8r +0)
sin 2 (8r 8)
r =....,..''
J. sin 2 (8r +8)
Finally, the amount of soLar radiation absorbed by the pond (qsolar) is expressed as:
qsolar =/(1 p')Apond
(36)
(37)
(38)
.,
.:2.
where I is the solar radiation flux incident on the pond surface (here, the total reflectance
is approximated by the beam reflectance) and Apond is the area of the pond surface. The
71
model also accepts solar radiation in the form of beam (Ib) and diffu e (Id) component , in
which case I is computed from:
I =I b cos (} +I d
3.4.1.2. Thermal Radiation Heat Transfer
(39)
This heat transfer mechanism accounts for heat tran fer at the pond surface due to
thermal or longwave radiation. This model uses a linearized radiation coefficient (hr )
defined as:
(310)
where c is the emissivity coefficient of the pond water, ais the StefanBoltzmann
constant, Tpond is the pond temperature in absolute units, and Tsky is the sky temperature in
absolute units. Tsky is computed from relationship given by Bli s (1961) which relate
sky temperature to the dew point temperature (Tdp) and the dry bulb temperatur (Tr/b):
I
(
T Td J4 slcy = Tdb 0.8 +p
250
(311 )

where all temperatures are in absolute units. The model uses the TRNSYS psychrometric
subroutine to compute Tdp if it is unknown. Tdp is computed from either of the following
pairs of state properties: (1) dry bulb temperature (Tdb) and wet bulb temperature (Twb) or
(2) Tdb and relative humidity. The thennal radiation heat transfer (qthennaD is then
computed by:
72
(312)
3.4.1.3. Convection Heat Transfer at the Pond Surface
This mechanism accounts for heat transfer at the pond surface due to free and
forced convection. Several empirical formulations exist for determining the convection
coefficient for different geometries. For a pond surface, correlations for a horizontal flat
plate are the most applicable.
In free convection heat transfer, the Nusselt Number (Nu) is often correlated to
the Rayleigh Number (Ra), which describes the relative magnitude of the buoyancy and
viscous forces in the fluid:
..
•)... )
)
Ra = gf3(!lT)L3
va
(313)
where g is the acceleration due to gravity, ais the thermal diffusivity of air, f3 i the
volumetric thermal expansion coefficient of air, v is the kinematic vi cosity of air, L1T
refers to the temperature difference between the pond and the air, and Lithe
characteristic length. The thermal properties a, fl, and v are aU evaluated at the film
temperature which can be approximated as the average of the air and pond temperatures.
In the model, the thermal properties of air are computed at each time step using
correlations given by Irvine and Liley (1984). For horizontal flat plates, the characteristic
length (L) can be defined as the ratio of the area (A) to the perimeter (P) (Incropera and
DeWitt, 1996):
A
L=
P
7
(314)
In external free convection flows over a horizontal flat plate, the critical Rayleigh
Number is about 107
. Therefore, two empirical relations for the Nu elt number ar used
in the model as described by Incropera and DeWitt (1996) for free convection from the
upper surface of a heated plate or the lower surface of a cooled plate:
I 
(104 < Ra < 107
Nu =0.54Ra 4  laminar flow) (315a)
and
I 
Nu=0.15Ra 3 (l07 > Ra > 1011  turbulent flow) (315b)
The convection coefficient (he) for free convection can then be determined from:
I
I
I.
I~•
)
),
..
Nu k
h =C
L
(316)
where k is the thermal conductivity of air evaluated at the film temperature a with the
other thermal properties described above and L is the characteristic length de cribed by
Equation 314.
In forced convection heat transfer, Nu is a function of the Reynolds (Re) and
Prandtl (Pr) Numbers. The Reynolds number is described as:
)1
Re= vL
v
74
(317)
where v is the wind speed, v is the kinematic viscosity of air, and the characteri tic length
(L) is now defined by the length dimension over which the wind blows and i detennined
from the pond orientation (degrees from north) and wind direction. The Prandtl Numb r
is defined as:
Cl'fl Pr=
k
(318)
where cp is the specific heat capacity of air, fl is the dynamic visco ity of air and k i the
thermal conductivity of air, all evaluated at the air film temperature.
For external forced convection over a flat plate (i.e. the pond surface), the critical
Reynolds Number is approximately 105 (Incropera and DeWitt, 1996). Therefore, two
empirical relations for the Nusselt number are used in the model as de crib d by
Incropera and DeWitt (1996) for forced convection over a flat plate:
to
)
to
),
)
J
I I
Nu =O.664Re 2 Pr 3 (laminar flow regime)
and
(319a)
., I
Nu = 0.037 Re 5 Pr 3 (mixed and turbulent flow) C319b)
The convection coefficient (he) for forced convection can then be determined by Equation
316 with the characteristic length value described by Equation 314 for a horizontal flat
plate.
~~ .d
75
Finally, the convection heat transfer at the pond urface (q nile tio ) i computed
by:
qeotlveclirm = heApOlld (Tair  Tpolld )
I I
(320)
where Tair is the ambient air temperature and he is taken as the maximum of the fr e
convection coefficient and the forced convection coefficient. This practice of choosing
the larger of the free and forced convection coefficients is recommended by Duffie and
Beckman (1991) and McAdams (1954) and is used in the absence of additional
experimental evidence regarding combined free and forced convection.
3.4.1.4. Heat Transfer to the Ground
This heat transfer mechanism accounts for heat conduction to/from the oi I or
rock in contact with the sides and the bottom of the pond. Th.is mechani m of heat _.._... 
transfer is hi hly sitespecific and complexand depends on man factors uch as . 
soil/rock thermal properties,. climatic factors, pond geometry, and thermall.oading
history. In this model, a semiempirical approach developed by Hull et a!. (1984) was
chosen to determine heat losses/gains from the bottom and sides of the pond. Hull et al.
(1984) used a threedimensional numerical code to compute steadystate ground heat
losses from solar ponds of varying sizes, geometries, and sidewall insulation type.
Hull et al. (1984) expresses ground heat losses from any pond as a function of the
pond area, pond perimeter, the ground thermal conductivity (kgruund), and the distance
II
,~
'... :.
.,...
c•.t:'
..
)
7
from the pond bottom to a constant temperature sink.. For practical urpo e th on tant
temperature sink can be taken as the groundwater table (Ki bore and Jo hi, 19 4). For a
rectangular pond with vertical side walls, a heat transfer coefficient for ground heat
transfer (Ugroufld) can be computed from:
Uground = o.999(d kgro",~ J+ 1.37(kgroundPpOlld J
groundwater d pond Apond
(321)
where kground is the thermal conductivity of the ground, dgroundwlIler is the depth to the water
table or the constant source/sink from the ground surface, dprmd is the pond depth, and
Ppond is the pond perimeter. The conduction heat transfer between the ground and the
pond is then given by:
...
I
,I'.
qgrou/ltl =Ugrollnd A,IOI'" (TgroundlValer  T"ond ) (322)
... :5
t is recognized that the above conduction heat transfer model is a relatively
simple representation of the true transient behavior of heat transfer in the ground.
However, ground heat conduction is a relatively minor proces affecting the overall heat
transfer within the pond as compared to other processes.
3.4.1.5. Heat Transfer Due to Ground Water Seepage
This heat transfer mcchanism accounts for inflows and outflows of ground water
to the pond. Although ground water contributions may not be expected in shallow heat
rejecter ponds, this heat transfer mechanism can be used to. a_.c_countfo~ ot.~~!JnflQ.W and ... ..
outflows, such as makeup water or rain water.
. ._.
I
:~
t
77
Tbe volumetric groundwater flow rate (Q) is computed by Darcy Law:
Q = Ki(Ppond (d pond  d gmundwalu) + Apond ) (323)
where K is the hydraulic conductivity of the soil/rock surrounding the pond and i i the
hydraulic gradient. The heat transfer contribution from ground water (qgroundwalu) i then
given by:
...
qgfOUndwUI" =Qpcp (Tgroundwat"  T pond ) (324)
where p and cp represent the density and specific heat capacity of ground water. These
properties of ground water are computed from relationships given in the Handbook of
Chemistry and Physics (CRC, 1980).
3.1.4.6. Heat Transfer Due to Evaporation
This heat transfer mechanism is the most important one contributing to pond
cooling. This model uses the jfactor analogy to compute the mass transfer of
evaporating water (m: )at the pond surface:
).. I..
•
... J·•
II
t:
t
)
t
(325)
where hd is the mass transfer coefficient, Wair is the humidity ratio of the ambient air, and
wsuifrepresents the humidity ratio of saturated air at the pond surface. The wsuifterm i
7
computed using the TRNSYS psychrometric subroutine by etting T db and Twb equal to
the pond temperature. The Wair term may also be computed u ing the TRNSYS
psychrometric subroutine depending on what two state properties of the air are known.
The model accepts the following pairs of state properties for the calculation of Wair if it is
unknown: (1) Tdb and TWb' (2) T db and relative humidity or (3) Tdb and Tdp. The mas
transfer coefficient (hd) is defined using the ChiltonColburn analogy as:
(326) ....
where he is the convection coefficient defined previously, cp is the specific heat capacity
of the air evaluated at the pondair film temperature, and Le is the Lewis number. Le is
computed as:
)·•
..
I.... •;.
a
Le=
DAB
(327) ....
where a is the thermal diffusivity of the air evaluated at the pondair film temperature
and DAB represents the binary diffusion coefficient. The thermal propertie (a and cp) of
air are computed at each time step using correlations given by Irvine and Liley (1984).
DAB is computed after Mills (1995) who references Marrero and Mason (1972):
1.87xl OloT2.072
DAB =
~lir
(280K < T < 450K) (328)
7
where T refers to the pondair film temperature in ab olute units and Pair i the
atmospheric pressure in atmospheres.
The heat transfer due to evaporation (qevaporalion) is then computed by:
q~vaporalion = h/8 Apondm: (329)
where hfg is the latent heat of vaporization and is computed at each time tep from the
relationship given by Irvine and Liley (1984).
3.4.1.7. Heat Transfer Due to the Heat Exchange Fluid
Heat transfer due to the heat exchange fluid represents the pond thermal load.
~.._~
This model has been developed to account for water or antifreeze as the heat exchange
fluid. The thermal properties of the fluid are computed at each time step from
correlations given in the Handbook of Chemistry and Physics (CRC, 1980) for water and
from correlations given by Wadivkar (1997) for an antifreeze solution. The thermal
properties are computed a~~w~e_~~iQJ~J!l.Per~!.~_T, u!d). This temperature is
computed as the average of the inlet and outlet temperatures at the given time step. Since
the outlet temperature at any current time step is not known, the previous converged
value is used as an initial guess and calculation of Tfluid is iterative. Solution of the pond
temperature is also an iterative procedure as discussed below.
The heat transfer due to the heat exchange fluid ({jj1uid) is computed by:
).
...
I
;
:
qjluid =UApipe (TjlUid  T pond )(Ncircuit) (330
where UApipe is the overall heat transfer coefficient for the pipe expressed in terms of
inside pipe area and Ncircui/ refers to the number of flow circuit (i.e. the number 0
spools) installed in the pond. Thus, Equation 330 is based on the assumption that one
spool is one flow circuit and that the flow rate is divided evenly between the circuit in a
parallel arrangement. The term UApipe is expressed in terms of the inside pipe area as:
(331 ) •
)
where rj is the inner pipe radius, Lspoo{ is the length of one spool or circuit, and J:R,
re resen_ts.t_he comosite thermal resistance which is defined by the following resistance
network:
•
•
:
LRt =Ri + Rpipe + R" + if (332)
where Ri is the thennal resistance due to fluid flow through the pipe, Rpipe i the pipe
thermal resistance, R, is the thermal resistance at external pipe surface, andifrepre ent a
fouling factor at both the inner and outer pipe walls. The resistance terms are defined as
follows (in terms of inner pipe radius):
1 R=i
h.'
r
(333)
and
r. (r J Rpip~ = 'In .....!!....
kpip~ ri
(334)
(335
where hi is the convection coefficient due to fluid flow through the pip , kpipe i the
thermal conductivity of the pipe material, ho is the convection coefficient at the out r
surface of the pipe, and ri and r o are the inner and outer radii of the pipe, respectively.
The above convection coefficients are determined u jng correlation for the
Nusselt number in flow through a horizontal cylinder, 'nc 0 s ecific correlation exi t  
'.
for a sl" A constant heat flux at the pipe surface is assumed.
..
In~~~aQ§J~r due to internal flow, Nu is a function of the Reynold
and Prandtl numbers. Determination of Re is described in Equation 317. For thi ca
v represents the mass flow rate of the heat exchange fluid, v represent the kin matic
viscosity of the heat exchange fluid, and the characteristic length (L) is the inner pipe
diameter. A Reynolds number of 2000 is a umed to be critical. For laminar, fullydeveloped
flow in the pipe (Re<2000), the Nusselt Number is a constant equal to 4.36
(Incropera and DeWitt, 1996, Equation 853). For turbulent flow, the DittusBoelter
relation is used to compute the Nus elt Number:
•
:
·
4
Nu; =O.023Re 5 prA' (336)
2
where Pr is defined by Equation 318 for f1, cpo and k representing the thermal prop rtie
of the heat exchange fluid evaluated at its average temperature. The value of the
exponent x in Equation 336 is dependent upon whether the entering fluid i being heat d
or cooled; x is equal to 0.3 if the entering fluid is greater than the pond temperature and x
is equal to 0.4 if the entering fluid is less than the pond temperature. The inside pipe
convection coefficient can then be determined by using Equation 316 where Nu is NUj, k
is the thermal conductivity of the heat transfer fluid, and L is the characteristic length
which is defined in this case as the inner diameter of the pipe (Incropera and DeWitt,
1996).
Infreecon_vection at the external pipe surface, the Nusselt Number is a function .._ _ _._ . ..
of the Rayleigh Number. Ra is computed using Equation 313 wnere ~ /3, and
v represent the thermal diffusivity of water, the volumetric thermal expansion coefficient
of water, and the kinematic viscosity of water, respectively, all evaluated at the pip film
temperature, which is approximated as the average of the pipe surface and pond
temperatures at the given time step. The term ,1T refers to the temperature difference
between the pipe surface and the pond temperatures. Nu for free convection from a
horizontal cylinder is defined as (Churchill and Chu (1975) ):
2
..
)

:
...
I
O.387Ra 6
Nu o = 0.60+ 8
(1 +(O.559/Pr)':r (337)
t t
where Pr is defined by Equation 318 for Jl, cp, and k representing the thermal propertie
of the pond water evaluated at the pipe film temperature. Now the external pip
convection coefficient can be determined by u ing Equation 316 where Nu is Nuo k i
the thenna! conductivity of the pond water evaluated at tbe pipe film temperature, and L
is the characteristic length which is defined now as the outer pipe diameter.
The outlet fluid temperature (Tout) is computed from an overall energy balance on
the pipe:
T  T _ qcirclJit
out  fluid 2'
mcl'
(338)
where rh is the mass flow rate of the heat exchange fluid per flow circuit cp is the specific
heat capacity of the heat exchange fluid, and qcircuil is the heat rejected/extracted by one
flow circuit. This outlet temperature is u ed to compute the average fluid temperature at
the next iteration as described above.
3.4.1.8. Solving the Overall Energy Balance Equation
The differential equation describing the overall energy balance on the pond
(Equation 32) is rearranged in the form:
:.
(339)

4
where T represents the pond temperature, XI contain all term of Equation 32 th
multiply T, and X2 contains all terms of Equation 32 that are independent of T.
Equation 339 is a linear firstorder ordinary differential equation which i solved
at each time step using the exponential function ~ an integr~~in~ factor. The olution i
. b ~
gIven y the TRNSYS differential equation solver subroutine as:'
. 
where Tt is the average pond temperature at the new time step and T,..1I is the average
pond temperature at the previous time step.
Many of the quantities in the heat transfer equations described above require that
the average pond temperature at the current time step be known. Thu , the actual pond
temperature is arrived at iteratively. A convergence criterion for the pond temperature of
3.4.2. Computer Implementation
Thc component configuration for the pond model is hown in Figure 32. A
C9J!!P'~ model was also developed which manipu~es any weather dat~Jleeded.ior the
pond model. The weather component model makes use of the TRNSYS psychrometric   
subroutine to compute moist air properties given two known state properties. The two

•
..

state properties are dry bulb temperature and one of wet bulb temperature, relative
humidity, or dew point temperature. The weather component model al 0 compute the
sky temperature, the solar radiation on a horizontal surface, and the olar incidence angle.
A computer algorithm is shown in Figure 33 in the fonn of a flow chart.
air sky wind solar angle inlet mass
temperature temperature direction of incidence
inlet
flow rate
1 humidity1wind 1solar 1 fluid 1 ratio speed radiation temperature
~ ~ + +
POND MODEL PARAMETERS:
1. initial pond temperature
3. pond length
5. pond depth
7. extinction coefficient for water
9. spool length
11. pipe thermal conductivity
13. fluid type (water or antifreeze)
15. ground thermal conducti vity
17. ground water or far field temperature
19. hydraulic gradient
2. pond orientation from north
4. pond width
6. emissivity coefficient
8. number of spools or coil
10. pipe outer diameter
12. pipe wall thicknes
14. antifreeze concentration if used
16. fouling factor
18. soil hydraulic conductivity
20. depth to water table
l 1 1 1
pond outlet fluid mass flow heat
temperature temperature rate rejected/ab orbed
Figure 32. Pond model component configuration.
Wealher dam from
compOllenl model
Auid temperature
and flow rate from
upstream component
y
Yes Set pond tempernlu Ie
equal 10 rUlal value
or previou lime step
Compute the average temperature of the
heat exchange fluid
Call TRNSYS differential equation solver
to compute the average pond temperature
Compute new outlet fluid temperature
No
Figure 33. Pond model computer algorithm.
87
3.5. Results and Discussion
3.5.1. Model Comparison to Experimental Results with No Heat Rejection
The first step in the ~odel verification pr.~~ was to compare the model pond
temperatures to measured pond temperatures during time when no heat was being
rejected to the ponds. This comparison allowed a validity check of the simulation of the 
several environmental heat transfer mechanisms occurrin within the ponds. Simulated
and actual pond average hourly temperatures are shown in Figure 34 for an 8da eriod
in July 1998 when no heat was rejected to the ponds. Therefore, in these case, the model
is driven by weather data input only.
A review of the plots in Figure 34 shows that the model temperatures compare
favorably to the measured pond temperatures. The simulated temperature are within 3°F
 .
(1.67 °C) of the observed temperatures throughout the test period. The difference
. '._.
between the average simulated pond temperature and the average o~_~~~~ pond
temperature for theentiretest_per.iod is 1.93°F (1.07°C) for the 2feet deep pond and 1.55°F (O.86°C) for the 3.5feet deep pond. Shallow ground water was not encountered at
r 
the site and therefore ground water contributions were not considered. .!'.?uling of the
heat exchanger pipe was also not considered.

110
105
100
~
~ 95
!! .a
I!! 90
~
E 85
~
80
75
(a)
.. ji fI a
I II I II "
j
.r. If
y, .~\ II rJ.."' I~.. : ,., ~
It .. ; ; of' I: 1\ ,.. ,..,
'1\ ; 1\ ; \'1 I '.\ . :\ .. • ° .. 1\ . • \, . I'"
!II \/
\ It I~ If,
~ I~ 1'1
43.3
40,6
37,8
U
35.0 ~
! .a 32,2 I!!
Q)
Do
29.4 E
~
26,7
23.9
m ~.1
21Jul 22Jul 23Jul 24Jul 25Jul 26Jul 27Jul 28Jul 29Jul 3O..Jul
I Air =Actual Pond· M~el Pond'
(b)
110
105
100
G:'
~ 95
! .a I.'.l.l 90
~
E 85
{!!.
80
75
._~..._~._ .........._ ,.._..._ ~
.II , fI a
JI ~
J II
j
~ Ii 1/ If f'\ 11'0
It : \ .. f.\ :I..
,
: r;;.
.. ,'" I' l' ° ,. . ... .. : ..
n
... ~
1\ 1\, 'II I'
I ~ l'j IV
43.3
40.6
37.8
G'
35.0 ~
!
32.2 .~
i
29.4 E
~
26.7
23.9
70 21.1
21..Jul 22..Jul 23Jul 24..Jul 25..Jul 26..Jul 27..Jul 28..Jul 29Jul 3OJul
AirActual Pond······· Model Pond
Figure 34. Comparison of observed and simulated average pond temperature
with no heat rejection in the (a) 2feet (0,61 m ) deep pond and (b)
3.5feet (1.07 m) deep pond.

8
3.5.2. Model Comparison to Experimental Results with Heat Rejection
Heat rejection to the ponds was simulated over a 25day period from November
12 to December 7, 1998. Input data to the model consisted of weather data as de cribed
previously in addition to measured slinky heat exchanger supply water temperatures and
flow rates on 6rninutely time intervals. The model performance was evaluated by
comparing (I) the simulated to the observed return temperature of the heat exchange fluid
and (2) the simulated cumulative heat rejected to the ponds to the measured water
heating element and pump power input. These comparisons are shown in Figures 35 and
36 respectively. As with the previous comparisons, ground water contribution and
fouling of the heat exchanger pipe were not considered.
A review of the temperature plots in Figure 35 shows that model fluid return
temperatures compare favorably to the observed fluid return temperatures. The average
observed and modeled fluid return temperatures over the test period i.n the 2Jeet (0.61
meter) deep pond were 70.5°F (21.4°C) and 70.2°F (21.2°C), respectively, and in the 3.5
feet (I.07meter) deep pond were 69.2OP (20.7°C) and 70AoF (2] .3°C), respectively. The
deeper pond has slightly larger differences between modeled and observed fluid return
temperatures. The error is small, however, and is probably acceptable for purposes of
simulating hybrid GSHP systems; even a 20P (1.1 ]°C) error in return fluid temperature
from the pond will cause only a slight difference in modeled heat pump performance.
/1
/
'.
Il ~
. ~ i " .'rI
i. 14
If I i~' ~
J l , IV ., .'
1\ ~
~ 1
\
,
76
75
74
73
li:'72 o
; 71
~ 70
~ 69
& 68
E 67
~ 66
65
64
63
62
12Nov 17Nov 22Nov
(8)
27Nov 2Dec
24.4
23.9
23.3
22.8
22.2 Y
21.7 
CD
21.1 ~
20.6 1ii
20.0 t.
19.4 E
18.9 ~
18.3
17.8
17.2
16.7
7Dec
o
E !ctual Return· M:>~,~h