The main purpose of this paper is to present a robust forward model for simulating extraction and storage of thermal energy in an aquifer. The model is a local three-dimensional finite element model with boundary conditions derived from an analytic large-scale model based on the regional water balance. Numerical investigations and thermo-hydraulic evaluation of a typical dipole injection/extraction system are presented. Most of the simulation results are focused on the spatio-temporal extension of the hot water plume close to the injection well where the main challenges occur with respect to numerical stability. Because the (aquifer thermal energy storage system is located close to the groundwater divide, the energy recovery is less sensitive to the well configuration with respect to the groundwater flow direction.

## INTRODUCTION

The demand for energy is omnipresent, and there is a pressing need for alternatives to fossil carbon energy for environmental and economical reasons. Deep geothermal energy has been exploited commercially for a long time, but shallow aquifers also offer some promising opportunities for energy storage and extraction. One option, explored in this study, is to employ heat pump technology. Ground source heat pumps transfer heat to or from the ground and use it for heating or cooling of buildings. They are also commonly referred to as geothermal heat pumps, groundwater heat pumps, geo-exchange systems, and earth energy systems. In this paper, the term aquifer thermal energy storage (ATES) has been used.

ATES is a system which utilizes inter-seasonal heat storage. This involves storage of excess energy from summer for use in winter heating applications, and the storage of cooling potential from winter for free cooling in summer (Figure 1). For typical summer conditions, low-temperature water from a cold well is pumped through a simple heat exchanger and used for cooling. This increases the temperature of the production water to . Before injection of production water into a hot well, temperature is further increased to . In winter, the process is reversed. Warm water is pumped from the hot well and sent through the heating system to pre-heat the buildings’ air intake. In transferring thermal energy to the air, water becomes cooler. This cooler water is returned to the cold well. The performance of the ATES system depends on the difference between temperature demand ( or ) and temperature of pumping water ( or ). The smaller the temperature difference , the less energy is required in the energy plant.

Researchers have long highlighted numerical modeling for analysis and optimization of ATES systems. A three-dimensional numerical model for simulating coupled groundwater flow and thermal energy transport was developed and validated by Molson *et al*. (1992). Kim *et al*. (2010) performed three-dimensional thermo-hydraulic modeling of a cyclic mode (pumping and injection wells are switched by the season) in order to identify the thermal interference by the distance between boreholes, the hydraulic conductivity, and the injection/pumping rate on the system performance. Their studies, however, derived conclusions based on the simulation results for ATES systems without considering the effects of regional groundwater flow. It is well known that groundwater flow may have a significant impact on the system performance (Nagano *et al*. 2002). Numerical simulation of the continuous operation of an ATES system under regional groundwater flow was carried out by Lee (2011). He emphasized the influence of regional groundwater flow on the performance of the system under various operation scenarios. The present work extends previously reported research to ATES systems in large-scale groundwater flow by implementing a numerical model that applies local boundary conditions obtained from observations of the regional water balance. A study of a shallow aquifer exposed to natural groundwater flow when the direction of the regional groundwater flow is perpendicular to or along the line joining the production and injection wells will be presented. In such natural systems energy may be leaking out from the well field. Optimizing the well geometry and production volumes will improve the performance of the ATES system, and the forward model we propose is a first step to implement the essential physics. Efficiency of the ATES at Oslo Airport under the regional groundwater flow is also analyzed. This study is mainly focused on the temperature anomaly around the hot water well. The most significant problems during production occur in the high-temperature domain of the aquifer. The problems are related to chemical precipitation, and bacteria growth reduces the pumping capacity (Harder 1919). The spatio-temporal extension of the hot water plume is therefore important for two reasons. First, it defines the energy storage efficiency of the ATES system which is decisive for the economical outcome; second, it defines what could be denoted the thermal contamination of the aquifer.

The ATES system models have been simulated by the commercial finite element software COMSOL Multiphysics (COMSOL 2012). The technology offers a high-computational efficiency in the solution of three-dimensional thermal transport. The software allows multiple source/sink wells, adaptive grid geometry, and a coupling of the transient groundwater flow equation with the thermal transport equation for the aquifer system.

## MATHEMATICAL MODEL

The thermo-hydraulic analysis requires calculation of groundwater and heat transport in an aquifer consisting of a solid porous matrix with pores saturated with water . The flows are governed by a set of coupled partial differential equations describing the mass and heat energy balance.

### Groundwater flow

*et al*. 1992) Here, is the specific discharge or the Darcy velocity, the intrinsic permeability tensor,

*z*the elevation of the piezometric head relative to a datum level,

*p*the fluid's pressure, water mass density,

*g*the acceleration of gravity, and the dynamic viscosity of the water. Furthermore, the parameter is the hydraulic conductivity and is the hydraulic head.

*et al*. 1992) ([

*T*] = °C) Assuming that the solid matrix as well as the water is incompressible, conservation of mass combined with Darcy's law results in the Poisson equation where is a water source/sink density that may change as a function of time according to the pumping/injection history.

### Heat flow

*c*is the specific heat, is the thermal conductivity, and is any heat source/sink. When the solid and fluid temperatures are the same, the aquifer's heat capacity, depending on the porosity

*n*, is The advective heat flux is (Nield & Bejan 2013). The heat conductivity has contributions both from molecular diffusion, often expressed as , and mechanical dispersion, which in its simplest form may be expressed as , where is the dispersivity length. The heat diffusivity is thus A similar expression for the heat diffusivity, , is used in the COMSOL software for stabilizing numerical diffusion, where the -parameter is of the order of the numerical cell diameter, and is a dimensionless tuning parameter.

### Energy efficiency

*E*in the well over a time interval may be written as where

*T*is the water temperature at the well, the ambient groundwater temperature, and

*Q*(m

^{3}/s) the water injection or extraction rate. The energy recovery, measured in relative terms is

*et al*. 1982). Similarly, for ATES systems having more than one injection/extraction well, the system efficiency may be obtained by dividing the sum of extracted energy from all wells by the sum of injected energy to all wells.

## THE ATES SYSTEM AT OSLO AIRPORT

### Location and size of the ATES system

The ATES system at Oslo Airport is located at Gardermoen, on the largest precipitation-fed aquifer of mainland Norway. The aquifer is part of a glaciofluvial delta structure (Tuttle 1997). Sediments consist of different layers of clay in the deeper parts, and sand and gravel in the upper parts of the delta.

A local domain of the aquifer is utilized as a source for heating and cooling of the indoor space at Oslo Airport. The ATES system consists of 18 wells, nine warm and nine cold, each with a diameter of 450 m and drilled down to 45 m below the surface. The cold wells are located 180 m east of the warm wells. Each well is supplied with its own groundwater pump and injection tube. Warm and cold water are extracted from or injected to the aquifer from the wells approximately at depths from 20 to 45 m below the surface level. The Gardermoen delta and the location of the geothermal energy wells are shown in Figure 2. The ATES system is located close to the regional groundwater divide.

### Flow parameters and experimental data

Initial physical parameters for the Gardermoen aquifer (Table 1) were obtained from Goshu & Omre (2003) and Tuttle (1997). The wells in the ATES system were tested after construction with respect to maximum pumping yield, . As part of the maintenance procedure, a similar pumping test was repeated for each well, . All pumping tests were performed after a temporary shut up of the ATES system in order to obtain a resting groundwater head in the aquifer. Thus, initial conditions were close to the natural groundwater flow for all wells. In addition to monitoring , the drawdown of groundwater head, , was monitored at the pumping wells. The purpose of the test was to identify reduction of pumping yield, . Reduction of the pumping yield is due to chemical precipitation and bacterial growth inside and around the screens (filter tips) of the wells. Observations from pumping tests of two wells in the period 20–27 March 2012 are shown in Figure 3. This result will be used later for validation of the numerical experiments.

Property | Symbol | Value |
---|---|---|

Porosity | 0.1507 | |

Horizontal permeability | 1.08 × 10^{−11} m^{2} | |

Vertical permeability | 4.64 × 10^{−12} m^{2} | |

Density of fluid | 1,000 kg/m^{3} | |

Density of aquifer | 2,630 kg/m^{3} | |

Specific heat of fluid | 4,200 J/kgK | |

Specific heat of solid | 800 J/kgK | |

Thermal conductivity of fluid | 0.6 W/mK | |

Thermal conductivity of solid | 2.0 W/mK |

Property | Symbol | Value |
---|---|---|

Porosity | 0.1507 | |

Horizontal permeability | 1.08 × 10^{−11} m^{2} | |

Vertical permeability | 4.64 × 10^{−12} m^{2} | |

Density of fluid | 1,000 kg/m^{3} | |

Density of aquifer | 2,630 kg/m^{3} | |

Specific heat of fluid | 4,200 J/kgK | |

Specific heat of solid | 800 J/kgK | |

Thermal conductivity of fluid | 0.6 W/mK | |

Thermal conductivity of solid | 2.0 W/mK |

### Model setup

To demonstrate the essential physics, we decided to simplify the simulations to a pure dipole system with one injection and one extraction well. Pumping and injection volumes are chosen close to the real production volumes. The boundary conditions are derived from the analytic large-scale model of Kitterød (2004) and consistent with the regional water balance.

*x*,

*y*, and

*z*directions, respectively (Figure 4). The domain is oriented according to the regional model such that the boundary , denoted by , is placed at the groundwater divide. The natural groundwater flow is in the -direction, and thus parallel to the boundaries and . Boundary is assumed to have constant head (see derivation below). There is no water flux at the bottom of the aquifer but the aquifer is recharged by infiltration,

*N*(mm/day) from the top boundary . Thus, the groundwater flow before pumping/injection from the wells is given by Equation (4) with the boundary conditions

The temperature of the injected hot water through is set to *T _{hw}* = 30 °C. A steady-state recharge, mm/day = m/sec, at the top boundary has also been applied (Kitterød 2004).

*H*is the thickness of the aquifer. Under the condition, and , the solution of Equation (13) is Applying Equation (14), with a reference head taken from a well observation, the calculated hydraulic head on the boundary is .

The setup in COMSOL is as follows: to avoid numerical oscillations at the temperature front, some artificial isotropic diffusion has been introduced. This method adds an artificial diffusion term (, where is a tuning parameter, *h* is the mesh element size and is the convective velocity vector) to the physical diffusion coefficient (COMSOL 2012). For this study, has been used.

For the coupled flow and temperature simulations, the COMSOL 4.3a adaptive mesh refinement options have been applied. The simulations apply about 167,000 tetrahedral elements with a total of 273,500 degrees of freedom.

### Parameter evaluation

To investigate the effect of coupling temperature with the density and viscosity, we have simulated the flow (with closed wells) for two different temperatures and compared with the analytical solution. Equations (2) and (3) are inserted into the equation for hydraulic conductivity, and the result is shown in Figure 5. Figure 5 also shows the hydraulic head in the model when infiltration of 30 °C water has been recharged from the top surface for 1 year with temperature dependent density and viscosity. Simulated results are fully consistent with the analytic solutions when *T* is constant, and demonstrate the impact of numerical coupling between temperature and Darcy flow.

Our aim has been to simulate the ATES system as it was originally. Thus, only monitoring data with minimum chemical precipitation or bacterial growth in the well screens were relevant, so all drawdown observations from the pumping test with a pumping yield less than 0.75 compared to the original pumping yield have been removed. Hence, only the four wells with are included in the simulations presented in Figure 6. To obtain a simulated hydraulic head consistent with the observed drawdown, the permeability was increased by a factor of 1.75 compared to Goshu & Omre (2003). As illustrated in Figure 6, the flow model does not reproduce exactly the monitoring data, but the simulated drawdown is consistent with the range of the observed drawdowns.

## NUMERICAL RESULTS

As mentioned above, the simulations are limited to a single dipole well system, and the first-model configurations deal with the well orientations relative to the regional groundwater flow. Three cases are considered, of which one corresponds to the well orientation at the Oslo Airport ATES system. There have been some uncertainties about the influence of the temperature dependence on density and viscosity, and that was investigated before proceeding to further simulations. If the dependence proved to be minor, the rest of the configurations can be simulated by constant density and viscosity. This also helps the code to run faster. The study continues with investigations of the temperature behavior at four locations in the field, again considering the effect from the orientation of the wells and the regional groundwater flow. A final configuration provides spatio-temporal simulations of the temperature propagation with a cyclic pumping/injection operation.

Steady-state flow for three different well configurations has been used to illustrate energy storage efficiency and risk for heat leakage out of the system (Figures 7(a)–7(c)). The temperature at four different points (A, B, C, and D), 14 m above the bottom of the model, has been extracted for further inspection. For Figure 7(a), the regional groundwater flow is perpendicular to the line joining the two wells. Since recharge from precipitation is relatively small compared to the injection/pumping rates, a large fraction of the velocity field is directed toward the pumping well in this case. In Figures 7(b) and 7(c), pumping and injection wells are oriented along the groundwater flow direction. If injection of hot water takes place according to Figure 7(b), and injection of cold water according to Figure 7(c), this well configuration should minimize energy leakage out of the ATES system. At the same time, this well configuration increases the risk of affecting the temperatures around the downstream well. Thus, Figure 7(a) can be considered as a trade-off between risks for heat leakage out of the system, and cross-contamination of temperatures.

Temperature at the observation points (locations A–D in Figure 7) were simulated for well configuration cases (a)–(c) (Figures 8 and 9). For simplicity, all configurations had injection temperature equal to 30 °C. Figure 8 shows the effect of temperature dependence on density and viscosity. Since the effect turns out to be relatively minor, further runs apply temperature independent parameters. The breakthrough curves in Figure 9 indicate minor temperature difference at the four points. For cases (b) and (c), temperature development in point A and B is almost identical, while the asymmetrical flow field for case (a) is evident in observation points A and B. The breakthrough curve in point C is sharper for all cases compared to location D, and the steady-state temperatures are higher. These results reflect the fact that observation point C is located close to the pumping well while point D is located most distant from the pumping well.

The following simulation results indicate the spatio-temporal temperature development for well configuration in Figure 7(a). After 6 months of continuous injection of hot water, the temperature anomalies around the injection well are almost symmetric, indicating minor influence of infiltration velocities from precipitation (Figure 10). In Figure 11, the pumping was reversed, i.e., there is pumping of hot water and injection of cold water. Here, the impact of the regional groundwater flow (in *x*-direction) is clearly captured in the spatial temperature distribution after 3–6 months. The temperature anomalies after a full injection/pumping cycle indicate the storage efficiency of the ATES system.

## DISCUSSION

After 6 months of constant pumping/injection (28 m^{3}/hr), the simulations indicate no risk of breakthrough of hot water in the cold water pumping well (Figure 10). Close to the injection well there is almost no difference in temperature between upstream and downstream locations with respect to the background groundwater flow (Figure 9). Apparently, since the infiltration rate is relatively small compared to the injection/pumping rate, the temperature difference between the observation points in the model is insignificant. In the case study, the total infiltration rate into the model domain is about 14 m^{3}/hr, which is about 50% of the pumping/injection rate.

The profile in Figure 9 shows the temperature at the four different points for all three cases. In all cases, the point in the direction of the original groundwater flow with respect to the injection well has a higher temperature than the others. However, since the infiltration rate is small compared to the pumping/injection rate, the natural groundwater flow, which is perpendicular to the line joining the two wells, has only a small impact on the flow propagation front.

The current model may be applied, with appropriate modifications, to several other scenarios such as the temperature distribution in the aquifer and hot water contamination of the cold well after many years of production.

To minimize the advective loss of heat energy, the ATES system is located close to the groundwater divide, which gives very low (natural) groundwater flow velocities. Thus, the final aspect of storage efficiency discussed here is related to the risk of thermal cross-contamination of injection/pumping wells. This problem is addressed by running a half-year continuous hot water injection. Half a year is the most likely maximum continuous injection period, because after half a year there will be a demand for heating and therefore pumping of hot water. After 6 months of continuous injection, there is no indication of thermal cross-contamination. The simulations after repeated pumping/injection cycles show that there will be a minor temperature increase around the hot water well before the next hot water injection takes place. This temperature anomaly will be transported either to the cold water well or out of the model domain. After one pumping/injection cycle the temperature increase is only 1–2 °C, which indicates the significant energy recovery of the Oslo Airport ATES system.

The storage efficiency depends on the Peclet number which expresses the ratio between the advective and the diffusion/dispersion transport component. Sharp fronts give higher energy recovery than dispersed fronts. The importance of mechanical dispersion for temperature transport is a subject of debate (Anderson 2005). Some authors (Bear 1972; Ingebritsen & Sanford 1999; Hopmans *et al*. 2002) argue that mechanical dispersion is small relative to the molecular diffusion for heat transport, while others (Smith & Chapman 1983; de Marsily 1986) argue that mechanical dispersion should be taken into account. It is well documented that mechanical dispersion increases with transport distance. In our case, where transport of heat is limited to a radial distance of about 50 m, mechanical dispersion may be suppressed to a minimum. In the present study case only thermal diffusion is considered, but to obtain numerical stable solutions in COMSOL, some artificial dispersion is introduced. The Peclet number was between 100 and 200 close to the wells, and from 20 to 30 far from the wells or if the ATES system is at rest. Figure 12 shows the average temperature of the extracted water. The temperature of the water extracted in case (c) is less than the others after a few days. Owing to minor impact of the regional groundwater flow, the energy recovery factors (Equation (10)) for the three different ATES configurations were very similar (Figure 12). The theoretical recovery factors increase as a function of increasing number of injection/pumping cycles, from 0.75 to 0.85 (see Figure 13). It should be emphasized that these results indicate an upper theoretical limit for the energy efficiency of the Gardermoen ATES system. In reality, the energy recovery declined after a few years due to chemical precipitation and iron-bacteria in the aquifer. The numerical model presented in this study may serve as a starting point for further studies of coupling chemical and biological processes to temperature transport in the aquifer.

### Limitations of the study

The simulation results are highly affected by the boundary conditions. In this study, no temperature flux at the boundaries were allowed. Some initial simulations, not presented in this study, indicate clearly that the thermal insulation of the aquifer is crucial for the energy recovery (Doughty *et al*. 1982).

The variation in density and viscosity across the ATES system in this case study had only a minor effect on the propagation of the temperature anomalies. Density differences cause buoyancy effects, thus producing additional vertical fluid motion. In fact, density gradients can introduce gravitational instabilities giving rise to circulation of water in the aquifer, while the viscosity only alters the resistance to groundwater flow (Ma & Zheng 2009). With the temperature range in this case study from 4 to 30 °C, the effect on density is less than 1%, while the differences in viscosity is about 50% (cf. Equations (2) and (3)). The impact on vertical velocities due to changes in density is also reduced due to the vertical/horizontal permeability anisotropy (Table 1). Thus, the importance of viscosity in this case is one to two orders of magnitude more important than changes in density.

The Oslo Airport ATES system is operated under cyclic flow regimes that may be reversed several times during a week and with variable pumping/injection rates. Rapid fluctuations and variable pumping/injection rates create dispersed temperature anomalies around the wells while long production periods with constant pumping/injection rates are expected to create sharper temperature fronts.

## CONCLUSION

This paper briefly reviews the differential equations for ATES and shows how boundary conditions for the equations may be obtained from regional data. By using estimated flow parameters from previous studies of the Gardermoen aquifer, we simulate drawdowns in the same range as recorded during operational well tests. Numerical simulations have been carried out to see the flow and thermal behavior of an ATES system with two wells under cyclic operation at the Gardermoen aquifer. The energy recovery of an ATES system depends on the natural groundwater flow. Recovery of hot water increases if dispersion is small compared to advective flow velocity. Repeated injection and pumping cycles increase the theoretical energy recovery. Because the ATES system is located close to the groundwater divide, the natural groundwater flow is minor. For that reason, the configuration of the injection and pumping wells with respect to the direction of the natural groundwater flow has only minor impact.

The simulation results indicate no risk for breakthrough of hot water from the injection to the pumping well during 6 months of steady-state production for the three different well geometries analyzed in this study.

From a hydrological point of view, the ATES flow model should have boundaries coinciding with the aquifer, but this requirement would make the ATES model too big for the computer. Thus, for feasible simulations, there is always a trade-off between computational efficiency and the size of the local flow domain. In the ATES system at Gardermoen, after a continuous injection of hot water for 6 months, the cold water reservoir was never affected by hot water. It is important to evaluate the efficiency of the ATES system with respect to the regional groundwater flow. In this case, the energy recovery was not sensitive to the direction of the groundwater flow.

## ACKNOWLEDGEMENTS

The authors want to thank Dr Kevin J. Tuttle, Norconsult, and staff at Avinor, Oslo Airport for inspiring discussions and for providing the data used in this study. The authors also express their appreciation to Dr Bertil Nistad for his assistance on the COMSOL software. We also thank the unknown referees for their many helpful comments.