Dimensionless Thermal Efficiency Analysis for Aquifer Thermal Energy Storage

Seasonal warm and cold water storage in groundwater aquifers is a cost‐effective renewable energy technology for indoor heating and cooling. Simple dimensionless analytical solutions for the thermal recovery efficiency of Aquifer Thermal Energy Storage (ATES) systems are derived, subject to heat losses caused by thermal diffusion and mechanical dispersion. The analytical solutions pertain to transient pumping rates and storage durations, and multiple cycles of operation, and are applicable to various well configurations and thermal plume geometries. Heat losses to confining layers, its implications for optimizing plume geometries and aspect ratios, and heat spreading due to free convection are also discussed. This provides a general tool for broad and rapid assessment of aquifers and ATES systems. We show that if heat exchange with the confining layers is negligible, the thermal recovery efficiency of thermal plumes with cylindrical geometry is independent of the aquifer porosity and heat capacity C0. Therefore, if mechanical dispersion is negligible as is often the case, the only aquifer property that affects the recovery efficiency is the aquifer thermal conductivity. The field‐scale aquifer thermal conductivity λ can therefore be inferred from the recovery efficiency of a push‐pull heat recovery test. Remarkably, an increase in C0 could either increase, decrease, or not affect the recovery efficiency, depending on the thermal plume geometry. Hence, the analytical solutions reveal that the recovery efficiency is affected by complex interactions between the thermal plume geometry and aquifer properties. Approximate analytical solutions for subsurface temperature profiles over an entire ATES cycle are also derived.

TANG AND RIJNAARTS 10.1029/2023WR035797 2 of 28 and free convection driven by density gradients between the warm and cold groundwater (Sheldon et al., 2021).Most studies on the thermal recovery efficiency of ATES invoke numerical modeling to capture these relevant processes, but such models require detailed data on (heterogeneous) aquifer composition and properties, and water injection and extraction rate timeseries, which may be unavailable or very costly to obtain from field characterization, while also requiring intensive model construction and simulation time.Numerical studies therefore yield precise characterizations of highly specific scenarios but few analytically generalizable insights, at large computational costs.For this reason, we investigate the potential of deriving simple, dimensionless, and easily evaluable analytical descriptions of the thermal recovery efficiency of ATES systems, that apply under transient pumping rates.
Analytical solutions are of scientific interest because they reveal, more generally than empirical numerical simulations, the relationships between system parameters (e.g., aquifer properties, flow field geometry, and operational choices such as the injection rate and volume) and outcomes of interest such as the thermal recovery efficiency.Furthermore, simple dimensionless analytical solutions are easy to understand and fast to evaluate while being sufficiently accurate, which makes them useful for governance and policymaking (Marazuela et al., 2022;Pophillat, Attard, et al., 2020;Regnier et al., 2023;Serianz et al., 2022).Such analytical solutions may also be used to replace or complement computationally intensive numerical simulations, or to rapidly identify meaningful regions of parameter space for more detailed numerical investigation, especially in urban settings where multiple ATES systems are located within a small area, which is a scenario that might otherwise require highly complex numerical models (Attard et al., 2020;Burns et al., 2020;Lyden et al., 2022;Previati et al., 2022;Walch et al., 2021).Simple analytical solutions are also useful for interpreting (thermal) tracer tests to characterize aquifer properties (Del Val et al., 2021;Park & Lee, 2021;Schroth & Istok, 2005; Tang & van der Zee, 2021a), as the alternatives of solving inverse problems with numerical models or exact analytical solutions with integral transforms, are highly computationally intensive (Jung & Pruess, 2012;Regnier et al., 2022).
In this study, we present simple dimensionless analytical solutions for the recovery efficiency F of ATES systems with various well configurations and thermal plume geometries.The analytical solutions consider primarily heat spreading due to thermal diffusion and mechanical dispersion, and can be modified to consider the effects of heat exchange with aquitards, aspect ratio optimization, and free convection due to fluid density gradients.The analytical solutions apply to scenarios with three distinct phases in the ATES operational cycle: an injection phase, storage phase, and extraction phase.An analytical method to calculate the recovery efficiency of transient flow scenarios with time-varying injection and extraction rates, and scenarios with multiple cycles of ATES operation, is also developed.Approximate solutions for subsurface temperature profiles across the ATES operational cycle are also derived.This analysis yields new insight into the sensitivity of ATES recovery efficiency to various aquifer properties and ATES operational parameters.

Problem Statement
Our analysis of the thermal efficiency involves a gradual build-up of complexity.In Sections 2.2-2.4,we analyze the thermal recovery efficiency when heat losses are caused solely by thermal diffusion, for transient flow scenarios involving a distinct storage phase.In Section 2.5, the effects of mechanical dispersion are integrated into the analytical solutions.In Section 2.6, the recovery efficiency after multiple cycles of ATES operation is evaluated.In Section 2.7, subsurface temperature profiles arising from ATES operation are derived.
We also discuss the effects of heat exchange with confining layers (Section 3.2), and the effects of free convection arising from density gradients (Section 3.4).In practice, the fraction of injected heat lost to the confining layers is typically smaller than 0.04 (Shi et al., 2023;Sommer et al., 2015), which is small compared to the magnitude of heat lost laterally within the aquifer, the main subject of this study.Negligible density-driven convection occurs for most ATES systems, as 99% of ATES systems involve small temperature differentials (Fleuchaus et al., 2018).As these two processes have minor effects for most ATES systems in practice, we account for them in a relatively approximate manner, which yields simple and generalizable insight.
The initial temperature of the aquifer is assumed to be spatially uniform, and the injection temperature is assumed constant.When background flow in an aquifer is negligible compared to well-induced flow, a typical fully penetrating ATES well produces a cylindrical flow field and correspondingly, a cylindrical stored thermal plume, which is TANG AND RIJNAARTS 10.1029/2023WR035797 3 of 28 a two-dimensional radial plume.The concept of a d-dimensional radial flow field and thermal plume is not limited to d = 2.A short and partially penetrating well in a thick aquifer is akin to a point source in three-dimensional space, resulting in a spherical plume (Schroth & Istok, 2005) which is a three-dimensional radial shape.In contrast, a row of fully penetrating wells gives rise to a planar plume (Molinari & Peaudeceff, 1977;Sauty, 1977), a one-dimensional radial shape.These plume geometries are illustrated in Figure 1.
In such radial flow fields, heat transport in the aquifer is governed by the following advection-dispersion equation (ADE) (Tang & van der Zee, 2021a) where t is time, and r is radial position, c is the temperature difference between the injected water and the aquifer (positive for warm water storage, and negative for cold water storage), and k is the thermal diffusion coefficient.
Local thermal equilibrium between the solid and liquid phases is assumed, as is common practice in the ATES literature (e.g., Lin et al., 2019;Pophillat, Attard, et al., 2020).Therefore, the advection velocity of a moving heat front v(r, t) in a d-dimensional flow field is given by (2) where θ is the aquifer porosity, and δ is the thermal retardation factor, Q is the constant injection (positive) or extraction (negative) rate, Γ is the gamma function, and the term in square brackets is the surface area of a d-dimensional sphere of unit radius.We have also used the relation , where C 0 is the aquifer heat capacity, ρ w is the density of water, and s w is the specific heat of water (Molson et al., 1992).Note that the thermal advection velocity depends only on one aquifer property: the aquifer heat capacity, but not the aquifer porosity, because   =  0

𝜌𝜌 𝑤𝑤 𝑠𝑠 𝑤𝑤
. A summary of mathematical symbols and their units is provided in Table 1.

Analysis of the Recovery Efficiency
In an ATES operational cycle, warm water is injected into the aquifer at a rate of Q in and a duration of T in , then stored without any forced flow for a period of T st .Finally, the plume of stored water is extracted at a rate of Q ex and a duration of T ex .Assume that the temperature of the injected water and the background temperature are constant, and that Q in = −Q ex and T in = T ex .The recovery efficiency of the ATES system is given by where M is the fraction of injected thermal energy lost from the stored thermal plume.Now we impose an additional restriction,  st = 0 , and name this problem Scenario 1 (injection immediately followed by extraction).An approximate solution for   in this scenario is (Tang & van der Zee, 2021a) (5) The plume radius under constant injection rates is    10.1029/2023WR035797 6 of 28 where V(t) is the injected volume of water at time t.Substituting Equation 3, T in = T ex , and the maximally achieved plume radius R T = R(T in ) into Equation 5yields The goal now is to derive M when the restriction of T st = 0 is relaxed.At some time since injection t, the heat flux out of the stored plume at the interface between the plume and aquifer (Philip, 1964) is where A(t) is the plume surface area.The temperature profile is closely approximated by where the term ω(t) characterizes the width of the diffuse region of the plume and the integral in ω(t) is an integral across the history of R(t) (Gelhar & Collins, 1971;Tang & van der Zee, 2021a).Substituting    |= , A(t) ∝ R(t) d−1 , and Equation 11 into ϕ(t) yields The width of the diffuse part of the plume ω is a function of time due to the ongoing process of thermal diffusion, and is also a function of plume radius as a plume with a larger radius has a larger differential volume    ∝  −1 ∝  over which to distribute the diffuse part of the plume.Hence, ω is a function of A(t), which means that ω is sensitive to transient pumping rates at the well.However, the functional relationship between ω and A(t) may, on average across an entire ATES operational cycle, diminish under the following conditions.If the injection and extraction rate time series and volume are mirrored (i.e., Q ex (t) = −Q in (T in − t)), then during the injection phase, ω is smaller when A is small and ω is larger when A is large, whereas during the extraction phase the opposite applies (i.e., ω is larger when A is small).Therefore, for a mirrored cycle, on average ω may be roughly independent of A(t), which means that on average According to the analysis of the previous paragraph, the total heat lost from the plume ψ(t) ∝ 〈ϕ(t)〉 is approximately proportional to where the integral is taken over the history of A 2 (t).For scenario 1, this yields Now instead consider a scenario where a storage phase of  st follows an injection phase of  in , and also assume that no heat losses occur during the injection and extraction phases.We refer to this case as Scenario 2.Then, Equating Equations 14 and 15 reveals that they are identical when This result suggests that regarding the magnitude of heat mass lost from the stored plume, time spent in the injection or extraction phase counts as   3−2 times time spent in the storage phase.This is a key result that underlies many of our subsequent analyses in this study.Substituting Equation 16 into Equation 7yields for Scenario 2.
For comparison with the approximate Equation 17, exact analytical solutions can be found for M in Scenario 2. These are where the subscripts 1D, 2D, and 3D refer to planar, cylindrical, and spherical plumes respectively, and I n is the nth modified Bessel function of the first kind.Appendix A details the derivations of these three solutions.The agreement between Equations 17 and 18 is excellent for all M ≤ 0.3, regardless of the values of d or the other parameters (Figure 2a).The relative overprediction (Equation 17/ Equation 18) due to the approximation (Equation 17) rises slowly as M increases: from around 1% at M ∼ 0.3, to merely 10% at M ∼ 0.6 (cylindrical and spherical plumes) and at M ∼ 0.5 (planar plumes) (Figure 2b).We note that most ATES systems in practice have F ≥ 0.7 (i.e., M ≤ 0.3) (e.g., De Schepper et al., 2019;Mahon et al., 2022;Shi et al., 2023), well within the accurate range of Equation 17.Therefore, we later use Equation 17for further analyses involving heat losses to the confining layers, where M ≤ 0.3 is assumed (Section 3.2).

Recovery Efficiency of a Three-Phase Cycle
The scenario of interest in this section is that of a three-phase cycle where T in , T st , T ex are all non-zero.Now we pose the following empirical hypothesis: M can be accurately predicted by defining an effective time parameter T f , and substituting that for T st in the exact analytical solutions (Equation 18) and using the appropriate equation for the plume geometry in question (planar, cylindrical, spherical).This hypothesis is tested with the following 10.1029/2023WR035797 8 of 28

Water Resources Research
TANG AND RIJNAARTS 10.1029/2023WR035797 8 of 28 example.Considering that we have found that time spent in the injection or extraction phase counts as   3−2 times time spent in the storage phase, then for this scenario an obvious choice of T f is To verify the hypothesis, we performed full numerical simulations of the ADE with distinct injection, storage, and extraction periods are performed in MATLAB using the pdepe function.Verification of the solution was  18to predict M yields an excellent agreement with the numerical simulations, for both cylindrical plumes and spherical plumes (Figure 3a), with an average relative error of 0.004 and average absolute error of 0.001.
Through Equation 18, F = 1 − M is completely characterized by two parameters: T f and and can be interpreted as a characteristic time scale of thermal diffusion.It can be written in terms of aquifer properties and operational parameters, such that where λ is the thermal conductivity of the aquifer.Here, we have used  =   0 (Molson et al., 1992).It is evident that the recovery efficiency, through its dependence on , is a function of two aquifer properties: λ and C 0 .The density ρ w and specific heat s w of water are known quantities.Although ρ w is temperature-dependent, the temperature difference between injected and native water is small for typical ATES systems, such that ρ w can be assumed constant (Bloemendal & Hartog, 2018).
Remarkably, for cylindrical plumes where d = 2, but not for other plume shapes, Equation 20 simplifies to The dependence on C 0 vanishes, and the recovery efficiency now depends on only a single aquifer parameter, λ.This allows the aquifer thermal conductivity λ, known to be difficult to measure in the field and laboratory (Dalla Santa et al., 2017;Markle et al., 2006), to be directly calculated from the recovery efficiency of a thermal tracer test.Once λ is known, the value of C 0 may also be inferred from thermal tracer tests using different plume geometries or other thermal tests.

Recovery Efficiency Under Periodic Flow Rates
Equation 13 may be used to transform any transient flow scenario into an equivalent three-phase scenario, whose recovery efficiency can then be calculated with Equation 18.However, the transformation may only be accurate within certain regions of parameter space.An example of this is illustrated in this section, where we discuss a practical application of Equation 13: the case of periodic pumping rates.ATES flow rates in practice tend to be annually periodic, perturbed by small amplitude high frequency (e.g., daily or hourly) noise (e.g., Boon et al., 2019;Fleuchaus et al., 2020;Sommer et al., 2014;Visser et al., 2015), as indoor heating and cooling demand is determined by the climate, which is itself annually periodic.Therefore, consider an injection-extraction scenario with a periodic flow rate at the well, such that corresponding roughly to periodic heating and cooling demands encountered in practice.The total volume injected in the periodic flow scenario is thus Here, the steady equivalent scenario is described by effective parameters of (Q eq , T eq ), which is to be solved for under the constraint that Q eq T eq = V in .Solving for (Q eq , T eq ) is done by equating the total heat loss ψ (Equation 13) of the transient scenario and its steady equivalent.For the transient case, the plume surface area A(t) in Equation 13 can be obtained by substituting the plume volume V(t) = ∫tQ(t)dt and Equation 6 into Equation 9.For cylindrical plumes, this yields  , in = ex = per .Numerical simulations confirm that this analysis is highly accurate in predicting M for scenarios with M ≤ 0.5 (Figure 3b), and that the accuracy of this method falls off gradually for scenarios with M > 0.5.As previously mentioned, most ATES systems are operated at thermal efficiencies well above 60% (i.e., M ≤ 0.4).Therefore, our analytical method for addressing transient flow will apply in most cases in practice.
The existing literature also supports these findings.Numerical simulations by Russo et al. (2014) with daily-averaged and monthly-averaged injection rate data resulted in highly similar temperatures as simulations with the actual hourly injection data and experimental measurements.Simulations of thermal plume development across a heating season using long-term averaged injection rates by Pophillat, Bayer, et al. (2020) resulted in similar thermal plume shapes as simulations with monthly injection rates, when background groundwater flow rates were small.Therefore, further analyses of ATES systems may be performed with the simplification of periodic injection-extraction rates to a constant value.

Mechanical Dispersion
In practice, aside from thermal diffusion, thermal spreading and heat losses may also occur due to other dispersive processes, such as free convection, macrodispersion, and mechanical dispersion.Sauty et al. (1978) showed that the effects of mechanical dispersion on the thermal efficiency can be modeled implicitly through the use of an effective thermal conductivity in the form of where λ α is an apparent thermal conductivity arising from the effects of mechanical dispersion.By performing extensive numerical simulations and fitting to the recovery efficiency, Doughty et al. (1982) found that applies to cylindrical plumes.This empirical expression was tested by Tsang et al. (1981), and was found to describe field data and further numerical simulations well.We divide this equation throughout by C 0 , and substitute k = λ/C 0 , which yields where k eff is the effective thermal diffusivity and k α is the apparent thermal diffusivity due to mechanical dispersion.
Now we derive k α analytically, through steps detailed in Appendix B. For cylindrical plumes, as an example, we obtain We performed numerical simulations of the ADE for cylindrical plumes that explicitly consider mechanical dispersion to verify the validity of using Equation 27to account for mechanical dispersion (Figure 4a).We consider scenarios with T in = 10, T ex = 10, and all combinations of k = {400, 1,000, 2,000}, α = {1, 10, 20, 50}, and T st = {0, 10, 20, 30} that fulfill the criterion     ≤ 2 .These combinations of parameter values resulted in scenarios with 0.55 ≤ F ≤ 0.9, in line with real ATES systems (e.g., De Schepper et al., 2019;Mahon et al., 2022;Shi et al., 2023).Figure 4a shows that modeling mechanical dispersion by substituting k eff = k + k α into Equation 18 leads to highly accurate predictions of the thermal recovery efficiency.27) and k eff = k + k α0 (Equation 30), as a function of k α0 and T st , when T in = 10 and k = 1,000.If k α0 or T st are small, then σ → 0 and k α → k α0 , which yields which is identical to the empirically-derived k eff (Equation 26), aside from the magnitude of the coefficient 2/3.The coefficient of 2/3 we analytically derived for the contribution of mechanical dispersion to k eff , is larger than the 0.3 found by Doughty et al. (1982) from numerical simulation data.This is because Doughty et al.'s numerical simulations also included the effects of thermal diffusion to the hydraulically impermeable confining layers, which decreases the relative contribution of mechanical dispersion to overall heat losses.We have therefore provided a physical explanation for Doughty et al.'s "effective thermal conductivity," which has till now been empirical in nature.As macrodispersion due to aquifer heterogeneity may be considered a form of mechanical dispersion with a space-dependent mechanical dispersion coefficient, it would also be possible to include the effects of macrodispersion into k eff , if the macrodispersion parameters are known (Tang & van der Zee, 2022).
For the same set of scenarios as in Figure 4a, we show in Figure 4b  become more accurate as k eff /k becomes smaller, and as T st becomes smaller.For the scenarios with T st = 10 (T st = 0), the mean and maximum relative errors of predicting M with k eff were 2% (1%) and 6% (1%) respectively, whereas the mean and maximum relative errors of ignoring mechanical dispersion were 14% (13%) and 26% (29%) respectively.Here, the relative error refers to |1 − M(analytical)/M(numerical)|.
The error associated with using the simplified Equation 30 instead of the full solution Equation 27 is illustrated in Figure 4c.Here, it is evident that the error is negligible for T st /T in ≤ 1 and k α0 /k < 0.2.Since mechanical dispersion is usually weak (i.e.,   0  ≤ 0.2 ) for ATES scenarios (Bloemendal & Hartog, 2018), and since T st ≤ T in typically (Fleuchaus et al., 2020), using Equation 30 to model mechanical diffusion is suitable for ATES systems in general.The validity of the simplified Equation 30, which we have shown, is important because unlike Equation 27, the simpler Equation 30 is independent of T st .This could mean one less parameter to consider in a sensitivity analysis or inverse problem, therefore reducing data and computational requirements.
In Sheldon et al. (2021), Equation 26(following Doughty et al.'s (1982) empirical observations) was applied to model mechanical dispersion in two scenarios with different relative storage durations: one with T in = T ex = T st , and another with T in = T ex = 2T st .However, Doughty et al. (1982) performed numerical simulations to fit Equation 26 only for the case of T in = T ex = T st .Therefore, the use of the empirical Equation 26 to model scenarios with T in = T ex = 2T st does not necessarily follow from Doughty et al.'s analysis.In fact, it is to be expected that the effect of mechanical dispersion should be weaker when T st is larger, as no mechanical dispersion occurs during the storage phase.In this study we have shown this to be true through analytical solutions, which implies that it is indeed not physically appropriate to use the same k α for various values of T st .However, the error of doing so is small when T st and k α are small.

Multiple Cycles
When ATES systems are operated for multiple cycles, the recovery efficiency increases with the number of cycles, according to the expression where n refers to the nth cycle, F n = 1 − M n is the recovery efficiency of the nth cycle, and z is an empirical parameter that is negatively correlated with F and analytically bounded in the range  1 2 ≤  ≤ 1 (Tang & van der Zee, 2021a).Extensive numerical simulations show that although the value of z is highly positively correlated with M (Tang & van der Zee, 2021a, 2022).However, the value of z is also affected by the duration of the rest phase, which refers to the period of time between an extraction phase and the injection phase of the subsequent cycle.Nevertheless, if the duration of the rest phase is fixed, then z increases as F decreases.In other words, an increase in  ∕ 2  will cause z to increase.An analysis of a limiting case with interesting behavior is given below.

Water Resources Research
TANG AND RIJNAARTS 10.1029/2023WR035797 14 of 28 Consider a four phase system, in which a rest phase of duration T rs follows the standard three-phase scenario discussed in Section 2.3.As an example, we will consider T in = T st = T ex = T rs here, which roughly divides the year into four periods.After each rest phase (e.g., spring), subsequent cycles begin once again with an injection phase (e.g., summer).Some pilot numerical simulations show that as long as the relative durations of the four phases do not change, then M and z remain constant as T in is varied.The reason that M remains constant follows analytically from Equation 17: M is a function of the dimensionless ratio   ∕ 2  , which is constant (only for cylindrical plumes) under changes in T in , if the ratio T in :T st :T ex :T rs is fixed.The pilot numerical simulations therefore show that the empirical parameter z also remains constant under changes in T in in such cases.This suggests that within each climate zone, where the ratio T in :T st :T ex :T rs would not be expected to vary significantly across multiple years, a unique functional relationship exists between z, M, and  ∕ 2  for the relevant T in :T st :T ex :T rs ratio.We performed a large set of numerical simulations that vary  ∕ 2  (by randomly varying k and Q in to yield 140 combinations) to obtain this unique functional relationship for the case of T in = T st = T ex = T rs .The results show that z is uniquely determined by M, with an empirical relationship that yields which applies to all scenarios with T in = T st = T ex = T rs (Figure 5).This reveals that a scenario with large M will also have an M n that decreases more slowly as the number of cycles grows, in a logarithmic sense, because This result emphasizes the need to optimally choose suitable aquifers and operational parameters for ATES at the outset, as the consequences of suboptimal choices stick even after many cycles of operation.Physically, this is because an aquifer with a low first-cycle recovery efficiency would also lose remnant heat more quickly during the rest phase and during subsequent cycles.

Well and Aquifer Temperature Profiles
For the purposes of making coarse estimates and sensitivity analyses of heat transport processes, approximate dimensionless analytical solutions for the temperature profile in the aquifer or at the wellbore may be useful (Attard et al., 2020;Guimerà et al., 2007;Pophillat, Attard, et al., 2020;Schroth & Istok, 2005).The temperature profiles resulting from ATES operation are also important for assessing whether anthropogenic heat changes in the subsurface result in environmental and ecological damage (Blum et al., 2021).From our result that time spent in an injection or extraction phase counts as   3−2 times time spent in a storage phase for heat loss, analytical solutions for the (a) temperature profile in the aquifer during the injection, storage, and extraction phases, (b) and the production temperature at the well, and (c) temperature evolution at the well during the storage phase may be derived, for all three radial plume geometries (planar, cylindrical, radial).Analytical solutions for these problems that consider successive injection, storage and extraction phases are new to the literature as far as we are aware.We provide analytical solutions for the temperature profile of (a) planar, (b) cylindrical, and (c) spherical plumes, that apply to the (a) injection phase, (b) storage phase after injection, (c) extraction phase after injection and storage.These yield 3 × 3 = 9 combinations of cases.Of these 9 combinations, as far as we are aware, simple analytical solutions that do not involve integral transforms or multidimensional numerical integration  (Brau et al., 2017).
The analytical solutions for the temperature are listed and derived in Appendix C. Figure 6 shows that these approximate analytical solutions may agree excellently with numerical simulations.However, additional testing shows that the agreement between the analytical and numerical solutions tend to worsen as the injection rate decreases, R T decreases, or k increases.This is evident in Figure 6c, where the analytical curves corresponding to large k have a large error at intermediate time scales, although at small and large time scales they converge to the "ground truth" numerical solutions.Since these temperature profiles are not the main focus of this study, we present these analytical approximations as-is without further discussion.A more in-depth investigation into the parameter ranges in which these solutions are valid could be a possible avenue for further research.
The fact that the analytical temperature evolution at the well during a storage phase (Figure 6c) always converges toward the numerical "ground truth" suggests that the solutions (Equations C9 and C10) may be used to solve inverse problems to infer the aquifer properties C 0 , λ (Schroth & Istok, 2005).Equation C9 for cylindrical plumes is independent of C 0 for the same reasons discussed in Section 2.3, which means that with cylindrical plumes, it is possible to estimate the thermal conductivity of the aquifer by simply observing the temperature evolution at the well.Furthermore, in Hermans et al. (2019), it was shown that partial (statistical) knowledge of aquifer properties could be used to statistically predict subsurface temperatures accurately during the extraction phase, but less accurately during the storage phase.This in turn means that storage phase temperatures are more sensitive to aquifer parameter values.Therefore, aquifer parameter values inversely fitted to observations of storage phase temperatures would be more accurate than those fitted to observations in other phases.
Recall that in the analytical solutions for the recovery efficiency discussed in Section 2.3, it was assumed that T ex = V in which leads to V ex = V in .If T ex ≠ T in , then the recovery efficiency can be calculated by integrating the approximate solutions for the production temperature (see Appendix C) over time (Doughty et al., 1982) Therefore, the analytical solutions for the production temperature allow for simple analytical sensitivity analyses of the effects of varying T ex , and other model parameters, on F.

Analytical Insights Into the Influence of Aquifer Properties on ATES Thermal Efficiency
This analytical result that C 0 is irrelevant to the recovery efficiency for cylindrical plumes (Equation 21), and that the only aquifer property that matters is λ, appears to be new to the literature, and agrees with the numerical simulation data of Shi et al. (2023).Shi et al. performed a numerical sensitivity analysis to relate the recovery efficiency to C 0 , and found that doubling or halving C 0 only changed M by ±0.01.The effect these authors found is not exactly zero, as they also simulated heat losses to the confining layers, which does depend on C 0 ; this will be discussed in Section 3.2.In their numerical sensitivity analysis, Shi et al. also found that varying the porosity θ between 0.1 and 0.3 only affected M by ±0.01, which we have showed in Equation 3 to be inconsequential for ATES thermal efficiency (see also Heldt et al., 2023).Once again, the effect Shi et al. found was not exactly null, due to heat losses to the confining layers.Therefore, the simulation data of Shi et al. (2023) confirm our findings that M is independent of C 0 for cylindrical plumes, and independent of θ for all studied plume geometries.
Equation 20 shows that for planar plumes, the recovery efficiency decreases as C 0 increases (   ∝  1∕2 0 ), whereas for spherical plumes, the recovery efficiency increases as ).Altogether, this analysis shows that the effects of the aquifer heat capacity C 0 on the thermal recovery efficiency is highly dependent on the geometry of the thermal plume, not only in magnitude, but also in a qualitative sense (i.e., positive or negative).Note that this also applies to the thermal retardation coefficient δ, as δ ∝ C 0 .This result also appears to be novel in the literature.Therefore, for problems involving heat storage and transport in porous media in general, it is not possible to generalize how varying C 0 will affect heat transport, without first considering the geometry of the thermal plume.
As the aquifer heat capacity C 0 is a volumetrically weighted average of the heat capacity of the fluid and solid phases (Molson et al., 1992), one could intentionally alter the heat capacity of either phase to achieve some

Water Resources Research
TANG AND RIJNAARTS 10.1029/2023WR035797 17 of 28 desired heat transport behavior.This also means that when calculating the thermal recovery efficiency of ATES, the chemical composition of the injected and native water may require increased scrutiny.For example, van Lopik et al. ( 2022) proposed an ATES system design in which saline hot water is injected into a freshwater aquifer using partially-penetrating wells (i.e., d ≈ 3 plumes).The advantage of their proposed system is that heat loss due to free convection is suppressed, because the increased salinity of the injected water (i.e., increased density) compensates for its increased temperature (i.e., decreased density).Although van Lopik et al.'s numerical simulations show that their proposed system is effective in counteracting heat losses from free convection, they did not take into account the decreased heat capacity of saline water.This effect may reduce the thermal efficiency of their proposed system, though the extent is likely small as   ∝  −1∕6 0 .

Heat Losses to Confining Layers
Commonly, the confining layers have the same thermal conductivity and specific heat as the aquifer matrix (Doughty et al., 1982).In this case, for a cylindrical plume, Rubinshtein (1959) shows that over the injection phase, the heat lost to the confining layers is exactly where H is the height of the aquifer.In the derivation of this solution, finite vertical heat diffusion within the aquifer, and finite horizontal diffusion in both the aquifer and confining layers was explicitly considered.Rubinshtein's solution has been shown by Spillette (1965) to agree well with full numerical simulations of the problem of warm water injection into a confined aquifer, and with experimental data.If we instead assume that heat losses to the confining layer occur only during the storage phase, then M b is approximately characterized by the solution of the one-dimensional heat equation (Equation 17with d = 1).Taking into account both upwards and downwards heat losses, It can be assumed that heat fluxes in the caprock and bedrock are primarily in the vertical direction (Li et al., 2010), if within the aquifer, advective heat fluxes are larger than dispersive fluxes (Yortsos & Gavalas, 1982).Tang and van der Zee (2021a) showed that this is strongly valid for scenarios with M ≤ 0.3, which applies to most ATES systems (e.g., De Schepper et al., 2019;Mahon et al., 2022;Shi et al., 2023).If this assumption holds, the diffuse outer edges do not play a significant role in thermal exchange with the confining layers, thus for planar and cylindrical plumes, A b = V/H are identical for a given aquifer thickness H, where A b is the interfacial area for thermal exchange between the aquifer and the confining layers.Hence, M b (t) for cylindrical and planar plumes can be considered identical, and Equations 36 and 37 may be applied for both cylindrical and planar plumes.In any case, in practice, the magnitude of heat losses to confining layers is 0 ≤ M b < 0.04, an order of magnitude lower than M, for both planar and cylindrical plumes (Shi et al., 2023;Sommer et al., 2015).Hence, an approximate characterization of M b is sufficient.
An extensive numerical sensitivity analysis of scenarios with T in = T st = T ex shows that during the extraction phase, the total heat loss to the confining layers during the injection phase is much larger than that of the storage phase for all scenarios (Shi et al., 2023).Furthermore, they also show that during the extraction phase, back-diffusion of heat from the confining layers to the aquifer causes a net heat gain in the stored thermal plume, instead of a net heat loss.In all of Shi et al.'s simulated cases, the net heat gain during the extraction phase is around 50%-80% of the total heat loss during the storage phase.Of all parameters considered in Shi et al.'s sensitivity analysis, the aquifer-aquitard heat fluxes were most sensitive to the aquifer height, yet the 50%-80% recapture ratio hardly changed.This implies that heat losses to the confining layers during the storage phase are small, and much of the heat lost during the storage phase is recaptured during the extraction phase.Hence, heat losses during the injection phase determine to a large extent the overall net heat losses to the confining layers of an ATES system.
For a fixed injected volume of water, the thermal efficiency of an ATES system is optimized at a certain optimal aspect ratio     of plume radius to height.For fully penetrating wells (cylindrical plumes) or rows of wells (planar plumes), the plume height is fixed by the aquifer thickness, which is determined by the local hydrogeological TANG AND RIJNAARTS 10.1029/2023WR035797 18 of 28

Water Resources Research
TANG AND RIJNAARTS 10.1029/2023WR035797 18 of 28 properties of the subsurface.Therefore, as H is not a degree of freedom, optimizing the aspect ratio for a fixed injected volume of water essentially means optimizing R T .For point sources (spherical plumes), the aspect ratio is not a meaningful concept, as spherical plumes arise only when no thermal interaction with the confining layers occur.The optimal plume aspect ratio is a subject of interest in the design of ATES systems (Collignon et al., 2020;Doughty et al., 1982;Shi et al., 2023) and the evaluation of an aquifer's suitability for ATES implementation (Bloemendal & Hartog, 2018;Gao et al., 2019).Given that net heat loss to the confining layers occurs primarily during the injection phase, we use Equation 36 to calculate the optimal aspect ratio.This is in contrast to using Equation 37, as was previously done by Doughty et al. (1982).
For a cylindrical plume, for the injection phase only, = 0 to obtain the optimal aspect ratio, and using √ 2 8 ≈ 0.53 .Repeating this for the storage phase yields a previously known result     = 1 2 (Doughty et al., 1982).Hence, the optimal aspect ratio for cylindrical plumes depends on the relative durations of the injection and storage phases, and falls within the narrow range  0.5 ≤    ≤ 0.53 .As total net heat losses are mostly determined by heat losses during the injection phase, the actual optimal aspect ratio should be closer to the injection phase result (0.53) than to the previously available storage phase result (0.5).The difference between 0.53 and 0.5 is small, so it does not really matter which value is used.However, this is not the case for planar plumes.
The optimal aspect ratio is now derived for a planar plume, where V = 2RH.The optimal aspect ratio for the injection phase corresponds to = 0 , which yields     = 3 4 .During the storage phase, a similar derivation yields an optimal aspect ratio equal to     = 1 2 .Therefore, the optimal aspect ratio for planar plumes falls within the relatively broad range 4 , but is closer to the injection phase result  3 4 as most net heat losses to the confining layers occurs during the injection phase.This suggests that the optimal aspect ratio of a planar plume is significantly wider (i.e., larger R) than that of a cylindrical plume.As aquifer thicknesses (i.e., H) are fixed, this is another factor to consider in the design and placement of ATES wells, especially in densely built-up areas.
If the aquifer and confining layer thermal conductivities and heat capacities are different, the aspect ratio can be estimated by multiplying the reference ratios (0.53 for cylindrical plumes, 0.75 for planar plumes) by (Doughty et al., 1982) 1 + where λ b and C b are the thermal conductivity and the heat capacity of the confining layers.Using the values of λ and C 0 Shi et al. (2023) simulate in their sensitivity analysis (their Sections 4.3 and 4.4 respectively), we estimate that the optimal aspect ratio of their simulations will range between  0.54 ≤    ≤ 0.84 .By inspecting their numerical simulation data, Shi et al. concluded that the optimal aspect ratio is between 0.5 and 1, for cylindrical plumes, in agreement with our analysis.Similarly, using this method we predict that the numerical simulations of Doughty et al. (1982) will yield optimal aspect ratios between  0.52 ≤    ≤ 0.93 , in agreement with their data.Therefore, we have provided analytical insights into the optimal aspect ratio, which has previously been studied using numerical simulations.

Optimizing Thermal Plume Geometry and Well Configuration
The influence of plume geometry on ATES thermal efficiency has been previously studied.Sommer et al. (2015) showed, using numerical simulations, that ATES plume zonation by alternating planar warm and cold plumes ("lane pattern") leads to a higher recovery efficiency than alternating cylindrical warm and cold plumes ("checkerboard pattern").Using the analytical methods introduced in this study, we can arrive at a similar conclusion.
To be specific, consider a scenario where the position of wells in a densely utilized aquifer is predetermined, with a fixed total injected volume V in , as in the scenarios studied by Sommer et al. (2015).If the subsurface is densely packed with thermal plumes, Sommer et al.'s lane pattern would be analogous to d = 1 and their checkerboard pattern would be analogous to d = 2, in the context of this study.Since V in is fixed, it is not necessary to compare vertical heat losses, as M b are identical for both patterns.From Equation 17, we find that horizontal heat losses Such findings on the effect of dimensionality in densely utilized aquifers is also supported by other literature.Duijff et al. (2021) found that when the plumes of multiple warm (or cold) wells are combined into a single large plume, the overall recovery efficiency of the cluster of ATES systems will increase.Therefore, the notion of reducing heat losses by using more compact plume shapes is already present in the literature, arising through empirical observations from numerical simulations and field data.The relation between plume compactness and dimensionality, and the magnitude of this effect, are analytically quantified for the first time in this study.
The findings of this study suggest that planar plumes are more efficient than cylindrical plumes, which are typically used.First, a larger space-use efficiency is achievable with planar plumes than cylindrical plumes, because three-dimensional space cannot be efficiently packed with cylinders.Furthermore, the same heat storage is achievable with a larger distance between opposite-type wells as planar plumes have a larger optimal aspect ratio, allowing each individual ATES operator to use fewer wells to store the same amount of water, which lowers capital costs.In addition, in some aquifers background groundwater flow may be responsible for a large yet difficult to estimate effect on the thermal efficiency of ATES systems (Bloemendal & Hartog, 2018;Pophillat, Attard, et al., 2020).Another advantage of planar plumes is therefore that is rows of wells of each type may be intentionally oriented parallel to the mean background flow direction to improve recovery efficiencies.This also allows the effects of regional groundwater flow on recovery and mutual interference, if applicable, to be more easily calculated using simple geometric methods (e.g., Bloemendal & Hartog, 2018).
Urban planning policies that allow the positions of ATES systems to self-organize due to market and social factors, such as proposed by Bloemendal et al. (2014), may not lead to optimal use of subsurface space, because the tendency appears to be for wells of each type to cluster in a radial rather than linear fashion (Beernink et al., 2022).History shows that decentralized urban expansion tends strongly toward concentric radial shapes, regardless of initial configurations, whereas linear shapes may arise and persist only under strict and highly centralized control (Batty, 2022).The findings of Beernink et al. (2022) suggest that the same appears to be true for ATES well positioning, but in any case, ATES wells cannot be placed beneath existing buildings, so ATES well positioning is necessarily determined by the present urban landscape.Beernink et al.'s (2022) study of free-market determined ATES positioning show that when the minimum permissible distance between same-type wells (warm vs. cold) increases, it becomes more likely for planar plumes to arise.However, Beernink et al. (2022) also shows that when this minimum permissible distance is decreased, denser well placement is encouraged, which would lead to the desirable outcome of more intensive ATES implementation.Therefore, allowing decentralized market-driven ATES positioning would lead to an inevitable tradeoff between thermal efficiency and ATES adoption.This is strongly in contrast to the recommended alternative where the subsurface is centrally planned and zoned into warm and cold lanes that allow planar plumes to be implemented.In this case, the thermal efficiency increases with ATES adoption, as the thermal efficiency intrinsically increases with the stored volume (Equation 17), even if some negative thermal interference occurs (Sommer et al., 2015).In other words, market-driven ATES positioning tends toward local (sub-)optima (in the language of mathematical optimization) in the reduction of non-renewable energy dependence.A higher-level global optimum is more likely achievable under strictly controlled ATES planning that zones the subsurface into distinct hot and cold lanes.

Effects of Free Convection on High Temperature ATES Systems
For high temperature ATES (HT-ATES), thermal efficiencies are generally lower than for regular low temperature ATES systems.This is due to heat transport by free convection, arising from large density differences between the warm and cold water plumes.Schout et al. (2014) showed that the thermal recovery efficiency for HT-ATES is (in terms of the parameters used in this study) given by where B and F emp are empirical parameters, G is the hydraulic permeability of the aquifer, and Δγ is the difference in specific gravity between the injected water and native groundwater.Also, Δγ is equal to the temperature difference between injected and native groundwater multiplied by α f , the coefficient of thermal expansion of water.Schout et al. (2014) and Sheldon et al. (2021) showed that this empirical equation is flexible: it is valid to describe the recovery efficiency after multiple operational cycles, that also include storage and rest periods.In the context of our study, this implies that .If H is large, such that heat losses to the confining layers are negligible, then F emp → F, where F = 1 − M is the thermal efficiency in that we have analytically derived in Equation 18.Our analytical solution of F derived in this study therefore provides a physical basis for directly calculating F emp when heat losses to the confining layers is small, without the need for empirical fitting.
It has also been previously reported that when density-driven flow is significant, the effects of free convection on heat spreading can be accounted for with an analogous "thermal diffusivity" parameter (Jalbert et al., 2000) FC = 1 5 Δ (40) Note that the original formulation in Jalbert et al. (related to saltwater transport) has a denominator of θ, which vanishes in our Equation 40 as the porosity and thermal retardation factor cancel out for heat convection, as previously discussed.We therefore propose a modified thermal diffusivity k HT to characterize the thermal efficiency of high-temperature ATES: If the effective thermal diffusivity is calculated using Equation 41, then the thermal recovery efficiency could be estimated without having to first fit the empirical parameter B from field or numerical data.A flaw of the analytical Equation 41 is that it does not consider heat losses to the confining layers, and interaction effects between free convection and thermal diffusion, whereas this is accounted for through the empirical coefficient B. Since free convection tilts the plume front, it has the effect of changing the plume-aquifer and plume-aquitard interfacial areas.We thus hypothesize that using k HT to model the additional heat losses due to free convection should be most accurate when the thermal recovery efficiency is least affected by perturbations to the interfacial areas: the discussion in Section 3.2 suggests that for cylindrical plumes, this is when the optimal aspect ratio  = 5, = 0, emp = exp − 2.17 − 0.12 , = − 1.99 2 + 0.608 − 0.00673 in = st = ex = 90 d, = 3 W∕mK, 0 = 2.21 ⋅ 10 −6 J∕m 3 K, = 2 ⋅ 10 −6 K −1 , = 5, = 0, emp = exp − 2.17 − 0.12 , = − 1.99 2 + 0.608 − 0.00673 ).We perform the comparison for two values of the hydraulic conductivity that differ by one order of magnitude: G = 1 • 10 −4 m/s and G = 1 • 10 −5 m/s, native groundwater temperature of 10 K, and injected water temperatures of 10-60 K. Sheldon et al.'s empirical equation is taken to be the "ground truth" for this comparison, as their predictions yield a coefficient of determination (r-squared) of 0.88 against a large data set of numerical simulations.
Figure 7 shows the resulting plots of M HT − M for both methods.The accuracy of the analytical method is good, especially considering that it was not necessary to fit an empirical parameter to a-priori data, and that the functional dependence of M HT − M on GΔγ is recovered well.
Since we have shown that Equation 39 and our analytical method has a similar functional dependence on GΔγ, it means that when the hydraulic conductivity is anisotropic, G can be replaced by the geometric mean of the horizontal and vertical hydraulic conductivities in our analytical Equation 41, because Schout et al. (2014) and Sheldon et al. (2021) have shown that this is valid for Equation 39.Furthermore, Sheldon et al. (2021) remark that the empirical parameter B also depends significantly on the cycle duration, mechanical dispersivity, and cycle number, and that B needs to be re-fitted if these parameters vary.In contrast, these aspects can be integrated into our analytical method directly through the procedures outlined in Sections 2.2-2.6.
Therefore, Equation 41 is useful for understanding the effects of free convection on HT-ATES, especially after further research reveals under what exact circumstances k HT is valid.This will require significant work due to the large number of parameters associated with heat losses due to free convection, and its interactions with heat exchange with the confining layers.For example, we also attempted to compare the analytical and empirical methods for scenarios that do not fulfill     = , and found that the agreement between both methods was not as good on average.However, it is possible that other unexplored regions of parameter space will yield a better agreement for such scenarios.Such an exhaustive search in parameter space was not possible in this study as the values of the empirical parameter B is only reported in the literature for a limited number of scenarios.

Conclusions
The thermal recovery efficiency of an ATES system is given by substituting Equation 19 into the appropriate form (depending on plume geometry) of Equation 18.If applicable, transient pumping is considered by using the transform procedure outlined in Sections 2.2-2.4.Other processes that affect heat spreading can be integrated into Equation 18 through the methods described in Section 2.5 (mechanical dispersion), Section 2.6 (multiple cycles of operation), Section 3.2 (heat losses to confining layers), and Section 3.4 (free convection).For scenarios with high recovery efficiencies F ≥ 0.7, the simpler Equation 17 can be used instead of Equation 18.All three phases of ATES operation: injection, storage and extraction, are fully accounted for.The simplicity of the analytical solutions allow for extensive physical insight into ATES system performance, that were previously observed through empirical methods.They also allow radial heat losses to be analytically compared against vertical heat losses to the confining layers, to derive optimal thermal plume aspect ratios.The synthesis of this study suggests that planar plumes are optimal for both thermal and spatial efficiency, especially in regions where the subsurface space is densely utilized, in agreement with the literature.However, foresight and planning is necessary for this to materialize.The key innovations of this study are: 1. Dimensionless analytical solutions of the thermal recovery efficiency for a complete injection-storage-extraction cycle, for planar, cylindrical and planar plumes.The combination of all three phases in an analytical solution is new to the literature.2. Analytical solutions for the thermal recovery efficiency that combine the effects of thermal diffusion and mechanical dispersion, which is new to the literature.3.An analytical method to transform transient pumping scenarios to steady pumping scenarios, which can then be characterized using these analytical solutions.Analytical methods to characterize the effects of transient flow on thermal diffusion for ATES are new to the literature.This allows for an analytical evaluation of ATES systems with periodic (e.g., sinusoidal) pumping rates, which is applicable in practice due to the periodic nature of indoor heating and cooling demands.4. The analytical result that if heat spreading is fully controlled by thermal diffusion, then the recovery efficiency of an ATES cycle is independent of the aquifer heat capacity, if and only if the thermal plume has a cylindrical geometry.This allows the aquifer thermal conductivity to be inferred from a simple push-pull thermal tracer test.The thermal conductivity can also be inferred by monitoring the temperature at the well during the storage phase.5.The dependence of the thermal recovery efficiency on the heat capacity of the aquifer is complex, both in terms of magnitude and directionality (i.e., positive or negative), because it is determined by the geometry of the stored heat plume.Hence, the suitability of a particular aquifer for ATES implementation is dependent on the design and configuration of the wells, which affects the plume geometry.

Figure 2 .
Figure 2. (a) A comparison of the general approximation of M, Equation 17, against the exact solutions for a storage phase with a planar (Equation18a), cylindrical (Equation18b), and spherical plume (Equation18c).The black dotted line demarcates M = 0.3, beneath which the approximation is essentially identical to the exact solutions.(b) Error of using Equation 17 to approximate Equation18.The black dotted line demarcates 10% relative error.

Figure 3 .T
Figure 3. Analytical predictions against numerical predictions for M in (a) three-phase scenarios, and (b) periodic flow rate scenarios.

Figure 4 .
Figure 4. (a) Accuracy of substituting k eff = k + k α (Equation 27) into Equation 18 to analytically predict M subject to mechanical dispersion in cylindrical plumes, compared to explicit numerical simulations.(b) The open circles show the accuracy when the simplified analytical method  eff =  + 0 =  + 2 3 (    in ) is used, and the solid crosses show the prediction error if mechanical dispersion is altogether ignored, for the same set of scenarios as in (a).(c) Contour plot of the magnitude of the difference in M calculated with k eff = k + k α (Equation27) and k eff = k + k α0 (Equation30), as a function of k α0 and T st , when T in = 10 and k = 1,000.

Figure 5 .
Figure 5. Empirical relationship between z and M discussed in Section 2.6.Each point marks one of the 140 simulated scenarios with varying Q in and k.

Figure 6 .
Figure 6.Analytical solutions (markers) and numerical simulations (lines) describing cylindrical (d = 2) and spherical (d = 3) plume scenarios with T in = T st = T ex = 10.(a) Temperature profiles (b) Production temperature at the well during the extraction phase.Dimensionless parameters used are T in = 10, T ex = 10,        0 = 200  ⋅ 10 3(−1) , k = 800, unless otherwise specified.(c) Analytical solutions (dashed lines) and numerical simulations (solid lines) of the temperature at the well for a cylindrical plume during a storage phase of duration T st = 40, preceded by an injection phase of duration T in = 10.From top to bottom in (c), each line corresponds to an increasing value of k, ranging from k = 100 to k = 3,000.

M
increase as d increases.Thermal interference between opposite-type (warm-cold) wells, if taken into account as Sommer et al. did, will only exacerbate the difference in recovery efficiency between planar and cylindrical plumes.Therefore, using lane patterns, which have lower dimensional radial geometry, would lead to decreased heat loss fractions under otherwise identical parameters.