Numerical simulation of aquifer thermal energy storage using surface-based geologic modelling and dynamic mesh optimisation

Aquifer thermal energy storage (ATES) has significant potential to provide largescale seasonal cooling and heating in the built environment, offering a low-carbon alternative to fossil fuels. To deliver safe and sustainable ATES deployments, accurate numerical modelling tools must be used to predict flow and heat transport in the targeted aquifers. This paper presents a simulation methodology for ATES based on surface-based geologic modelling (SBGM) and dynamic mesh optimisation (DMO). DMO has been previously applied in other fields of computational fluid dynamics to reduce the cost of numerical simulations. DMO allows the resolution of the mesh to vary during a simulation to satisfy a user-defined solution precision for selected fields, refining where the solution fields are complex and coarsening elsewhere. SBGM allows accurate representation of complex geological heterogeneity and efficient application of DMO. The paper reports the first systematic convergence study for ATES simulations, and demonstrates the application of these methods in two ATES scenarios: a homogeneous aquifer, and a realistic heterogeneous fluvial aquifer containing meandering, channelised sand bodies separated by mudstones. It is demonstrated that DMO reduces the required number of mesh elements by a factor of up to 22 and simulation time by a factor of up to 15, whilst maintaining the same accuracy as an equivalent fixed mesh. DMO offers significant potential to reduce the computational cost of ATES simulations in both homogeneous and heterogeneous aquifers.


Introduction
Heating and cooling of buildings account for half of the world's energy consumption, of which 75% is produced from fossil fuels (Fleuchaus et al. 2018). Cooling demand is forecast to increase in the coming decades as a consequence of climate change (Stephenson et al. 2019). To reduce carbon emissions and transition to a net zero society, it is essential to adopt low carbon technologies for the heating and cooling of buildings. Aquifer thermal energy storage (ATES) has significant potential to supply large-scale space heating and cooling, providing a sustainable alternative to fossil fuels (Bloemendal et al. 2015). The operating principle of ATES is illustrated in Fig. 1.
Aquifer thermal energy storage installations are associated with a range of induced thermal, hydraulic, chemical and mechanical impacts within the storage aquifer which must be taken into account when planning, installing and operating ATES systems (Bauer et al. 2013;Fleuchaus et al. 2018). To develop safe and sustainable installations, numerical modelling must be used to predict these processes and induced impacts.
Numerical codes typically used to simulate ATES can be broadly separated into two categories: those that use k-orthogonal structured or unstructured meshes and those that use nonk-orthogonal unstructured meshes. In the case of k-orthogonal meshes, the governing equations are discretised using the finite difference method (FDM) or the finite volume method (FVM). Both these schemes require k-orthogonal meshes to accurately calculate fluxes between cells using the two-point flux approximation. Simulators of this type that have been previously applied to ATES include MODFLOW-MT3DMS (Bloemendal and Olsthoorn 2018), UTCHEM (Lee 2011), TOUGH2 (Gao et al. 2017) and HSTwin (Sommer et al. 2014). The use of k-orthogonal grids, usually pillar grids where pillars extend from the base to the top of the models, can be computationally expensive-a large number of cells or elements may be required to accurately represent geological heterogeneity, as all layers in the model have the same areal grid resolution. This can be particularly expensive when modelling large domains in aquifers that are heterogeneous across a range of lengthscales. The use of k-orthogonal grids also causes problems when modelling heterogeneities which do not follow the grid orientation, introducing stair-stepping effects and corresponding artefacts in the numerical solution. Such heterogeneities include common geological features such as dipping faults, meandering channels or pinched-out layers (Jackson et al. 2015).
Simulators which use non-k-orthogonal meshes can alleviate these issues. In these simulators, the governing equations are discretised using the finite element method (FEM) or hybrid FVM-FEM schemes that ensure local mass conservation. Common FEM-based simulators applied to ATES include COMSOL (Gao et al. 2017), FEFLOW (Winterleitner et al. 2018) and OpenGeoSys (Kolditz et al. 2012;Wang et al. 2021) and the hybrid FEM-FVM simulator CSMP (Yapparova et al. 2014). Due to the flexibility in shape and dimension of the mesh elements, geologic heterogeneity can be more accurately represented with fewer elements; however, FEM-based simulators may experience slow convergence or fail to converge if the discretisation creates poor quality mesh elements with large internal angles approaching 180°(e.g. Babuska and Aziz 1976). Many geological heterogeneities have a high aspect ratio (e.g. Massart et al. 2016;Salinas et al. 2017); thus, poor quality elements are created during discretisation unless a very fine resolution is used. In either case, the computational cost is high. As described later, this problem is alleviated here by using a robust Double Control Volume Finite Element Method (DCV-FEM; Salinas et al. 2017).
Current tools for simulation of ATES face a number of significant challenges. Geological heterogeneity can impact significantly on ATES efficiency and sustainability Bridger and Allen 2014;Winterleitner et al. 2018); thus, the complex architecture of aquifers targeted for ATES must be captured and preserved during numerical simulation. Moreover, in this study, it is reported that a high mesh resolution is required to accurately predict key ATES metrics such as well head temperature and thermal recovery efficiency. Together, these requirements can impose very high and infeasible computational costs, particularly when modelling realistic large-scale ATES systems. The high cost of simulations leads to unrealistic assumptions being made concerning well spacing and flowrates, to artificially reduce the cost of the simulations (Possemiers et al. 2015).
In this paper, a new approach to simulate ATES is presented to address the issues outlined in the preceding. The first application of dynamic mesh optimisation (DMO) to numerical simulation of ATES is reported. DMO, which has been used in a range of fields of computational fluid dynamics (Alauzet and Loseille 2016), allows high-fidelity simulations to be obtained at a reduced computational cost because the geometry of the mesh adapts over time to minimize a userspecified error metric for solution fields of interest such as pressure or temperature. DMO allows the mesh to refine in areas where the solution field is complex and coarsen elsewhere. Fig. 1 Working principle of ATES, showing ATES operation in summer for direct cooling of buildings and storage of warm water for the following winter, and ATES operation in winter for heating of buildings coupled with a heat pump (HP) and storage of cool water for the following summer. Modified from Pellegrini et al. (2019) Application of DMO in aquifer simulations is nontrivial, as the underlying geological heterogeneity must be preserved when the mesh changes. Conventional aquifer modelling approaches typically use a fine grid or mesh resolution to represent geological heterogeneity, in which case DMO is prohibitively expensive because the fine-scale model properties such as permeability must be up-, cross-or down-scaled onto the new mesh each time the mesh changes (Renard and de Marsily 1997;Rühaak et al. 2015). An alternative approach to represent geological heterogeneity that allows efficient application of DMO is to use surface-based geological modelling (SBGM; Sech et al. 2009;Hassanpour et al. 2013;Jackson et al. 2015;Massart et al. 2016). In SBGM, all heterogeneities of interest are represented by surfaces; properties within the domains defined by the surfaces are assumed constant or are modelled with simple trends (Jackson et al. 2014). SBGM allows efficient application of DMO because the mesh adapts within the domains defined by the surfaces.
The numerical methods reported are implemented in the Imperial College Finite Element Reservoir SimulaTor (IC-FERST). IC-FERST is a tetrahedral, mesh-adaptive, parallel reservoir simulator that implements a double control volume finite element method (DCVFEM) to solve the governing equations. The open-source ICFERST code is available under the terms of the Affero General Purpose License (AGPL v3.0). IC-FERST has been applied previously in multiphase reservoir flows (Jackson et al. 2015;Teoh et al. 2021), salt transport in coastal aquifers (A. Hamzehloo, Imperial College London, unpublished data, 2022) and geothermal reservoir modelling ). The methods have been validated previously for an analogous problem to ATES of singlephase flow and heat transport against the one-dimensional (1D) analytical solution for the parabolic advection-diffusion equation .
In this paper, the first systematic convergence study for ATES simulation is reported, demonstrating that a high mesh resolution is required to properly capture key ATES metrics and that this imposes high computational cost when using a conventional fixed mesh. SBGM and DMO are then applied to a simulation of ATES under cyclic operation in two different scenarios: (1) a homogeneous aquifer, and (2) a highly realistic and heterogeneous fluvial aquifer comprising meandering channelised sandbodies interbedded with low permeability mudstones and siltstones. The latter test case is representative of target aquifers for ATES such as the Sherwood sandstone which underlies large parts of the UK and comprises fluvial sandstone bodies interbedded with mudstones, with spatially variable sandstone:mudstone ratios. Both test cases are simulated using fixed meshes and DMO, and it was shown that significant computational speed-up was achieved through use of DMO. The first application of DMO for ATES simulations is reported and it is demonstrated that this method offers significant advantages in accuracy and cost compared to conventional methods to simulate realistic ATES scenarios.

Governing equations
The governing equations solved here for thermal storage problems are the mass conservation, momentum and heat transport equations . Symbols are listed in the Appendix. Flow through a porous medium is described by Darcy's law: where ∇ is the grad operator, u is the Darcy velocity, p is the pressure, ρ f is the fluid density, μ is the fluid viscosity, K is the permeability tensor and g is the acceleration due to gravity. The continuity equation for a compressible flow, assuming an incompressible porous medium, is defined as: where s c is a source term allowing coupling between the reservoir and wells, as described later in the article. In this study, the density of the fluid (water) is assumed to be constant, consistent with previous numerical simulations of lowtemperature ATES (Bridger and Allen 2010;Sommer et al. 2014;Possemiers et al. 2015). Given the small variations in temperature induced in the aquifer, the variations in density can be neglected. Finally, heat transport is described by where subscripts p and f denote the porous medium and fluid, respectively, T the temperature (assumed to be in thermal equilibrium across the fluid and porous medium), ξ the total volumetric heat capacity, κ the total thermal conductivity, φ the porosity, C P the heat capacity at constant pressure and s t a source term. In this approach, thermal dispersion was not considered; it is often omitted, as it is typically small compared to thermal conduction (Bear 1972;Woodbury and Smith 1985;Hidalgo et al. 2009). Previous studies of ATES have also neglected the impact of thermal dispersion (Bakr et al. 2013).
The governing equations are discretised using the DCVFEM presented in Salinas et al. (2018). Contrary to the classic Control Volume Finite Element Method (CVFEM), in the DCVFEM, pressure is discretised control volume-wise, whereas all other variables are discretised following the classic CVFEM. The DCVFEM shows better convergence than CVFEM schemes as well as greater robustness against poor quality elements, which is important when solving porous media problems where high aspect ratios and complex geometries can lead to highly distorted meshes and elements with large internal angles (Salinas et al. 2018).
The DCVFEM is based on the tetrahedral element pair P n (DG) -P n + 1 (CV), where n is the discretisation order, CV denotes control volume, and DG denotes discontinuous Galerkin. The particular element pair used in this paper is P 0 (DG) -P 1 (CV). Velocity has a discontinuous zero-order finite element representation (constant across elements), whilst temperature and pressure, which are defined control volume wise, have a continuous discretisation of order one (linear across elements). Time is discretised using the adaptive implicit Θ-method, where Θ varies between 0.5 (Crank-Nicholson) and 1 (implicit Euler), following a total variation diminishing approach. Equations (1)-(6) are solved using a Picard-iterative nonlinear solver, with a specified convergence criterion of 1% (i.e. the solution is assumed to be converged when the change in the infinite norm of the normalised difference between two nonlinear iterations is less than 1%). Solutions of Eqs. (1)-(6) in IC-FERST were validated in Salinas et al. (2021) against the 1D analytical solution for the parabolic advection-diffusion equation given by Siemieniuch and Gladwell (1978).
Wells are modelled as 1D lines. In order to make a connection between the aquifer and the wells, a dual domain approach is used: the first domain represents the aquifer, denoted by superscript j = 1, and the second one represents the wells, denoted by j = 2. Both domains share the same mesh, allowing for consistency when using DMO.
To model flow within the wells, Eqs. (1), (3) and (4) are modified  to represent flow within a pipe. Equation (1) is modified to obtain where In Eq. (8), λ is the Fanning friction factor for pipes of general roughness (Fanning 1882): λ ¼ −3:6 log 10 6:9 Re þ a=r i 3:7 10 where Re is the Reynolds number inside the pipe, and a and r i are the roughness and inner radius of the pipe respectively. In Eqs.
(3) and (4), coupling between the two domains is ensured at the shared nodes along the well through corresponding source terms. In Eq. (3), the source term represents the mass flowrate between both domains: where ρ f is the fluid density in either the aquifer or the well domain based on the upwind direction, L se is the screen length or portion of the well which is open to flow and connected to the node, and γ is defined using the Peaceman correction (Peaceman 1978). The element-wise discretisation of the Peaceman correction γ e is given by in which K ne is the permeability of the element normal to the well, S is the so-called 'skin' representing enhanced or reduced rock permeability adjacent to the well, r w is the outer radius of the well and r e is defined as r e = 0.14dn (Peaceman 1978), where dn is the length of the element normal to the well trajectory (Fig. 3c). The Peaceman correction was derived for regular hexahedral meshes, but it was demonstrated in Salinas et al. (2021) that the correction can also be applied to unstructured tetrahedral meshes so long as dn is approximately uniform at a given coupled well/reservoir node. As described later in section 'Surface-based geologic modelling', the mesh is constructed in such a way that dn satisfies this requirement.
Here, S = 0 is considered. For Eq. (4), the source term describes the heat transfer between the wells and the aquifer. If the well is open to flow, heat transfer is dominated by advection and the source term is given by: where the well is closed and heat transfer is dominated by conduction. The resulting source term is obtained using the equation of heat loss through a pipe (Willhite 1967): where k w is the effective thermal conductivity of the casing and surrounding cement, and L we the portion of the well which is closed to flow and connected to the node.

Dynamic mesh optimisation
Dynamic mesh optimisation enables the mesh to refine in parts of the domain where the solution is complex and coarsen elsewhere during simulation. The basic principles of adaptivity are summarised here; a more complete explanation can be found in Pain et al. (2001). Variables that are defined control volume-wise, are interpolated from the control volume space to the finite element space. Following this initial step, the mesh resolution and aspect ratio are varied in an iterative process to minimize an error metric. The error metric is approximated from the Hessian of the solution field as well as the precision required by the user. An approximation of the error in the numerical solution is calculated based on Cea's lemma, which states that the interpolation error between a given smooth field and its linear interpolation over a given fixed mesh is bounded by a function of the Hessian matrix H. The mesh adaptivity routine uses a functional I to measure the solution quality, which depends on the interpolation error bound (Pain et al. 2001) and where E is a matrix composed of entries E ij , v i are the vectors describing the element edge lengths on the finite element mesh, ϵ is a normalisation constant, H is the Hessian matrix of the specified solution field δ is the dimension of the problem, η is the desired precision for the solution field and ι is the degree of the polynomial of the norm used, which is equal to the discretisation order of the solution field(s) that the mesh is adapting to; here, ι = 1. During the process of adapting the mesh, the functional I is minimised through an iterative process. The mesh optimisation process is also constrained by user inputs such as minimum and maximum edge lengths, maximum aspect ratio and maximum number of nodes (Pain et al. 2001). Mesh optimisation is performed using hr-adaptivity, where h-adaptivity refers to adding/removing elements and r-adaptivity to repositioning of nodes. DMO implements a series of operations during the optimisation: node repositioning, element splitting by node insertion, coarsening by node deletion and face-edge swapping (Fig. 2). These are carried out on all the mesh elements, and the mesh is modified only if these operations lead to a smaller solution error. The final step of the optimisation process is the interpolation of the solution fields from the old mesh to the new mesh, using the supermesh technique (Farrell and Maddison 2011) and a conservative, bounded interpolation method (Adam et al. 2016). In the simulations presented here, the mesh is optimised every timestep for the pressure and temperature fields. The target precision for temperature is 0.1 K and the relative precision for pressure is 0.01.

Surface-based geologic modelling
Surface-based geologic modelling (SBGM) uses surfaces to represent key geological features such as stratigraphic surfaces, lithology boundaries or faults (Pyrcz et al. 2005;Jackson et al. 2014;Jacquemyn et al. 2019). Here, parametric Non-Uniform Rational B-Spline (NURBS) surfaces are used which are widespread in a range of computer-aided design applications, implemented for geological modelling as described by Jacquemyn et al. (2019). Once surfaces are created, Boolean operations are applied to surface intersections to construct 3D watertight domains defined by the surfaces (Jacquemyn et al. 2019). In this paper, the domains represent different lithologies and petrophysical properties are assigned per lithology. Domains that represent a particular lithology are assigned the same value. Properties in each domain are therefore constant.
The SBGM approach that is implemented here does not require models to be created on a predefined mesh, which yields a number of advantages for efficient application of DMO. Once the watertight domains have been created, an initial 3D mesh is generated, defining the accuracy with which geological heterogeneity will be modelled in the simulation. From this mesh, 'master nodes' are identified that represent the discretized surfaces (Jackson et al. 2015). These nodes cannot be modified when the mesh is optimised, thus ensuring that an accurate representation of geological heterogeneity is maintained during application of DMO. When the mesh is adapted, petrophysical properties do not need to be interpolated to the new mesh; therefore, reducing the computational effort associated with DMO. Defining properties per domain rather than per element also means that properties do not need to be upscaled or downscaled following changes in mesh resolution.
The geological models used in this study are shown in Figs. 3 and 4 and corresponding fluid (f) and porous media (p) properties given in Table 1. Model A (Fig. 3) represents a homogeneous sandstone aquifer of thickness 100 m, which extends from a depth of 50-150 m, confined above and below by mudstone layers. Model B (Fig. 4) represents a 30-m-thick fluvial aquifer which comprises meandering channelised sandbodies interbedded with low permeability overbank deposits of fine-grained sandstones, siltstones and mudstones. This aquifer architecture is representative of several target aquifers for ATES projects such as the upper part of the Triassic Sherwood sandstone in the UK (Gluyas et al. 2018). The channelised sandbodies erode into each other, leading to complex connectivity of permeable sandstone and, therefore, complex flow paths (Fig. 4c,d). Similar to model A, the under-and overburden consist of mudstones. Because realistic channel meandering is not easily captured by a parametric equation, creation of meandering channels here deviates slightly from the cross-section and extrusion approach described in Jacquemyn et al. (2019). Instead, channel thalwegs and channel sides are created directly from the meandering trajectory. The process is summarised in Fig. 5.
Meandering trajectories follow work by Howard and Knutson (1984) and Sylvester et al. (2019) that describe how local and upstream curvature control lateral migration and river meandering, based on natural examples. Here, each trajectory starts off as an arbitrary, undulating sinusoidal curve that is iteratively updated by its predicted lateral migration rate A cutaway of this model is also shown (b); the wells and associated triangular element sleeves are also displayed. b) A zoom from the plot focussing on the sleeves and the associated computational mesh, is shown (c). Annotations for dn and L se , L we are provided to clarify the meaning of these variables required in Eqs. (10), (11) and (13). As shown here, the use of sleeves constrains the length of the normal for each element dn to be approximately uniform around the well at a given location and, therefore, the correct application of the Peaceman correction for 250 iterations or until a specific sinuosity is achieved. The resulting trajectory forms the basis of the channel shape. The thalweg of each channel is created by translating the trajectory vertically downwards by the chosen channel thickness. Copies of the thalweg trajectory are translated perpendicular to its course to create the left and right channel banks. Local Fig. 4 Surface-based geologic model B, comprising meandering channelised sandbodies with different petrophysical properties (domains 4,5, and 6) interbedded with overbank deposits of fine-grained sandstones, siltstones and mudstones-domain 3; removed for clarity (b); the aquifer is bounded by mudstones; see domain 1 (a). The aquifer representation shown (b) is therefore embedded within domain 1, indicated by the white box (a), which shows a cut away of the full model. A cross-section (c-b) exposes the channelised sandbodies and their complex connectivity. This type of deposit can commonly be observed at outcrop (d): uninterpreted (upper) and interpreted (lower) photographs of cliff-face exposures in which channelised sandbodies (red) form vertical ledges and surrounding mudstones (green) form partly vegetated slopes. Here, NTG denotes 'net-to-gross' where 'net' denotes the sandstone and 'gross' denotes the entire interval. Example taken from the Cretaceous Blackhawk Formation exposed in the Wasatch Plateau, Utah, USA -modified from (Flood and Hampson 2015). a-c Numbers and colours denote domains with different petrophysical properties (Table 1) Table 1 Rock and fluid properties (Allen et al. 1997;Bricker et al. 2012;Boon et al. 2021;Parkes et al. 2021). Values of permeability K specified are isotropic

Model/ domain
Rock type curvature is used to control the amount of translation to either side, thus enabling asymmetrical channels. The channel basal surface is then created by fitting a NURBS surface to the riverbanks and thalweg. The channel top surface is formed by fitting a surface to the riverbank trajectories. Combining basal and top surfaces forms the required set of surfaces to perform Boolean operations as described in Jacquemyn et al. (2019). Given that both these surfaces share the same edge geometry, they can be joined together forming a 'watertight' connection. In the model used in this study, three different channel average widths are defined which correspond to different depositional ages. Channel width ranges from 10 to 70 m and channel thickness from 5 to 10 m. The deeper set of channels (channel sand 1) correspond to the oldest deposits, and are therefore assigned lower values of permeability of order 10 −14 m 2 , whilst the most recent, shallower channels (channel sand 3) are assigned higher values of the order 10 −12 m 2 (Table 1). Both models A and B include a well doublet. Wells are vertical and represented by 1D lines. A 3D triangular sleeve, split into three sections is created, with the intersection of the sections corresponding to the well trajectory (e.g. Fig. 3b). The sleeve ensures that the trajectory of the well is preserved when DMO is applied; moreover, it guarantees that dn, the length of the normal to the well trajectory for each element, is approximately uniform and isotropic around the well (Fig. 3c), which ensures the Peaceman correction is valid ).

Numerical simulations
Model setup In both models tested, cyclic operation of ATES was simulated where injection of warm and cool water alternates via a well doublet. Cool water was always injected and produced via the cool well, whereas warm water was always injected and produced via the warm well (Fig. 1). The initial temperature at the top of the domains was set to 283 K and the temperature increased with depth following a typical geothermal gradient of 3 K/100 m (Busby 2010). Injection occurred during 6-month intervals at a fixed temperature and volume rate q in via the injection well, whereas production occurred via the production well at constant pressure P out , defined at the top of the well. After 6 months, the boundary conditions were exchanged between wells. The injected volumes in the cold and warm wells were the same, to maintain a balanced thermal load. This was achieved by using the same injection flowrate for both cold and warm wells over a fixed 6-month period. Boundary conditions and additional parameters used in the simulations are given in Table 2. Unless otherwise stated, all other model boundaries were sealed (no mass flow) and insulated (no heat flow). The systems modelled here are comparable to current small-scale ATES deployments: given the simulated pumping/injection rate of 14 m 3 h −1 (Table 2), the peak  (Howard and Knutson 1984). b Left and right channel banks (blue) are created by translating the final trajectory perpendicular to its local tangent. The amount of translation takes into account curvature to enable channel asymmetry. c The final trajectory is translated vertically by the channel thickness, and a surface is fitted to the three curves (white surface). d The top surface is created by fitting a surface to left and right channel bank curves which, with the base surface, forms a watertight volume heating and cooling capacity achieved is 160 kW. In model A, conductive heat exchange between the well and the mudstone in the overburden was included but was neglected in model B.

Model A: mesh convergence study
To determine the required mesh resolution for a converged solution, a fixed mesh convergence study was performed for the homogeneous model (model A). For all meshes, far away from the region of interest, a coarse element size was specified: 200 m in the x,y-directions and 50 m in the z-direction. For the region of interest around the wells, the area which needs to be refined based on the thermal radius is estimated. The thermal radius provides a theoretical prediction of the area around the well, in a homogeneous reservoir, where the temperature field will be affected by the injected water, and is defined as (Doughty et al. 1982) where V is the volume of water injected over a cycle and ξ the volumetric heat capacity of the reservoir (Eq. 5). The thermal radius as defined in Eq. (16) does not account for conduction, dispersion, vertical flow or geologic heterogeneity and typically represents a minimum estimate of the actual thermal radius (Sommer et al. 2015). For the case simulated here, the thermal radius is R th = 16.4 m. A high mesh resolution was therefore specified for an area extending from the top surface to 200 m depth, and of 30 m radius around each well, which corresponds approximately to 2R th . A convergence study was also performed using DMO. Table 3 summarises the different models and their associated mesh parameters. All simulations were run on a single core for fair comparison of computational cost. Parallel implementation is discussed later. Injection and production were simulated for a period of 5.5 years. To quantify the effect of mesh refinement, several metrics were compared: the average production well head temperature during the last production period, and the thermal recovery factor for the last production/injection cycle. The thermal recovery factor T r measures the efficiency of heat/ cool recovery from each well and is defined as : where T n is the average initial undisturbed aquifer temperature over the well screen length. The average initial temperature is used to account for the geothermal gradient. These metrics for the final production/injection cycle were chosen because they capture errors that accumulate during the simulation. Moreover, the average well head temperature reflects the temperature of the water produced along the whole well screen and over the whole production period, and thus is representative of a large portion of the aquifer rather than a single location at a single time. Both metrics are important from an operational point of view. A converged solution was identified when the variation between simulations with increasing resolution varied by less than 0.1°C for average temperature and by less than 1% for the thermal recovery factor.

Model B: mesh parameters
For model B, simulations were run with fixed and adaptive meshes for two different mesh resolutions: a coarse resolution of 10 × 10 × 2 m in the x-, y-, and z-directions respectively and a fine resolution of 3 × 3 × 2 m (B.1-B.4 in Table 3), based on the results of the mesh refinement analysis described in section 'Model A: homogeneous aquifer'. Given the high aspect ratio of the channels in this model and the potential effect of vertical flow, finer mesh resolutions of 1 and 0.5 m in the vertical direction were also tested (B.5-B.6). To compare these simulations, the same metrics were used as for the mesh convergence study: the average well head temperature during the last production period and the thermal recovery factor for the last injection/production cycle. Between simulations B.4-B.6, the average well head temperatures varied by less than 0.1°C and the thermal recovery factors by less than 1%. Therefore, simulation B.4 (resolution 3 × 3 × 2 m) was taken to correspond to the converged solution. Given the complex geological heterogeneity of this model, predicting the effective thermal radius is challenging: the high permeability channels lead to migration of the thermal plumes far beyond the theoretical thermal radius (Eq. 14). Therefore, to ensure that a sufficient portion of the model was refined to the chosen fixed mesh resolution, a fine mesh was used over the entire domain representing the channels and interbedded mudstones (Fig. 4b). Because conductive heat loss in the overburden was neglected, the mesh in the over-and under-burden was kept coarse, with maximum edge lengths of 500 × 500 × 250 m in the x-, y-and z-directions. The fine fixed mesh model contained 3.5 × 10 6 elements (Table 3).

Model A: homogeneous aquifer
As described in Section 'Model A: mesh convergence study', a mesh convergence study was performed for both fixed and adaptive meshes. As the mesh was refined, the simulated temperatures at the warm well increased and those at the cold well decreased, whilst the thermal recovery factor for both wells increased, until a converged solution was obtained (Fig. 6). In simulations with fixed meshes, a converged solution was identified when the number of elements reached 7.6 × 10 5 . In the case of adaptive meshing, a converged solution was reached with 9.1 × 10 4 elements, or ×8.4 fewer elements (Table 3). The adaptive and fixed mesh simulations converged to the same solution for all four metrics, and within the precision outlined in the preceding section.
Reasonable predictions of well head temperatures were obtained at a lower fixed mesh resolution (within 0.5°C of the converged solution), but the thermal recovery factor was then significantly underestimated; for example, in the case of simulation A.3 (fixed resolution of 15 × 15 × 5 m; see Fig. 6), the average warm wellhead temperature was within 0.5°C of the converged value, but the thermal recovery factor was 10% lower than the converged value. When developing ATES systems, high mesh resolution is Having demonstrated the importance of a high mesh resolution, the cost of simulations for adaptive and fixed meshing in a homogeneous aquifer was then compared. A comparison of the fixed (A.6) and adaptive (A.14) simulations with the highest mesh resolution after 2.5 years of simulation time is shown in Fig. 7. The temperature isolines for both simulations are observed to be very similar, consistent with the convergence study results ( Fig. 6) but, in the case with DMO, use of a coarse mesh away from regions with steep temperature and/or pressure gradients allowed results to be obtained using much fewer elements (Fig. 7c).
The reduction in required elements for a converged solution with DMO leads to a significant computational speed-up compared to the equivalent fixed mesh (Fig. 8). For the Fig. 6 Average produced wellhead temperatures (a,b,e,f) and thermal recovery factor (c,d,g,h) for the warm (red) and cold (blue) wells respectively for the last injection/production cycle for homogeneous model A and the different mesh resolutions detailed in Table 3. The dashed lines indicate the converged solutions. The green and white stars indicate simulations at which convergence is reached for the adaptive and fixed meshes respectively. Number of elements refers to average number of elements in the adaptive simulations Fig. 7 Vertical section through wells showing the temperature field after 2.5 years for the converged solution obtained in the homogeneous model A using a adaptive mesh (simulation case A.14 in Table 3) and the corresponding mesh (c). b Plot shows the temperature isolines from 8-16°C for the fixed mesh simulation (A.6) in green and for the adaptive mesh simulation (A.14) in black. Both simulations have comparable mesh resolutions of 3 × 3 × 2 m, in the x-, y-, and z-directions, respectively. An animation of simulation A.14 is shown in the electronic supplementary material (ESM1) coarser meshes, which did not yield converged solutions for the fixed mesh case, the computational cost of fixed and adaptive mesh simulations was similar. However, as the resolution increased, the ratio of elements and CPU times for fixed versus adaptive cases considerably increased, showing increasing speedup using DMO. The converged solution was obtained in 7.17 h using DMO, ×4.1 faster than with the fixed mesh.
Having demonstrated the accuracy and efficiency of using DMO compared to a fixed mesh, the converged results obtained from the simulation with DMO (A.14 in Table 3) were studied further. To assess the longer-scale behaviour, the simulation was repeated over a period of 20 years. The well head temperatures and thermal recovery factors over time are shown in Fig. 9. Over the initial production/injection cycles, the temperatures at the end of the production stage at the cold well rapidly decreased, and rapidly increased at the warm well. In the warm well, with every cycle, some heat was not recovered, leading to a gradual warming of the aquifer around the well. Due to this increase in temperature, the range of temperatures that were recovered during the following winters was reduced. An analogous process occurred in the cold well, leading to progressive cooling of the aquifer around the well. This local warming and cooling yields an increase in the thermal recovery of the system, as has been observed in previous studies (Bakr et al. 2013;Sommer et al. 2013;Collignon et al. 2020). The initial thermal recovery was approximately 0.6; it then increased and stabilised to a value of approximately 0.86 after 20 years. The warm and cold wells had approximately the same thermal recovery factor, because the difference between the injection temperatures and the average aquifer temperature was the same in both cases. The system tends towards a steady state; due to its symmetry, the cold and warm water injection balances. Fig. 8 Comparison of a number of elements and element ratio and b CPU time for simulations using fixed and adaptive meshes of equivalent resolution for homogeneous model A (simulations A.1-6 and A.8-12, A.14 in Table 3). Resolutions are given separately for the x-,y-and zdirections: for example, 3/2 m refers to 3 m in the horizontal (x,y) direction and 2 m in the vertical direction. The simulations run with the resolution required for convergence (3 × 3 × 2 m) are highlighted with white stripes. For the adaptive simulations, the number of elements corresponds to the average number of elements Fig. 9 a Well-head temperatures for warm and cold wells and b corresponding thermal recovery factor for converged adaptive simulation A.14 of homogeneous model A (Table 3) Model B: meandering channelised sandbodies To demonstrate DMO in a realistic heterogeneous model, the approach used here was applied to model B representing a fluvial aquifer with complex sandbody connectivity. As shown in Fig. 10, results of the simulations with adaptive meshing closely matched those obtained using the equivalent fixed mesh for the fine resolution required to obtain converged solutions. Similar to the homogeneous model results, errors in the predicted well temperature were small when using coarser meshes, but the thermal recovery factors differed by approximately 5%. When developing an ATES system, fine mesh resolution is essential to obtain accurate predictions of system efficiency.
For both coarse and fine resolutions, the number of elements required when using adaptive meshing was much smaller than for the corresponding fixed mesh (Fig. 11a). With DMO, there is no need to predict a priori the effective thermal radius, as the mesh automatically refines where required. With DMO, the average number of elements during the simulation was 1.6 × 10 5 , ×22 fewer than for the equivalent fixed mesh, leading to a ×15 reduction in simulation time, whilst maintaining the same solution accuracy (Fig. 11b). Notably, the number of elements in the fine adaptive case was smaller even than the number of elements required for the fixed coarse mesh. It can be concluded that DMO yields significant advantage in reducing computational cost in realistic models of highly heterogeneous aquifers which can be targets for ATES.
Focusing on the results from the simulation with a fine adaptive mesh (B.4 in Table 3), the water advanced fastest through the high permeability channels (Fig. 12), whilst heat was exchanged with the lower permeability mudstones by conduction (Fig. 12c,d). During the period of cold water injection, the mesh refined around the cold water front and coarsened at the location of the warm well; the opposite process occurred during warm water injection. As expected, the effective thermal radius of the system was much larger than the theoretical thermal radius R th = 16.4 m for both wells (Eq. 14), owing to the complex geologic heterogeneity present in the aquifer. The geometry of the plumes was also markedly different from the simple cylinder assumed in predictions of thermal radius. Here, the complex plume geometry leads to interaction between the cold and hot plumes. These results are discussed further in the next section. Fig. 10 a Well-head temperatures over a period of 5.5 years; b close-up of the warm well-head temperatures for the first injection/production cycle, and c thermal recovery factors of warm (red) and cold (blue) wells for the adaptive and fixed mesh simulations B.1-4 of heterogeneous model B (Table 3)

Discussion
In this paper, it was demonstrated that a high mesh resolution is essential for modelling ATES and that using DMO can provide significant advantages. DMO significantly reduces the computational cost of simulations, whilst maintaining the same solution accuracy compared to a fixed mesh. To accurately predict thermal recovery factor, high mesh resolution is essential. Coupling DMO with SBGM allows modelling of realistic and highly heterogeneous aquifers for ATES Corresponding CPU time for the adaptive and fixed mesh simulations and fixed/adaptive CPU time ratio Fig. 12 Plan view of the temperature field for the heterogeneous model B after a 2 years of storage, at the end of the injection cycle for the warm well (red) and b 2.5 years of storage at the end of the injection cycle for the cold well (blue), or simulation B.4 (Table 3). Due to the geothermal gradient (3°C/100 m), the deeper channels are warmer than cooler, shallower channels and so show warmer colour. Mudstones are omitted for visualisation purposes. The locations of the cold and warm wells are indicated by the blue and red triangles respectively. c Plot shows the temperature fields after 2.5 years for the channel domains with a mesh overlay, and for a horizontal cross section at the midpoint through the channel stack, showing heat transfer to the lower permeability mudstones. d Plot shows a cross section of the temperature field after 2.5 years through both wells. An animation of simulation B.4 is shown in ESM2 applications at minimal computational cost. DMO is particularly advantageous in models of heterogeneous aquifers, as the associated thermal plumes are more complex and difficult to predict; fine fixed meshes must therefore be used over a larger model volume to ensure the plumes are captured. All simulation times here were reported for a single core; however, DMO has also been implemented in a parallel computational framework.
Given the computational speed up from DMO, further sensitivities and uncertainties can be studied at lower cost compared to using a fixed mesh. Here higher injection/production rates were tested in the homogeneous model A (Table 3), consistent with many commercial ATES systems that employ higher rates to deliver larger heating and cooling capacities. Higher flow rates cause an increase in the thermal radius (Eq. 15) and therefore require a larger model domain with a fine fixed mesh which is known to increase computational cost (Possemiers et al. 2015). The effect of higher flowrates was tested with a larger well spacing of 280 m, to avoid interactions between thermal plumes, recognizing that the theoretical thermal radius for the highest rate tested corresponds to 45 m. Results are shown in Fig. 13 and also include the results from section 'Model A: homogeneous aquifer' for comparison, obtained at the lower flowrate and well spacing of 140 m ( Table 2). Simulations with DMO were run for each flowrate tested; equivalent fixed mesh simulations were run for the lowest and highest flowrates only and with a refined mesh over a radius of 100 m from each well in the higher rate case.
Increasing flowrate led to an increase in thermal recovery factor, because advection of heat became more dominant compared to conduction. Given thermal dispersion is not included in the model, the higher flowrates do not lead to increased thermal losses in the aquifer. In fixed mesh simulations, increasing flowrate leads to a significant increase in computational cost because of the larger, fine-mesh domain required to capture the thermal plumes: CPU time is increased by 155 h (Fig. 13b). In contrast, increasing flowrate in DMO simulations increased CPU time by only 11 h between the lowest and highest rates tested (Fig. 13b). This demonstrates the practical advantage of DMO, allowing modelling of larger systems and higher flowrates with a limited increase in cost. DMO also removes the requirement to estimate a priori the size of the domain in which a fine fixed mesh is required to capture the thermal plumes.
As shown in section 'Model B: meandering channelised sandbodies', geological heterogeneity leads to more complex thermal plumes compared to a homogeneous aquifer. If the warm and cool wells are placed too close together, the plumes interact, leading to a decrease in the thermal recovery factor. This is what was observed in the heterogeneous model B. The cold and hot plumes migrate beyond the predicted thermal radius R th = 16.4 m, as the injected water moves preferentially through the connected, high permeability channels, leading to interaction of the plumes (Fig. 12), despite the well spacing (85 m) being much larger than the theoretical thermal radius. After 5.5 years of operation, the average thermal recovery factor for the warm and cold wells remains low and is approx-imately~63%. Geological heterogeneity was found to have a critical impact on the efficiency of ATES and therefore, its effect on flow and heat transport must be properly captured in simulation models when planning an ATES deployment, consistent with previous studies (e.g. Bridger and Allen 2014;Possemiers et al. 2015;Winterleitner et al. 2018).
To test the relationship between well spacing and thermal recovery factor, the well spacing was increased in the heterogeneous model B from 85 to 215 m (Fig. 14). Increasing the spacing did not lead to any significant increase in computational cost of the simulations with DMO despite the larger region of temperature change in the reservoir. Simulation B.7 for a larger spacing (resolution of 3 and 2 m in the The larger well spacing (B.7) yields a smaller range of produced temperatures, and a lower rate of cooling/heating during the storage period (Fig. 15a). The thermal recovery factor increased with the larger well spacing (Fig. 15b). After 5.5 years of operation, the average thermal recovery factor for both wells reached 70%, compared to 63% with the original well spacing. This highlights the importance of accurately predicting flow and heat transport in heterogeneous aquifers; numerical modelling tools which can capture complex geological heterogeneity and the resulting heat transport and fluid flow are essential to make informed decisions regarding well placement to maximize thermal recovery.
By reducing the computational cost of simulations, use of DMO can play a key role in assessing uncertainty associated with aquifer heterogeneity, and in optimizing well placement. The impact of geological heterogeneity on ATES deployments is recognized Bridger and Allen 2014;Possemiers et al. 2015;Winterleitner et al. 2018); the coupled approach of DMO and SBGM demonstrated here allows efficient investigation of uncertain aquifer heterogeneity on thermal recovery. Compared to conventional grid-based methods, SBGM models are cheap to create and manipulate; a large number of probabilistic realisations can rapidly be generated and simulations performed for each realisation. In systems with complex, uncertain geological heterogeneity, the flow patterns are complex and hard to predict. The theoretical thermal radius is no longer valid, as shown here. This was demonstrated for only one type of geological deposit, but the theoretical thermal radius will not be valid for a range of heterogeneous reservoir types; for example, aquifers containing mudstone aquitards or high permeability fractures. The use of DMO can therefore reduce the computational cost of simulations, by removing the need to mesh a large domain at high resolution, the required size of which is unknown at the start of the simulation.
Dynamically adapting meshes allow us to capture complex features with high accuracy, and at a reduced cost compared to structured, fixed meshes. In the example of Model B, surfacebased modelling and unstructured meshes are key to preserving the connectivity between channelised sandbodies; a fixed, structured mesh requires a far greater number of elements to ensure that connectivity is maintained (Jacquemyn et al. 2019). Previous studies of other subsurface applications have demonstrated the benefits of SBGM and DMO in capturing  and preserving a range of geological features such as inclined, intersecting faults, growth faults, folded strata and complex stratigraphy, or complex intersecting fractures (Jackson et al. 2015;Jacquemyn et al. 2019;Salinas et al. 2021). In the models presented here, smaller-scale variability within the channelised sandbodies is neglected, but this is not essential to the method: it is possible to include, for example, heterogeneous channel-fill deposits (Jacquemyn et al. 2019). SBGM does not include the cell-to-cell scale variability often observed in pixel-or grid-based models, but this has been shown to have a negligible impact on flow (Osman et al. 2021). It is capturing and preserving the spatially correlated heterogeneity and connectivity of geological flow zones and flow barriers that is key to predicting fluid flow and heat transport accurately, and therefore the behaviour of ATES systems.
The computational speedup obtained with DMO offers a number of practical advantages for the understanding and development of ATES systems-for example, it can be highly beneficial when optimising large-scale ATES deployments where a number of doublet wells are considered. Interference between wells and well doublets can significantly affect thermal recovery. Identifying an optimal layout for wells is essential when planning large-scale projects to ensure their sustainability. DMO allows plumes associated with each well to be tracked at high resolution and at reduced cost compared to a fixed mesh. This also allows a (much) larger number of simulations to be run, to optimise well parameters such as spacing, depth and orientation.
The reduced run times obtained using DMO can also facilitate the inclusion of additional physical or chemical processes (Yekta et al. 2021). The addition of contaminant transport or geochemical modelling, for example, could be of particular interest for some ATES deployments which target freshwater aquifers, in order to understand the impact on groundwater quality of long-term ATES operation.

Conclusions
Accurate numerical modelling tools are essential to the safe and sustainable development of ATES systems. It was demonstrated that high mesh resolution is required in numerical simulations of ATES systems to obtain fully converged solutions and accurately predict key ATES performance metrics such as well head temperature and thermal recovery efficiency. Moreover, complex geological heterogeneity was found to lead to a larger thermal radius than predicted and may mean that the well spacing required for efficient ATES operation is large, so numerical simulations require a large model domain. Together, these requirements mean that ATES simulations have high associated computational cost when using conventional fixed mesh simulation approaches. Motivated by these findings, a different approach to simulate ATES was applied, using surface-based geologic modelling, unstructured tetrahedral meshes, a double control volume finite element method, and dynamic mesh optimisation (DMO). In this approach, it was shown that DMO allows flow and heat transport to be simulated with the same high accuracy as high resolution fixed meshes, but at significantly reduced computational cost.
In the examples described, the number of elements required in simulations using DMO was up to ×22 smaller than in the equivalent fixed mesh, and the associated reduction in simulation time was up to ×15. The approach offers significant advantages compared to conventional methods in modelling realistic large-scale ATES scenarios, in assessing the impact of uncertain geologic heterogeneity on ATES operation and efficiency, and in optimising individual and multiple ATES deployments.