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.)