## Abstract

Aquifer thermal energy storage (ATES) systems offer reduced energy costs, lower carbon emissions, and increased energy resilience. The feasibility, however, depends on several factors and usually require optimization. We study an ATES system with injection and extraction wells (cf. graphical abstract). The purpose of the investigation was to calculate the recovery factor of an ATES system with a cyclic repetition of injection and pumping. In the paper, we discuss analytical and numerical radial solutions of differential equations for heat transport in water-saturated porous media. A similar solution was obtained for a 2-D-horizontal confined aquifer with a constant radial flow. Numerical solutions were derived by using a high-resolution Lagrangian approach suppressing spurious oscillations and artificial dispersion. The numerical solution and the analytical solutions give consistent results and match each other well. The solutions describe instantaneous and delayed heat transfer between fluid and solid, as well as time-varying water flow. In hydrological terms, these solutions are relevant for a wide range of problems where groundwater reservoirs are utilized for extraction and storage (namely, irrigation; water supply; geothermal extraction).

## HIGHLIGHTS

Analytical and numerical solution of symmetric aquifer thermal energy storage is analyzed.

Numerical solution for coupled non- equilibrium temperature for liquid and solid phase is presented.

Energy efficiency of a simplified aquifer thermal energy system for cyclic pumping and injection system is calculated for different grain size as a function of time.

## INTRODUCTION

Heat transfer in porous media has received considerable attention and has been the topic of a number of investigations during the last decade (Bejan & Kraus 2003). A driving force for research on this subject is engineering applications, such as geothermal systems (Ganguly & Kumar 2014), heat exchangers (Diao *et al.* 2004), thermal insulation (Kim *et al.* 2010; Birhanu *et al.* 2015b), safety issues regarding storage of nuclear waste (Sun *et al.* 2010), flow in porous media is in the context of nuclear waste repositories (Medici & West 2022), and nested shallow geothermal systems (García-Gil *et al.* 2020).

In addition to equations for the fluid flow, the mathematical model of heat transfer in porous media is given by second-order partial differential equations for heat energy conservation and flow in the model domain. Two kinds of models can be applied to investigate thermal characteristics of conduction and advection within a porous medium, namely, a thermal equilibrium model and a thermal non-equilibrium model (Yang *et al.* 2011). The difference between the two models is the thermal coupling between the liquid and the solid phases. For the equilibrium model the coupling is modeled as an instantaneous heat transfer. This assumption is close to the reality for homogeneous aquifers with solid particles of minor size (diameter mm). For the non-equilibrium model, there is a time delay attached to the heat transfer between the two phases. In the literature of transport in porous media, this model is usually called the double porosity model which may also be expanded to a dual permeability model (Stadler *et al.* 2012). The dual permeability model refers to a conceptual framework used in hydrogeology to describe subsurface fluid flow in heterogeneous porous media. This model accounts for two distinct types of permeabilities within the subsurface: matrix permeability, which characterizes fluid flow through the bulk of the rock matrix, and fracture permeability, which pertains to fluid flow through interconnected fractures or fissures within the rock. Fractures play a significant role in the dual permeability model by providing preferential pathways for fluid movement, often exhibiting significantly higher permeability compared to the surrounding matrix. These fractures can significantly influence the overall flow behavior and fluid transport in the subsurface, enabling faster movement of fluids and altering the distribution of groundwater, solutes, and contaminants. In essence, fractures are integral components in the dual permeability model, introducing complexity and heterogeneity that can have substantial implications for various subsurface processes, including groundwater flow, contaminant transport, and resource extraction (Aguilar-López *et al.* 2020; Agarwal *et al.* 2021).

Researchers have highlighted different analytical and numerical methods to find the solution for this model based on different physical phenomena. Bear & Cheng (2010) presented the solution of heat transfer in porous media for 1-D by using the similarity solution method. A 2-dimensional (2-D) numerical model for heat transport in a heterogeneous porous aquifer thermal energy storage (ATES) system is presented by Ganguly *et al.* (2014). They considered the transient heat transport phenomenon in a heterogeneous porous aquifer due to hot water injection and validated the solution analytically.

Underground thermal energy storage (TES) systems can be broadly categorized into four main types, each with its unique characteristics and applications. These categories are tank TES, pit TES, borehole TES, and aquifer TES (also known as ATES). Tank TES and pit TES are constructed independently of the geological conditions of the aquifer and remain unaffected by groundwater flow. This is achieved by incorporating an insulation layer on the outer surface of the container or structure. In contrast, borehole TES operates as a closed-loop system, storing energy by utilizing ground heat exchangers within the aquifer. This stored energy can then be supplied to buildings for either heating or cooling purposes (Lin *et al.* 2019). ATES is an example of technology where subsurface storage and transport of heat are used to save energy. ATES systems may utilize inter-seasonal heat storage, which means storage of excess energy from summer that is used in winter time for heating purposes. For cooling purposes low-temperature water is extracted from cold wells and heated water is injected into hot water wells (Birhanu *et al.* 2015b). Thus, ATES installations actively store cooled and heated groundwater in the ground from respective heating and cooling mode cycles (Dickinson *et al.* 2009). An ATES system involves both the flow of water and heat transport. In order to predict the performance and efficiency of an ATES system, one possibility is to run detailed numerical simulations, and researchers have long highlighted numerical modeling for the analysis and optimization of ATES systems (Lee 2011). In some cases, the injection and pumping wells may be simplified to a classical dipole geometry. In this case, the flow field may be simplified to axis symmetry where the flow velocities are governed by the injection/pumping rate and the aquifer porosity. In the present study, the advantage of a simplified flow field is utilized, and the transport equation is solved using simple analytical and numerical methods to assess the energy efficiency of an idealized ATES system. Numerical solutions to transport problems are usually affected by artifacts, and because all mathematical models are simplifications of reality, boundary conditions should be specified with great precautions. Birhanu *et al.* (2015a) showed that the most evident boundary condition of temperature at the top of an unconfined aquifer gave unphysical energy efficiency for simulation experiments of a real ATES system. In this study, analytical solutions are employed to assess the quality of numerical simulations through the execution of (simple) numerical experiments. This study aims to bridge the existing research gap by addressing the complexities associated with radially symmetric ATES problems. Our primary objectives are as follows: (i) to examine both analytical and numerical radial solutions for the differential equations governing heat transport within water-saturated porous media; (ii) To obtain numerical solutions by employing a high-resolution Lagrangian approach that effectively eliminates unwanted oscillations and artificial dispersion. (iii) To determine the recovery factor for an ATES system through a cyclic process involving injection and pumping. An experiment that is conducted involves an idealized ATES production sequence, wherein hot water is repeatedly injected and pumped within a confined aquifer. The performance of the alternative solutions is quantified by a recovery (or efficiency) factor. To fulfill these objectives, our study employs a comprehensive research methodology encompassing analytical derivations, numerical simulations, and practical case studies. Existing analytical solutions will be modified and expanded upon while implementing numerical models that incorporate the radially symmetric characteristics of aquifers. The combination of theory and application will provide valuable insights into the optimization and management of thermal energy storage systems in such aquifers.

## MECHANISMS AND EQUATIONS

*et al.*1992; Nield & Bejan 2013):

**k**is the intrinsic permeability tensor,

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

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

*g*is the acceleration of gravity, and is the dynamic viscosity of water. Furthermore, is the hydraulic conductivity and is the hydraulic head. The volume average velocity differs from the velocity of the water in the pores, the so-called seepage velocity, , where

*n*is the (effective) porosity (Kangas & Lund 1994). The water density and, more pronounced, viscosity varies with temperature. However, in this context, we will assume that the flow field represented by

*q*is unaffected by temperature changes. Assuming that the solid matrix as well as the water are incompressible, mass conservation combined with Darcy's law leads to the Poisson equation for the hydraulic head:where is a source/sink term.

*c*is the specific heat, and subscripts

*w*and

*s*refer to water and solid. At a local temperature equilibrium, where , the heat content may be expressed as , where (see Kangas & Lund 1994; Hecht-Méndez

*et al.*2010):see Kangas & Lund (1994) and Hecht-Méndez

*et al.*(2010). In the following, we will use the convention for . The water flow causes the advection of heat:whereas conduction/diffusion of heat takes place both in the solid and the liquid:and is the heat diffusion coefficients. If the two media are at a local thermal equilibrium, the volume average diffusive heat flux may be expressed by:where is a bulk aquifer heat diffusion coefficient (see Kangas & Lund 1994; Nield & Bejan 2013).see Kangas & Lund (1994) and Nield & Bejan (2013). Other expressions for , e.g. porosity-weighted geometric and harmonic means are also discussed in the literature (Nield & Bejan 2013). In addition, the heterogeneity of the pores induces a certain amount of thermal dispersion, parametrized, in its simplest form as:Here, is the thermal dispersivity length, and the total diffusion flux becomes (Bakr

*et al.*2013).

*h*is a heat transfer coefficient (Kreith

*et al.*2011; Nield & Bejan 2013). The coefficient varies with temperature and the flow, in particular, for large flows. Following the discussion in Nield & Bejan (2013), may be expressed as: , where is the surface area of the water/solid interface per volume unit. For low Reynolds numbers, may be expressed as:

A similar time scale may actually be derived from the heating of spheres discussed in Gockenbach & Schmidtke (2009).

We observe that when is kept constant and diffusion is neglected, the natural time scale (inverse rate of change) in Equation (19) is essentially . For less than about a millimeter, the thermal equilibrium is virtually spontaneous and we may assume that and are equal.

**q**are assumed to be independent of

*T*, then dividing through with in Equation (22) leads to:

**q**is caused by water injected or pumped with a discharge rate from a well located at the origin. When , water is injected from the well into the aquifer, causing a flow away from the well. During pumping, and the flow is directed toward the well. Utilizing the symmetric geometry of the aquifer near the well, the flow is in which the discharge velocity where:

*W*and height

*H*are constants characteristic of the aquifer. In this case, Equation (23) becomes:

## ANALYTICAL SOLUTIONS

Analytical solutions express the relation between the principal variables involved directly. This provides basic insight into the problem without any further numerical evaluation. Besides this intuitive conceptual advantage, the classical application of analytical solutions is practical parameter estimation and sensitivity analysis. In practical hydrology, physical parameters like permeability or storage capacity are not always accessible for direct measurements. Response functions (namely, hydraulic head and temperature) on the other hand, are usually more easy to monitor. In such cases, the physical parameters are estimated by solving the inverse problem. Before the advent of computer technology, this was done by dimensionless solutions of the analytical expression, which provided tables or so-called type-curves. The inverse problem was solved by curve fitting of the empirical data to the analytical type curve. Today, analytical solutions are applied in a similar way, but the curve fitting is substituted by numerical perturbation of the involved variables. Sensitivity analysis is another example where the application of analytical solutions is convenient. After the estimation of optimal parameters, the relative impact of the uncertainties might be evaluated by simple numerical perturbation of the involved parameters. Here, in this context, the motivation for using analytical expressions of the transport equation is due to the numerical challenge of solving the transport equation if advective flow is dominant. In such cases, numerical solutions are prone to numerical dispersion. Even though analytical solutions simplify real transport, numerical artifacts do not affect the solutions. Therefore, by using the same simplified velocity field for both the analytical and the numerical solutions, the performance of the numerical algorithm can be evaluated directly by using the analytical solutions as benchmarks.

If the diffusion term is negligible, Equation (25) becomes a simple hyperbolic equation which, for any initial temperature distribution, , has the solution . In particular, for the conditions in Equation (28), the hyperbolic solution is the moving front where is the complementary Heaviside function ( for , for ). Since is typically about twice as large as , the ratio depends on the porosity *n* and varies between 1 and 2. The temperature front thus moves significantly faster than the discharge velocity , but slower than the average seepage velocity, .

### The 1-dimensional (1-D) case

*x*rather than . In this case, a well-known similarity variable is , which, inserted into the equation results in:with the general solution . However, no solution from this collection satisfies the boundary condition . Nevertheless, the similarity solution of the closely related problem satisfying the initial values , is a good approximation:

### The 2-D radial symmetric case

The spherical symmetry 3-D-case is easily seen to have intrinsic scales for *r* and *t* involving and , and no simple similarity solution exists. The scaled equation may however be transformed to a 1-D heat equation with a space-dependent diffusion coefficient (Philip 1994). Actually, the numerical algorithm below applies a similar transformation.

## COMPUTATIONAL PROCEDURE

We shall now consider a numerical algorithm for solving Equations (19) and (20), or the equilibrium model Equation (23) in the case of symmetric geometry. It turns out that the numerical solutions to these problems are nontrivial. They are typically advection-dominated, and we have already seen in the previous section that the temperature profile is a sharp front moving away from the source. In the radial and spherical case, the flow becomes very large close to the origin, leading to an almost hyperbolic equation in this region. Advection-dominated problems are notoriously difficult to solve numerically. Popular schemes like central differencing schemes result in unstable or spurious oscillatory solutions. Upwind discretization for the advection term avoids oscillations, but does create artificial diffusion, leading to a smoothed temperature front when applied on a coarse grid. Several other methods have been proposed and discussed in the literature, see e.g. Strikwerda (2004) and LeVeque (1992).

A computational procedure utilizing the special structure of Equation (25) or Equations (26) and (27) is proposed. The procedure is composed of the well-known Lagrangian approach combined with a coordinate transformation. The idea will first be explained for the equilibrium case Equation (25).

*t*, and for the space variable . A reasonable relation between the two scales is , is the distance a hyperbolic temperature front moves in time for a constant . The temperature

*T*will typically be scaled with the temperature of the injected water. With these scales, we get the dimensionless equation:where ,. As a curiosity, notice that in the case of a constant , by choosing and such that , the parameters are completely absorbed, in the sense that also become equal to one. This is possible for and , but not for . However, we will not pursue this curious issue any further here.

*T*with respect to time becomes:which, inserted into Equation (41) gives:

*v*is constant, in which case the characteristics are just straight lines.

In order to be able to resolve a sharp front, the characteristics used in the discretization can be concentrated around it.

**Example 1.**

*The algorithm described above is used for solving*Equation (40)

*for*

*, using*

*,*

*and*

*. The exact solution*

*is given by*Equation (39)

*. The problem is first solved by the algorithm described above. The initial computational*

*domain is*

*where*

*is*

*chosen sufficiently large to avoid any influence from the boundaries. The concentration of*

*characteristics around the front is achieved by using*:

*where*

*p*

*is a positive*

*integer, the higher*

*p*

*the stronger concentration. In our experiments, we have used*

*,*

*and*.

*For comparisons,* Equation (40) *is also solved by a standard difference scheme with constant**stepsize. In this case, the advection term is approximated with an upwind scheme**.**For the diffusion term a central difference scheme is applied. The spatial gridsize is*.

*The problem is solved for different values of**, and**the results are shown in*Figure 4 *together with the exact solutions given by* Equation (39)*. For**, the**diffusion is large and the artificial diffusion of the upwind method is insignificant. For smaller values**of**,**the front remains sharp, and the effect of the artificial diffusion of the upwind method**becomes quite pronounced. The Lagrangian approach preserves the sharp temperature**front*.

### The non-equilibrium case

**Example 2.**

*Consider the nondimensional equations*

*Equations (48) and (49)*

*, using*

*and*

*parameters*:

*and a time dependent flow,*

*.*

*For*

*, both*

*injection,*

*,*

*and pumping,*

*,*

*are demonstrated. The equations are solved with the Lagrangian approach, using a central*

*difference approximation for the diffusion terms, and a downwind scheme*:

*for the advection term. We have used a spatial stepsize*

*, varying*

*from*

*to*

*and*

*concentrated around the temperature front of the water. The semidiscretized system is solved by*

*MATLABs ODE15s*.

*The result of the simulation is presented in*Figure 5. *We can clearly see how the hot water plume develops with time**. The**concentration of characteristics moves with the temperature front, which remains sharp. At the**same time, there is a heat exchange from water to solid, and the temperature of the water behind**the front is reduced at the same time as that of the solid temperature increases. The diffusion is**almost negligible most of the time, but will cause a slight smoothing of the front. When**the**front approaches the well with increasing velocity, and the smoothing of the front becomes more**pronounced. This example illustrates the numerical challenge, namely the coupling of advective and**diffusive physics. Usually, the numerical solutions of such problems suffer from artifacts**such as numerical diffusion and/or oscillations, but in this case, it is possible to suppress the**numerical artifacts to negligible levels*.

## CASE STUDY: TEMPERATURE PROFILE NEAR A WELL IN AN ATES SYSTEM

In this section, we consider the temperature propagation around a hot well in an ATES system. The physical and thermal properties of the Gardermoen aquifer were obtained from Goshu & Omre (2003) and Tuttle (1997), and are given in Table 1. The ATES system typically operates in one of the two modes: Injecting hot water during daytime and extracting it during nighttime, resulting in a full operational cycle of 24 h. Alternatively, the warm water is injected during the summer months, and extracted in the winter, giving a cycle of one year. In reality, a combination of these is used, but in our study, we only consider the first modi, assuming the injection and extraction periods to be of equal lengths, 12 h.

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

Porosity | 0.1507 | |

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

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

Specific heat of fluid | 4,200 J/kg K | |

Specific heat of solid | 800 J/kg K | |

Thermal conductivity of fluid | 0.6 W/mK | |

Thermal conductivity of solid | 2.0 W/mK | |

Injection/pumping rate | 28 m^{3}/h | |

Temperature of the injected water | 30 °C | |

Aquifer initial temperature | 4 °C | |

Aquifer height | 24.4 m |

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

Porosity | 0.1507 | |

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

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

Specific heat of fluid | 4,200 J/kg K | |

Specific heat of solid | 800 J/kg K | |

Thermal conductivity of fluid | 0.6 W/mK | |

Thermal conductivity of solid | 2.0 W/mK | |

Injection/pumping rate | 28 m^{3}/h | |

Temperature of the injected water | 30 °C | |

Aquifer initial temperature | 4 °C | |

Aquifer height | 24.4 m |

. | 500 mm . | 100 mm . | 10 mm . | 5 mm . | 1 mm . |
---|---|---|---|---|---|

49.3 | |||||

10.3 h | 25 min | 15 s | 3.8 s | 0.15 s |

. | 500 mm . | 100 mm . | 10 mm . | 5 mm . | 1 mm . |
---|---|---|---|---|---|

49.3 | |||||

10.3 h | 25 min | 15 s | 3.8 s | 0.15 s |

The equations are solved by the numerical approach outlined in Section 4.

### Transient injection phase

The similarities of the top left and the lower right plots are depicted in Figure 6. This is due to the fact that the thermal transfer coefficient *h* given by Equation (14) is proportional to . For the 2-D radial symmetry case, the two solutions may be proved to be identical up to scalings of *t* and .

### Energy efficiency

*E*in the well over a time interval is given by:where is the water temperature at the well. The efficiency can be measured in terms of the energy recovery factor given by (Doughty

*et al.*1982):Clearly, if the injected water has a constant temperature and the injection rate

*Q*is constant, then . During pumping, the water temperature at the well will vary depending on the dispersion of temperature in the aquifer, which includes natural heterogeneity and the heat transfer between the liquid and the solid phases. A numerical simulation of the temperature of the water and solid at the well over five consecutive cycles is given in Figure 7, as well as the recovery rate over these cycles.

We observe that the heat exchange has a significant impact on the efficiency rate for mm, otherwise not. We notice that the efficiency recovery rate based on this simplified model corresponds well with the rates achieved for an ATES system in the same aquifer presented in Birhanu *et al.* (2015a).

## CONCLUSION

This paper has briefly reviewed the differential equations for heat transport in water-saturated porous media, and presented numerical and analytical solutions for radially symmetric flows. In particular, a simple similarity solution was obtained for the heat transfer in a 2 horizontal confined aquifer in local fluid/solid thermal equilibrium. For a time-varying groundwater flow and different fluid and solid temperatures, that is, the non-equilibrium or delayed case, solutions have to be obtained numerically. The numerical algorithms have been based on a semi-discrete Lagrangian formulation.

The major conclusions of this study are summarized as follows: (i) the numerical models have enabled us to consider the primary purpose of this investigation, namely to calculate the recovery factor of a one-well ATES system with a cyclic repetition of injection and pumping; (ii) The performance is found to rely on the overall cycle duration in relation to the timescale for heat transfer between the fluid and solid. The latter may be linked to the typical grain size as shown in Table 2. For a total cycle of length 24 h, referring to Figure 7, the performance is seen to be virtually independent of the grain size as long as is less than about 100 mm ( min), but significantly affected for mm ( h). In the latter case, the efficiency is also significantly reduced; (iii) The analytic and numerical solutions should provide a consistent tool for the understanding of water and solid temperatures near wells with radial flows.

## ACKNOWLEDGEMENTS

The authors thank the NOMA-Mathematical and Statistical Modeling project along with the Norwegian Educational Loan Fund (Quota program) for the financial support to conduct this research.

## DATA AVAILABILITY STATEMENT

All relevant data are included in the paper or its Supplementary Information.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

*A Posterior Inverse Model for Porosity and Hydraulic Conductivity in A Ground- Water Aquifer*

*Ph.D. thesis*

*Sedimentological and Hydrogeological Characterisation of A Raised ice-Contact Delta – the Preboreal Delta-Complex at Gardermoen, Southeastern Norway*