ATES systems in aquifers with high ambient groundwater flow velocity

Aquifer thermal energy storage (ATES) is a technology with worldwide potential to provide sustainable space heating and cooling using groundwater stored at different temperatures. In areas with high ambient groundwater flow velocity (g25 m/y) thermal energy losses by displacement of groundwater may be prevented by application of multiple doublets. In such configurations two or more warm and two or more cold wells are aligned in the direction of the ambient groundwater flow. By controlling the infiltration and extraction rates of the upstream and downstream wells, the advection by ambient groundwater flow can be compensated by storing thermal energy through the upstream well, while re-extracting it from the downstream well. This study uses analytical and numerical tools and a case study to analyze the relevant processes, and provides guidelines for well placement and an operation strategy for ATES wells in aquifers with considerable groundwater flow. The size of the thermal radius relative to ambient groundwater flow velocity is an important metric. With multiple wells to counteract groundwater flow, this ratio affects the pumping scheme of these wells. The optimal distance between them is around 0.4 times the distance traveled by the groundwater in one year. A larger distance negatively affects the efficiency during the first years of operation.


Introduction
Many governments and companies set targets to reduce greenhouse gas (GHG) emissions (EU, 2010;UN, 2015;SER, 2013;Ministry-of-Economic-affairs, 2016). To meet these goals, the heating and cooling demand in the built environment is important because it consumes about 40% of the total fossil energy worldwide (Jong, 2016;EIA, 2009;RHC, 2013). Aquifer Thermal Energy Storage (ATES) systems help reduce energy use by providing seasonal storage and recovery of heat, which allows sustainable space heating and cooling for buildings. Buildings in moderate climates generally require cooling in summer and heating in winter. Where aquifers of sufficient capacity exist, the seasonal heat and cooling discrepancy can be overcome by seasonal thermal energy storage and recovery in the subsurface (Bloemendal et al., 2015). Thus, in winter a building is heated by means of a heat pump that extracts heat from warm groundwater that was stored in the previous summer. While delivering its heat to the building, the heat pump simultaneously cools this groundwater, which is re-injected into the subsurface with a second well, the "cold" well. During the summer, the flow is reversed, and then cold water is extracted and used to directly cool the building (generally bypassing the heat pump). While cooling the building, the groundwater is warmed up in the heat exchanger and immediately injected into the other well, the "warm" well.
This system thus balances out seasonal discrepancies in the supply and demand of heating and cooling. The warm and cold ATES wells can be separated horizontally; each pair thus formed is then called a doublet. The well screens can also be installed vertically in a single borehole, forming a pair called a monowell (see Fig. 1).

Problem statement
Several studies have shown how the wells of ATES systems can best be designed (Bloemendal and Hartog, 2018;Doughty et al., 1982;NVOE, 2006). In aquifers with a high ambient groundwater flow velocity, losses of thermal energy caused by groundwater advection can be limited by choosing a shorter screen length, as was shown by Bloemendal and Hartog (2018). However, in many cases this strategy is neither possible nor desirable, because short screens limit the capacity of the wells and increase their thermal radius, which precludes optimal use of available aquifer space. In practice, under high ambient groundwater flow conditions the so-called recirculation system (see Fig. 1) is often applied. Compared to the normal ATES systems, these systems have a smaller temperature difference between the warm and cold well, a lower efficiency and a large downstream thermal plume, which may affect other ATES systems or groundwater uses. Bloemendal and Hartog (Bloemendal and Hartog, 2018) showed that at groundwater flow velocities of > 25 m/y, heat losses due to the groundwater displacement become considerable, relative to the conduction losses. In areas with a high ambient groundwater flow velocity, thermal losses due to groundwater displacement can be prevented by installing multiple doublets, where at least two wells of the same type (warm-cold) are aligned in the direction of the ambient groundwater flow. By injecting the yearly storage volume in the upstream well and extracting it from the downstream well in the next season, the ambient groundwater flow is counteracted, resulting in higher recovery efficiency. This principle is schematically represented in Fig. 2. The optimal design and operation strategy of such an ATES system depends on several conditions; -The actual ambient groundwater flow velocity and direction.
Field estimates of groundwater flow velocity and its direction always come with an uncertainty, because limited groundwater head measurement locations are available. The same is true for aquifer samples to identify the spatial variation of horizontal anisotropy and hydraulic conductivity. Reducing the uncertainty of the flow velocity, its direction and its range is expensive. NB. In this study, variations in groundwater flow velocity and direction are not considered. In the Dutch aquifers used for ATES such variations rarely occur. -The possibility to place wells at a required distance and in line with the groundwater flow.
In urban areas, well placement is often limited by buildings and by infrastructure in the shallow subsurface, which may lead to wells at suboptimal locations (e.g. not precisely in line with the groundwater flow, and/or at larger or smaller mutual distance than optimal). -The possibility to control the wells for counteracting the groundwater flow.
The basic scheme for counteracting the groundwater flow depicted in Fig. 2 is to install two warm and two cold wells, in order to infiltrate in the upstream well during one season, and extract from the downstream well in the next. However, for ATES systems with multiple warm and cold wells, the full pumping capacity of both warm and cold wells is required during the warmest and coldest days of the year to provide peak heating and cooling capacity to the building. Under such conditions, when heating or cooling demand exceeds the capacity of one well, both warm and cold wells are needed to supply energy (so that all wells are involved, for either pumping or injection)regardless of any desired groundwater flow counteraction.

Goal and approach
The goal of this paper is to identify under which conditions individual ATES systems can counteract ambient groundwater flow by using multiple wells of the same type aligned with the flow direction. The approach to meet this goal is to first identify the dominating variables and control possibilities, which is done in Section 2 by describing the analytical relations, working conditions and numerical simulation tools. During extraction and infiltration at high groundwater flow velocity, the flow lines from and to the wells cannot be treated as radial, as would be valid for ATES systems in aquifers with low flow velocity (Bloemendal and Hartog, 2018 Recirculation. Recirculation systems always use the same wells for extraction and infiltration; water is extracted from the upstream well and injected into the downstream well (sometimes also referred to as "pump and dump" systems), the arrow represents the ambient groundwater flow direction.
The processes and the dependency of the energy efficiency on dominant parameters are described, as well as the controls that help optimize design and operational strategy. Secondly, results of the basic and detailed conditions specific for ATES systems in practice are analyzed in Section 3 by numerical simulations, to assess the required control under more realistic conditions encountered in practice. This is done with a general simulation following the working conditions identified in Section 2, as well as for a case study.

Working conditions
The range of working conditions that are analyzed in this paper are derived from the characteristics of existing ATES systems in the Netherlands presented by Bloemendal and Hartog (Bloemendal and Hartog, 2018). The average groundwater flow velocities of interest in the Dutch aquifers range from 25 to 200 m/y. Distances between ATES wells of the opposite type range from several meters to building plot sizes (regularly around 50-150 m), so well distances of the same type can be in the same range or even smaller. ATES system storage volumes (V) range from 50,000 m 3 /y to 500,000 m 3 /y for doublet systems, with well screen lengths (L) between 5 and 150 m.

Loss of heat due to displacement by ambient groundwater flow
Significant ambient groundwater flow is known to occur at ATES sites (Bonte et al., 2013a;Groot, 2013;Hartog et al., 2013), which leads to displacement of the injected volumes (Bonte et al., 2013b;Bear and Jacobs, 1965). This may considerably reduce the thermal energy recovery efficiency of ATES systems as ambient groundwater flow (u) contributes to thermal losses by displacing the injected water during storage. The heat transport velocity (u * ) is retarded with respect to ambient groundwater flow (Doughty et al., 1982;Hecht-Mendez et al., 2010) due to heat storage in the aquifer solids. The thermal retardation (R) depends on porosity (n) and the ratio between volumetric heat capacities of water (c w ) and aquifer (c aq , with c aq = nc w + (1-n)c s and c s the solids volumetric heat capacity), following: Resulting in a heat transport velocity at approximately 50% of the groundwater flow velocity (u). Under conditions of ambient groundwater flow, thermal energy stored in an aquifer will thus be displaced and can only be partly (Bear and Jacobs, 1965) recovered.

Retrieving heat from the downstream well
The heat transport follows the same rules as water transport, only at a lower rate due to thermal retardation. Disregarding conduction, dispersion and diffusion then allows the use of water particle tracking to assess the distance at which wells should be placed to capture all heat in the downstream well that was infiltrated by the upstream well.
The approach to counteract the effect of ambient groundwater flow (u) on heat displacement in the subsurface with multiple wells was also used for ATES systems by Groot, (Groot, 2013). This work focused on a system with two wells aligned along the groundwater flow, in which the upstream well injects a volume of water that is recovered by the downstream well after a possible intermediate storage period in which no pumping takes place. It follows from Fig. 2 that the distance D between the upstream and downstream well should equal the distance by which the heat is displaced by the groundwater, during the time between infiltration of the first warm/cold water, and extraction of the last (this period is a year under theoretical operational conditions, where storage volume is completely extracted and the pumping pattern has a symmetrical shape for infiltration and extraction). At this distance, any particle that leaves the injection well under an angle with the ambient flow will be captured at the downstream well at an intermediate time. Hence, for the analysis of the ideal placement of the two wells, it suffices to consider only the flow along the straight line through the two wells, which are aligned with the ambient groundwater flow. The specific discharge q x [m/d] by a well with flow Q [m 3 / d] (injection positive), in an aquifer of thickness L [m], with a uniform ambient flow with ambient groundwater flow q = u n [m/d] along a line through the well parallel to the ambient flow, can be expressed as: During the storage period Q = 0, the water moves entirely with the ambient flow, while during infiltration and extraction the combined effect of the wells and ambient flow influences the flow of the water and heat. During injection, water particles cannot move beyond the upstream stagnation point (x s , where q x = 0) 1 of the injection well; during extraction, water particles beyond the downstream stagnation point of the extraction well cannot be retrieved. Therefore, the specific discharge can be rewritten as: or, equivalently: 1 With u = dx/dt, the obtained ordinary differential equation can be integrated and solved with initial location x 0 at t = 0, to find the position of x as a function of time; The solution is only valid when the argument of the logarithm is positive, i.e. when x and x 0 are both on the same side of the stagnation point x s . The optimal distance between the well depends on the stagnation point, which in turn depends on the well discharge and screen length. For a given yearly storage volume, V [m 3 /y] and discharge Q [m 3 /h], the required distance between the infiltrating and extracting well also depends on the time it takes the well to inject and extract this volume and, of course, on the length of the storage period. An added complexity is that ATES wells only operate at maximum capacity Q during the coldest and warmest days, which is a limited time. During low to moderate energy demand, the wells do not pump at full capacity, which shifts the position of the two stagnation points closer to the wells. To fix the locations of the stagnation points it is assumed that ATES wells operate at a constant capacity, like is often done in ATES simulation (Bonte et al., 2013a;Sommer et al., 2013). This approach allows to iteratively identifying the optimal distance D between the upstream infiltrating and downstream extraction well. The iteration is done with the Newton-Raphson method (Thomas and Finney, 1984) The velocity of a water particle traveled to/from any well (i) can be expressed by: With x-x i the distance along the x-axis from the upstream well and y-y i along the y-axis, Q is the flow (injection positive, extraction negative), u x/y the ambient groundwater velocity in the x and y direction, n porosity, and L aquifer thickness. This expression can be solved by numerical integration, allowing tracking of water particles which are respectively infiltrated and extracted by an upstream and downstream well, with known ambient groundwater flow velocity and direction. The numerical solution is found by implementing the Euler method for temporal discretization using 4-order Runge-Kutta (Boyce and DiPrima, 1977).

Numerical model → MODFLOW
Next to analytical evaluation, the combined effect of hydrodynamic dispersion, conduction and advection requires numerical simulations that simultaneously solve the groundwater flow and transport equations. The numerical simulations apply the assessment framework and model of Bloemendal and Hartog (Bloemendal and Hartog, 2018). The simulation results are assessed relative to those for a single well. The simulations are carried out for the warm wells only, assumed to also be representative for cold wells. This is valid when there is no thermal interaction between warm and cold wells, and when the differences in both density and viscosity have a negligible effect on the behavior of the stored heat in the aquifer. These assumptions were validated by Bloemendal and Hartog (Bloemendal and Hartog, 2018).
As losses due to conduction, dispersion and displacement occur simultaneously, MODFLOW (Harbaugh et al., 2000) simulations are used to evaluate their combined effect on recovery efficiency. For the simulation of ambient groundwater flow and heat transport under various ATES conditions, a geohydrological MODFLOW model (Harbaugh et al., 2000) is coupled to the transport code MT3DMS (Hecht-Mendez et al., 2010;Zheng and Wang, 1999). These model codes use finite differences methods to solve the groundwater and (heat) transport equations. The use of groundwater wells for injection and extraction can thus be simulated along with the distribution of groundwater temperature, as was done in previous ATES studies e.g. (Caljé, 2010;Bonte, 2013;Visser et al., 2015). In the different modeling scenarios, the storage volume is varied between 12,000 and 300,000 m 3 with flow rates proportionally ranging from 8 to 200 m 3 / hour, screen lengths between 10 and 105 m, and ambient groundwater flow velocities between 0 and 50 m/y. This follows characteristics from Dutch practice, as will be introduced in the next section. Density differences are neglected, as this is considered a valid assumption (Caljé, 2010) for the considered ATES systems that operate within a limited temperature range (< 25°C). The parameter values of the model are given in Table 1, and the model was discretized as follows: -Model layers; the storage aquifer is confined by two 10 m thick clay layers. The storage aquifer is divided in 3 layers: 5 m thick upper and lower layers, and a middle layer for which thickness is set according to the required screen length of the modeled scenario. -The spatial discretization used in the horizontal direction is 5 × 5 m at well locations, gradually increasing to 100 × 100 m at the borders of the model. A sufficiently large model domain size of 6 × 6 km was used to prevent boundary conditions affecting (< 1%) simulation results. The gradually increasing cell size with distance from the wells results in a cell size of 15 m at 200 m of the well. This discretization is well within the minimum level of detail to model the temperature field around ATES wells, as was identified by Sommer (Sommer et al., 2014). -A temporal discretization of one week is used, which is sufficiently small to account for the seasonal operation pattern, and resulting in a courant number smaller than 0.5 within the area around the wells where the processes of interest occur. The simulation has a horizon of 10 years, which is sufficiently long to achieve stabilized yearly recovery efficiencies.
The PCG2 package is used for solving the groundwater flow, and the MOC is used for the advection package simulating the heat transport with a courant number of 1. To set the desired ambient groundwater flow velocity for the different scenarios simulated, the constant hydraulic head boundaries were used to set the required hydraulic gradient. An ATES doublet is placed in the aquifer with a well distance of five times the maximum thermal radius of the wells to avoid mutual interaction between the warm and cold storage volumes. In scenarios with groundwater flow, the ATES wells are oriented perpendicular to the flow direction.
Variations in weather conditions and building use affect the energy demand profile of ATES systems accordingly, which is of importance for the actual value of the thermal efficiency. Efficiencies are determined under 12 varying scenarios for both a weather-dependent and the regular energy demand profile, showing that the efficiencies of the corresponding conditions differ. However, they show the same relation according to the changes in conditions; the two sets of simulation results show a Pearson correlation coefficient of 0.97. Based on this evaluation, all simulations are done with one basic energy demand profile. To allow comparison with the analytical solutions, the constant storage volume energy demand pattern will also be used: heat injection, storage, extraction and again storage during 13 weeks each, as is commonly done in other ATES research (e.g. (Sommer et al., 2014;Zuurbier et al., 2013)).

Definition of thermal recovery efficiency for ATES systems
The thermal energy stored in an ATES system can have a positive and negative temperature difference between the infiltrated water and the surrounding ambient groundwater, for either heating or cooling purposes (Fig. 1). In this study, the thermal energy stored is referred to as heat or thermal energy; however, all the results discussed equally apply to storage of cold water used for cooling. As in other ATES studies (Doughty et al., 1982;, the recovery efficiency (η th ) of an ATES well is defined as the amount of injected thermal energy that is recovered after the injected volume has been extracted. For this ratio between extracted and infiltrated thermal energy (E out /E in ), the total infiltrated and extracted thermal energy is calculated as the cumulated product of the infiltrated and extracted volume, with the difference of infiltration and extraction temperatures (ΔT = T in -T out ) for a given time horizon (which is usually one or multiple storage cycles): With Q being the well discharge during time step t, and ΔT the weighted average temperature difference between extraction and injection. Injected thermal energy that is lost beyond the volume to be extracted is considered lost as it will not be recovered. To allow unambiguous comparison of the results, the simulations in this study are carried out with constant yearly storage and extraction volumes (V in = V out ).

Distance between the wells
3.1.1. Analytical/particle tracking results This analysis was carried out for several different groundwater flow velocities and pumping schemes, with the results given in Table 2. The results in Table 2 distinguish between the displacement during infiltration, storage and extraction.
The first and the last pumping schemes in Table 2(A & D) are beyond the situations in which ATES systems operate. During the night and in weekends, climate systems of buildings are often idle; during spring and autumn, energy demand may be so small that ATES systems are not running all the time. The ratio between storage and infiltration/ extraction depends on the building and its operation, as well as on climatic conditions. ATES systems typically have a limited operation at nighttime, when buildings are not used. During daytime, ATES systems may also often be inactive when energy demand is small, as buffer tanks prevent the ATES wells from alternating and pumping at low capacity. It is therefore safe to assume that the storage time is between one half to one third of a year (B-C scenarios in Table 2). The results in Table 2 show that during pumping, the required distance between the wells to counteract the groundwater flow is invariably about 2/3 of the expected groundwater flow during the pumping period (underlined values in Table 2). The groundwater is also moved when the ATES system is inactive, which then suggests that the ratio between storage and pumping time determines the required distance between the wells. The well distance between the upstream infiltrating and downstream extracting wells should be around 80% of the yearly groundwater flow, following the results for the B and C scenarios (bold) in Table 2. Corrected for the thermal retardation of 2, this becomes around 40% of u. Fig. 3 shows that all water particles are captured by the downstream well for the four scenarios indicated in Table 2.
The explanation for this constant 2/3, is that the average distance that all the particles travel in the direction of the groundwater flow is always 2/3 of the groundwater flow during the infiltration/extraction period. For example the B case in Fig. 3: the upstream distance of the particle on the xaxis is −57.4m, the downstream distance 90.7m, so the average distance is 16.6m, the infiltration time is 90 days, so the ambient groundwater displacement 24.6 m → 16.6/24.6 = 0.67. For example the D case in Fig. 3: upstream distance is −44 m, the downstream distance 110m, so the average distance is 33m, the infiltration time is 180 days, so the ambient groundwater displacement 49.3 m → 33/49.5 = 0.67.
The explanation and implications of these results are schematically depicted in Fig. 4. The pattern of flowlines of the upstream infiltration well during injection, are mirrored for the extraction well during extraction. The stored thermal energy is also displaced when the ATES system is inactive (e.g. during nighttime). This then results in a required distance of: With a and b corresponding to the distances indicated in Fig. 4, and t storage the time during which the ATES system is not operating. An important note is that this distance does not depend on the thermal radius, which is usually used as a metric to identify well distance for ATES wells (NVOE,   Table 2. Groundwater flow direction is to the right and velocity is 100 m/y. Fig. 4. Schematic representation of required distance between upstream and downstream well, for continuous infiltration followed by continuous extraction, and same infiltration and extraction volumes (V in = V out ).
2006; Li, 2014). This is caused by the fact that the upstream and downstream distances (a and b) do not depend on storage/ extraction volume and screen length, but on groundwater flow velocity and pumping time.

Basic numerical simulations of well distance to counteract groundwater flow
The optimal strategy for the optimal distance between the upstream and downstream well under different ambient groundwater velocities is to place them at a mutual distance of around 40% of yearly groundwater flow (0.4 ut with t = 1 year), in the direction of groundwater flow. However, due to existing infrastructure, such optimal placement is often difficult to achieve in most (dense) urban settings. In this section, numerical simulations are used to assess the effects of well distances which are either smaller or larger than optimal under ambient flow conditions. Fig. 5 shows the thermal well efficiency as a function of well distance (as a multiplier of the yearly groundwater velocity u), averaged over a 10year simulation period (lines) and as achieved during the last simulated year (markers). The simulations were done for a range of specific-storage volumes 2 and a range of ambient groundwater flow velocities encountered in practice by Bloemendal and Hartog (Bloemendal and Hartog, 2018). As with normal/single well storage, both the largest storage volumes and smallest groundwater flow velocity result in the highest efficiency. Efficiency is also less sensitive with the largest storage volumes and at the lowest groundwater velocity, compared to small ATES systems in aquifers with high groundwater flow velocity.
The calculated efficiencies confirm that the distance between the wells needs to be around the 0.4 u identified in previous subsection. Counteracting the groundwater flow also works when the wells are placed at larger mutual distance, as can be seen by the high and constant last year efficiency at D > 0.4 u (markers). But this has a negative effect on the total efficiency, as during the first year(s) of operation the heat does not yet reach the downstream well; this is reflected in the decrease of overall efficiency at larger mutual distance of the wells. This is illustrated in Fig. 6, which shows the temperature for the downstream well during a 20-year simulation period for different well distances (note that u = 100 m/y), where the reference refers to the case with a single well (distance is zero).
Using a downstream well is always efficient in situations with high ambient groundwater flow. For larger distances, e.g. 2.5 u in Fig. 6, it takes about 10 years to reach the highest efficiency. The efficiency can be derived from the extraction temperature of the downstream well always being lower than the 15°C injected. However, the last year efficiencies (markers) show a flat optimum at D > 0.4 u in Fig. 5, which indicates that efficiency hardly decreases when well distance increases. Only in Fig. 5D subplot is a small decrease visible. The limited magnitude of these losses is partly caused by the model set-up for these basic simulations. The simulation of specific-storage volumes (V per meter screen length) allowed using thin model layers (which saved computational resources), but in turn results in relative larger losses to confining layers. Simulations with conditions closer to practice show a steeper decline in last year efficiency, which is to be expected due to the increased area over which conduction losses occur in the aquifer (i.e. across the circumference area).
The groundwater flow causes the stored heat in the aquifer to have a stretched shape, which in turn results in an increasing total area over which losses occur. The A/V-ratio can, therefore, be used as an indicator for the increase in losses with increasing well distance. 3 The A/ 2 To allow for easy comparison between ATES systems sizes, the specific storage volume is used, which is defined as the storage volume per meter well screen (V/L). For the Dutch situation, the specific storage volume varies between 2000 to 10,000 m 3 /m/y. 3 A/V-ratio is a metric for the expected conduction losses (Bloemendal and Hartog, 2018). The yearly storage volume must be applied here, and not the volume that is actually stored over multiple cycles in the aquifer (both are the same for the cylindrical case (Bloemendal and Hartog, 2018)). Applying the volume of the stretched cylinder for V results in a smaller and constant A/V-ratio of 1/R th + 2/L at increasing well distance, instead of the 2/R th + 2/L for cylindrical shaped storage volume (Bloemendal and Hartog, 2018).
V-ratio increases linearly with increasing well distance; the simulation results indicate that the optimal well distance has a relatively flat optimum for D > 0.2 u, but it is best to try to stay as close to the recommended 0.4 u as possible.

Well screen length
The results in Fig. 5 were obtained for a range of specific storage volume of ATES systems. In practice, aquifer thickness and storage volume, which both affect the thermal radius, will vary from place to place. Therefore, it is important to obtain insight in how mutual distance, thermal radius and screen length affect the heat losses. Bloemendal and Hartog (Bloemendal and Hartog, 2018) and Doughty et al. (Doughty et al., 1982) used the geometric ATES well property L/R th (i.e. screen length over thermal radius) to assess conduction losses of ATES wells. The thermal radius (R th ) is defined as: The size of the thermal cylinder thus depends on the storage volume (V), screen length (L, for a fully screened aquifer), porosity (n) and water and aquifer heat capacity (c w , c aq ). Fig. 7 shows the results of the last year's efficiencies for scenarios as a function of L/Rth. The differences in last year's efficiency among the distance scenarios are small at least for D > 0.2 u. In addition, although heat losses increase with well distance, their magnitude remains relatively small as was already noticed in Fig. 5. Fig. 7 shows that changing the wells' screen length has a large influence on the efficiency. Fig. 7 shows that in the reference case (D = 0) the L/R th -ratio has its optimum around 0.1, confirming that with singlewell operation in high ambient groundwater flow, the thermal radius needs to be large (Bloemendal and Hartog, 2018). The results for (D > 0.2 u) show that to counteract ambient groundwater flow, the screen lengths need to be smaller compared to the 1-4 L/R th bandwidth identified by Doughty et al. (Doughty et al., 1982). The required well distance has a flat optimum around 0.5 L/R th . This is similar to the values found by Bloemendal and Hartog (Bloemendal and Hartog, 2018), which they show in their Fig. 15 for higher groundwater velocities with single wells.

Pumping scheme of the upstream and the downstream well
This subsection discusses how the pumping schedule affects the overall efficiency. The basic scheme for counteracting the groundwater flow analyzed in previous sections is to only infiltrate in the upstream well and only extract from the downstream well. However, as discussed in the introduction, ATES systems with multiple wells have to use both wells during the coldest and warmest day to provide heating and cooling to the building. To allow for flexibility during operation, different pumping schemes for upstream and downstream wells are analyzed next, to identify how such schemes may affect optimal well distance and/or screen length.

Predefined pumping schemes
Four different pumping schemes were evaluated: -The "1/0" scheme, in which total storage volume is infiltrated upstream and extracted downstream -The "0.9/0.1" scheme, in which the second well is only used to cover peaks, so that 90% of the storage volume is infiltrated upstream and extracted downstream -An intermediate scheme (0.75/0.25) in which the second well is used more frequently. In this Scheme 75% of the volume is infiltrated upstream and extracted downstream -The usual operation mode for a two-doublet ATES system in a nonflowing aquifer, which stores and extracts 50% of the storage volume in each well For ambient groundwater velocities of 50 and 100 m/y, Fig. 8 shows the recovery efficiency achieved in the last year of the 10-year simulation period for different distances and pumping schedules. This shows that the pumping scheme is more important with higher ambient groundwater velocity. At a small mutual distance between the two warm wells, it is best to stay close to the 1/0 pumping scheme defined above, but the efficiencies for a pumping regime of 0.9/0.1 or 0.75/ 0.25 are almost the same. At D > 0.4 u (and R th < D), the heat infiltrated in the upstream well does not reach the downstream well within one storage cycle. Fig. 9 shows that some (25% of V) infiltration in the downstream wells helps improve the downstream well's efficiency during the first years of operation, resulting in higher extraction temperatures. Fig. 6. Extraction temperatures of the downstream well for different well distances. u = 100m/y, constant injection temperature of 15°C, sinusoidal injectionextraction of volume of 250,000 m 3 in an aquifer of L = 40 m with thermal retardation of 2. Fig. 7. Recovery efficiency as a function of geometric properties of the storage volumes at different mutual well distances, for 100 m/y groundwater and a yearly storage volume of 250,000 m 3 /y. Only the screen length of the well is changed within the different simulation distance scenarios.

Automatic control of the applied pumping scheme
During operation in practice, the downstream well is also needed for injection and the upstream well for extraction of heat, in order to cover maximum demands. Because the pumping scheme has an impact on recovery efficiency, it is evaluated in more detail. However, identifying the optimal pumping strategy is complex: the number of possible pumping strategies is practically infinite, and under practical working conditions the required pumping capacity cannot be predicted, as it follows the building's energy demand − which in turn depends on highly uncertain factors such as weather conditions and the use of the building. As a result of differences in heating and cooling demand, the pumping scheme may commonly be asymmetrical.
The objective for identifying the optimal pumping scheme is to avoid the "escape" of heat away from the downstream well. Therefore, in the MODFLOW model, a temperature monitoring point downstream of the downstream well is used to automatically control the pumping scheme. The 0.5/0.5 pumping schedule is applied to the wells as a basic pumping scheme, which is then altered by the controlling temperature measured at the monitoring point at run-time. When the temperature in the monitoring point increases above a given threshold, pumping at the next time step is set to either a higher percentage of infiltration in the upstream well, or to a higher percentage of extraction in the downstream well, depending on whether the wells are in heating or cooling mode. The downstream distance of this monitoring point from the downstream well should not be farther than the downstream stagnation point. It must also be within the maximum thermal radius of the downstream well to provide rapid feedback.
Simulations were done with this control scheme for varying storage volumes (125,000; 250,000 and 500,000 m 3 /y) and varying groundwater flow velocities (50; 75; 100; 150 and 200 m/y)). The optimized pumping schemes obtained from these simulations are presented in Fig. 10; they show that schedules vary from 0.7/0.3 to 0.9/0.1. The temperature threshold or trigger temperature of the temperature monitoring point had virtually no effect on the obtained pumping schemes. The following observations and conclusion appear from the results in Fig. 10; -The larger the ambient groundwater flow velocity is, the closer the pumping scheme needs to be to the 1/0 scheme. This is logical, as the efficiency is more sensitive to the groundwater flow at higher velocity values. -The larger the storage volumes/thermal radii, the closer the pumping scheme needs to be to the 1/0 scheme. This is counter-  intuitive, as efficiency is generally less sensitive to disturbances at larger storage volumes. This is caused by the fact that at larger storage volumes, the temperature threshold in the monitoring point is exceeded earlier. This indicates that the location of the monitoring point should, therefore, consider the expected/required amount of time the ATES system needs to run on both wells. -In all situations, the pattern between infiltration and extraction is more or less symmetric, with a light preference to more infiltration upstream. This can be explained by the fact that during infiltration, the temperature threshold will be reached sooner in the monitoring point.
This type of automatic control works well because the groundwater flow responds virtually instantaneously to the control signal, which prevents oscillations.

Case study
A case study is used to illustrate how the separately discussed aspects come together in practice. The city of Apeldoorn, in the east of The Netherlands, is situated on inclined 4 aquifers with large head gradients, in combination with high permeable (k h = 30 m/d) sandy aquifers. The groundwater flow velocity in the local aquifer used for ATES in Apeldoorn is about 35 m/y. At a site with multiple office buildings, a joint ATES system with two doublets and a total seasonal storage capacity of 425,000 m 3 is situated in this aquifer. The wells have a screen length of 40 m and a capacity of 175 m 3 /hr each; the mutual distance between the wells is 200 and 150 m for the warm and cold wells respectively (5.7 u and 4.3 u), and both sets are aligned in the (expected) direction of the groundwater flow (Groot, 2013). To prevent thermal interaction between the warm and cold wells, the sets of warm and cold wells are separated over 200 m perpendicularly to the groundwater flow direction. All wells are required to supply thermal energy to the different buildings during the warmest and coldest periods of the year.
The simulation results of this case in Fig. 11A shows that it takes 15 years before the extraction temperature of the downstream well stays above 13°C (with an infiltration temperature of 15°C), for different pumping schemes. An optimal pumping scheme of 0.65/0.35 is identified using downstream temperature monitoring to control the pumping scheme for the case study, Fig. 11C and D. The relatively large mutual distance between the wells causes much lower efficiency compared to the values indicated in Fig. 5C and D. The overall efficiency is of course low, due to the many years it takes before the heat reaches the downstream well. The temperature distribution in Fig. 11B and D also explains why the last year efficiency is much lower; the surface area of the "footprint" of the warm water between the wells increases considerably, compared to the situation where wells are located closer to each other. This results in conduction losses at the boundary of the warm zone in the aquifer and to confining layers, lowering efficiency. Contrary to the findings in the previous subsection, at very large mutual distance (> 2.5 u) the increased size of the warm zone in the aquifer, therefore, results in such large conduction losses, that a pumping scheme closer to 0.5/0.5 works better.
The wells were not placed closer together in this case as A) this research on identifying optimal distance was carried out after the ATES system in Apeldoorn was constructed, but mainly as B) the aquifer is phreatic, in which a small mutual distance of the wells would result in groundwater head changes that could affect the flora in the adjacent nature reserve. This illustrates that other design aspects may also affect choices in well design/location, which, sometimes unavoidably, result in a suboptimal well lay-out − as is the case at this ATES system. However, the principles elaborated in this section do allow for technical optimization with respect to the geohydrological conditions in Apeldoorn.

Discussion
The analysis in this paper was done with a regular (sine-shaped) energy demand profile of the building and pre-defined pumping schemes. However, in practice, the operation of the wells is not as straightforward as presented here. The automatic adjustment of the pumping scheme at runtime based on downstream temperature measurements showed to be an efficient way to optimize the pumping scheme. At the cost of an extra monitoring well with a temperature sensor and transmitter, such a direct control is applicable in practice and will also reduce heat loss (heat pollution) downstream. An important aspect to consider is that the location of the monitoring point should consider the expected volume to be pumped in the downstream well and the stagnation point of the well at low flow rates. Making smart choices during well installation may also help improve efficiency. With the presently applied sine-shaped energy demand profile and a combined pumping capacity of 100% for both wells, around 30% of the required volume has to be pumped at the "wrong" well. When installing wells with a capacity of 75% of the maximum required capacity instead of 50%, the percentage of pumping at the wrong well can be reduced to only around 10%. Especially for larger systems, and ATES systems in areas with groundwater flow velocity larger than 100 m/y, this considerably improves the recovery efficiency.
In cases where flow direction and velocity change considerably, it is difficult to anticipate on with ATES well design and/or operation.
Despite the legal periodic requirement for energy balance, ATES systems never have exactly the same heating and cooling demand. This affects how efficiently the displacement of heat by groundwater can be compensated. However, under single doublet and no ambient flow conditions, an imbalance would affect operation similarly; more water is extracted from one of the well types, while thermal energy remains in the subsurface at the other well. This then results in a higher efficiency in one, and lower efficiency in the other well. When energy balance is not periodically met, one well type will completely recapture the infiltrated thermal energy, while the other will have a thermal plume downstream of the downstream well. In the condition with high ambient groundwater flow, such situations may even be preferable (when there are no other interests downstream), since the heat accumulation in one of the wells will not affect the other well in the long run.
A detail that was not discussed in this study is the effect of simultaneously injecting and extracting groundwater in the other well pair, regarding the distribution of the warm and cold water in the aquifer. Theoretically, the crosswise injection and extraction pattern, as schematically depicted in Fig. 2, means that warm and cold groundwater is also pulled towards the downstream well of the other well pair. As such, collecting all the warm and cold groundwater requires the two downstream wells to be closer together than the upstream wells, as exaggerated and schematically represented in Fig. 12A. It is difficult to generalize the extent to which this mechanism affects the distribution of warm and cold groundwater, as it depends on local conditions and the position of the wells with respect to each other. For the averaged sized doublet, this appears to be negligible, as the change in efficiency of the warm well is smaller than 1% compared to the situation where the cold wells were not present in the model. This is also shown by the temperature distribution in Fig. 12B, where the wells are aligned in the direction of groundwater flow. This confirms the validity of the simulation approach in which only one of the well types was incorporated in the model.
At large distances between the wells, a deviation in flow direction of a few degrees may already result in a situation in which the downstream well can only partly recover the heat. Horizontal anisotropy may often cause the groundwater flow direction to be misinterpreted.
Several exploratory simulations were carried out to assess this effect. The results show that a misinterpretation of the flow direction of 25°r educes efficiency by 7%, 17% and 29% for well distances of 0.4, 1 and 2 times the groundwater flow velocity, respectively. This confirms that at small distances between the wells, the recovery efficiency is less  sensitive to the ambient groundwater flow direction.
The cost of additional wells may render the application of multiple wells unfeasible, at least for small ATES systems or buildings. For the smallest buildings, this may thus not be affordable. It is therefore wise to develop collective ATES systems in areas with high ambient groundwater flow. For aquifers with sufficient thickness, using two monowells aligned in the direction of the groundwater velocity may be a cost-effective substitute for two doublets.

Conclusions
This paper provides guidelines for the optimal well placement and operation of ATES wells to counteract ambient groundwater flow. This can be done by using multiple warm and cold wells aligned in the direction of the groundwater flow, and is a suitable alternative for the recirculation systems (Fig. 1) which are often applied in areas with high ambient groundwater flow velocity. The optimal distance between the wells is 0.4 times the yearly groundwater flow velocity (u). The well distance should at least be 0.2 u, while a larger mutual distance has a considerable negative effect on thermal efficiency during the first years of operation. Also, at larger distances, the heat losses to confining layers and surrounding aquifer increase due to the increase of overall area over which losses occur. Therefore, the well distance should be kept as close to 0.4 u as possible. Another advantage of a small mutual distance is that a deviation of the expected direction of the groundwater flow does not have a large impact on the efficiency of the system. Screen length should be chosen shorter compared to when there is no groundwater flow: the ratio between screen length and expected thermal radius should be around 0.5, in order to minimize conduction losses.
While using two wells to counteract the groundwater flow, it is always best to infiltrate upstream and extract downstream. However, during energy demand peaks in mid-winter and mid-summer, it may be needed to also use the upstream well for extraction and the downstream well for infiltration. The highest efficiency is achieved when the upstream well is used preferentially for infiltration and the downstream well for extraction, while using the other wells as little as possible during periods when the required discharge is larger than a single well can deliver.
At low groundwater flow velocities and small ATES systems, the efficiency is less sensitive to the proposed optimal distance and to the pumping strategy. The higher the groundwater flow and the larger the thermal radius of the ATES system, the more important it is to stay as close to the 0.4 u distance, and to the pumping scheme where infiltration is in the upstream well and extraction in the downstream well.
The used code can be send to you on request to the corresponding author, running the code requires flopy: https://modflowpy.github.io/flopydoc/ introduction.html.