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

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

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.

A thermo-hydraulic analysis requires the calculation of simultaneous water and heat transport in an aquifer consisting of a solid porous medium () with pores filled with water (). The water flow depends on the properties of the water as well as the solid, and the gradient of the hydraulic head, as stated in Darcy's law (Molson et al. 1992; Nield & Bejan 2013):
(1)
Here, is the specific discharge or the Darcy velocity, 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:
(2)
where is a source/sink term.
The heat energy content per aquifer volume unit may be written:
(3)
where 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):
(4)
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:
(5)
whereas conduction/diffusion of heat takes place both in the solid and the liquid:
(6)
(7)
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:
(8)
where is a bulk aquifer heat diffusion coefficient (see Kangas & Lund 1994; Nield & Bejan 2013).
(9)
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:
(10)
Here, is the thermal dispersivity length, and the total diffusion flux becomes (Bakr et al. 2013).
If a thermally insulated aquifer with initial temperature and , energy conservation implies that the system attains an equilibrium temperature equal to the weighted mean:
(11)
How fast the temperature equilibrium between water and solids is reached depends on the efficiency of the energy exchange between the two media. It turns out to be reasonable to express the heat exchange per time and volume unit as:
(12)
where 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:
(13)
Here, is the size of the grains making up the solid. solid. The expression for h then becomes:
(14)
A rough estimate of the time scale toward thermal equilibrium may be obtained from the energy exchange per time unit at the start of the heating, , compared to the required amount of energy to be transferred, :
(15)
The time scale is thus only dependent on basic material constants and the grain size. With typical values for rock:
(16)

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

For an elemental aquifer volume R with boundary , the solid's integral conservation law reads:
(17)
Note that since the solid is not moving in this case, there is no advective heat for the solid phase. Similarly, the integral conservation form for heat in the water is:
(18)
The differential forms of the conservation laws with the assumptions above become:
(19)
(20)

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.

For the case where , we obtain by adding Equations (17) and (18):
(21)
and the corresponding differential form:
(22)
If the parameters like , and the flow q are assumed to be independent of T, then dividing through with in Equation (22) leads to:
(23)
In our model, the flow 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:
(24)
The width W and height H are constants characteristic of the aquifer. In this case, Equation (23) becomes:
(25)
Using the same symmetry considerations in the non-equilibrium case, Equations (19) and (20) become:
(26)
(27)

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.

The formation of a hot water plume in a local thermal equilibrium aquifer will be examined, which is generated by a constant hot water source located at the origin. Consider Equation (25) for convenience normalized such that and with initial and boundary conditions:
(28)

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

When Equation (25) becomes:
(29)
where we choose the more standard spatial variable x rather than . In this case, a well-known similarity variable is , which, inserted into the equation results in:
(30)
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:
(31)
A modification of this solution, satisfying all conditions in Equation (28) exactly has been derived from Ogata & Banks (1961); see also Bear & Cheng (2008, Eq. 6.4.30):
(32)
The similarity solution Equation (31) and the exact solution Equation (32) are presented in Figure 1 together with the hyperbolic front solution . As expected, the similarity solution in Equation (31) does not satisfy the boundary conditions at . Still, as tends to 0, the solution approaches the hyperbolic front, and Equation (31) becomes a very good approximation.
Figure 1

The 1-D similarity solution Equation (31) and the exact solution Equation (32) compared with the hyperbolic front for two different values of the diffusion coefficient .

Figure 1

The 1-D similarity solution Equation (31) and the exact solution Equation (32) compared with the hyperbolic front for two different values of the diffusion coefficient .

Close modal

### The 2-D radial symmetric case

For a 2-D problem, assuming radial symmetry, Equation (25) becomes:
(33)
which may be rewritten as:
(34)
Again, it turns out that assuming the similarity variable we obtain an equation:
(35)
with general solution:
(36)
The solutions may be written in terms of the incomplete -function, defined as:
(37)
The radial 2-D similarity solution becomes:
(38)
or:
(39)
and the solution is shown is Figure 2.
Figure 2

The exact similarity solution Equation (39) compared with the hyperbolic front for two different values of the diffusion coefficient in the 2-D radial symmetric case.

Figure 2

The exact similarity solution Equation (39) compared with the hyperbolic front for two different values of the diffusion coefficient in the 2-D radial symmetric case.

Close modal

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.

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

Before discussing the numerical scheme, it is convenient to say something about the scaling of the problem. Consider a time scale for the time 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:
(40)
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.
To handle the singularities in origin in the radial and spherical cases, we introduce the transformation so that , valid for all dimensions . In this case, Equation (40) becomes:
(41)
Notice that for the , elucidating the hyperbolic nature of the problem close to the origin. The numerical difficulties of hyperbolic problems can be resolved by using a Lagrangian method: Given a path in the plane. The solution along this path is and the total derivative of T with respect to time becomes:
(42)
which, inserted into Equation (41) gives:
(43)
If we let the path satisfy , the advection term is eliminated. In fact, the paths are the characteristics of the hyperbolic equation we obtain for . As a result, Equation (41) can be solved as a system of differential equations:
The first equation is an ordinary differential equation, whereas the second one is a heat equation with a space-dependent diffusion coefficient. This can be discretized in space by some appropriate finite difference schemes, e.g.:
(44)
(45)
with , and initial values and . The procedure is significantly simplified if v is constant, in which case the characteristics are just straight lines.
The spatial domain can be extended to by defining for . In this case, we may solve Equation (41) with the boundary conditions:
(46)
When water is injected, and the temperature of the water at the well is . This is realized by choosing whenever and for . Here, is either the initial time or a switching time, that is whenever changes from negative to positive (from pumping to injection). The procedure is illustrated for the injection phase in Figure 3.
Figure 3

The extended domain and the characteristic lines in the injection phase.

Figure 3

The extended domain and the characteristic lines in the injection phase.

Close modal
Figure 4

Comparison of the numerical solutions based on the Lagrangian method, and a classical upwind scheme of the 2-D radial symmetric case (Equation (40)) for different values of the diffusion parameter . The Lagrangian solution overlaps the exact solution for all three values of .

Figure 4

Comparison of the numerical solutions based on the Lagrangian method, and a classical upwind scheme of the 2-D radial symmetric case (Equation (40)) for different values of the diffusion parameter . The Lagrangian solution overlaps the exact solution for all three values of .

Close modal
Figure 5

Temperature distribution of water (left) and solid (right) in the aquifer as a function of time t and radial distance r for . Hot water is injected for and pumped for .

Figure 5

Temperature distribution of water (left) and solid (right) in the aquifer as a function of time t and radial distance r for . Hot water is injected for and pumped for .

Close modal

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 solutionis given by Equation (39). The problem is first solved by the algorithm described above. The initial computationaldomain iswhereischosen sufficiently large to avoid any influence from the boundaries. The concentration ofcharacteristics around the front is achieved by using:
(47)
wherepis a positiveinteger, the higherpthe stronger concentration. In our experiments, we have used,and.

For comparisons, Equation (40) is also solved by a standard difference scheme with constantstepsize. 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, andthe results are shown inFigure 4 together with the exact solutions given by Equation (39). For, thediffusion is large and the artificial diffusion of the upwind method is insignificant. For smaller valuesof,the front remains sharp, and the effect of the artificial diffusion of the upwind methodbecomes quite pronounced. The Lagrangian approach preserves the sharp temperaturefront.

### The non-equilibrium case

The non-equilibrium case, Equations (26) and (27) are now being considered, which, in its scaled form, is presented as follows:
(48)
(49)
where:
(50)
and:
(51)
Again, for a given time scale it is appropriate to choose a spatial scale such that is of the size of 1. The boundary conditions in the injection case () is:
(52)
with if the temperature is scaled.
By applying the transformation and using the Lagrangian approach, we can find the solutions and on the characteristics from:
(53)
(54)
(55)
where is the velocity of the temperature front, typically . The formulation is used to construct a spatial grid which is dense around a sharp solution profile, and moves with it. This is illustrated in the following example.
Example 2.Consider the nondimensional equationsEquations (48) and (49), usingandparameters:
(56)
and a time dependent flow,.For, bothinjection,,and pumping,,are demonstrated. The equations are solved with the Lagrangian approach, using a centraldifference approximation for the diffusion terms, and a downwind scheme:
(57)
for the advection term. We have used a spatial stepsize, varyingfromtoandconcentrated around the temperature front of the water. The semidiscretized system is solved byMATLABs ODE15s.

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

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.

Table 1

Physical and thermal properties of fluid and aquifer for the thermo-hydraulic modeling of the Gardermoen aquifer

PropertySymbolValue
Porosity  0.1507
Density of fluid  1,000 kg/m3
Density of aquifer  2,630 kg/m3
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 m3/h
Temperature of the injected water  30 °C
Aquifer initial temperature  4 °C
Aquifer height  24.4 m
PropertySymbolValue
Porosity  0.1507
Density of fluid  1,000 kg/m3
Density of aquifer  2,630 kg/m3
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 m3/h
Temperature of the injected water  30 °C
Aquifer initial temperature  4 °C
Aquifer height  24.4 m

The Gardermoen aquifer is a delta structure deposited in a glacio-fluvial/glacio-marine environment during the last deglaciation of the Scandinavian crust (approximately. 10,000 B.P., Tuttle (1997)). The river discharge and the sediment load from the melting glacier were significant, which explains the wide range of grain sizes of the sediments from boulders ( mm, at the proximal side of the delta), to fine sand and silt ( mm) at the distal side of the delta. The ATES system for this case study was located in the delta foresets with homogeneous fine sand, but it is interesting to compare the energy efficiency of this ATES with alternative locations. We, therefore, let the sediments vary form mm, which corresponds to a foreset location, to mm mimicking a location close to the glacial portals. In that case, the aquifer permeability would have been better, but to keep the experiment as simple as possible, we let the pumping rate and the porosity be the same for all grain sizes. In Table 2, the value of the heat transfer coefficient , Equation (14), and the time scale toward thermal equilibrium , Equation (16), for different particle sizes are shown. So, we can conclude that within the time scales given by the injection/extraction periods, there is almost thermal equilibrium for realistic particle sizes. It is still of interest to see what happens in the initial injection phase before thermal equilibrium is established, so we will solve Equations (26) and (27). In a horizontally confined aquifer, we can assume radial symmetry in the vicinity of a well, so . Initial and boundary conditions for the first injection phase are:
(58)
(59)
Table 2

The heat transfer coefficient (Equation (14)) and the estimated time scale (Equation (16)) toward thermal equilibrium for different particle size

500 mm100 mm10 mm5 mm1 mm
49.3
10.3 h 25 min 15 s 3.8 s 0.15 s
500 mm100 mm10 mm5 mm1 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

In Figure 6, we present the temperature profiles for the first few seconds of the injection period. The first rows shows the situation for particle size mm. Thermal equilibrium happens almost immediately in this case, but the energy transfer still has an effect in the sense that the temperature front becomes more smooth. Also notice that after 0.15 s, the solid temperature at the wall has reached about two-thirds of the water temperature, while the water is almost cooled at the front. This is consistent with the fact that , thus we expect the water to cool down approximately three times as fast as the solid heats up. The lower row gives the same profile for mm and mm, and as expected, the heat exchange is significantly slower in these cases. As a consequence, the width of the front increases.
Figure 6

The temperature profile of the solid and the water for different values of (mm) in 2-D radial flow near a well. The solid lines indicate the temperature of the water and the broken lines for solid temperature. The upper row emphasizes the water and solid temperature profile of the same particle size with different timescales while the bottom row emphasizes the water and solid temperature profile of different particle sizes over the same timescale.

Figure 6

The temperature profile of the solid and the water for different values of (mm) in 2-D radial flow near a well. The solid lines indicate the temperature of the water and the broken lines for solid temperature. The upper row emphasizes the water and solid temperature profile of the same particle size with different timescales while the bottom row emphasizes the water and solid temperature profile of different particle sizes over the same timescale.

Close modal

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

In this part, we study the energy recovery from an ATES well based on 12-h injection and extraction periods. In general, the energy transfer E in the well over a time interval is given by:
(60)
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):
(61)
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.
Figure 7

The temperature of the water and the solid at the well for five consecutive cycles of 24 h for aquifers of different particle sizes. In the lower right corner, the corresponding recovery rates.

Figure 7

The temperature of the water and the solid at the well for five consecutive cycles of 24 h for aquifers of different particle sizes. In the lower right corner, the corresponding recovery rates.

Close modal

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

The temperature profile over one cycle (injection and pumping) is given for the two extreme cases mm and mm as shown in Figure 8.
Figure 8

Temperature profiles of the water (solid lines) and the solid (dashed lines) for different particle sizes. The left column shows the temperature profiles during injection, the right during pumping.

Figure 8

Temperature profiles of the water (solid lines) and the solid (dashed lines) for different particle sizes. The left column shows the temperature profiles during injection, the right during pumping.

Close modal

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.

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.

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

The authors declare there is no conflict.

Agarwal
R.
,
M. P.
&
Agarwal
R. P.
2021
Fractional flow equation in fractured aquifer using dual permeability model with non-singular kernel
.
Arabian Journal of Mathematics
10
,
1
9
.
Aguilar-López
J. P.
,
Bogaard
T.
&
Gerke
H. H.
2020
Dual-permeability model improvements for representation of preferential flow in fractured clays
.
Water Resources Research
56
(
8
),
e2020WR027304
.
Bakr
M.
,
van Oostrom
N.
&
Sommer
W.
2013
Efficiency of and interference among multiple aquifer thermal energy storage systems; a Dutch case study
.
Renewable Energy
60
,
53
62
.
Bear
J.
&
Cheng
A.
2008
Modeling Groundwater Flow and Contaminant Transport
, Vol.
23
.
Springer Verlag, Dordrecht
, pp.
89
103
.
Bejan
A.
&
Kraus
A. D.
2003
Heat Transfer Handbook
, Vol.
1
.
John Wiley & Sons
,
Hoboken, NJ, USA
.
Birhanu
Z.
,
Kitterød
N.-O.
,
H.
&
Kværnø
A.
2015a
Temperature boundary conditions for ATES systems
. In:
AIP Conference Proceedings
, Vol.
1648
.
AIP Publishing
, p.
030032
.
Birhanu
Z. K.
,
Kitterød
N.-O.
,
H. E.
&
Kværnø
A.
2015b
Numerical modeling of aquifer thermal energy efficiency under regional groundwater flow: A case study at Oslo airport
.
Hydrology Research
46
(
5
),
721
734
.
Diao
N.
,
Li
Q.
&
Fang
Z.
2004
Heat transfer in ground heat exchangers with groundwater advection
.
International Journal of Thermal Sciences
43
(
12
),
1203
1211
.
Dickinson
J.
,
Buik
N.
,
Matthews
M.
&
Snijders
A.
2009
Aquifer thermal energy storage: Theoretical and operational analysis
.
Geotechnique
59
(
3
),
249
260
.
Doughty
C.
,
Hellström
G.
,
Tsang
C. F.
&
Claesson
J.
1982
A dimensionless parameter approach to the thermal behavior of an aquifer thermal energy storage system
.
Water Resources Research
18
(
3
),
571
587
.
Ganguly
S.
,
Seetha
N.
&
Kumar
M. M.
2014
Numerical modeling and analytical validation for the movement of thermal front in a hetrogeneous aquifer thermal energy storage system
.
International Journal of Numerical Analysis And Modeling, Series B
4
,
413
424
.
García-Gil
A.
,
Mejías Moreno
M.
,
Garrido Schneider
E.
,
Marazuela
M. Á.
,
Abesser
C.
,
Mateo Lázaro
J.
&
Sánchez Navarro
J. Á
.
2020
Nested shallow geothermal systems
.
Sustainability
12
(
12
),
5152
.
Gockenbach
M.
&
Schmidtke
K.
2009
Newton's law of heating and the heat equation
.
Involve, A Journal of Mathematics
2
(
4
),
419
437
.
Goshu
A.
&
Omre
H.
2003
A Posterior Inverse Model for Porosity and Hydraulic Conductivity in A Ground- Water Aquifer
.
Ph.D. thesis
,
Norwegian University of Science
,
Trondheim
.
Hecht-Méndez
J.
,
Molina-Giraldo
N.
,
Blum
P.
&
Bayer
P.
2010
Evaluating MT3DMS for heat transport simulation of closed geothermal systems
.
Ground Water
48
(
5
),
741
756
.
Kangas
M.
&
Lund
P.
1994
Modeling and simulation of aquifer storage energy systems
.
Solar Energy
53
(
3
),
237
247
.
Kim
J.
,
Lee
Y.
,
Yoon
W.
,
Jeon
J.
,
Koo
M.
&
Keehm
Y.
2010
Numerical modeling of aquifer thermal energy storage system
.
Energy
35
(
12
),
4955
4965
.
Kreith
F.
,
Manglik
R.
&
Bohn
M.
2010
Principles of Heat Transfer
.
Cengage Learning, Inc.
,
Stamford, CT, USA
.
Lee
K.
2011
Numerical simulation on the continuous operation of an aquifer thermal energy storage system under regional groundwater flow
.
Energy Sources, Part A: Recovery, Utilization, and Environmental Effects
33
(
11
),
1018
1027
.
LeVeque
R. J.
1992
Numerical Methods for Conservation Laws
, 2nd edn.
Springer
,
Basel
.
Medici
G.
&
West
L.
2022
Review of groundwater flow and contaminant transport modelling approaches for the sherwood sandstone aquifer, UK; insights from analogous successions worldwide
.
Quarterly Journal of Engineering Geology and Hydrogeology
55
(
4
),
qjegh2021–176
.
Molson
J.
,
Frind
E.
&
Palmer
C.
1992
Thermal energy storage in an unconfined aquifer: 2. Model development, validation, and application
.
Water Resources Research
28
(
10
),
2857
2867
.
Nield
D. A.
&
Bejan
A.
2013
Convection in Porous Media
, 4th edn.
Springer
,
New York
.
Ogata
A.
&
Banks
R. B.
1961
A Solution of the Differential Equation of Longitudinal Dispersion in Porous Media
.
Tech. rep., US Government Printing Office
,
Washington, DC
,
USA
.
Philip
J.
1994
Some exact solutions of convection-diffusion and diffusion equations
.
Water Resources Research
30
(
12
),
3545
3551
.
L.
,
Hinkelmann
R.
&
Helmig
R.
2012
Modeling macroporous soils with a two-phase dual-permeability model
.
Transport in Porous Media
95
(
3
),
585
601
.
Strikwerda
J. C.
2004
Finite Difference Schemes and Partial Differential Equations
.
Society for Industrial and Applied Mathematics
,
.
Sun
Y.
,
Buscheck
T.
,
Lee
K.
,
Hao
Y.
&
James
S.
2010
Modeling thermal-hydrologic processes for a heated fractured rock system: Impact of a capillary-pressure maximum
.
Transport in Porous Media
83
(
3
),
501
523
.
Tuttle
K.
1997
Sedimentological and Hydrogeological Characterisation of A Raised ice-Contact Delta – the Preboreal Delta-Complex at Gardermoen, Southeastern Norway
.
Ph.D. thesis
,
University of Oslo
,
Oslo
,
Norway
.
Yang
C.
,
Kuwahara
F.
,
Liu
W.
&
Nakayama
A.
2011
Thermal non-equilibrium forced convective flow in an annulus filled with a porous medium
.
The Open Transport Phenomena Journal
3
,
31
39
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).