Predicting the risk of saltwater contamination of freshwater aquifers during aquifer thermal energy storage

Aquifer thermal energy storage (ATES) is an underground thermal energy storage technology with a large potential to decarbonise the heating and cooling of buildings. ATES installations typically store thermal energy in aquifers that are also exploited for potable water, so a major consideration during development is ensuring that system operation will not lead to groundwater pollution. In this study, the risk of contamination due to upconing of a shallow freshwater/saltwater interface during ATES operation is investigated. Fluid flow, and heat and salt (chloride ion) transport are simulated in a homogeneous aquifer during ATES operation via a well doublet. The impact of geological, hydrological and operational parameters is investigated in a sensitivity analysis. Two new dimensionless numbers are proposed to characterise salt upconing and redistribution during ATES operation and provide a close match to simulated concentrations: CR,w characterises the contamination risk at the ATES installation, and CR,d characterises the risk at locations downstream of the ATES installation with respect to background groundwater flow. ATES systems with CR,w and CR,d < 10 introduce low risk of contamination in a homogenous aquifer, with chloride concentration at, and downstream of, the ATES system, remaining below the World Health Organisation’s advised limit. ATES installations with CR,w and CR,d > 10 cause a rapid increase in aquifer chloride concentration. The results are used to estimate an exclusion distance beyond which ATES system operation will not cause contamination in a homogenous aquifer. The dimensionless parameters proposed allow rapid assessment of the potential for saltwater contamination during ATES operation.


Introduction
More than 75% of heating and cooling demand worldwide is met by combustion of fossil fuels, representing approximately 50% of the world's total primary energy consumption (Fleuchaus et al. 2018). Aquifer thermal energy storage (ATES) is a technology that could be used at a large scale to provide low-carbon heating and cooling of buildings by storing, seasonally, hot and cold water in aquifers (Bloemendal et al. 2015). Most installations are currently located in the Netherlands (approx. 90% of all ATES deployments) but there is a large potential for the deployment of ATES systems worldwide (Bloemendal et al. 2015;Fleuchaus et al. 2018;Pelligrini et al. 2019).
Aquifers utilised for ATES are often also sources of potable water (Bonte et al. 2011b). One major concern when deploying ATES systems is, therefore, ensuring that their operation does not cause chemical or biological pollution. As the number of ATES systems grows worldwide, this concern is increasing amongst water companies and environmental agencies (Anibas et al. 2016). Groundwater pollution can take many forms, including contamination by sewage, salt, fertiliser, microbes and bacteria. Heat may also be considered a pollutant (Blum et al. 2021). Bonte et al. (2011a) found that treatment of contaminated freshwater is energy-intensive and is likely to negate the energy savings from ATES operations. To deploy ATES systems safely and sustainably, it is therefore essential to avoid groundwater contamination. Previous research on the impact of ATES on groundwater quality has largely focused on general reviews, specific site investigations and laboratory experiments, and has mostly considered the impact of temperature changes on reaction kinetics, metal and trace element concentrations (Bonte et al. 2013;Lüders et al. 2020) and mineral dissolution and precipitation (Jesußek et al. 2013;Bonte et al. 2013). Overall, at the low temperatures (<25 °C) typical of ATES systems used for heating and cooling, these studies found that ATES operation had a limited impact on groundwater quality.
One of the main risks to groundwater quality is mixing (Bonte et al. 2011b). Natural gradients in pH, redox potential or salinity are common in aquifers (Possemiers et al. 2014) and mixing of groundwater sourced at different depths can have a major influence on its chemistry and quality. During ATES operation, groundwater is produced through the well screen(s) and effectively mixed as it flows upwards through the abstracting borehole(s), through the surface pipework and heat exchanger(s), and back into the aquifer via the injecting borehole(s). This process leads to mixing of deeper and shallower groundwater, homogenising any natural gradients that were initially present. Groundwater in many aquifers is fresh near the surface, but becomes increasingly saline with depth; often, the transition between freshwater and saltwater takes place over just a few metres, allowing an interface between freshwater and saltwater to be defined (e.g. Fig. 1; MacAllister et al. 2018). Upconing of saltwater to freshwater abstraction boreholes is a well-known contamination risk (Griffiths et al. 2003;Cai et al. 2014).
Previous reviews on the impact of ATES on groundwater quality have identified mixing of fresh and saltwater as a potential risk (Bonte et al. 2011b;Possemiers et al. 2014). The World Health Organisation (WHO) currently recommends that sodium (Na) and chloride (Cl) concentrations in drinking water should be below 200 and 250 mg/L, respectively. These recommendations have been widely adopted, for example, in the EU and the USA. Changes in salt concentration can also affect mineral equilibria (Possemiers et al. 2014), leading to operational issues such as corrosion and clogging (Bloemendal et al. 2015).
Saltwater upconing toward abstraction boreholes during abstraction has been characterised previously using simple mathematical models (Dagan and Bear 1968;Garabedian 2013), laboratory experiments (Goswami and Clement 2007;Werner et al. 2009) as well as numerical models (Saeed et al. 2002;Goswami and Clement 2007;Jakovovic et al. 2011;Cai et al. 2014). Most studies have focused on groundwater abstraction and have shown that upconing is more likely to occur in low-permeability aquifers in which boreholes are operated at high abstraction rates, and with a small standoff (vertical separation) between the lowermost part of the well screen and initial saline interface; upconing is suppressed if the aquifer has lower vertical than horizontal permeability. However, ATES operation is characterised by both abstraction and reinjection of water and, as is reported in this study, the response of the saline interface to cyclic abstraction and injection via a typical ATES doublet differs from upconing during abstraction.
The aim of this study is to characterise the risk of groundwater contamination caused by upconing of a freshwater/saltwater interface during ATES operation. Groundwater flow, and heat and salt (chloride) transport are modelled using a high-resolution numerical simulation Fig. 1 Chloride concentration with depth, recorded in an observation borehole in the Triassic Sherwood Sandstone Aquifer at Chat Moss near Manchester, UK. Red-dashed line indicates the WHO limit of 8.86 mol/m 3 . Modified from Griffiths et al. (2003) methodology based on dynamic mesh optimisation (DMO), previously applied to ATES modelling in Regnier et al. (2022). The method is applied to a three-dimensional (3D) geological model of a homogeneous aquifer with a saltwater interface situated at shallow depth. ATES operation is simulated via a well doublet over fixed 6-month cycles. A sensitivity analysis is conducted to investigate the risk of contamination for a realistic range of geological, hydrological and operational parameters. Two new dimensionless numbers (C R,w and C R,d ) are proposed which can be used to characterise upconing and redistribution of salt during ATES operation. The parameters can be used to identify high and low-risk regimes for ATES operation, allowing rapid assessment of the contamination risk associated with a given ATES system. They can also be used to define zones of exclusion around ATES wells to avoid contamination of freshwater abstraction wells. Note that this study considers only balanced ATES systems in which the same volumes of groundwater are injected and produced at both wells, because (1) operating balanced systems is a legal requirement in some countries like the Netherlands (Schuppler et al. 2019) and (2) maintaining thermal balance is key to sustainable ATES operation (Pellegrini et al. 2019).

Model setup
The model used in this study represents a homogeneous sandstone aquifer which is top-confined by a 50-m-thick mudstone layer (Fig. 2a). The aquifer extends from a depth of 50 m to the bottom of the model (z = 500 m). Petrophysical and fluid properties are given in Table 1.
The initial geothermal gradient is 3 °C/100 m and the initial chloride concentration assumes a saline interface is present at depth, as observed in the Triassic Sherwood Sandstone Aquifer, at the Chat Moss observation borehole near the city of Manchester, UK (Griffiths et al. 2003). The chloride ion concentration in this example varies from 0.48 mol/m 3 in the upper, freshwater portion of the aquifer, to 291 mol/m 3 in the lower, saltwater portion with a sharp freshwater/saltwater interface at 190 m depth. The freshwater concentration of chloride is in the range typically measured for potable groundwater (Kristiansen et al. 2011;Amani et al. 2015). Similar sharp freshwater/saltwater interfaces are observed in many aquifers (e.g. Voss andSouza 1987, Saeed et al. 2002;MacAllister et al. 2018). The concentration of Fig. 2 a Geological model consisting of low-permeability mudstone overburden (brown) and sandstone aquifer (blue). The domain comprises two wells, a warm well and cold well, which extend from the top of the model to a depth of 150 m. b The wells are indicated by the thick red and blue lines, respectively, showing a cross-section of the initial chloride concentration along the solid white line. The initial concentration gradient with depth is based on the measured data shown in Fig. 1 Table 1 Rock and fluid properties (Wilson and Luheshi (1987) The groundwater density is the reference value for freshwater ρ 0 . Where more than one value is given, * indicates the value used for the example 'high-risk' case and * * for the example 'low-risk' case. Porosity, permeability and permeability anisotropy were varied during the sensitivity analysis. The values in brackets represent the range that was explored 1 * , 0.5 * * (0.01-1) chloride ions is modelled here as a proxy for salinity and this study tested a range of values for the concentration of the saline portion of the aquifer (Table 2). A linear equation of state is used to model the concentration and temperature dependence of water density: where ρ 0 = 1,000 kg/m 3 is the reference density at the freshwater reference concentration C 0 = 0.48 mol/m 3 , reference temperature T 0 = 298 K, α = 6.89 × 10 −5 m 3 /mol is the concentration density coefficient and β = 2.41 × 10 −4 K −1 is the thermal expansion coefficient.
The ATES system comprises a single well doublet. The wells are separated by 150 m and extend to a depth of 150 m; the well spacing is chosen to avoid thermal interference during operation (Bloemendal and Hartog 2018). Injection/ abstraction cycles are fixed to a duration of 6 months and the flowrate during warm and cool storage is uniform to ensure a balanced system (Bloemendal and Hartog 2018;Pelligrini et al. 2019). In 'summer' periods, cold water is pumped at a constant rate q p from the cold well, whilst warm water is reinjected at a fixed rate q in into the warm well. Conversely, in 'winter' periods, warm water is pumped from the warm well, whilst cold water is reinjected into the cold well. It is assumed that water is always reinjected at a fixed temperature T in,c and T in,w for the cold and warm wells respectively, having delivered cooling or heating to the associated building(s). However, the chloride ion concentration of the injected water varies because it is assumed there is no treatment of the abstracted water. In the simulations, at each timestep, the chloride ion concentration of the reinjected water is equal to the concentration of the water that was abstracted during the previous timestep.
Constant, horizontal groundwater flow u w is included using a velocity boundary condition along the left face as shown in Fig. 2a. The well doublet is oriented with respect to groundwater flow such that the thermal plume from one well is not advected toward the other, because this could lead to thermal interference and a significant reduction in storage efficiency (Bloemendal and Hartog 2018). The domain modelled extends significantly downstream of the wells, because the groundwater flow advects the plumes of injected water during ATES operation. The wells are therefore placed at a distance 1.5 km from the right boundary of the domain to avoid boundary effects; a domain of twice the size was also tested, and simulated temperatures and concentrations were found to vary by less than 3%, confirming the chosen domain size is sufficiently large.

Sensitivity analysis
To gain a comprehensive understanding of the contamination risk associated with ATES operation, a sensitivity analysis was conducted with different combinations of input parameters describing geological and hydrological properties, and operational controls. The parameters varied and the range of values tested is reported in Tables 1 and 2. The effect of larger well spacing was also investigated for one simulation case where the thermal radius was larger than 0.5× the original 150-m well spacing. As this was not found to have a significant effect on monitored concentrations, it therefore was not further tested. The sensitivity analysis comprised two sets of simulations. First, to understand the effect of individual parameters, simulations were run, in which each parameter was varied whilst maintaining the others fixed to a base value. Second, to ensure that the parameter space was appropriately explored and to understand the combined effect of varying multiple properties, additional simulations were run, with parameter combinations chosen using the maximin Latin Hypercube sampling method. Latin Hypercube sampling methods are often used to design sensitivity analyses, allowing for an efficient selection of model inputs and exploration of the parameter space (McKay et al. 1979;Seaholm et al. 1988). The design was created using the python package Python-DOE.
Based on these parameters, two dimensionless numbers are proposed to characterise the risk of saltwater contamination, designed by inspection of the system and simulated cases. The first characterises the risk of saline contamination at the ATES system wellheads, and is defined as: where R H is the hydraulic radius, given by (Doughty et al. 1982): and V is the total volume injected during one storage cycle (all units are SI).
The second dimensionless parameter characterises the risk downstream of the ATES installation at a distance d m from the system and a time t m since the onset of system operation: (2) This parameter should be used to assess the risk of contamination at distances d m > R H . For each simulation, the dimensionless parameters were computed to allow comparison across different simulation cases. The simulations run in the sensitivity analysis are characterised by values of C R,w ranging from 0.07 to 248, spanning 5 orders of magnitude and encompassing the parameter range expected in typical ATES deployments.
To demonstrate that, C R,w and C R,d are correlated with concentration during ATES operation, simulated chloride concentrations were plotted against the corresponding value of the dimensionless parameters and the coefficient of determination R 2 was used to quantify the quality of the fit (Bryhn and Dimberg 2011).

Classifying contamination risk
Freshwater is considered here to be contaminated where the chloride ion concentration exceeds the WHO recommended limit of 250 mg/L or 8.86 mol/m 3 . As previously mentioned, this limit is followed by the EU and US.
To assess the risk, simulated chloride concentrations at the wellhead and downstream of the ATES system were recorded. The system is classified as low risk if concentrations remain below the WHO limit. The system is classified as medium risk if the concentration exceeds the WHO limit only at the ATES system wellheads. The system is classified as high risk if the downstream concentration exceeds the WHO limit, because contamination in this case could affect neighbouring installations such as freshwater abstraction wells.

Modelling approach
The open-source IC-FERST numerical simulator was used to model combined groundwater flow, heat transport and contaminant transport. IC-FERST is a 3D tetrahedral, mesh-adaptive, parallel reservoir simulator (Jackson et al. 2015). ICFERST has been previously validated and used to model single phase fluid flow and heat transport in the context of geothermal energy production ) and ATES operation (Regnier et al. 2022). It has also been validated for component transport and applied to the problem of saline intrusion in coastal aquifers (Hamzehloo et al. 2022). An additional validation of the method is given in the Appendix, against the analytical equation for upconing of a saltwater/freshwater interface (Dagan and Bear 1968) which is of particular relevance to this study.
The governing equations are discretised using a double control volume finite element method and solved using a fixed-point Picard-iterative nonlinear approach (Salinas et al. 2018. The convergence criterion for the nonlinear solver was set here to 0.03 (i.e. the solution was assumed to be converged when the change in the infinite norm of the normalised difference between two nonlinear iterations is less than 3%). The adaptive implicit theta method is used to discretise time, with theta varying between 0.5 (Crank-Nicholson) and 1 (implicit Euler), based on a total variation diminishing approach . A high order method is used to calculate interface fluxes, using a normalised variable diagram ).
An unstructured tetrahedral adaptive mesh is used to discretise space, which adapts to solution fields during a simulation using dynamic mesh optimisation (DMO). The process of DMO has been described previously . In this study, the mesh adapted to the pressure, temperature and concentration fields. The maximum element edge lengths allowed in the domain were 500 m in the x and y directions, and 100 m in the z direction. The minimum edge lengths were 3 m in the x and y directions, and 0.5 m in the z direction. An adaptive time-stepping method was used, which is based on the convergence rate of the nonlinear solver: the faster the nonlinear solver converges, the larger the next timestep used. No minimum timestep was specified, but the maximum timestep was limited to 7.6 days.

Results
Two example cases from the ensemble of models run in the sensitivity study are first reported to illustrate end-member examples of ATES systems with a 'high' and 'low' risk of causing salt contamination of the freshwater aquifer.

'High risk' example case
First, the near-well effects on the saltwater/freshwater interface during ATES operation are shown. The chloride concentration fields at different snapshots in time are shown in Fig. 3a-c. There are several key processes that can be observed in the near-well region. As expected, water abstraction leads to upconing of the saltwater interface below the abstracting well, leading to intrusion of water with a chloride concentration higher than the WHO limit. However, unlike simple groundwater abstraction, downconing of the saltwater interface below the injecting well is also observed. Production of higher salinity water into the abstracting well leads to reinjection of water with a higher chloride concentration at the injecting well along the entire screen length. Alternating abstraction and injection at each well creates a plume of higher salinity water around both wells, which is transported laterally away from the wells and diffuses upwards into the overlying low-permeability mudstone layer. Background groundwater flow transports the plume of higher salinity water downstream (Fig. 4a) and, in this case, after 25 years of operation, the chloride concentration of the groundwater exceeds the WHO limit to a distance of 230 m downstream from the ATES wells.  Figure 5a shows the concentrations simulated at the wellheads and at offset monitoring points downstream from the wells. Upcoming of saline water soon after the system starts operating causes a rapid increase in the wellhead concentrations, exceeding the WHO limit. The concentrations monitored downstream show a more gradual increase, with the closest monitoring points increasing most rapidly. After 25 years of operation, the monitored concentration at a distance of 120 m downstream of the ATES wells exceeds the WHO limit; however, at 360 and 480 m, the concentrations remain far below the limit. Overall, in this case, the risk of freshwater contamination is high. The chloride concentration is above the WHO limit at the wellhead and downstream of the ATES system.

'Low-risk' example case
The chloride concentration fields at different snapshots in time for an example 'low-risk' case are shown in Fig. 3d-f. Compared to the example 'high-risk' case, upconing and downconing of the saltwater interface is significantly reduced. The concentration isolines from 8.86 to 100 mol/m 3 remain broadly flat, which is a consequence of the combined effect of the lower abstraction and injection flowrates in this case, and anisotropy of the aquifer permeability with low vertical permeability suppressing upconing. The limited movement of the saltwater interface means that the concentration of the groundwater which is circulated through the ATES system never exceeds the WHO limit (Figs. 3d-f and 5b); moreover, the downstream monitored concentrations show a gradual increasing trend but always remain well below the limit. This example ATES system therefore exhibits a low risk of salt contamination of the groundwater during operation.

Dimensionless analysis
The processes previously described in the two example cases are generally representative for long-term ATES operation in homogeneous aquifers; however, as illustrated by the differences between the example high-and low-risk cases, the magnitude of upconing and downconing of the interface, and the resulting risk of contamination, will inevitably change for different geological and operational parameters. A sensitivity analysis is therefore essential to understand the effect of these parameters on the risk of salt contamination in homogeneous aquifers.
Following the sensitivity analysis described in section 'Methodology', the dimensionless parameter C R,w was plotted against monitored wellhead concentrations at early and late times (Fig. 6). There is a clear trend, with increasing values of C R,w corresponding to systems with increased salt concentration at the wellhead (Fig. 6a). Statistically meaningful second order polynomial trends can be identified, associated with an R 2 value of 0.87 for both early and late times. For systems characterised by C R,w < 10, the chloride concentration remains approximately constant and below the WHO limit, irrespective of the value of C R,w . These systems are therefore at low risk of saltwater contamination during ATES operation, whereas for C R,w > 10, there is a gradual increase in the concentration at the wellhead with increasing C R,w . Above this threshold, the risk of saltwater contamination at the ATES wellheads markedly increases, with simulated concentrations above the WHO limit.
Chloride concentration in the aquifer downstream of the wells was also plotted as a function of C R,d (Fig. 6b). A second-order polynomial trend is observed between C R,d and concentration, characterised by an R 2 value of 0.67. For C R,d < 10, simulated concentrations remain below the WHO limit. For larger values of C R,d > 10, increasing C R,d corresponds to increased concentration, beyond the WHO limit. Note that simulated concentrations at monitoring distances of 360 and 480 m always remained below the WHO limit, a result that is discussed further in a later section.

Sensitivity study
The sensitivity study allows the impact of varying individual reservoir and operational parameters on salt concentration at the wellheads and downstream of the ATES installation to be quantified (Fig. 7).
For both wellhead and downstream concentrations, varying the standoff d has the most significant impact on concentration, with decreasing standoff leading to increased concentration. The injection/abstraction flowrate Q has the next largest impact, with increased flowrate leading to increased concentration. Upconing of a saline interface during abstraction is known to be sensitive to these parameters (e.g. Dagan and Bear 1968;Garabedian 2013;Werner et al. 2009;Jakovovic et al. 2011). The other parameters varied in the sensitivity analysis have a comparably smaller impact on salt concentration, and the rank order of significance is different for concentration at the ATES wellheads (Fig. 7a) and downstream of the ATES installation (Fig. 7b). Salt concentration is lower for smaller values of k z /k x and longer screen length L s , consistent with previous upconing studies; salt concentration is also lower for higher background groundwater flow |u w |. In this base case scenario, permeability k x has almost no impact on salt concentration downstream of the installation and only a small impact on salt concentration at the wellheads, with higher permeability corresponding to lower concentration. Perhaps surprisingly, the salinity of the intruding saltwater C i was found to have negligible impact on salt concentration at the wellheads or downstream of the installation, a consequence of the dependence of density on b Ratio of simulated chloride concentrations to the WHO limit after 15 years, recorded downstream from the ATES system, at 100 m depth, for different values of C R,d with d m = 120, 360 and 480 m evaluated at t = 10 and 15 years. The black-dashed line is defined by a = 1.097 × 10 −5 , b = 1.29 × 10 −2 and c = 0.246. The limit (at which C/C L = 1) is indicated by the red-dashed line salinity (Eq. 1). More saline water increases contamination, but also has higher density, suppressing upconing of the saltwater/freshwater interface. The concentration of the intruding saltwater was, therefore, omitted from the dimensionless parameters C R,w and C R,d .

Risk of saltwater contamination at the ATES system wellheads
Despite the high risk of saline contamination for some systems, there is a large parameter space within which ATES systems can be designed and operated with low risk of upconing causing contamination at the wellheads (Fig. 6a): the risk of upconing increases as the value of the dimensionless parameter C R,w increases above the threshold value of 10. Inspection of Eq. (2), which defines C R,w , shows that contamination risk increases with decreasing standoff d (varying as d -6 ), increasing hydraulic radius R H (a function of screen length L s , varying as √ 1∕L s ) and increasing stored volume V (where R H and V are both functions of the injection/abstraction flowrate Q, varying as Q 2.5 ), consistent with the results of the sensitivity analysis (Fig. 7a). Contamination risk also increases with increasing permeability ratio k z /k x (varying as k z /k x ), decreasing background aquifer flow |u w | (varying as |u w | -5 ), and decreasing permeability k z (varying as 1/k z ), but the results of the sensitivity analysis show that the impact of varying these parameters over the range reasonable for ATES installations is comparably smaller than varying d and Q. These are the most important operational parameters controlling the risk of upconing and saline contamination of the ATES system.
A sustainable ATES deployment must operate at low risk of contamination over long timescales. For a given value of C R,w , the concentration at the wellhead is typically similar or lower after 15 years of operation, compared to the value after 0.5 years of operation (Fig. 6a). In particular, systems at high risk of contamination, with large values of C R,w > 10, show a large decrease in concentration at the wellhead after 15 years of operation, which is significant, as it suggests that contamination risk does not increase with time. The cyclical nature of injection and abstraction at each well during ATES operation reduces the long-term risk of contamination because upconing in each well during abstraction is at least partly balanced by downconing during injection across each seasonal cycle. As this differs from upconing during groundwater abstraction, in which the cone always moves towards the well, the risk of contamination only increases with time.

Risk of saltwater contamination downstream of the ATES installation
There is also a large parameter space within which ATES systems can be designed and operated with low risk of upconing causing contamination downstream of the system (Fig. 6b). Contamination risk increases as the value of the dimensionless parameter C R,d increases above the threshold value of 10. Inspection of Eq. (4), which defines C R,d , shows that downstream contamination risk increases with decreasing standoff d, increasing hydraulic radius R H (again, a function of screen length L s ) and increasing stored volume V (where R H and V are again both functions of the injection/abstraction flowrate Q), consistent with the results of the sensitivity analysis (Fig. 7b). Contamination risk also increases with increasing permeability ratio k z /k x and decreasing permeability k z , but the results of the sensitivity analysis show that the impact of varying these parameters over the range reasonable for ATES Fig. 7 Impact associated with each variable on simulated chloride concentrations: a at the wellheads after 0.5 years of operation; and b 120 m downstream of the ATES system after 15 years of operation. Each bar extends over the range corresponding to the minimum and maximum value of the sampled variable in the sensitivity analysis, where each parameter was varied individually 1 3 installations is comparably smaller than varying d and Q. These are the most important operational parameters controlling the risk of saline contamination both at, and downstream of, the ATES system.
Comparing C R,w (Eq. 2) to C R,d (Eq. 4), it is evident that some controls on contamination at the ATES wellheads are broadly the same as those on downstream contamination; however, C R,d has an increased dependence on the hydraulic radius R H (varying as R H 2 as opposed to R H for C R,w ). Larger values of hydraulic radius R H are associated with a larger spread of the saline plume around the ATES system and, therefore, increasing risk of downstream contamination. C R,d also includes a dependence on monitoring time (varying as t m 0.5 ) and monitoring distance (varying as d m -1.5 ). As expected, downstream contamination risk increases with time and decreases with increasing distance; however, simulation results show that concentration remains below the WHO limit at sufficient distance from the ATES system, even at late times. Concentration profiles along section C-C′ in Fig. 2, at a depth of 100 m, after 15 years of operation, are shown in Fig. 8a. Consistent with the results shown in Fig. 6b, simulations with the largest values of C R,d are associated with the highest risk of downstream contamination, with concentrations higher than the WHO limit; however, in all cases, concentrations remain below the WHO limit at distances greater than 500 m from the ATES system. To determine whether this holds for longer-term ATES operation, eight cases with values of C R,d > 10 were simulated for 25 years (Fig. 8b). At this later time, the saline plumes have been further advected downstream, but concentrations remain below the WHO limit beyond 500 m from the ATES system. The reason is that the plumes spread laterally as well as downstream, such that the concentration at late times remains below the WHO limit, even downstream of ATES systems with high concentration at the wellheads. This distance can be compared with the hydraulic radius which, for ATES systems with C R,d (120 m, 15 years) > 10, has an average value of ca. 74 m and, for the highest risk case, a value of ca. 96 m.
Finally, a key difference between C R,w and C R,d is the effect of groundwater velocity u w , which appears in both the numerator and denominator of C R,d , but appears only in the denominator of C R,w . Increasing groundwater velocity only decreases contamination risk at the wellhead, but also Fig. 8 Ratio of simulated chloride concentrations to the WHO limit C L (250 mg/L = 8.86 mol/ m 3 ) at a depth of 100 m along line C-C′ in Fig. 2 after a 15 years for all simulated cases and b 25 years for selected cases with high downstream contamination risk. a Curve colour denotes the corresponding value of C R,d (120 m, 15 years); blue curves correspond to low-risk cases with C R,d (120 m, 15 years) < 10; red curves correspond to high-risk cases. b Solid lines show results after 15 years; dashed lines show the corresponding results after 25 years for eight cases values of C R,d > 10. The corresponding exclusion distances for each case, calculated using Eq. (5) based on a 25-year operation, are indicated by dots of corresponding colours. The WHO limit (at which C/C L = 1) is indicated by the black-dashed line contributes to increased risk of contamination downstream due to advection of the plume.; however, cases with high risk of downstream contamination are not correlated with the highest groundwater velocities and vice-versa ( Fig. 9). At the ATES wellheads, increasing groundwater velocity always correlates to a decrease in concentration, due to reduced upconing. Conversely, downstream of the ATES system, concentration rapidly increases and then slowly decreases with increasing groundwater velocity over the broad range of values tested ( Table 2). None of the ATES cases tested here with the largest values of groundwater velocity are characterised by a value C R,d (120 m, 15 yrs) > 10. Cases with low groundwater velocity show the highest salt concentrations around the ATES wells; however, the contamination plume does not spread far downstream, so there is a low risk of downstream contamination. In cases with higher groundwater velocities, the saline plume is advected further downstream, but concentrations remain low due to the reduced upconing.

Application of the dimensionless parameters
The dimensionless numbers proposed in this study allow rapid characterisation of upconing of saltwater during ATES systems operation in a homogeneous aquifer, based on key aquifer and operational parameters. A threshold for these numbers can be identified below which ATES systems show low risk of contamination. The dimensionless numbers can therefore be used as a tool when developing ATES systems to make a rapid assessment of the risk of saltwater contamination assuming a homogeneous aquifer, and during design of an ATES system to understand how specific parameters can be optimised to reduce the risk of contamination. For systems with values of C R,w and C R,d > 10, operational parameters such as the depth of the well completion, the screen length and the abstraction flow rate can be adjusted to reduce risk. If possible, the location of the ATES system can also be adjusted to minimise the risk caused by an unfavourable geological/hydrological setting.
In this study, the case of saltwater contamination was explored, but a similar methodology could also be applied to other forms of contamination. Gradients in aquifers are a common occurrence; the dimensionless number proposed here could be used to assess the risk of mixing of redox conditions, heavy metals or other pollutants forming an interface below an ATES system. It should also be noted that, in this study, the risk of contamination arises only through mixing between freshwater and saltwater. However, saline intrusion can cause other processes that further modify groundwater quality such as changes in microbial activity, mineral dissolution, metal mobilisation and ion exchange (Barker et al. 1998;Weston et al. 2006;Eagling et al. 2013). These additional interactions between the solute and the aquifer can affect spreading of the contaminant plume and, therefore, the overall contamination risk. For a more complete assessment of the risk of groundwater contamination associated with saline intrusion, these processes could be integrated in the simulation methodology; for example, by coupling the simulator with geochemical package such as PHREEQC (Yekta et al. 2021). This would then allow additional parameters to be taken into account in the dimensionless numbers introduced in this study.

Exclusion distance for ATES installations
For ATES installations with C R,d > 10, downstream contamination occurs, but concentrations above the WHO limit are restricted to within 500 m of the system. The area of influence for a particular installation, however, will vary depending on the parameters which describe the system. Using Eq. (4) and the observed threshold C R,d > 10 for which contamination is expected, this area of influence can be estimated. An exclusion distance d ex is defined for any given ATES system: To ensure these are installed at a necessary distance from existing groundwater usage-(for example, the drawdown region around a freshwater abstraction borehole) this distance could be used to define an exclusion distance around planned ATES systems at risk of upconing. The exclusion distances for the cases shown in Fig. 8b are shown, corresponding to 25 years of system operation. In all cases, concentration levels beyond the corresponding exclusion distance remains below the WHO limit after 25 years. It should be noted that abstraction boreholes operating above a saline interface are always at risk of contamination caused by local upconing, irrespective of any additional risk introduced by ATES system operation.

Impact of geological heterogeneity
The results presented here assume a homogeneous aquifer. Heterogeneity is known to have a significant impact on flow in the aquifer and therefore will affect salt transport Bridger and Allen 2014;Possemiers et al. 2015;Winterleitner et al. 2018). High-permeability layers will transport saline water further away from the wells, reducing the effective screen length and increasing the effective hydraulic radius (Bridger and Allen 2014), whilst lower permeability features will tend to limit the upwards movement of saltwater. Geological heterogeneity is complex and diverse; it is therefore not possible to generalize the risk of saltwater contamination for ATES systems irrespective of geological heterogeneity. However, one approach to increase the scope of applicability of the dimensionless numbers introduced here is to define an effective hydraulic radius and effective permeability ratio that account for the heterogeneity. If these effective properties can be estimated empirically or numerically, they could then be used to calculate C R,w and C R,d to provide a preliminary prediction of contamination risk.

Contamination risk for multi-well systems
In this study, the cases considered included only a single well doublet. ATES developments regularly consist of multiple well doublets or are located in the vicinity of other ATES systems (Bakr et al. 2013;Duijff et al. 2021). There are likely to be additional complexities arising when considering largescale ATES developments; however, the dimensionless number defined in this study may still be used to provide a first assessment of the risk of saltwater contamination. If each well pair is designed and operated to have values of C R,w and C R,d < 10, the overall risk of the system will likely remain low.

Conclusions
ATES systems are installed and operated in aquifers which also supply potable water and a key contamination risk is mixing of saltwater that underlies the freshwater aquifer during system operation. Here, ATES operation via a well doublet in a homogeneous sandstone aquifer is modelled to quantify the contamination risk caused by upconing of a saltwater/freshwater interface. Nondimensional numbers C R,w and C R,d are defined as a function of key reservoir and operational parameters and are strongly correlated with wellhead and downstream salt concentrations, respectively. For values of C R,w and C R,d < 10, the risk of contamination by upconing saltwater is low. Above this threshold value, there is a marked increase in the risk of saltwater contamination, although downstream contamination is restricted to 500 m of the ATES system in all cases, even at late times, owing to the lateral as well as downstream spread of saltwater. High groundwater flowrates are observed to be associated with a lower risk of downstream contamination, because upconing beneath the ATES system is suppressed. The standoff between the ATES wells and the saline interface, and the ATES abstraction/injection flowrates, are the most important operational parameters controlling the risk of saline contamination both at, and downstream of, the ATES system. Contamination risk is reduced with larger standoff and lower flowrates. The dimensionless numbers defined here can be used to assess the potential risk of saltwater contamination resulting from ATES system operation, and also yield an estimated exclusion distance around ATES installations in aquifers at risk of saltwater upconing to protect other groundwater uses (such as freshwater abstraction). The findings presented here were obtained assuming homogenous aquifer behaviour; future work will investigate the impact of geological heterogeneity.

Appendix 1: Model validation-Comparison against benchmark solution
The numerical methods implemented in the IC-FERST code used here have been previously validated for single-phase flow and heat transport against the one-dimensional analytical solution for the parabolic advection-diffusion equation given by Siemieniuch and Gladwell (1978). The methods have also been verified previously against the Henry problem (Henry 1964), a standard benchmark problem for density-dependent flow and transport models. Results for these validations are reported in Salinas et al. (2021) and Hamzehloo et al. (2022).
In this study, an additional validation of the method is provided which includes well boundary conditions. The model validation reported here was previously presented as a conference paper at the 2022 European Geothermal Congress. The numerical model presented here is validated against the analytical expression describing interface upconing under a pumping well presented in Dagan and Bear (1968).
The rise of the interface Z as a function of time t and distance x from the pumping well is given by: with subscripts 1,2 denoting the fresh and saltwater phases respectively, μ the fluid viscosity, Q the well discharge rate, ∆ρ the fluid density difference, g the acceleration due to gravity, k x ,k z the horizontal and vertical permeabilities of the aquifer, ϕ the porosity of the aquifer and d the distance between the bottom of the well and the saltwater interface at t = 0.
The analytical expression assumes a two-dimensional (2D) model, a homogeneous aquifer and an abrupt interface between two incompressible fluids. A 2D computational domain of size 1,000 m and 1,000 m was defined with a single abstraction well extending from z = 0 m to z = 150 m, 50 m above the initial interface. The freshwater and saltwater phases are modelled as a single phase, with different initial concentrations of salt. In order to fulfill the assumption of an abrupt interface between the two phases, the solute diffusivity and dispersivity were set to zero. The mesh around the interface was also refined to a resolution of 5 mm in the z direction to ensure a sharp initial contact. The parameters used for the validation are given in Table 3. Table 3 Inputs and simulation parameters for validation of the numerical modelling method against the analytical solution given by Dagan and Bear (1968) Parameter (unit) Value Q (m 2 /s) 9.8 × 10 −5 -3.9 × 10 −4 d (m) 50 k x (m 2 ) 10 −14 k y (m 2 ) 10 −14 ϕ (-) 0.2 ρ 1 (kg/m 3 ) 1,000 ρ 2 (kg/m 3 ) 1,030.5 μ 1 (Pa.s −1 ) 0.001 μ 2 (Pa.s −1 ) 0.001 g (m/s 2 ) 9.81 ∆ t (s) 5,000 Min edge length, x-direction (m) 0.005 Min edge length, y-direction (m) 1 Previous experimental studies of saline upconing have been compared to the analytical solution derived in Dagan and Bear (1968) by comparing the rise of the apex of the interface, which is located directly below the pumping well (Werner et al. 2009). Here, the rise of the apex was compared between the analytical solution and numerical simulations run in IC-FERST. Simulations were run for a range of flowrates and a good agreement was obtained Fig. 10c,d.

Fig. 10
Numerical method validation for density-dependent flow to a vertical abstraction well, against the analytical expression for saline upconing described by Dagan and Bear (1968). a Example of saline upconing experiment run by Werner et al. (2009), with Z(x,t) denoting the analytical expression from Dagan and Bear (1968) describing the interface rise as a function of distance from the well (x) and time (t) and b comparison of interface apex rise for a range of experiments against analytical expression (solid lines). c Model setup in IC-FERST for simulation of interface upconing towards an abstraction well, with d the initial distance between the well and the interface at t = 0 s and d comparison of interface apex rise for numerical simulations (circles) and analytical expression (solid lines) for a range of pumping flowrates