WATER RESOURCES RESEARCH, VOL. 18, NO. 3, PAGES 571-587, JUNE 1982
A DimensionlessParameter Approach to the Thermal Behavior
of an Aquifer Thermal Energy Storage System
CHRISTINEDOUGHTY,GORANHELLSTR6M,
1 AND CHIN FU TSANG
Earth Sciences Division, Lawrence Berkeley Laboratory, University of California
Berkeley, California 94720
JOHAN CLAESSON
Department of Mathematical Physics, Lund Institute of Technology, Lund, Sweden
To predict the performanceof an aquifer thermal energy storagesystem, an understandingof the
system'shydrothermalbehavior is needed.One possibilityis to run a detailed numerical simulationof
the system. However, for a single-wellsystem in which fluid flow is limited to steady radial flow, a
characterization
scheme
basedon a setof fourdimensionless
parameter
groups
allowsproduction
temperaturesand energyrecoveryfactorsto be read from graphs.The assumptionof radial fluid flow
is valid when buoyancyflow can be neglectedand a well is fully screenedin a horizontalaquifer which
is confinedabove and below by impermeablelayers. Criteria for little buoyancy flow include a low
permeability or vertically stratified aquifer, a small temperaturedifference between injected and
ambient water.
and short cycle
leneth
.,
,..,
The ha.•ie enerov
tranannrt
•,.•
.......
r
........
layer system with steady radial fluid flow in the aquifer are nondimensionalizedto derive the key
parameter groups. Next a numerical model which calculatesthe heat transfer in the aquifer and
confininglayers for an injection-storage-production
cycle is run for a rangeof valuesof thesegroups.
The calculatedproductiontemperaturesand energyrecoveryfactorsare then presentedgraphicallyas
a function of the parametergroups.Comparisonsbetween resultsof field experimentsand recovery
factors read from the graphsshow good agreement.
1.
INTRODUCTION
The concept of heat storage in an aquifer has attracted
increasinginterest during the last 5 years [Lawrence Berkeley Laboratory, 1978;Tsangeta!., 1980].The basicidea is to
inject hot water into an aquifer during periodsof low energy
demand and then, when energy demand is high, to extract
the water and use the heat. The energy supply may come
from a variety of sources, such as solar energy, industrial
waste heat, or power plant cogeneration.An aquifer naturally providesa large volume of water as the storagemediumat
a relatively low cost. Water is extracted from the aquifer
through a supplywell and, after beingheated, is injectedinto
a storage well in the same aquifer. After the hot water has
been producedfrom the storagewell and the heat used, the
water is reinjected into the supply well, thereby creating a
closed system with little net groundwater consumption.
Aquifer thermal energy storage(ATES) is consideredto be
one of the most promisingand cost effective alternativesfor
low-temperature heat storage on a large scale. Numerous
theoretical investigationsestimatethat as much as 80-90%
of the injected energy may be recovered during a seasonal
[Molz eta!., 1981; Yokoyarna eta!., 1978] or unconfined
aquifers [Mathey, 1977; Fabris and Gringarten, 1977; Iris,
1979; Reddell eta!., 1978] have been completed. Most of
these aquifersare rather shallowand of high permeability.
Experiments involving the storage of cold water for air
conditioning have also been done [Yokoyarna et al., 1978;
Reddell eta!., 1978]. The highest injection temperature
(55øC)andlargestinjectionvolume(55,000m3) experiment
to date has been carried out by Auburn University (Alabama) where two consecutive6-monthcyclesyieldedrecoveries of 66 and 76% of the injected energy. These are
promising values consideringthe relatively small injection
volume. However, many problems, both theoretical and
practical, need to be solved in order to make it possibleto
assessthe feasibility of hot water storagein a specificaquifer
at a given site. These problems relate to legal and environmental issues,water chemistry, soil mechanics,and thermohydraulics.
This study considersthe thermal behavior around a storage well in the casewhen buoyancyeffectscan be neglected.
A dimensionlessformulation of the energy transport equaATES cycl• [Tsangeta!., 1977;Fabris et al., 1977]. tions for the aquifer system is presented, and the key
Recently, Sauty et al. [1980a] have studied the thermal dimensionlessparametersare discussed.A numericalmodel
using a steadyflow field is used to simulateheat transportin
behavior of a simplifiedsingle-wellATES systemin terms of
the aquifer and confininglayers aroundthe injection/produca set of dimensionlessgroups.
A few heat storage field experimentsin either confined tion well during the ATES cycle. The key parameters are
varied in order to understandtheir influenceon the percent
of the injected energy that can be recovered and the temperature of the extracted water. The results are presented
•Alsoat Departmentof Mathematical
Physics,LundInstituteof
graphically. Finally, some comparisonswith field experiTechnology, Lund, Sweden.
Thispaperis not subjectto U.S. Copyright.Published
in 1982by ments are given to illustrate the use of the dimensionless
groups and graphs. This study follows similar lines to that of
the American GeophysicalUnion.
Sauty et al. [1980a], but considers a more general system,
Paper number 2W0171.
,
and extends the results in different
571
areas.
572
DOUGHTY ET AL' AQUIFER THERMAL BEHAVIOR
during the productionperiod. The energyrecovery factor e
for each cycle is defined as the ratio betweenthe produced
and injected energy when equal volumes of water are
injected into and produced from the aquifer. The energy
content of the water is defined using the original ambient
temperature of the aquifer, To, as a reference. The recovery
To
Xc,Cc
ConfiningLayer
:
O
:
,
Xe,Ce,d)l,d.L
Aquifer
factor
'
'
R
•c,Cc
e is
'
e=
Confining Layer
(CwTe - CwTo)Qedt
fti+ts+tp
J ti+ ts
Fig. 1. Schematicdrawingof the ATES system.
2.
2.1.
CONCEPTUAL
This expressioncan be written as
MODEL
Definitions
The ATES system considered (Figure 1) consistsof a
singleinjection/productionwell that fully penetratesa laterally infinite, horizontal aquifer of uniform thickness H.
Resultsare also applicablefor a multiple-wellsystemwhere
well spacing is large enough so that the thermal behavior
around the storage well is not significantly affected by
neighboringwells. This has been studiedto someextent by
Tsang et al. [1978]. Under the single-wellidealization,there
is radial symmetrywith respectto the well. Furthermore,the
aquifer is assumedto be homogeneouswith bulk thermal
conductivity ha, bulk heat capacity per unit volume (solid
plus fluid) Ca, and dispersivelengthsd• and d• (see notation
list for definitionsof all symbols).It is boundedabove and
below by impermeableconfininglayers which may be of
arbitrary thicknessand composition.This studyconsidersa
systemwith a cap rock of thicknessD and an infinitelythick
bedrock. Both cap rock and bedrock are homogeneousand
have bulk thermalconductivityhcandheatcapacityper unit
volumeCc. The heat capacityper unit volumeof water is Cw.
All materialpropertiesare assumedto be temperatureindependent. Initially, the entire systemis at constanttemperature To; above the cap rock the temperatureis kept at To
throughoutthe cycle.
The length of the ATES cycle is tc. The duration of
e=
Tl-
To
(3)
where•p denotes
theaverage
temperature
of theproduced
water during the production period. A dimensionlesstemperature T' is defined as
T' =
T-To
Tl-
To
(4)
The expressionfor the recovery factor (2) then becomes
= t.'
(5)
where•p' is theaverage
dimensionless
temperature
of the
produced water.
In most applicationsthere will be a cut-off temperature,
below which energy is not usable, which will probably be
higher than To. The recovery factor eref calculated with
respect to a dimensionlessreference (cut-off) temperature,
rref' , is
eref= (e -- rref')/(1 - rref')
(6)
is purely radial.
At the end of injection, the volume of injectedwater,
If the temperature of the produced water falls below the
reference temperature during a part of the production period, (6) is not correct becausethe productionwill stop when
the reference temperatureis reachednot when the produced
volume is equal to the injected volume. In such a case, (2),
with To replaced by Tref, may be used to calculate the
percent of energy above Trefthat is recovered. However,
this percent is not a recovery factor as we have defined it
becausethe produced and injected volumesare not equal.
Throughoutthis paper, the injectedenergyis referred to as
heat; however, all the resultsdiscussedalso apply to chilled
occupies
a regionin theaquiferVw= cb•rRw2H,
where•bis
water storage.
aquifer porosity and Rw is the maximum radial extent of the
injected water. It is convenientto defineanothervolume, the
thermal volume V • (Cw/Ca)Vw.The thermal volume is the
The influenceof regionalgroundwaterflow on the thermal
recovery of the storage system has not been studied. However, if this influence may becomelarge, the storageregion
in the aquifer will have to be protected, for example, by
using boundary wells [Tsang and Witherspoon,1975; Molz
and Bell, 1977; Whiteheadand Langhetee, 1978].
injection,storage,production,
andrestperiodsare ti, ts,tp,
and t•, respectively.The temperatureof the injectedwater is
constant, T•. The magnitude of the volumetric fluid flow
rate, Q, is kept constant during injection and production
periods,and the injectedandproducedvolumesare takento
be equal,thatis, Qiti-- Qptp• Vw.In the aquifer,fluidflow
cylindrical volume in the aquifer (solid plus fluid) which
would have, at constant temperature T•, thermal energy
equal to the total heat energy of the injected fluid. The
thermalvolumemayalsobe writtenas V = •rR2H,whereR
is the thermal radius. Hence, R is defined as
R = (V/qtI-1)
1/2= (Vwqw/qaqtI-1)
1/2
(1)
An essentialresult of our calculationsis the time-varying
temperatureTp of the water extractedfrom the aquifer
2.2.
Buoyancy Flow
In this study the only fluid flow is a steadyradial flow, and
fluid density is assumedto be constant. Hence buoyancy
flow induced by the density difference between injected
warm water and the cold water in the aquifer is neglected.
This is somewhatjustifiedfor low permeabilityor vertically
DOUGHTYET AL.: AQUIFERTHERMALBEHAVIOR
stratified aquifers, small temperature difference between
injected and original waters, or short cycle length. Under
theseconditions,buoyancyflow can be neglected.Additionally, for geometriesin which the thermal radius is much
greater than the aquifer thickness,a moderatetilting of the
thermal front caused by buoyancy flow may have only a
573
Further, if we define
Pe = QCw/2•rhaH
(12)
as the Peclet number, then usingthe dimensionlesstemperature T' (4), (7)-(10) become
li2T'
small effect on the overall thermal behavior of the system. A
1 tiT'
li2T'
ha liT'
lir'2+ r' lir'+ &,2- hclit'
theory proposedby Hellstr6m et al. [1979]givesa formula
z'>
H
-2L
(13)
for the characteristic time constant, to, for the buoyancy
tilting rate of the hot-cold water interface. This formula,
which is given in Appendix A, may be used to determine
caseswhere buoyancyflow may be neglected.On the other
hand, the conclusionspresentedin this paper may still be
applicablein a qualitativesense,even in caseswherebuoy-
lieT,
1 liT'
lieT,
1 liT'
Ca liT'
(14)
T'lz,=(m•+o= T'lz,=(mm-o
ancy flow is significant.
liT'
2.3. DimensionlessFormulation and Definition
of Dimensionless Parameters
t
The thermalbehaviorof the aquiferstoragesystemmay be
expressedin dimensionlessform. We assume,as a reason-
H
lir'• + r' lir'+ li'-•7•-Per'lir'- Cclit' z'<•2L
(15)
(16)
z'=(H/2L)+O hc •Z'
z' =(H/2L)-O
Basedon the form of theseequations,thetemperatureat any
point may now be written as
able approximationto the seasonalvariationof supplyand
demand of energy, that the injection, storage,production,
and rest periodsare of equaldurationsothat ti -- ts = tp = tr
= td4. In this sectionwe shallalsoassumethe caprock to be
The temperatureTp'of waterproduced
fromtheaquiferwill
of infinite thicknessand the dispersivelengthsto be zero.
be an average of the vertical temperaturedistributionin the
The effects of unequallengthperiods, a finite cap rock and
aquifer at the well, which is located at r' = 0. It is given by
velocity-dependentdispersionwill be treated as special
cases (sections3.2.1-3.2.3) and are not consideredin the
' =
T z',t' Pe,
dz'
following formulation.The systemis symmetricabout the
midplaneof the aquifer, z = 0. Hence only the regionz >- 0
need be considered.
The factor H/2L can be written in terms of C,/Cc, ha/hc,and
The temperature field in the confininglayers is governed a parameter, A, introducedby Fabris et al. [1977], which is
defined as
by the ordinary heat conductionequation, namely,
=
r',z', ,Pe,cc,hc,2L
Tp • j-w2L
li•'T
hc•
1 lit
+ her
lir
li:T
lit
z>--2
(7)
1 lit
Xa-• -[-Xar
lir
li2T CwQ1 lit
-[-ha liz2
2•rH r lir'
lit
= Ca•lit
(8)
t qahaAI (20)
r•'= r•' t',Pe,
cc,hc,
The recovery factor (5) is the averagetemperatureof the
produced water during the production period. The dimensionlessgroups controllingthe recovery factor are summarized as follows:
H
Pe = QCw/2•rhaH
= CaR2/2hati
The last term on the left siderepresentsthe convectiveheat
transport in the aquifer due to pumpingat the well. Temperature
and heat flow
must be continuous
at the interface
Ca/Cc
T[z=(m2)+
o = T[z=(m2)_
o
= ha
z=(H/2)+O
(9)
(10)
z=(H/2)-O
Let us choosethe following dimensionlessparameters:
t' = t/ti
r' = r/L
(21)
A = Ca2U2/Cchcti
between aquifer and confininglayer, i.e.,
hc
(19)
This implies that the temperatureof the producedwater is
Assumingan incompressibleradial flow in the aquifer, the
heat balance equation in the aquifer may be written as
li2T
Ca'ha'2L
A = Ca2U2/Cchcti
H
+ hc•=
Cc-liz2
lit
'
(17)
z' = z/L
where
ha
hcti,•
1/2
(11)
the number of cycles
2.4.
Steady Flow Model
Detailed numerical models which solve the coupled mass
and energy transfer equations of fluid flow in a porous
medium have been successfully used to match analytical
results [Tsang et al., 1977]as well as field data [Tsang et al.,
1980; Papadopulos and Larson, 1978] for ATES systems.
Generally, these modelsare expensiveand time consuming
to run. For the purposeof gaininga better understandingof
the processescontrolling heat loss in an ATES cycle, it is
often
desirable
to do a series
of simulations
in which
574
DOUGHTY ET AL..' AQUIFER THERMAL BEHAVIOR
Zl
perature T• is assigned to the first cell in each row. This
translation every time step, ti/M, simulatesa constantvolumetric fluid flow rate at the well:
Qi=CaTrR2H/Cwtim3/s
Z$
(22)
and a horizontal Darcy velocity:
o(r) = Qi/2,rHr m/s
(23)
at radius r in the aquifer. When t = ti, the temperature field
has been translated M times, and the thermal front radius is
R.
Heat transfer by convection is accountedfor by translation of the aquifer temperaturefield every time step ti/M.
Heat transfer by conductionis described by the ordinary
heat equation:
C•
&
= V (XVT) = -V.q
(24)
where q is the heat flow per unit area. The integral form of
(24) is solved numerically for each cell in the mesh during
every timestep At. For the (rn, n)th cell in the mesh it is
written
Radial
distance
Fig. 2. Scale drawing of part of the equal volume mesh. The
mesh may extend farther vertically than is shown to simulate
infinitely thick confininglayers.
dtfferentparametersare systematicallyvaried. To this end,
the present study uses a simplified, but fast, numerical
model of massand energytransfer in an aquifer system.This
model, the steady flow model (SFM) was developedby Lund
University [Hellstr6rn and Claesson,1978]and modifiedat
Lawrence Berkeley Laboratory for our present study.
Rather than solvingthe masstransfer equationto obtain a
fluid flow field as a function of time, a steady horizontal flow
field is prescribed in the aquifer, leading to an energy
transfer equation given by (8).
In general, the numerical simulation of a combined convection and conduction equation such as (8) introduces a
mesh-dependentnumerical dispersion,which causesspurious thermal front smearing.The SFM avoids this effect by
simulatingthe flow field in a specialway describedbelow. In
contrast, when the coupled equations for mass and heat
transport are solvedwith temperaturedependentproperties,
as in most detailed numerical models, it is not possible to
avoid a certain amount of numerical dispersion.
The flow field is simulated in the following way. The
thermal front radius at the end of the injectionperiod, R, the
number of mesh cells in each row between r = 0 and R, M,
and the duration of the injection period, ti, are given as input
to the program. A mesh is constructedin which all the cells
of a horizontal row have equal volumes. Owing to the
cylindrical symmetry, the radial dimension of the cells
decreases
as the radial distance from the well to the cell
increases,as shownin Figure 2. The lengthof the rnthcell is
fvc 15Tdv=
- fur,
, V'qdV=
- faq'•da
m,n
,n
m,n
(25)
where the right-hand side describesthe heat transfer to all
neighboringcells. Using the explicit finite differenceapproximation, (25) becomes
n[rm'n(t
+At)
- rm
'n(t)]
At
Cm,nVm, '
= qrm'n2•rRm(Zn+•- Zn)- qrm+l'n2•rRm+•(Zn+•-- Zn)
+ qzm'n•'(Rm+l2 - Rm2) - qzm'n+l•r (Rm+l
2 - Rm2) (26)
where
qr
m'n
=(Trnl,n
-Tm
n)
(Rm
-Rm-1
+Rm+l
--Rm,)-I
'
2km-l,n
2}•m,n
(27)
Zn
--Zn_.._E1
Zn
+l+Zn
)-1
qzm'n
-' (Trn'n-l
- Tm'n)
'•mn-1
+ •,,n
In the above equations,Cm,n,Tm,n,•,m,n,and qm,nrepresent
the average value of C, T, h, and q, respectively, over the
(rn, n)th cell which has volume Vm,n.The parameterh may
include dispersioneffects, as describedin section3.2.1. The
time step At is chosenso as to ensure numerical stability of
the solution.
During the storageand rest period no translationof the
temperature field occurs, and heat transfer is purely by
conduction.During the productionperiod, the convectionis
treated as during the injection period. The length of the
Rm+l- Rm,whereRm= [(m - 1)/M]I/2R.Thenthevolume productionperiod,t•,,is givenas input.For everytime step
of the rnthcell, whichis proportional
to Rm+l2 - Rm2, is tp/M the temperaturedistributionis shiftedonecell toward
independent of rn. During the injection period, whenever
time t is equal to ti/M, 2ti/M, 3ti/M, ß ß ß , ti, the temperature
distribution in the aquifer is translatedhorizontally one cell
away from the well, and the user-specifiedinjection tern-
the well, and the temperaturesfrom all the cells directly
adjacent to the well are weighted according to their heat
capacity and averagedto give the productiontemperature.If
t•, is differentfrom ti, the fluidflow rate will be adjustedso
DOUGHTY ET AL.' AQUIFER THERMAL BEHAVIOR
575
TEMPERATURE
FIELDS
SIM.ULATED BY STEADY FLOW MODEL
With conduction
and convection
With no conduction
t = 30 days
During injection
-2O
-20
-3O
-30
--•••99
,6
,I
-40
••_••_
,4
½3 -50
-60
-70
......
0
t = 90 days
End of injection
I0
-70
20
30
40
50
60
0
I0
20
30
40
50
60
-IO
-20
-50
-40
-50
-70
0
t = 180 days
Endof storage
I0
20
30
40
50
60
-2O
-•øo ,'o •'o •'o .;o ;o ;o
-20
-50
-40
:?
-50
0
I•
2'0 3'0 4'0 50 60
0
I0
20
•'0
40
50
•0
,
-10
t = 270 days
End of production
-2O
-2O
-7O
-70
o
•o
2o
30
40
50
60
0
I0
20
50
40
50
60
Radial distance (m)
ß =1
•1
Fig. 3. Temperaturedistributionsat various times during the first cycle generatedby the SFM with and without
conduction.
that the volumes of water injected and producedwill remain
3.
RESULTS
fixed, i.e., Qiti = Qptp.
Figure 3 shows the temperature distributions at various
times during the first cycle as generatedby the SFM with
and without conduction to illustrate the superpositionof
conductionand convection.A typical meshconsistsof about
1000 cells. Nearly 300 different caseshave been simulated,
with M ranging from 10 to 40. The computer time required
for a typical annual cycle is about 15 s on a CDC 7600.
The resultsof the largenumberof differentATES systems
that have been simulatedusingthe steadyflow model are
presentedin graphicalform. Section3.1 showsthe dependenceof the thermalbehavioron the dimensionless
groups
derived
in section
2.3aswellasthedependence
onsome
of
the individual parametersthat make up the dimensionless
groups.Section3.2 discusses
the effectof a velocity-depen-
576
DOUGHTY ET AL.: AQUIFER THERMAL BEHAVIOR
RECOVERY
FACTOR vs CAPACITY RATI 0
1,0
Pe, A
200,800
-0,8
200,20
0,6
_
5,800
5,20
cycles is shown in Figure 6 for different combinationsof Pe,
A, and ha/hc.For values of Pe larger than 200, production
temperature shows little dependenceon Pe.
To demonstratethe effect of a cut-offtemperature,consider the casewith ha/hc- 1, Pe = 20, A = 50. From Figure 5,
the first cycle recovery factor is 0.59. From Figure 6, the
final dimensionlessproduction temperature for the first
cycle is about 0.3. For an application that can only use
energy above Tref' = 0.25, equation (6) gives Eref-' 0.45. If
Tref'were greater than 0.3, equation (6) would underestimate
Eref.
0,2
-
)to
--'---- Fifthcycle
• =I
I
First
cycle
a number of cases have been studied in more detail.
I
Co/Cc
Fig. 4.
First and fifth cycle recoveryfactorsas a functionof heat
capacity ratio Ca/Cc, for various values of Pe and A.
dent dispersion, unequal length periods, finite cap rock,
long-term behavior, and multilayer flow.
3.1.
3.1.2. Results in terms of individual parameters. To
showthe explicitdependenceon certainphysicalparameters
Dependence on Parameters of the Dimensionless
Formulation
In section 2.3 it was shown that the recovery factor is a
function of five dimensionlessparameters. One of these is
the number of cycles. The recovery factor increasesas the
number of cycles increases, with the increment decreasing
for later cycles. The results given in the following sections
(3.1.1 and 3.1.2) are for the first and fifth cycles. Long-term
effects are discussed in section 3.2.3.
Another parameter, the ratio between heat capacity in the
aquifer and heat capacity in the confining layers, Ca/Cc,
varies within a small range and is mainly determinedby the
water content (or porosity) of the two layers. The correspondingly small variation of the recovery factor with respect to the capacity ratio, which is shown in Figure 4, is
largest when A and ha/hcare small. The resultswhich follow
are given for Ca/Ccequal to 1.25, i.e., when the aquifer has a
higher water content (porosity) that the confininglayers.
3.1.1. Results in terms of dimensionlessgroups. Of the
five dimensionlessparametersintroducedin section2.3, the
effects of three of them, Pe, A, and ha/hc,are more critical
and remain to be examined. The range of Pe and A is chosen
with respect to conditions met in seasonalor daily storage.
The ratio between thermal conductivity in the aquifer and
the confininglayers, ha/hc,dependslargelyon the magnitude
of dispersive effects caused by the flow in the aquifer.
Dispersion is discussedin section 3.2.2.
Recovery factor: Figure 5 shows the energy recovery
factor as a function of Pe and A for the first and fifth cycles.
Three different values of the thermal conductivity ratio, ha/
hc, have been used. The sensitivity of the recovery factor to
ha/heis most pronouncedfor small values of the Pe and A
numbers. Additionally, for a given ha/hc the large initial
increase in recovery factor with increasesof Pe and A is
followed by a more gradual increase. Hence the recovery
factor is most sensitiveto a changein the parametersha/he,
Pe, and A at small values of Pe and A.
Production temperatures: The temperatureof the water
extracted during the production period of the first and fifth
Reference case: The following parameters are used in
the reference case: aquifer thickness, 50 m; injection vol-
ume, 60,000m3; an annualcycle;no dispersion;
and an
infinitely thick cap rock.
These parameters are summarizedin Table 1, which also
shows the values of the dimensionlessgroupsformed from
these parameters: Pe = 39.6, A = 396.4, ha/hc= 1, Ca/Cc =
1.25.
The thermalvolumeis V = 98,000m3. Unlessotherwise
noted, reference case parameters are always used.
The recovery factors for the first five cycles are (1) s =
0.77, (2) s = 0.81, (3) e = 0.83, (4) e = 0.84, and (5) e = 0.85.
Figure 7 showsthe temperature of extracted water during
the production period.
As an illustration of the size of the reference case, an
average single family house in the United States requires
about 15,000 kWh of energy for heatingduring 1 year. Let us
assumethat the storage system operateswith a temperature
difference of 25øC. If the recovery factor is about 0.8, and f
half of the total energy requirement is met with stored
energy, the reference aquifer should sufiScefor 180 houses.
Volume: The size of the heated region in the aquifer is a
fundamentalparameterof the storagesystem.Figure 8 gives
the recovery factor as a function of cycle numberfor several
thermal volumes with the shapeof the thermal volume kept
the samein all cases,a cylinder with an aspectratio, H/R, of
1. The amount of stored heat is proportionalto the thermal
volume; the heat losses to surroundingmaterial take place
through the surfate of this volume. Therefore the relative
heat loss is roughly proportional to the surfaceto volume
ratio, which is 2/R + 2/H. This ratio decreases and thus
becomes more favorable
as the thermal volume increases.
The temperature of the extracted water during the first
production period is also shown in Figure 8. In the case
whereV = 3,100,000m3theinitialslopeof thecurvereflects
the vertical heat loss to the confining layers. The faster
decreasein the temperature curve after 60 days is due to the
radial heat loss in the aquifer. When the volume becomes
much smaller, the radial heat loss affectsthe temperatureof
the extracted water during the whole productionperiod.
Shape: Although the recovery factor increaseswith the
thermal volume, there is one optimal aspectratio, or shape,
for which the recovery factor attains a maximal value for
each volume. Figure 9 showsrecovery factor as a function of
aspect ratio for the first five cycles. The curves have a rather
fiat maximum at an aspect ratio of 1.5. The first cycle
production temperature for different values of the aspect
ratio are alsogiven in Figure 9. If the thermalconductivityof
DOUGHTYET AL..' AQUIFER THERMALBEHAVIOR
I
[
I
[
o
.
I I i,I '1
,
I
,
I
577
578
DOUGHTY ET AL.' AQUIFER THERMAL BEHAVIOR
o
the confining layers is decreased, the maximal recovery
factor is obtainedfor a somewhatsmalleraspectratio. This
is shown in Figure 10. The differencebetween the recovery
factors in these cases is, of course, most pronounced for
small aspectratios, when the area facingthe confininglayers
is relatively large.
An approximate analytical expression for the optimal
aspect ratio, derived in Appendix B, is given by
HIR = 2 [(xq)f/XaCa]
1/2
(28)
[
(29)
where
2
(•C)f
= l/(XaCa)l/2
+ l/(XcCc)l/2
Thus the variation of the optimal aspectratio with thermal
properties is slow.
In practice, the choice of aquifer may be limited, thus
fixing the aquiferthicknessand the thermalpropertiesof the
system. Also, seasonal energy supply and demand may
determine the length of cycle time periods. In this case, the
volume of injectedwater is the only parameterwith which to
optimize the recovery factor. As the thermal volume becomes larger, the recovery factor initially increasesrapidly,
then levels off. The initial rapid increaseoccursbefore the
thermal radius attains the value which yields the optimal
aspect ratio for that aquifer thickness.
When the volume of the injectedwater is limited and small
compared to the thicknessof the aquifer, it may be advantageousto use a well that penetratesonly a part of the aquifer
in order to get a more compact shapeof the heated region.
However, the use of partial penetration may lead to increasedmixing of hot and cold water in the aquifer, which
will lower the recovery factor.
If the quantityof energyto inject is fixed, onepossibilityis
to inject a larger volume of water at a lower temperature.
This will increase the recovery factor, provided that the
energy content is calculated using the original ambient
temperature as a reference. It may reduce problemsrelated
to thermal stratification and water chemistry. However, a
decrease
in energyquality(temperature)
is unacceptable
in
many applications.
Aquifer thermal conductivity: Heat loss in the aquifer is
due to the mixing of cold and warm water (dispersion,
discussedin section 3.2.1), and heat conduction.The recovery factor for the first five cycles is given for a number of
different thermal conductivityvaluesin Figure 11. Figure 11
also shows the correspondingfirst cycle production temperature and the recovery factor as a function of the thermal
conductivity for two caseswith an aspectratio of 1.
The rangeof thermal conductivity shownin theseplots, up
to )ta = 20 J/m s K is larger than would be likely to be
measured on a laboratory sample of aquifer material. However, as will be discussedin the next section,dispersionmay
contribute to a large effective thermal conductivity in the
aquifer.
Figure 11 shows dependenceon aquifer thermal conductivity for cases with aspect ratio, H/R, near the optimal
value. For a system with H/R << 1, vertical heat losseswill
dominate, and recovery factor will be nearly independentof
aquifer thermal conductivity.
DOUGHTY ET AL' AQUIFER THERMAL BEHAVIOR
579
REFERENCE
CASE
PRODUCTION
TEMPERATURE
1,0
0.8
cycle
0,6
4
5
Ii
0,4
0,2
0
,
20
40
60
80
Time (days)
Fig. 7. Productiontemperature
versustimefor thefirstfivecycles
of the reference
case.
3.2. Dependenceon Factors Not Accountedfor in
the Dimensionless Formulation
This section deals with several factors which influence the
behavior of an ATES system, which are not included in
equations (7) and (8).
3.2.1. Velocity-dependentdispersion. During periods
whenthe water flowsthroughthe aquiferthereis, in addition
to ordinaryheat conduction,a dispersionof heat due to the
velocitydistributionacrosseachflow channel,the irregular-
ity of the pore system,and large-scaleaquiferheterogeneities. Accordingto the theoryfor dispersionof a nonadsorbenttracerin uniformporousmedia[Scheidegger,1960],the
dispersion
is proportional
to Ivlm, wherevv is the pore
õ
velocity.The valueof m rangesfrom 1 to 2. When m is 1, the
moleculartransversediffusionbetweenadjacentstreamlines
can be neglected,and when m is 2, the transversediffusionis
important.The transverse
diffusionbecomes
moreimportant
as pore velocity decreases.The thermalconductivityof the
stagnantliquid-solidmixtureandthe heatdispersionmay be
combinedto form an effectivethermalconductivity.GenerTHERMAL
VOLUME
VARIATION
A RecoveryFactor
B Production
Temperature
1,0
•
1,0
0,8
0,8
0,6
O•.•',
0,4 •
0,2
o
A
B
C
,
I
2
;5
4
O
oCycl•
First
5
20
Cycle
R(m)
H(m)Pe fl.
A 100 I00
Fig. 8.
40
60
80
Time (days)
V(m
3)
634 1585 3,100,000
B
50
50
159 396
390,000
C
20
20
25,4 63,4
25,000
D
I0
I0
6,3
15,8
3,100
Vw(m
3)
1,900,000
240,000
15,000
1,900
Recovery factor for different thermal volumes for the
first five cyclesand first cycle productiontemperatureversustime
for different thermal volumes.
580
DOUGHTYET AL.: AQUIFERTHERMALBEHAVIOR
ASPECT RATIO VARIATION
i
i
i
i
i
AQUIFER
i
THERMAL
1,0,
Cycle
0,8
VARIATION
i
i
2.5
5,
0,8
I0,
15,
0,6•
V: I00,000 m•
1,0
, , X'
o,
0.8•
0,6
0,4
CONDUCTIVITY
i
20.
0.4
0,4
0.2
0.2
0,2
H/R
,iI
.25 .5
-i
I
2
I
4
F•rst
Cycle
I0
•
o
i
Log H/R
0
' 2'o' 4'o ' 6• ' 8•0
Time (doys)
Cycle
FIRST CYCLE PRODUCTION TEMPERATURES
H/R_<!
H/R>!
i
i
t
H R (m)
1.0
0.8
0.8
•
0.5
50,50
0,6
0.2
""•
'"
'"
'"
0.2
Fifth cycle
First cycle
0,2
o
2b
n'o
6b
eb
o
2o
4o
6o
so
Time (days)
Fig. 9. Recoveryfactor as a functionof aspectratio for the first
five cycles and first cycle productiontemperatureversustime for
20,20
0
0
i
5
i
I0
•
i
15
i
20
Xo(m--•)
Fig. 11. Recovery factor for different values of the aquifer
different aspect ratios.
thermal conductivityfor the first five cycles,first cycle production
temperature versus time for different values of the aquifer thermal
ally, the effectivethermalconductivity
is a tensor.However, conductivity, and first and fifth cycle recoveryfactorsas a function
we will assumethat for our meshdesignthe off-diagonal of the aquiferthermalconductivity.
elementsarezeroandtheeffectivethermalconductivity
has
differentvaluesparallelto andperpendicular
to the direction
of fluidflow. The effectivethermalconductivity
is writtenin
termsof the Darcyvelocity,whichis theproductof thepore
velocity and the aquiferporosity.When rn = 1, we have
•a,,--ha-• d,lIolCw
(30)
•a -- ha-• d-LIolCw
First,we consider
thecasewhere•a is linearlydependent
and when rn = 2,
•a,,--ha-• dll*Iv2Cw
(31)
•a -- ha-• d_L*Iol2 Cw
RECOVERY FACTORvs ASPECT RATIO
ConfiningLoyerConductivity
Voriotion
3,5
O,4-
First Cycle
V= IO0,000m •
0.2-
0
H/R
_
.I
.25
,5
I
2
4
I0
I
I
,
I
•
I
•
on v (equation (30)). Figure 12 shows recovery factor and
production temperature as d• is varied. The dependenceof
the recovery factor on d• appearsto be very similar to that
for X (Figure 11). The thermal conductivity used to obtain
Figure 11 can be consideredto be a scalar effective conductivity of the form
•a = ha Jr-•d
I'0
i' , , , , , , _
-
The parameters d[[, d-L, d]•*, d-L* are properties of the
aquifer; d• and d-Lare often referred to as dispersionlengths.
Laboratory experimentshave shown that d-Lis an order of
magnitude smaller than d• for most field samples.We have
found that the recovery factor and production temperature
exhibit only a weak dependenceon the value of d-Land d-L*.
Log H/R
Fig. 10. First cycle recoveryfactor as a functionof aspectratio for
different values of confininglayer conductivity.
(32)
where ha is a constantadditionto the thermal conductivityto
account for dispersion. To examine further the effect of the
velocitydependence
in •a, Figure13 showsthe firstcycle
production temperature for three cases which all have a
recovery factor of 0.60 but which use three different formulas to describethe dispersion(equations(30), (31), and (32)).
Clearly, the productiontemperaturecurvesare very similar.
For later cycles the recovery factors and production temperature curves begin to diverge slightly. The similarity
among the curves in Figure 13 implies that a constant
effective conductivity of the form of (32) could be used to
analyze a system's thermal behavior instead of the more
complicated velocity dependentform of (30) and (31) if a
general equation relating the different forms could be found.
This would have the advantageof allowing the incorporation
of dispersioninto the parametersPe and ha/hc.
DOUGHTY ET AL.' AQUIFER THERMAL BEHAVIOR
PRODUCTION
In the field, the tracer dispersionlengths may be determined by a tracer injection test. Experimental data show that
tracer and heat dispersion lengths (d•, dñ) are practically
TEMPERATURES
FOR VARIOUS DISPERSION
FORMULATIONS
identical for flow in a uniform material [Bear, 1972]. Hence
the correlation needed is that between d• and dñ and Xd. For
,
I
A
,6
the parameter ranges consideredin this study, we obtained
the empirical relation
•a= ha+ 0.3d•RCa
ti
581
0,8-
-
(33)
between the scalar effective thermal conductivity, which is
used for the whole ATES cycle, and the dispersionlength d•.
Equation (33) was generatedby runningthe SFM for a range
of the parameters R, H, ti, ha, Ca, using (30) and (32) to
describe dispersion effects and correlating cases with equal
first cycle recovery factors. The ratio d•/d• = 10 was used.
Different values of d• do not change(33) appreciably.
In analyzing the thermal behavior of an ATES system, (33)
0.6
_c•
0,4- ...,
o,2
B, Xo=
I vl Cw
C. =2,5+0,22-107[v[
maybe usedto determine
•'afromd•, R, Ca,andti. Then
may be used instead of ha in Pe and ha/hcto account for
dispersion.
The effect of the addition of dispersion to the thermal
conductivity may be observedby examiningFigures 5 and 6.
For example, considera casewith Pe = 50, A = 50, and
Xc = 1; an increase of effective thermal conductivity by a
factor 10 to include dispersion changes the value of the
J
- A,)to= 10,2msK
0
•
20
•
40
•
60
80
Time (doys)
Fig. 13. First cycle productiontemperatureversustime for different dispersion formulations.
parameters
to Pe = 5, and•a/kc= 10.Thisreduces
thefirst duction, and rest periods. Figure 14 showsthe first cycle
cycle recovery factor from 0.67 to 0.40. The decrease in
production temperaturesis also substantial.
3.2.2. Unequal lengthperiods. The recovery factor depends on the length of the cycle, tc, and the intracycle
scheduling: the relative duration of injection, storage, proDISPERSION
LENGTH
production,andrest periods.CaseC consistsof a hypothetical cycle with instantaneous
injectionand production.Here,
the durations of the injection and productionperiods are
zero, while the storageandrest periodsare eacha half cycle.
For a given tc, the recovery factor is higher for a shorter
VARIATION
I,O
0.8
d"O.
0,$
recovery factor as a function of tc for different schedulesfor
three thermal volumes.In caseA the fluid is injectedduring
the first half of the cycle and produced during the second
half. There is no storageor rest period. Case B refers to a
cycle which is subdivided into equal injection, storage,
storage period, ts.
The resultsof section3.1 were calculatedfor a cycle with
equallengths,ti, of the individualperiods.In order to usethe
0,6
0,4
RECOVERY
0,2
I,O
0,;:'
First Cycle
2
3
4
0
5
i
i
20
t
i
i
4•0 60
i
FACTOR
i
vs LENGTH
OF CYCLE
i
i
i
I
400
I
500
8•0
Time (doys)
Cycle
•
0,8
0,6
0,6
0.4
0,4
First cycle
o,2
0
0,2
J
Xo=2,5(•--•)
b- ti=ts=tp=tr='•-
dx=du
c- ts=tr=•- ti=tp=0
t•
5
I0
15
d. (m)
Fig. 12. Recovery factor for different values of the dispersion
length for the first five cycles, first cycle production temperature
versustime for different values of the dispersionlength, and first and
fifth cycle recovery factors as a function of dispersionlength.
O0
I
I00
I
200
300
I
600
tc (doys)
Fig. 14. First cycle recovery factor versus cycle length for
several injection-storage-production-restschedulesfor three thermal volumes.
582
DOUGHTY ET AL' AQUIFER THERMAL BEHAVIOR
FINITE THICKNESS CAPROCK EFFECT ON RECOVERY FACTOR
0,8
--(•a/•c)(•cti/Cc)
d
Becauseof the small interdependenceof radial and vertical
losses,D = o• is used for the plots in section3.1, and the
effect of a finite cap rock thickness is discussedhere
i•{ 1,0
0.6
0.5
0,4
separately.
0,2
0,2
(35)
0
0'1
d
0,8
1,0
0,6
0,4
0,5
0,2
0,2
1,0
0,8
The new parameterused to describethe finite caprock
effectis d = D/H. Figure 15showsthe recoveryfactorfor an
aquifer with a finite cap rock relative to that with an infinite
cap rock, s(d)/s(o•),as a functionof 1/A. The figuredisplays
resultsfor variousvaluesof d, for the first and fifth cycles,
and for Xa/Xc= 1, 2, and 10. Hence the recoveryfactorfor
the finite cap rock casescan be calculatedby usingthe
appropriate figure which gives the recovery factor for the
infinitecaprock (figure5) andthen multiplyingthe resultby
the number obtainedfrom Figure 15.
Severalobservations
may be madeaboutFigure15. For a
relatively thick cap rock (d > .5), vertical heat loss is
primarilydeterminedby the aquitardthermalproperties,but
d
0,6
LONG TERM
BEHAVIOR
Recovery Factor
0,4
0,2
0,8
0,1
0,2
0.3
0,4
0,5
I
_
0,6
Fig. 15. Finite thicknesscap rock effectfor the first and fifth cycle,
recovery factors for X,/Xc = 1, 2, and 10.
section3.1 figuresfor caseswhere the differentperiodshave
unequal durations, a time measure more suitable than ti is
required. An adequatemeasureis the time an injectedwater
particle resides in the aquifer, averagedover all particles in
the injection volume. This time, called r, is given by
r = «(ti+ tp)+ ts
H,R - 20,20 (m)
0,2
D =2m
I
i
1
i0
2o
3o
Cycle
Temperature Field
]
(34)
,
,
TI
0,1
provided that the flow is constant during the injection and
production periods. To allow for any scheduling,the time ti
is replacedby r/2 in theformulas
for Pe, A, •a, andthefinite
thickness cap rock effect (equations (21), (33), (36), and
(37)). The error made by usingthe recovery factor curves of
section 3.1 with r/2 is small. According to Figure 14, the
largest departurefrom case B causesa changein recovery
factor of less than 5%. Variation in schedulingaffects the
production temperature more strongly. When there is no
storage period, the production temperature begins at T•;
with a storage period, the production temperaturebegins
below T but decreasesmore slowly.
When the aquifer and confining layer thermal properties
are the same, the first cycle recovery factor for caseC can be
determined analytically. The recovery factor is given by the
formula for the mean temperature decline in a cylinder,
which is shownin Appendix C.
3.2.3. Finite thicknesscap rock. The vertical heat loss
from a shallow aquifer system may be substantiallyenhanceddue to the influenceof the groundsurface,actingas a
constant temperature boundary. A cap rock with thickness
D can be accounted
for in the dimensionless
formulation
4O
• 20 - -/9-•.......................
rr
,
I0t•''
I
I
I
i0
20
30
Cycle
Fig. 16. Recoveryfactor for 33 cyclesand the radial distanceto
discussedin section 2.3, by forming another dimensionless severalisothermsat the end of the first 33 cyclesfor infinitelythick
parameter:
cap rock and thin cap rock cases.
DOUGHTYET AL.' AQUIFERTHERMALBEHAVIOR
there is an increase in the effect of ha as d decreases.The
thicknessof the aquifer,H, appearsboth in d and in 1/A. The
parameter
d decreases
asl/H, and1/Adecreases
asl/H2;an
increasein H yields an increasein e(d)/e(m).The ratio e(d)/
e(oo)initially decreases as the number of cycles increases
then reachessteady state. The numberof cyclesrequiredto
reach steady state increasesas D increases.
For the parameterrangesconsideredin this study, if either
of the followingextreme conditionsholds,then the finite cap
rock effect will be small (lessthan a 5% reductionof the fifth
cycle recovery factor):
583
During each cycle, energyis lost from the injectedwater
to the surroundings.If the aquifer may be usedfor purposes
other than heat storage,the thermalpollutioncausedby the
residualheatmaybe a problem.Theradialdistance
fromthe
well to a certain isothermat the end of each storagecycle is
given in Figure 16 for the infiniteand finite cap rock cases.
The isothermsgive a roughoutlineof the radial distribution
of the residual heat just before the injection period of the
next cycle starts. As expected, in the thin cap rock case,
with heat lost to the surface,the radialextentof the residual
heat is lessened and can be seen approaching a constant
value.
Condition
1
The radial distance at which the temperature increase
d > 10
A > 2
(36)
equals10% of the temperaturedifferencebetweeninjected
and ambient water has not reached farther than 47 m after 33
Condition
2
A > 300
any d
(37)
cyclesfor the infinite cap rock case. At this time the rate at
which this distanceincreasesis only 0.4 m/yr. The residual
heat is spreadingout in the aquiferat a slowrate. However,
Whenthefirstcondition
holds,
thethickness
ofthecaprock in this examplethe thermal conductivityin the aquifer is
is so large that the effect of the surfaceis not appreciablyfelt
by the aquifer. When the secondconditionholds,the aquifer
is thick enough that the heat loss to the ground surface
results in only a minor change in the aquifer heat content.
only 2.5 J/m s K, whichis a ratherlow valueif dispersionis
significant.Of course, in these calculationswe have not
considered
thepossible
presence
offaultsor otherheterogeneitiesthat may serveas specialchannelsfor fluid flow nor
For intermediate
cases,Figures15providesa goodindica- has the effect of regionalflow been considered.
tion of the size of the finite cap rock effect. When the aquifer
and confining layer thermal properties are the same, the
meantemperaturedeclineformulafor a cylinderwith a finite
cap rock can be used to obtain an analytical solutionfor the
finite thickness cap rock effect on recovery factor for the
first cycle. This is shown in Appendix D.
3.2.4. Long-term effects. The transient energy loss for
each cycle decreasesas the number of cycles increases,so
that the recovery factor improves. Figure 16 shows the
recovery factors for 33 annual cycles for an aquifer with an
infinitely thick cap rock and with a relatively thin cap rock (d
= D/H = 0.1). The increase in recovery factor is quite rapid
during the first few cycles for both cases, but the thin cap
rock case approachessteady state faster.
I.O
I
I
I
I
I
I
i
I
_
3.2.5. Multiple layers with differentflow rates. Buoyancy flow inducedby the densitydifferencebetweenhot and
cold water is discussed
briefly in AppendixA. While the
steadyflow modeldoesnot solvethe coupledheat and mass
transfer equations,the final effect of buoyancyflow is a
tilted thermal front. When H/R < 4, the surface area of the
hot regionis increasedwhenthe thermalfront is tilted. If the
tilting is extreme,duringproduction,unheatednativewater
in the lower part of the aquiferwill be recoveredalongwith
heated water. Both these effects lower the recovery factor,
so it is desirable to examine them. A tilted front can be
generatedroughlyby multiplelayerswith differentspecified
horizontal
flowrates
within the SFM.
As an example of this technique,we considera threelayered case, where the flow rate during injection is enhancedin the upperlayer and reducedin the lower layer.
The flow rate during productionis the same in all layers.
Figure17shows
theproduction
temperature
forthreecases
0.8
with different combinations of flow rates in the layers, as
indicatedby the accompanyingfigures.The aquifer has a
0,6
thickness of 20 m. The total flow is the same in all cases.
Curve A representsthe casewhenthe flow rate is the same
in the three layers. The recoveryfactor is 0.65. In caseC,
where the tilting angleafter the injectionperiod is closeto
0.4
0.2
45ø, the recovery factor is reduced to 0.58.
H,R =20,20 (m)
The permeabilityof the aquifermustbe quitelow in order
that the tilting angledoesnot exceed45øduringthe storage
i
I
20
i
I
I
40
6
cycle [Hellstr6rn et al., 1979].
i0 I, 80I
Time (days)
•=0,65
4.
•=0,62
•=0.58
Q
1,8Q
O
Q
Q
A
Q
0.5Q
B
0,2Q
COMPARISONS WITH
FIELD
EXPERIMENTS
The use of the results presentedin section 3 is illustrated
by comparing these results with the recovery factors and
production temperaturesfor two recent ATES field experiments.
4.1.
Auburn
C
Fig. 17. First cycle productiontemperatureversustime for the
different combination of flow rates indicates schematically in the
figure for an aquifer of thickness20 m.
The Water
Resources
Research
Institute
at Auburn
Uni-
versity, Auburn, Alabama, hasperformeda two-cycle ATES
field experimentusinga 21-m-thickaquiferwith an ambient
584
DOUGHTY ET AL.' AQUIFER THERMAL BEHAVIOR
temperatureof 20øC[Molz et al., 1979, 1981].Althoughthe
top of the storageaquiferlies40 m belowthe groundsurface,
it is overlain by a 9-m-thick clay layer, above which is a
shallow aquifer, whose temperatureremained constant
throughoutthe experiment. The injection/production
well
penetrated the middle 9 m of the aquifer.
During each 6-month cycle, 55,000 m3 of water was
injected at 55øC,stored, and recovered.The recoveryfactors for the two cycles were 66 and 76%.
Independentmodelingwork [Tsanget al., 1980;Buscheck
et al., 1982], usinga coupledheat and masstransfernumerical model, has matched the simulated and observed aquifer
temperature fields at the end of the first injectionperiod. It
prescribedin the radial direction is used to run numerical
simulationsof ATES cycles. The net thermal behavior of the
storage system is given by the recovery factor and the timedependent temperature of the recovered water. The recovery factor is defined as the quotient between recoveredand
injected energy; the energy is calculated using the initial
groundwater temperature as a reference. The effects of
buoyancy flow are neglected in this study, thus the results
shouldbe applicableto a systemwith either a low permeability or vertically stratified aquifer, small temperaturedifference between injected and ambient water, short cycle
length, or small aquifer thickness. A series of graphs are
presented to show recovery factor and production tempera-
groups:
Pe, A, Xa/
indicatesthat an effectiveconductivity
•a = 2Xais appropri- tureasa functionof thekeydimensionless
ate for this experiment. A summary of the parameters
needed is shown in Table
1. Also shown are the dimension-
less groups used to predict recovery factor: Pe = 78.4, A =
80.3, ha/hc = 2, and d = 0.4.
The first cycle finite cap rock effect, shownin Figure 15, is
negligible.Hence the first cycle recovery factor is estimated
from Figure 5 to be 0.71. One of the effectsof the partially
penetrating injection/production well may be shown the
following way. From temperatureobservationsmade in the
aquifer, the thermal radius is inferred to have reachedabout
43 m, which implies the effectivethicknessof the aquifer was
16 m. Using R = 43 and H = 16givesPe = 101.9,A = 46.6, d
= 0.5. Using these valuesin Figures 15 and 5 yields e = 0.68.
The decreasefrom 0.71 to 0.68 is due to the worseningof the
aspect ratio. Experimental observationsalso indicated some
thermal front tilting. Accordingto equation(A1), the characteristic tilting time, to, is 28 days, to be comparedwith a •' of
110 days. Hence the observed first cycle recovery factor,
0.66, may be lower than that predicted due to buoyancy
Xc, and d. The definitions and significancesof these groups
are summarized in Table 1. The graphs are designed to
provide a quick method for characterizationof potential
ATES
sites.
Some specific conclusionsare summarizedbelow:
1. The storagevolume is a fundamentalparameterof the
system. The relative heat loss is roughly proportionalto the
surfaceto volume ratio. Even for an optimal aspectratio for
our reference case, a minimal injection volume of about
50,000m3 is requiredin orderto obtaina goodrecovery
factor (0.7) during the first annual cycle. For other cases,
similar minimal volumes can be estimated using the figures
in section3. The aquifer storageconceptfor a seasonalcycle
must be applied on a large scale in order to ensure a high
recovery.
2. The thermal behavior of the storage system is very
sensitive to the value of the effective thermal conductivity in
the aquifer. The effective thermal conductivity, which is
defined to include a contribution from dispersion,may be
quite large. The constantform of Xa may be usedin the
formulasfor Pe and•a/•.c.
flow.
3. After an initial transient period, the heat loss through
the upper surface of the aquifer is determined by the
The Bureau de Recherches G6ologiques et Mini/•res,
thicknessof the cap rock. The transientperiod may be quite
Orleans, France, has conducteda four-cycle experiment on
long. The effect of a finite cap rock on the fifth cycle
a 2.5-m-thick aquifer with ambient temperature 12.5øC at
recovery factor is less than 5% if d > 10 and A > 2. Further,
Bonnaud, France [Sauty et al., 1980b]. During each 12-day
the influence of the finite cap rock loss decreasesas the
cycle,490m3 of waterwasinjectedat 35øCandrecovered.
thickness of the aquifer increases.
The fourth cycle recovery factor was 67.7%. A summaryof
4. The cycle is divided into periodsof injection, storage,
neededparametersis shownin Table 1. Also shownare the
production, and rest in accordance with the supply and
dimensionlessgroupsused to predict recovery factor: Pe =
demand of heat. The relative duration of these periods
15.4,A = 25.2,•a/•.c= 13,andd = 1.6.Thecurvesfor •a/•.c
appears to have a fairly small influence on the recovery
= 10 are usedto predict the recoveryfactor. The fifth cycle
factor when the length of the cycle is constant.A parameter
cap rock effect, shown in Figure 15, is about 0.99. From
ß = (ti + tp)/2+ tsis definedasthe averagelengthof timean
Figure 5, the fifth cycle infinite cap rock recovery factor is
injected fluid particle spendsin the aquifer. For caseswith
0.65. Combiningthesetwo yields e = 0.64. The fourth cycle
periods of unequallength •/2 may be usedin place of ti in the
recovery factor would be slightly smaller. This underpreformulasfor Pe, A, and •a.
dicts the observedfourth cycle recoveryfactor of 0.677. The
5. The recovery factor increases with the number of
difference is due, at least in part, to the scheduling, as
storage
cycles. During each cycle heat is lost to the cold
discussed in section 3.2.2.
surroundings.The increasing amount of residual heat improves the performance of the storage systemsduring the
next cycle. However, the residual heat may be in conflict
5. SUMMARY AND CONCLUSIONS
with other usesof the aquifer. The thermal pollutionappears
The thermal behavior of an aquifer thermal energy storage to spreadout at a rather slow rate, when there is no regional
system which consistsof a single well in a laterally infinite flow in the aquifer.
6. The shapeof the heatedvolume shouldbe as compact
horizontal aquifer has been studied. The heat transfer equations are presentedin dimensionless
form in orderto identify as possible in order to minimize the heat loss. The recovery
the set of parameters which define the system's thermal factor has a rather flat maximum when the radius of the
heated volume is about 2/3 of the aquifer thickness.
behavior. A numerical model in which a steady fluid flow is
4.2
Bonnaud
DOUGHTY ET AL.'. AQUIFER THERMAL BEHAVIOR
APPENDIX
A: BUOYANCY
TILTING
The density differencebetween hot and cold water induces
a buoyancy flow of the hot water toward the upper part of
the aquifer. An order of magnitudeestimateof the tilting rate
of a more or less vertical front may be deduced from the
following idealized case. Consider a sharp, vertical thermal
front separating regions of temperature T• and To in a
confined aquifer of thickness H. The aquifer layer has an
infinite horizontal extension. The vertical permeability k'
may differ from the horizontalpermeabilityk. The tilting rate
of the front is given by a characteristictilting time to, which
multipliersis usedto find the aspectratio which maximizese
for a given volume. This optimal aspectratio is
H_2((hC)f•'/2
•-- •)kaea/
H
conductivity
•xa in theaquifer(section
3.2.1),theexpression
for the optimal aspectratio becomes
H/R = 2 ()ka/Xa)
1/2
Ca (Ixo+
(A1)
where It is fluid viscosityand p is fluid density.
The buoyancytilting of an initiallyverticalfront duringa
time to is about 60ø. The most important factors are the
temperature levels, which determine /x and p, and the
permeability. If the time of the cycle is smallerthan to, then
the tilting is expected to be moderate.
The time constantto was derivedfor the plane case, but
the magnitudedoes not changeappreciablyfor the radial
case. If the thermal front is diffuserather than sharp, the
tilting rate is slightly lowered. Furthermore, as the thermal
front tilts, the flow resistancein the hot part of the aquiferis
reducedbecauseof the lower viscosityof hot water. Forced
convection,then, givesan increaseof the tilting rate during
injectionperiodsand a decreaseduringproductionperiods
for hot water storage. For chilled water storage, forced
convection decreasesthe tilting rate during injection and
increasesit during production.
The derivation of to and further developmentsincluding
the forced convectioncontributionto tilting are given by
Hellstr6m et al. [1979]. Criteria for minimal tilting for
specificexamplesare also given.
APPENDIX
B: VARIATION
(B3)
For the common situation in which Ca and Cc are similar
and the laboratory measuredvalues of ha and hc are similar
but dispersion creates a large effective horizontal thermal
is
to=0.034
(kk,)l/•
Cw(Po
- P•)
585
APPENDIX
MEAN TEMPERATURE
(B4)
C:
DECLINE
IN A CYLINDER
Let us assumethat the storagevolume has the shapeof a
cylinder with radius R and height H. Initially, the temperature is T• throughoutthis volume and Toin the surrounding,
infinite medium. The whole system is assumed to have
homogeneousthermal properties, with thermal conductivity
h, heat capacity per unit volume C, and thermal diffusivity K
= h/C. The temperaturefield is governedby the ordinary
heat equation.
The mean temperaturein the storagevolume, Tin,at a time
t is given by a product solution [Claesson et al., 1980],
namely,
Tm(t)= To+ (T• - To)gm(Kt/R
2)fm(4•t/H2)
(Cl)
where the radial factor is
grn(X)= 1 - e-1/2X[Io(1/2x)
+ Ii(1/2x)]
and the vertical
(C2)
factor is
fm(Y)= err(1/V•) - (y/z')m (1 - e-vy) (C3)
Here I0 and I• are modified Bessel functions and erf is the
error function. The following approximationsmay be used:
OF THE OPTIMAL ASPECT RATIO
gm(X)• 1 - 2 (x/z')v2
Some insightinto the variation of the optimal aspectratio
with thermal propertiesmay be gained by consideringthe
following expression,which givesa roughapproximationfor
the first cycle recovery factor. The heat loss per unit area
across the plane interface between two semi-infinite media
originally at different temperaturesin a time t [Carslaw and
Jaeger, 1959] is multiplied by the surface area of a cylinder
of thicknessH and radiusR. This yields
x << 1
(C4)
MEAN TEMPERATURE DECLINE RECOVERY FACTOR
I
25
)to/)t.c=I Co/Cc=I
20
(B1)
15
where
()kC)f-'
,1/()kaqa)l/2
+ 1/()kcqc)l/2 (B2)
Thefactor((hC)f)
m is theharmonic
meanof theaquifer
and confininglayer (hC)v: values.The first term in the
expressionfor e in bracketsis proportionalto horizontalheat
loss, the second to vertical heat loss. This approximation
always underpredictsrecovery factor and is worse for small
recovery factors. For the reference case, it is 0.03 less than
the numerically simulated value of 0.77. For a recovery
factor of 0.48, it gives 0.34. The method of Lagrange
I
5
Fig. C1.
I
IO
I
15
2
I0
l
25
Analytical recovery factor from the temperaturedecline
of a cylinder.
586
DOUGHTYET AL.' AQUIFERTHERMALBEHAVIOR
NOTATION
gin(X)
• •XX1--
X)) 1
(C5)
fro(Y)• 1 - (y/z')m
y << 1
(C6)
volumetricheatcapacity,J/m3 K.
aquifervolumetric
heatcapacity,J/m3K, equal
fro(Y)
• (•y)•/2
1--
y>>1
(C7)
to &aCw + (1 - &a)Cr, where Cw and Cr are
water and rock volumetric heat capacities,respectively.
miTt,it surfaceareaof the (m, n)th meshelement,m2.
The recoveryfactorfor the firstcyclecanbe estimatedby
d
substitutingr for t in (C1), where the time constantr is
r-
(ti + tv)/2 + ts
equal to D/H.
first-order dispersion constants, dispersion
(C8)
The recovery factor becomes
dll*, d.*
t• = gm(n:r/R
2)f m(4Kr/H
2)
ratio of caprock thicknessto aquiferthickness,
D
(C9)
H
lengths, m.
second-orderdispersionconstants,s.
cap rock thickness, m.
aquifer thickness, m.
horizontal and vertical aquifer permeabilities,
Figure C1 showsthe recoveryfactor as a functionof the
m2/s.
parameters
R/(Kr)1/2= pe1/2
L
H/(trr/2)
1/2= •
M
APPENDIX D:
INFLUENCE OF A FINITE CAP ROCK
ON THE MEAN TEMPERATURE DECLINE
Q
Qi, Qp
fro', which givesthe verticaleffects,includesthe dependenceon the finitethicknessof the cap rock, that is,
frn'(d,Y)
=fm(y)+
X/-•y
ierfc
(2d+
1)
number of mesh elements in each row between
r -
Pe
The effect of a finite cap rock of thicknessD can be
analyzedusingthe methoddescribedin AppendixC andthe
superpositionprinciple[Claessonet al., 1980].The function
characteristic
length,m, equalto [(Xa/Xc)(Xcti/
Cc)]1/2.
F
0 and R.
Pecletnumber,equalto CaR2/2kati.
heatflow rate per unit area,J/m2 s.
volumetricfluid flow rate, m3/s.
injection and productionvolumetricfluid flow
rates,m3/s.
radial coordinate, m.
F•
dimensionlessr, equal to r/L.
R
thermal radius, m.
Rm
distanceto the inner edgeof the mth columnof
meshelements,m, equalto [(m - 1)/M]1/2R.
t
V•ieffc
X/-•y
]+ieffc
(D1)
2
dimensionless
time, equal to t/ti.
tc
lengthof onecycle,s, equalto ti + ts+ tp+ tr
where ti, ts, tp and tr are the lengthof the
injection,storage,production,andrestperiods,
whered = D/H, y = 4•r/H2,and• = X/C.Thefunction
fmis
given in AppendixC and ierfc denotesthe integralof the
complementary
error function.The functionfro(Y)givesthe
verticaleffectsfor an infinitecap rock. The functionfro'
replacesfm in (C9) in AppendixC. We get the following
to
At
T
expressionfor the recovery factor:
e = gm(•cr/R
2)fm'(D/H, 4Kr/H2)
(D2)
FigureD1 showse(d)/e(oo)
= fm'(d, Y)/fm(Y)as a function
of 1/A for several values of d.
time, s.
t'
respectively.
characteristictilting time, s.
time step for conduction,s.
temperature, K.
To
originalambienttemperature,K.
T•
injection temperature, K.
T'
dimensionless
temperature,equalto (T(r•-
Tv
•v
To)/
r0).
productiontemperature,K.
production
temperature
averaged
overproduction period, K.
MEAN TEMPERATURE
dimensionless
averageproductiontemperature,
DECLINE
FINITE THICKNESS CAPROCK EFFECT ON RECOVERY FACTOR
Xo/Xc=l
equal to e.
Tref
Co/Cc=l
reference(cut-off)temperature,K.
pore velocity, m/s.
1,0
steadyradial darcy velocity, m/s, equalto Q/
0,8
2 rrHr.
0,6
0,4
0.5
0,3
v
thermalvolume,m3, equalto rrR2H.
volumeof the (m, n)th meshelement,m3.
volumeof water injectedand produced,m3,
equalto Qiti = Qptp= (Ca/Cw)V.
z
vertical coordinate, m.
0,2
0,2
i 13,1
0
0:l
0,'2
013
0'.4
0',5
I
Fig. D1. Analyticalfinite thicknesscap rock effecton recovery
factor.
z'
Zn
•
dimensionless
z, equal to z/L.
verticaldistancefromthe top of the caprockto
the top of the nth row of meshelements,m.
recoveryfactor, ratio of producedto injected
energy,with energiesmeasuredrelativeto To.
DOUGHTYET AL..'AQUIFERTHERMALBEHAVIOR
Eref
recovery factor with energiesmeasuredrelative
to Tref.
thermaldiffusivity,m2/s,equalto h/C.
thermal conductivity, J/m s K.
effective aquifer thermal conductivity, of the
form•a = ha-{-dispersion
term,J/ms K, where
the dispersionterm may be ha const,in which
case•a,,= •al • •a, all
first-
ordervelocitydependent,
and
second-order
velocitydependent.
((}kC)f)
1/2 effective(hC)v2,harmonicmeanof aquiferand
confining
layervalues,
J/m2sv2K, equalto2/[1/
A
(haqa)
1/2-{-1/(hcqc)l/2].
Lambdanumber,equalto Ca2H2/hcCcti
.
viscosityof ambientandinjectedwater, kg/ms.
densityof ambient
andinjected
water,kg/m3.
porosity.
effective time a fluid particle spends in the
aquifer,s, equalto (ti + t•,)/2+ ts.
Subscripts
a
c
aquifer.
confininglayer.
radial.
vertical.
label for rnth column of elements in mesh.
label for nth row of elements in mesh.
label for the (rn, n)th elementin mesh(also used as
a superscript).
587
water flow, Theoretical Analysis and Computer Simulation of
Solid-Fluid Heat Storage Systemsin the Ground: Extraction of
Earth Heat, interim report part II, Dep. of Math. Phys. and Bldg.
Sci., Lund Inst. of Technol., Lund, Sweden, 1978.
Hellstr6m, G., C. F. Tsang, J. Claesson,Heat storagein aquifers:
Buoyancy flow and thermal stratification problems, report, Dep.
of Math. Phys., Lund Inst. of Technol., Lund, Sweden, 1979.
Iris, P., Exp6rimentation de stockagedes chaleur intersaisonnieren
nappe phr6atique Campuget (Gard) 1977-1978, report, Ecoles des
Mines de Paris Cent. d'Inf. G .., Fountainebleau, France, 1979.
Lawrevce Berkeley Laboratory, Proceedingsof Thermal Energy
Storagein Aquifers Workshop, Rep. LBL-8431, Berkeley, Calif.,
1978.
Mathey, B., Development and resorptionof a thermal disturbancein
a phreatic aquifer with natural convection, J. Hydrol., 34, 315333, 1977.
Molz, F. J., and L. C. Bell, Head gradient control in aquifers used
for fluid storage, Water Resour. Res., 13(4), 795-798, 1977.
Molz, F. J., A. P. Parr, P. F. Andersen, V. D. Lucido, and J. C.
Warman, Thermal energy storagein a confined aquifer: Experimental results, Water Resour. Res., 15(6), 1509-1514, 1979.
Molz, F. J., A. P. Parr, and P. F. Andersen, Thermal energy storage
in a confined aquifer: Second cycle, Water Resour. Res., 17(3),
611-645, 1981.
Papadopulos, S.S., and S. P. Larson, Aquifer storage of heated
water, II, Numerical simulation of field results, Ground Water,
16(4), 242-248, 1978.
Reddell, D. L., R. R. Davison, and W. B. Harris, Thermal storageof
cold water in groundwateraquifersfor coolingpurposes,Proceedings of Thermal Energy Storage in Aquifers Workshop, Rep.
LBL-8431, pp. 51-55, Lawrence Berkeley Lab., Berkeley, Calif.,
1978.
Sauty, J.P., A. C. Gringarten, A. Menjoz, and P. A. Landel,
Sensibleenergy storagein aquifers, 1, Theoreticalstudy, report,
Bur. de Rech. G6ol. et Minib•res, Serv. Geol. Natl., Orleans,
France, 1980a.
Sauty, J.P., A. C. Gringarten,H. Fabris,P. Thiery, A. Menjoz, and
P. A. Landel, Sensible energy storage in aquifers, 2, Field
experimentsand modelingresults,report, Bur. de Rech. G6ol. et
Acknowledgments. This work was supportedby the Assistant
Secretaryof Conservationand Solar Energy, Office of Advanced
Minib•res, Serv. Geol. Natl., Orleans, France, 1980b.
ConservationTechnology,Division of Thermal and Mechanical Scheidegger,A., The Physicsof Flow throughPorousMedia, 2nd
StorageSystemsof the U.S. Departmentof Energy undercontract
rev. ed., pp. 256-259, MacMillan, New York, 1960.
W-7405-ENG-48.This reportrepresentswork performedwithin the Tsang, C. F., and P. A. Witherspoon,An investigationof screening
seasonalthermalenergystorageprogrammanagedby PacificNorthgeothermalproductionwells from effectsof reinjection,Geotherwest Laboratory.On the Swedishside,the work hasbeensupported
mal Reservoir Engineering, Stanford Geotherrn. Program Rep.
by the National SwedishBoard for Energy SourceDevelopment
SGP-TR-12, pp. 62-64, StanfordUniv., Stanford,Calif., 1975.
(NE).
T sang, C. F., M. J. Lippmann, C. B. Goranson, and P. A.
Witherspoon,Numerical modelingof cyclic storageof hot water
in aquifers,Rep. LBL-5929, LawrenceBerkeley Lab., Berkeley,
REFERENCES
Calif., 1977.
Bear, J., Dynamics of Fluids in Porous Media, pp. 650-651,
Elsevier, New York, 1972.
Buscheck,T., C. Doughty,andC. F. Tsang,Aquiferthermalenergy
storage-Parameterstudy, report, in press, Lawrence Berkeley
Lab., Berkeley, Calif., 1982.
Carslaw, H. S., and J. C. Jaeger,Conductionof Heat in Solids,2nd
ed., pp. 87-89, Clarendon, Oxford, 1959.
Claesson,J., B. Eftring, andG. Hellstr6m,Temperaturedeclineof a
heated region in the ground, Rep. Lund-MPh-80122,Dep. of
Math. Phys., Lund Univ., Lund, Sweden, 1980.
Fabris, H., and A. C. Gringarten,Etudesexp6rimentales
de stockage et de transfert de chaleur en aquifb,
re, Summary of the
PrincipalScientificandTechnicalResultsof the NationalGeological Servicefor 1977,supplementto the bulletinof the Bureaude
RecherchesG6ologiqueset MiniSres, National GeologicalService, Orleans, France, 1977.
Fabris, H., A. C. Gringarten,P. A. Landel, M. L. Noyer, and J.P.
Sauty, Etude des possibilit6sdu stockaged'eau chaudeen aquifb•reprofond, ler rapport d'avancement,Rapp. BRGM 77 SGN
Tsang, C. F., T. Buscheck, D. Mangold, and M. Lippmann,
Mathematical modeling of thermal energy storage in aquifers,
Proceedingsof Thermal Energy Storagein AquifersWorkshop,
Rep. LBL-8431, pp. 37-45, LawrenceBerkeley Lab., Berkeley,
Calif., 1978.
Tsang, C. F., D. Hopkins, and G. Hellstr6m, Aquifer thermal
energystorage•A survey,Rep. LBL-10441,LawrenceBerkeley
Lab., Berkeley, Calif., 1980.
Tsang, C. F., T. Buscheck,and C. Doughty, ATES-•A numerical
simulation of Auburn University field experiments, Water Resour. Res., 17(3), 647-658, 1981.
Whitehead, W. R., and E. J. Langhetee,Use of boundingwells to
counteract the effects of preexisting groundwater movement,
Water Resour. Res., 14(2), 273-280, 1978.
Yokoyama, T., H. Unemiya, T. Teraoka, H. Watanabe,K. Katsuragi, and K. Kasamaru, Seasonalregenerationthrough undergroundstrata,Proceedingsof ThermalEnergyStoragein Aquifers
Workshop,Rep. LBL-8431, Lawrence Berkeley Lab., Berkeley,
Calif., 1978.
658 HYD, Bur. de Rech. Geol. et Minieres, Orleans, France,
1977.
Hellstr6m, G., and J. Claesson,Heat lossesand temperaturefields
for heatstorageaquifers:A computational
modelwith a simplified
(Received June 22, 1981;
revised January 10, 1982;
acceptedJanuary20, 1982.)