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
) 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):
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:
is a source/sink term.
, the heat content may be expressed as 
, where (see Kangas & Lund 1994; Hecht-Méndez et al. 2010):
for 
. The water flow causes the advection of heat:
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:
is a bulk aquifer heat diffusion coefficient (see Kangas & Lund 1994; Nield & Bejan 2013).
, 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:
is the thermal dispersivity length, and the total diffusion flux becomes 
(Bakr et al. 2013).
and 
, energy conservation implies that the system attains an equilibrium temperature 
equal to the weighted mean:
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:
is the size of the grains making up the solid. solid. The expression for h then becomes:
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, 
:
A similar time scale may actually be derived from the heating of spheres discussed in Gockenbach & Schmidtke (2009).
, the solid's integral conservation law reads:
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.
, and the flow q are assumed to be independent of T, then dividing through with 
in Equation (22) leads to:
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:
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.
and with initial and boundary conditions:
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
Equation (25) becomes:
. In this case, a well-known similarity variable is 
, which, inserted into the equation results in:
. 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:
. 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. 
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 
.
The 2-D radial symmetric case
we obtain an equation:
-function, defined as:
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.
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.
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).
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:
,
. 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.
so that 
, valid for all dimensions 
. In this case, Equation (40) becomes:
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:
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:
, and initial values 
and 
. The procedure is significantly simplified if v is constant, in which case the characteristics 
are just straight lines.
by defining 
for 
. In this case, we may solve Equation (41) with the boundary conditions:
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. 
The extended domain and the characteristic lines in the injection phase.
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 
.
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 .
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 
.
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 .
In order to be able to resolve a sharp front, the characteristics 
used in the discretization can be concentrated around it.
, using
,
and
. The exact solutionis given by Equation (39). The problem is first solved by the algorithm described above. The initial computationaldomain is
where
ischosen sufficiently large to avoid any influence from the boundaries. The concentration ofcharacteristics around the front is achieved by using:
,
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
it is appropriate to choose a spatial scale 
such that 
is of the size of 1. The boundary conditions in the injection case (
) is:
if the temperature is scaled.
and using the Lagrangian approach, we can find the solutions 
and 
on the characteristics 
from:
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.
andparameters:
.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:
, varyingfrom
to
andconcentrated 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. When
thefront 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.
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.
Physical and thermal properties of fluid and aquifer for the thermo-hydraulic modeling of the Gardermoen aquifer
| Property | Symbol | Value |
|---|---|---|
| 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 |
| Property | Symbol | Value |
|---|---|---|
| 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 |
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:
The equations are solved by the numerical approach outlined in Section 4.
Transient injection phase
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. 
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.
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.
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
is given by:
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):
. 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. 
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.
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.
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).
mm and 
mm as shown in 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.
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.
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.
