Assessing the potential of low-transmissivity aquifers for aquifer thermal energy storage systems: a case study in Flanders (Belgium)

The Member States of the European Union pledged to reduce greenhouse gas emissions by 80–95% by 2050. Shallow geothermal systems might substantially contribute by providing heating and cooling in a sustainable way through seasonally storing heat and cold in the shallow ground (<200 m). When the minimum yield associated with the installation of a cost-effective aquifer thermal energy storage (ATES) system cannot be met, borehole thermal energy storage, relying mostly on the thermal conductivity of the ground, is proposed. However, for large-scale applications, this requires the installation of hundreds of boreholes, which entails a large cost and high disturbance of the underground. In such cases, ATES systems can nevertheless become interesting. This paper presents a case study performed on a Ghent University campus (Belgium), where the feasibility of ATES in an area with a low transmissivity was determined. The maximum yield of the aquifer was estimated at 5 m3/h through pumping tests. Although this low yield was attributed to the fine grain size of the aquifer, membrane filtering index tests and long-term injection tests revealed that the clogging risk was limited. A groundwater model was used to optimize the well placement. It was shown that a well arrangement in a checkerboard pattern was most effective to optimize the hydraulic efficiency while maintaining the thermal recovery efficiency of the ATES system. Hence, for large-scale projects, efficient thermal energy storage can also be achieved using a (more cost-effective) ATES system even in low-permeability sediments.


Introduction
Shallow geothermal systems have proven to be locally available, green, and renewable alternatives to fossil fuels both for cooling in summer and heating in winter (Perego et al. 2020).On average 0.5 kg of CO 2 per m 3 of pumped water can be saved compared to conventional technologies (Fleuchaus et al. 2018).Implementing such systems in the building sector has the potential to significantly reduce global greenhouse gas emissions (European Commission 2012, 2019;Ramos-Escudero et al. 2021).
Both borehole thermal energy storage (BTES) systems and aquifer thermal energy storage (ATES) systems are shallow geothermal systems that make use of a heat pump to extract the heat out of the subsurface reservoir.
BTES systems are closed-loop systems that use ground heat exchangers in the subsurface i.e. long loops through which water, sometimes mixed with an antifreeze, circulates (Fig. 1).Their capacity is mostly dependent on the thermal conductivity of the ground and its capacity for thermal recharge (Bayer et al. 2012;Hecht-Méndez et al. 2013;Glassley 2015).In contrast, the potential for ATES systems is mostly dependent on the hydraulic conditions (hydraulic conductivity and hydraulic gradient) in the aquifer.These are open-loop systems which extract and inject groundwater from an aquifer through a well (Bloemendal et al. 2015).
For buildings with high energy demand, ATES systems are generally favoured over BTES systems as the costs of the drilling become more important compared to the other costs (pipework, controls, etc.) for similar performance.A BTES system for that kind of application would require the installation of tens to hundreds of boreholes.
In Flanders (Belgium), the potential of aquifer layers for installing ATES systems has been estimated based on their transmissivity (WTCB 2017).When the transmissivity is below 50 m 2 /day, it is deemed unsuitable, while above 250 m 2 /day, the potential is recognized (WTCB 1 3 2017).In between those values further investigations are recommended.The threshold on the transmissivity implicitly accounts for the fact that a minimum yield of about 10 m 3 /h is required to justify the investment costs (Bloemendal et al. 2015;Hermans et al. 2018).Following these recommendations, the implementation of ATES systems was deemed unfeasible in many areas (Fig. 2).In these areas, BTES systems, which are less dependent on the variability of subsurface properties, are seen as the only viable option, even to fulfil a high power demand.BTES systems in these circumstances nevertheless require significant investment costs and a large available surface area.
The interest in ATES systems is increasing, pushed by high energy prices and government incentives to invest in green energy.More specifically, geothermal energy is labelled the most attractive option for renewable energy production according to Batac et al. (2022).This will inevitably add additional stress to valuable subsurface systems, and specifically to the most promising reservoirs that already sustain a freshwater supply.Considering this, more complex reservoirs, such as low-transmissivity aquifers, will also be targeted as the technology for ATES systems evolves.Targeting low-transmissivity aquifers could also be beneficial to reduce the effect of buoyancy on the storage efficiency of ATES systems operating at a higher temperature (Beernink et al. 2022).In general ATES systems produce more energy per well resulting in fewer and less deep drillings compared to a BTES system.Therefore, especially when high investment costs for a BTES system are at play, it becomes nevertheless interesting to further explore the feasibility of ATES systems even in low-transmissivity aquifers.However, because of its dependence on aquifer properties, the implementation of an ATES system is more complex and more prone to failure.Hence, there is a need for methods to assess the feasibility of ATES in aquifers with low transmissivity.
This paper demonstrates the feasibility of an ATES system functioning at a low pumping rate (5 m 3 /h) on Campus Sterre, at Ghent University in Belgium.This was accomplished by first estimating the maximum yield of the aquifer through pumping tests.In addition to this, considering the fine grain size of the aquifer, membrane filtering index (MFI) tests and long-term injection tests were carried out to estimate the clogging risk.For the actual ATES system, interactions between the wells might result in excessive pressure differences potentially flooding the surface and damaging the wells and the confining layer.Interactions could also cause a short circuit between the warm and cold well areas (i.e. a thermal breakthrough).To limit these risks, a groundwater model calibrated with the pumping test data was used to optimize the well placement.

Setting of the study area
Ghent University aims to become CO 2 -neutral on the Faculty of Science campus (Campus Sterre) by 2050.Reaching this objective is challenging, and different sustainable, green, innovative technologies should be studied and evaluated.Based on an energy audit, the most sustainable alternative to cover the heating and cooling demand on this campus was to combine the residual heat from cooling the servers with a shallow geothermal system.The geothermal system would store the residual heat of the servers during the summer period when heating is not needed, so it can be released in winter in addition to the directly used residual heat from the servers.The power of the geothermal buffer was estimated at 0.63 MW.
A quick evaluation of the study area showed that, because of the absence of a thick productive aquifer at the campus, the transmissivity of the available aquifers might not be sufficient to reach a high enough pumping rate (Fig. 2).Alternatively, to cover the power demand in the study area, a BTES system of 175 boreholes of 100 m depth was proposed.As this would result in an investment cost of approximately € 3,500,000 and a large surface area occupied by the borehole heat exchangers, it was decided to determine the potential for ATES.
A description of the hydrogeology of the study area was made by Lebbe et al. (1992).Quaternary deposits, consisting of clays, silt, sand and gravel with a thickness of 9.5 m and variable hydraulic conductivity, can be found at the surface, whereas below, the Formation of Gentbrugge is considered a confining layer.It has an irregular extension: it thickens towards the northeast and disappears in the southwest of the study area.It is composed of silty clay or clayey silty glauconiferous very fine sand.Sand lenses with organic material and small pyritic concretions as well as layers of clayey sandy coarse silt might occur.
The main groundwater reservoir is situated below this formation.It can be subdivided into 6 units according to Lebbe et al. (1992) (Fig. 3), labelled Yd 6-Yd 1 from top to bottom.Yd 6 consists of slightly clayey glauconitic fine sand which might contain small shell fragments.The lithology of Yd 4 and Yd 2 resembles the one of Yd 6, without shell fragments.These three units are considered aquifer units and have a cumulative thickness of 18.5 m, an average hydraulic conductivity of 1.08 m/day, and hence an average transmissivity of about 20 m 2 /day.Yd 5, Yd 3 and Yd 1 contain more clay and are only semi-pervious: Yd 5 is very sandy clay, Yd 3 is sandy clay to clay, and Yd 1 is a sandy and silty clay with intercalations of thin clayey fine sand beds.This aquifer unit is bounded from -40 mTAW (Tweede Algemene Waterpassing, groundwater level in relation to sea level) by the clayey Formation of Kortrijk, which is up to 95 m thick.

Pumping and injection
The hydraulic efficiency of an ATES system relies on the hydraulic conductivity of the subsurface which governs the maximum pumping rate and the corresponding drawdown.Lebbe et al. (1992) conducted a triple pumping test on Campus Sterre to determine the hydraulic parameters of the different layers.Pumping wells were drilled in the pervious layers Yd 2, Yd 4, and Yd 6 (named PP 2, PP 4 and PP 6 respectively).Each pumping well was accompanied by three observation wells located at fixed distances from the pumping well (Fig. 4)-for detailed information on the triple pumping test see section S1 in the electronic supplementary material (ESM).For this study, additional field tests were carried out to: 1. Estimate the maximum pumping rate in the entire aquifer system.In the absence of a fully filtered pumping well, this was estimated in each pervious layer separately.The total maximum pumping rate can be estimated by summing the individual rates if all pervious layers can be considered fully confined and independent of each other.2. Estimate the maximum injection rate.For the injection test, a fully filtered well (PB JE) was used (Fig. 3).However, this well was constructed as a piezometer and has a limited inner diameter (63/57 mm).As such the well efficiency is not optimal, limiting the injection capacity.
Because the water level in the aquifer is shallow (~3 m below the surface), the injection might also cause water pressures above the ground level, potentially causing flooding and/or instability of the confining layer between Yd 6 and the top layer.3. Simulate the long-term stability of a well pair consisting of one injection and one pumping well mimicking the behaviour of an ATES system.4. Generate data sets for validation of the groundwater model.
In practice, the pumping and reinjection tests took place in the same wells that were used by Lebbe et al. (1992).In August 2021 the maximum pumping rate in each pumping well was estimated.To simulate the long-term stability of a well pair, PP 4 was selected as the pumping well and PB JE as the injection well.The test started in October 2021 with a constant flow rate of 3.8 m 3 /h based on the limited capacity of the injection well.Pressure transducers were installed in PP 2, PB 2.1, PP 4, PP 6, PB 6.1, PB 6.2 and PB JE to record the hydraulic head and temperature-for more details on these tests, see section S2 in the ESM.

Membrane filtering index
Because the aquifer contains fine particles, the (injection) wells could clog relatively rapidly.As a result, the injection pressure can increase over time and therefore decrease the injection capacity.For ATES systems, the wells are alternately injection or pumping wells according to the season.As such, if the capacity of the wells would decrease over time due to clogging, this would be detrimental to the overall efficiency and capacity of the ATES system (Jenne et al. 1992;De Zwart 2007).
Therefore, a membrane filtering index (MFI) test was carried out on-site to determine the clogging risk.It is a measure of the rate at which a filter paper (0.45 µm) becomes clogged under constant water pressure (2 bar) (Schippers andVerdouw 1979, 1980;Olsthoorn 1982).The index can be derived by plotting the ratio of the filtration time and the filtered sample volume (t/V) as a function of the total filtered volume (V).When the slope of the curve is less than 10, the water purity is considered acceptable for reinjection purposes for ATES systems.When the slope of the graph is less than 3, the water purity is considered excellent for reinjection (Schippers and Verdouw 1980;Olsthoorn 1982).The actual MFI test for this study was carried out twice from 11h45 to 13h00 on 28 October 2021, after the well was sufficiently developed and the sand content (>70 µm) of the pumped water was visually analysed (Aalten and Witteveen 2015).Next, a long-term pumping and reinjection test (25 October-20 December) served to verify whether there was a decrease in injection capacity with time.

Modeling approach
Because of the limited maximum yield of the aquifer, an ATES system with several well pairs each operating at a flow rate close to the maximum yield is required to fulfil the energy demand.The system must therefore be carefully designed to avoid hydraulic and thermal interaction between the wells.Creating a numerical model is a viable and indispensable tool to assess the feasibility of the project.It not only helps in understanding and predicting the behaviour of complex systems but it also helps to optimize the desired project implementation (Yapparova et al. 2014).
For this project, the freely available USGS MODFLOW 6 software was used (Langevin et al. 2017a, b) together with the software ModelMuse as a graphical user interface (Winston 2019).MODPATH was used to simulate advective transport (Pollock 2012) and MT3D-USGS was used to model the full transport processes (Bedekar et al. 2016).MODFLOW 6 uses the control-volume finite-difference method to solve the mathematical equation numerically (Langevin et al. 2017a).
Based on the hydrostratigraphic setting of the study area shown in Fig. 3 and the already calibrated model parameters by Lebbe et al. (1992), a 3D model of 5 × 5 km around Campus Sterre was made (Table 1).All layers, except for the Quaternary, were set to be confined in the model.A structural grid was used for the spatial discretization (DIS package) (Langevin et al. 2017a).To improve the solution but avoid an exaggerated computation time, the grid was only refined where a steep gradient is expected, i.e. around the pumping/injection wells.The largest grid size was 100 m and was decreased to 5 m, 0.5 m and finally approximating the drilling diameter of the wells (roughly 0.25 m).For heat transport, the grid had the same extent but the smallest grid size was set to 1 m around the well area to limit the computational time.The previously smaller cell size was needed for calibration, to assess the wellbore effect on the accuracy of the numerical model, and to limit numerical dispersion for the advective transport simulation.
For the bottom of the model, a no-flow boundary was set, as below Yd 1 the aquitard corresponding to the Kortrijk Formation is present.The northern, eastern, southern and western boundaries were set as constant head boundaries with a hydraulic head of + 10 mTAW for the calibration period.For the long-term simulations, a hydraulic gradient of 0.14%, deduced from monitoring wells located outside of the study area, was imposed by assigning a constant head to the northern and southern boundaries.At the start of each simulation, an initial head of + 6.97 mTAW was chosen.This is the hydraulic head that was measured in PP 4 before the installation of the diver in October 2021.A zero-dispersion/diffusion heat flux was imposed for the transport boundary conditions.Only transient simulations were used due to the interest in the evolution of the groundwater level/thermal storage with time.The boundary conditions remained constant during the simulation period and no atmospheric heat loss was considered.
Because of the similarity between solute and heat transport and because of the disregarding of density/viscosity effects, MT3D-USGS can be used to model heat transport processes (Zheng 2010;Hecht-Méndez et al. 2010;Sommer et al. 2013;Possemiers 2014).The heat transport by conduction and thermal dispersion is analogous to molecular diffusion and mechanical dispersion in the solute transport equation.Through the buffering effect of conduction, the value of thermal dispersion is small in comparison with mechanical dispersion and can therefore be neglected as a heat transport process (Hopmans et al. 2002;Vandenbohede et al. 2011).Next, the transport of heat by groundwater flow is analogous to the advection term in the solute transport equation (Zheng 2010).To implement this in MT3D-USGS, the (thermal) distribution coefficient (K d t , in m 3 /kg) and the molecular diffusion coefficient (D m t in m 2 /s) must be defined as follows (Zheng 2010): Table 2 provides an overview of the parameters.It is important to mention that water density was considered as a constant (1) (2) because the temperature changes will remain limited (<15 °C) so that free convection (driven by density differences) is negligible compared to forced convection (due to the imposed hydraulic gradient) (Zuurbier et al. 2013;Zeghici et al. 2015).
To solve the heat transport equation, the third-order TVD (total variation diminishing) method with a backward approximation was used (Zheng and Wang 1999).The method of characteristics (MOC) was also tested in this study, which simulates more gradually thermal diffusion but it did not change the conclusions.The generalized conjugate gradient (GCG, convergence criterion of 10 -10 ) solver was used to implicitly solve the dispersion, sink/source and reaction terms with a finite-difference method.Finally, a linear sorption isotherm was selected in the RCT package, which accounts for the heat transfer process between the fluid and the solid (conduction).This conduction results in a retardation of the movement of the warm/cold temperature plume in comparison to the average linear groundwater flow velocity (Zheng and Wang 1999;Vandenbohede et al. 2011).

ATES system design
To design an ATES system, first the energy demand of the building, both for heating and cooling, that must be fulfilled by geothermal energy must be determined.With this knowledge, the total volume of water to be extracted can be estimated based on the heat capacity of water.This is illustrated by the following equations, ignoring the coefficient of performance of the heat pump (Glassley 2015): where E (kWh) is the thermal energy that can be stored/ extracted from a given volume of water V (m 3 ), c w is the volumetric heat capacity of water (1.16 kWh/m 3 K), ∆T the temperature difference between the extracted and injected water (K), Q the total flow rate of the system (m 3 /h), t the time (h), and P the power (kW).As such, knowing the maximum pumping/injection rate for a single well pair, the required number of well pairs can be calculated.
Knowing the number of wells, optimal use of (sub)surface space and optimal thermal recovery(/storage) efficiency should be targeted.One of the parameters influencing this performance, besides the groundwater flow and pumping rate, is the well placement (Yapparova et al. 2014).Previous studies showed that storage efficiency decreases with decreasing distance between the warm and cold well areas and increasing hydraulic conductivity (Kim et al. 2010;Yapparova et al. 2014).As such, the relatively low hydraulic conductivity in the study area might be an advantage in this project by optimizing the usage of space without compromising on storage efficiency.The ATES storage efficiency decreases with increasing hydraulic conductivity because it increases the risk of a short circuit between the cold and warm well areas.This is also often called a thermal breakthrough (Kim et al. 2010;Gao et al. 2013;Yapparova et al. 2014;Bloemendal et al. 2018).In theory, it is sufficient to ensure that the distance between the cold and warm wells in the design of the system is sufficient.This safe distance can be estimated from the thermal radius of influence (R th ) (Bloemendal et al. 2018;Bloemendal and Olsthoorn 2018).The distance between wells of the opposite and same type should be 2.5•R th and 1•R th respectively (Bloemendal et al. 2018).To apply this formula in practice, first, the hydraulic radius of influence was estimated analytically using the Thiem-Dupuit method for a confined aquifer in steady-state (Dupuit 1863;Thiem 1906) where R is the hydraulic radius of influence (m), s the drawdown at a certain distance x (m) from the well due to pumping/injecting (m), K is the horizontal hydraulic conductivity of the medium (m/s), e the thickness of the groundwater reservoir (m), and Q the pumping/injection rate (m 3 /s).Finally, using the earlier defined thermal retardation factor, the thermal radius of influence (R th = hydraulic radius/ thermal retardation factor) can be estimated.
Following the guidelines of well placement for ATES systems of Bloemendal et al. (2018) will result in a clustered well placement with one warm and one cold cluster (each composed of several wells) which are alternately pumping or injecting.While this type of well placement might be a better option considering the thermal recovery efficiency, the opposite (i.e.alternating all injection and pumping wells) might be beneficial for the hydraulic head.The hydraulic efficiency is closely related to the superposition principle in a confined aquifer implying that the resulting drawdown at a location is the algebraic sum of the effect of multiple pumping/injection wells in the neighbourhood.On the one hand, it must be able to maintain a low enough pressure in the injection wells to avoid increasing the risk of flooding and damaging the confining clay layer.On the other hand, it must be able to sustain all pumping wells with a sufficiently high pumping rate as the pumping wells will operate close to their maximum capacity.An alternating well arrangement where the pumping and injection wells are placed close enough so they can hydraulically interact could in this way help to reduce excessive head changes.Both the clustered and alternating well arrangements were evaluated for this project.

Efficiency assessment framework
To objectively determine the feasibility of an ATES system in the low transmissivity aquifer the overall efficiency was evaluated using two key parameters which were each qualitatively and quantitatively assessed.The first parameter is the hydraulic efficiency of the system which relates to the risk of flooding of the surface or damaging the borehole vicinity due to excessive pressure differences in the wells (NVOE 2006).The hydraulic efficiency of the system was qualitatively assessed by comparing the hydraulic head after a pumping-injection cycle of the ATES system to the average natural groundwater level of about 7.25 mTAW in the study area.Quantitatively, the change in hydraulic head (in m) due to the operation of the ATES system was compared to the Dutch ATES design standards (NVOE 2006): where the filter top depth corresponds to the top of Yd 6 (9 m) (Table 1).When a suitable arrangement was found in terms of hydraulic efficiency, the transport simulation was done to evaluate the thermal recovery(/storage) efficiency of the ATES system (Fleuchaus et al. 2020).This is the second key evaluation parameter.The simulation covers a period of 20 years in which warm water injection is assumed to take place in summer (6 months) and cold water injection is assumed to take place in winter (6 months), both at the maximum rate.These long periods represent an extreme scenario.The temperature of the injected water was imposed in the model while the temperature of the extracted water varies according to heat transport processes.The injection temperature is normally approximately +5 °C for the warm well area and -5 °C for the cold well area (relative to the natural groundwater temperature of 13.8 °C); furthermore, a yearly balance between heating and cooling demand from the underground was assumed.
In the first instance, the thermal recovery efficiency was qualitatively assessed by confirming that there is no thermal (6) Maximum allowable hydraulic head change = 0.2 × f ilter top depth breakthrough between the warm and cold well area as this would mean that the outlet temperature of the pumped water will change with the inlet temperature of the injected water as the water mixes (Gao et al. 2013).Quantitatively, the thermal recovery efficiency η th was calculated for each cycle of 6 months (i.e. one season) as the percentage of thermal energy that can be extracted from the energy that was stored in the previous cycle (Duijff et al. 2021): where E ex and E in (kWh) are the extracted and injected energy, Q ex and Q in (m 3 /h) the total extraction and injec- tion flow rate of the system, c w the specific heat capacity of water (1.16 kWh/m 3 K), ΔT (°C) is the absolute tempera- ture difference between the injected/extracted water and the background temperature (13.8 °C) of the aquifer, and t (h) is time.

Maximum aquifer yield
The maximum pumping rate in the currently available wells Yd 4 and Yd 6 was estimated to be, respectively, 4 and 1 m 3 /h.According to Lebbe et al. (1992), a pumping rate of 1.66 m 3 /h could be reached in PP 2, while this study's equipment was limited to a rate of 1 m 3 /h.
During the pumping/injection test, the drop in water level due to pumping in Yd 4 could be observed in the pumping layer itself but also in Yd 6 (Fig. 5).In Yd 2, no drop in water level could be observed.This drop in water level illustrates that the semipervious layer Yd 5 does not prohibit the connection between Yd 6 and Yd 4 so the two layers cannot be considered as separate units.In contrast, the fact that no drop in water level could be observed in Yd 2 illustrates the confining nature of Yd 3.This is confirmed by the observed influence of the injection.Despite the fact that PB JE is filtered from the top of Yd 6 to the bottom of Yd 1 (Fig. 3), the influence of the injection was only clearly visible in Yd 2 (Fig. 5).
Considering the absence of a pumping well with filter screens in all three pervious layers, the maximum pumping rate in the layered aquifer had to be estimated by means of the principle of superposition.This principle assumes that each pervious layer is fully confined (thus isolated from each other).On the one hand, this assumption seems acceptable for Yd 2 and Yd 4. On the other hand, it was shown that a strong connection between Yd 4 and Yd 6 exists indicating that Yd 5 is not a good confining layer.Hence, the maximum pumping rate in Yd 6 and Yd 4 combined is most likely smaller than the sum of their individual rates.Assuming that (7) PP 6 can only account for 10% of its estimated maximum pumping rate (i.e.0.1 m 3 /h), a total maximum pumping rate of 5.76 m 3 /h was estimated in a fully screened well in the layered aquifer.It is however not recommended to exploit a well at the maximum rate, especially for a long period of time such as in an ATES system; therefore, the rate will be limited to 5 m 3 /h for the rest of the study.
The maximum injection rate in the fully penetrating well PB JE was estimated to be 3.8 m 3 /h with a water level reaching the surface but not flooding it (head increase of about 3 m).This indicates that the maximum allowable head change of 1.8 m (Eq.6) is quite strict as in practice no problem was encountered on the field with a change of 3 m.The injection well PB JE was the limiting factor during the field test as it was constructed as a piezometer with a small diameter, not as a pumping or injection well.A higher rate can probably be reached when using injection wells with a larger diameter and adapted filters; thus, most likely, the estimated maximum rate of 5 m 3 /h can also be used for reinjection.
The injection capacity was also determined by the MFI tests.No sand could be visibly observed in the mesh netting and the slope for the two tests was 5.2 and 8.3, respectively.The filters of both MFI tests were visibly still relatively clean, which is a good indicator that the low MFI values measured are reliable.During the long-term injection, no decrease of injection capacity with time was observed, confirming that the risk of clogging should be limited.
The estimated maximum flow rate of 5 m 3 /h could not be confirmed in the field because of the absence of a fully screened well with a large diameter and adapted filter.The impact of injection at this flow rate was therefore evaluated with the groundwater model (section 'Model validation').Operating an ATES system at a flow rate of 5 m 3 /h per well pair would be significantly below the usual pumping rate per well which is used in operating ATES systems (Table 3) and therefore challenges in terms of hydraulic efficiency are expected.

Model validation
The new model is based on an axisymmetric model made earlier by Lebbe et al. (1992) which put great effort into calibration with a triple pumping test.To ensure that the new model is representative of the layered aquifer, it was validated by first simulating the same triple pumping test of Lebbe et al. (1992) and comparing the results.A good agreement with the triple pumping test data was observed, indicating that the new numerical model is a good proxy for the model of Lebbe et al. (1992)-for detailed results see Fig. S1 in ESM.
Second, using the same model parameters, the simulated drawdowns (relative to the natural groundwater level) were compared to the drawdowns observed during the pumping and injection tests that were carried out for this project (Fig. 6).In general, when the injection is implemented, it becomes more difficult to simulate the observations.There is a good agreement in Fig. 6 for the observation wells (PB x.y); however, there remains a discrepancy for the pumping/injection wells (PP x and PB JE).When considering a logarithmic time scale, the observed and the simulated drawdowns are however relatively parallel, indicating they are characteristic of a reservoir with similar properties (Cooper and Jacob 1953).
The simulated pressure at the bottom of the well is too high compared to the observed one.This deviation can likely be attributed to the imperfect sealing of the semipervious layers during completion (Lebbe et al. 1992).The discrepancy at the injection well is larger, and this is attributed to the difficulty in modeling well behaviour.The injection well was explicitly represented as a fully filtered well in the layered aquifer in the model used for validation.This was simulated by setting the hydraulic conductivity very high at this location from Yd 1 to Yd 6 and setting the injection in Yd 4, which is different from the more straightforward representation of the pumping well that is only filtered in layer Yd 4.
Modeling a pumping/injection well is complex as the detailed set-up of the well (bentonite seal, gravel pack, inner tube) is almost impossible to implement explicitly in MOD-FLOW.The size of the discretization grid was reduced to a size close to the well diameter, approximating the dimensions of the well (Klepikova et al. 2016).The presence of water in the well instead of sediment was simulated by setting a very high value of the hydraulic conductivity (vertical and horizontal), decreasing the water pressure at the well when simulating injection.Nevertheless, the positive skin factor, which is a reduction of the permeability in the immediate vicinity of the well due to drilling/well completion/production processes, was not implemented in the model (Van Everdingen 1953).In general, by reducing the permeability, the drawdown in the well itself would be larger and the cone of depression would be steeper.As such, at the pumping well PP 4, a positive skin factor could be present, which causes the observed drawdown to be larger than the simulated one as observed in Fig. 6.At the injection well PB JE, a positive skin factor would cause the negative drawdown to be larger (in absolute values) than the simulated one; however, the opposite is observed (Fig. 6).Although a negative skin factor could be present, as a gravel pack is present over the whole thickness of the aquifer including around semipervious layers, it is believed that the discrepancy is related to an inadequate representation of the injection well in the model, which is corroborated by the higher discrepancy compared with PP 4. Since the model was initially calibrated by the pumping test as a homogenous aquifer, heterogeneity could also explain the difference: a higher hydraulic conductivity in the vicinity of PB JE would reduce the simulated pressure increase.However, including lateral heterogeneity at this stage would be highly speculative as the test zone is not representative of the whole modelled area.
The validation results show that field tests in this kind of circumstance are an added value.By comparing the field test results with the simulated values, a discrepancy between the model and reality at the location of the injection/pumping wells can be identified.This deviation is not present near the observation wells; thus, the reason likely lies in an inadequate representation of the pumping/injection wells.Therefore, the model was considered valid for simulating the ATES scenarios, noting that the simulated pressure is likely overestimated in both the pumping and injection wells.This should be taken into account when interpreting the model results.The results will be interpreted at a distance of 1 m from the centre of the wells.
To quantify the head changes that are expected at a flow rate of 5 m 3 /h, a short simulation was done using the groundwater model (Fig. 7).Combining these findings with the validation results shows that the head change at 1 m distance from the well simulated by the model is 3.45 m, which is likely an overestimation, as discussed earlier.According to Eq. ( 6), the maximum allowable head change in the wells would be 1.8 m (NVOE 2006).However, the long-term injection test has shown that such an increase in hydraulic head could be sustained by the well.Therefore, this head change is considered realistic when simulating multiple wells for the ATES system; however, this should be confirmed with a properly designed well at the rate of 5 m 3 /h.

Hydraulic efficiency
Considering the power requirement of 0.63 MW, ignoring in the first instance the coefficient of performance, and a standard temperature difference between the extracted and the injected water of 5 °C, the required total pumping rate was estimated to be 108.62 m 3 /h (Eq.4).Assuming a pumping rate of 5 m 3 /h per well pair, this results in the need for 22 well pairs and hence 44 wells.Next, the hydraulic and thermal radius of influence were estimated in unit Yd 4 as it is the most permeable.This was calculated analytically taking into account a drawdown of 1.56 m at a distance of 5 m from the well as was indicated by a simulation of 6 months, resulting in a radius of influence of 36 m (Eq.5).Accounting for the thermal retardation factor of 1.83 (Eq.3), the thermal radius of influence in Yd 4 is ~20 m, an estimation that was also confirmed by the model.Considering the guidelines for well placement drawn up by Bloemendal et al. (2018), first, a well arrangement in two PB JE PP 4 PB 6.1 Fig. 6 Comparison of the simulated to the observed drawdown at wells PP 4, PB JE and PB 6.1.The drawdown is positive when the water level decreases and negative when the water level increases.On the right side, this is plotted using a logarithmic time clusters (i.e. one group of cold wells and one group of warm wells) was simulated (Fig. 8), allowing the space occupied by the ATES system, using the existing buildings as a constraint to be minimized.However, this configuration yields a total absolute increase in water level in the injection cluster of ~53 m, while the drawdown resulting from a cluster of pumping wells is 46 m.Compared to the maximum allowable head change, those values are too high and not sustainable in practice.They result from the interactions between wells of the same type which are each operating close to the maximum flow rate of the aquifer.Other clustered well arrangements were tested, but the simulated (negative and positive) drawdown remained too high (Fig. 8); hence, increasing the risk of soil outbursts or flooding.This indicates that the hydraulic conditions, rather than the thermal conditions, constitute the main limiting factor in this case.
To limit these excessive head changes, a well configuration in lanes (alternating injection and pumping wells) was tested.This configuration should effectively limit the head changes by the principle of superposition as the lanes of pumping wells can counterbalance the influence that the lanes of injection wells have on the groundwater level.Taking into account the thermal radius of influence of 20 m and the guidelines for placement drawn up by Bloemendal et al. (2018), the distance between the lanes should be at least 50 m.It was shown that a thermal breakthrough between the warm and cold well area was established within the first 6 months of operation, implying that the distance between the lanes was too small.Consequently, the distance between the lanes was incrementally adjusted up to 90 m to avoid a thermal breakthrough.This configuration significantly reduced the drawdown and the injection pressure, but the maximum head change of 15 m was still too high (Fig. 9) because wells of the same nature within a lane are too close to each other, and therefore interact.
Since the hydraulic conditions are the limiting factor in such a low-transmissivity aquifer, a configuration where the wells were placed in a checkerboard pattern was tested to further limit the head changes.In this pattern, all injection and extraction wells are alternating, constituting the best option to limit the change in the hydraulic head in the aquifer, but it could reduce the thermal recovery efficiency.The wells were placed as far as possible from each other, within the available space of Campus Sterre, resulting in a well spacing of minimum 80 m, i.e. 2.2× the estimated hydraulic radius of influence.With this configuration (Fig. 10), the maximum hydraulic head while injecting remained limited (12.3 mTAW, or 2.3 m above the ground surface) as was the case for the single well in Fig. 7.This maximum head change of ~5 m is higher than is advised in the guidelines of NVOE ( 2006), but taking into account the results discussed

Thermal recovery efficiency
Both the checkerboard-type and lane-type arrangements were tested for 20 years to determine the thermal recovery efficiency of the ATES system.For the lane-type arrangement, the maximum temperature of the cold well area decreases with time to roughly 10.5 °C, while the minimum temperature of the warm well area increases to roughly 17 °C (Fig. 11).The minimum temperature difference between the warm well area and the cold well area also increases with time, from ~8.5-9.5 °C, indicating that the storage is efficient and no thermal breakthrough occurs.
For the latter, a difference between the seasons of roughly 0.5 °C can be observed.There is also a limited difference in temperature between the observation locations in the NE and SW, related to the inhomogeneity of the layer thickness.
Figure 12 also shows that for the checkerboard-type arrangement, the minimum temperature of the warm well area increases while the maximum temperature of the cold well area decreases.The system is therefore thermally efficient, although the temperature difference is, as expected, slightly lower than in the previous scenario.This is also indicated by the thermal recovery efficiency η th which increases with time and becomes close to 1 for both systems but is slightly lower for the checkerboard-type arrangement (Fig. 13).As the maximum temperature difference is roughly twice as large as the initially estimated 5 °C, the maximum produced power was also estimated to be about twice the power requirement of this project (2 × 0.63 MW).Next, the extracted thermal energy per season (integration over 6 months for heating or cooling) was calculated based on the temperature difference in Figs.11 and 12, and the extraction rate (Fig. 13).The initially estimated power demand for this project was based on a cumulative energy demand of 1.5 GWh with a peak cooling demand in summer.Figure 13 shows that the simulated energy output is 3× as large.However, as long winter and summer seasons (6 months) are simulated in which the system operates at its maximum flow rate and thus power, it must be emphasized that the simulated scenario represents the maximum upper limit for energy production.

Discussion
The field tests and initial modeling results have validated that the low-transmissivity aquifer located on Campus Sterre could only sustain a maximum pumping and injection rate of about 5 m 3 /h.This limiting extraction rate results in new constraints for the arrangement of the ATES wells.
Arranging the wells in two large clusters as was first done in Fig. 8 is not a feasible option for the future ATES project considering the excessive drawdowns.This results from the principle of superposition which implies that the resulting hydraulic head at a certain location in a confined aquifer is the combination of all influences (resulting from different pumping/injection wells) at that same location.When many wells are grouped into one cluster with an inter-well distance that is less than the hydraulic radius of influence, this will result in an excessive drawdown (positive or negative), especially in a configuration where the extraction rate is close to the maximum yield.Even though the model tends to overestimate the injection pressure, an increase in water level of minimum 31 m relative to the natural groundwater level is not acceptable.Using such a configuration, flooding or well collapse will likely occur at the injection wells and the relatively thin confining clay layer might not be able to withstand such high pressures.Also, the minimum drawdown of ~28 m at the pumping wells is not feasible.This would probably cause the aquifer to become partly unsaturated and consequently the necessary pumping rate of 5 m 3 /h per well could not be reached anymore.
Next, the arrangement of wells in several lanes showed that a distance of 90 m between the lanes was necessary to avoid a thermal breakthrough.This is significantly larger than the guidelines proposed by Bloemendal et al. (2018).Because of the superposition effect of multiple wells, the radii of influence which were deduced for a single well pair are not valid anymore for clusters or lanes.These configurations result in an increased gradient between warm and cold wells which speeds up groundwater flow.The lane arrangement showed an overall improvement in the calculated drawdown at the location of the injection wells due to the closeness of the pumping wells.However, it is also shown that the hydraulic head in the SW of the well area is significantly higher than in the NE, which might be the result of the fact that the warm lanes in the SW are not located in between two cold lanes, which counterbalances the increase in hydraulic head.
From a thermal point of view, no thermal breakthrough occurred if the distance between the wells was large enough.In this case, the safe distance was 4.5× larger than the thermal radius calculated for one well pair.As expected, the thermal efficiency increases with time as illustrated by the minimum temperature difference between warm and cold wells.A seasonal variation of about 0.5 °C remains within realistic limits.The small difference between the temperature at the observation locations in the NE and SW might be explained by the fact that in the NE the lanes consist of more wells close to each other, slightly decreasing the thermal recovery efficiency.The thermal recovery efficiency after 20 years of ATES simulation showed to be close to 1 and the simulated energy output is almost 3× as large as the initially estimated cumulative energy demand for the project.This relates to the temperature difference which increases with time.As the well number was calculated based on this temperature difference, the number of wells could be decreased based on the results.However, this would also impact the coefficient of performance of the heat pump (more energy required to reach this temperature).A better approach would be to decrease the flow rate per well pair as it is close to the critical flow rate of the aquifer and the hydraulic circumstances are the limiting factor in the design.Important to emphasize is that the calculations represent an upper limit for energy production and that the efficiency of the heat exchange system is not considered in this study.In practice, the efficiency of the system is often ~50%, which is generally included in the calculation of the coefficient of performance which was, however, not done for this project.
These investigations demonstrate that the injection pressure is the main limiting factor in the design of this ATES system and in general this should be taken into account when looking to optimize well placement in lowtransmissivity aquifers.It is not surprising, as the pumping/injection rate (5 m 3 /h) is close (86%) to the estimated maximum rate for this aquifer (about 5.8 m 3 /h).Consequently, a well arrangement in a checkerboard pattern is the best approach to limit the increase in hydraulic head in injection wells.However, this pattern utilizes all the available space and is hence less cost-efficient when considering the piping and installation costs.In this specific case, the heat network connecting all the buildings to the servers would anyways require a vast number of piped connections, including those for a BTES system.A simplified economic analysis showed that the total cost of installation of an ATES system would be about 50% the cost of the foreseen BTES system.So even though more piped connections (€ 56/m) are required for the checkerboardtype arrangement, the ATES system still remains far more cost efficient.In terms of thermal recovery efficiency, the checkerboard-type arrangement might also increase the risk of a thermal breakthrough as cold and warm wells are alternating; however, as the injection pressure is more limited, it reduces the effective thermal radius.Such a configuration can easily be adapted to accommodate the buildings and optimize the usage of the available space.Because of the limited hydraulic interaction, the distance between the wells might be decreased to ~60 m without significantly increasing the risk of thermal breakthrough.
The proposed design of the well placement strongly depends on the estimation of the increase in hydraulic head related to the injection at an estimated rate of 5 m 3 /h by a numerical model.In this case, the simulated increase in head was overestimated by the model, but this observation is likely dependent on local conditions, including well completion.Since the design, and therefore the risk of failure of such an ATES system is strongly dependent on the results of the numerical model, the authors suggest that more field investigations should be performed before envisaging ATES systems in low-transmissivity aquifer.To limit investigation costs, a two-phase experimental plan is likely the most efficient.First, the maximum pumping and injection rate should be estimated, together with the risk of clogging.This can be best done using a pair of wells and performing a step-wise pumping test, and a long-term (at least a few months) pumping and injection test.In that way, the viability of a single pair of wells can be validated.Based on the energy demand, the number of pairs should then be estimated, and the possibility to install them based on an estimated hydraulic radius of influence should be checked.It is recommended to consider wells at a distance of at least 2× the hydraulic radius of influence.At this stage, a groundwater flow model can also be developed to verify the efficiency of such a system.Second, it is recommended to install a second pair of wells, to verify in field conditions the interactions between the wells.This will ensure that the estimated maximum rates can be sustained when several wells are operated together and would reduce the dependency on model results for the design of the ATES system.
Although this study demonstrated an ATES system should be sustainable on Campus Sterre, there are still some challenges and uncertainties which need to be further investigated: 1.The pumping and injection rate of 5 m 3 /h should be confirmed using a newly constructed well pair with optimally placed filters in the three pervious layers.These wells would have a large diameter and should be thoroughly developed to reduce well losses.2. To improve the model and optimize the efficiency of the system, heat losses at the surface should be introduced in the model by adapting the boundary conditions for heat transport.3. A detailed sensitivity analysis should be carried out for the porosity and thermal parameters.Heat tracer experiments could also be carried out to validate the porosity and thermal parameters (e.g.Wildemeersch et al. 2014).
Next to this, a thermal response test can be carried out to validate the thermal conductivity of the subsurface.4. Homogeneous parameters within each layer were used for the groundwater models.This is likely an oversimplification and should be investigated.Heterogeneity is known to influence the overall efficiency of ATES systems (Possemiers et al. 2015;Sommer et al. 2013;Hermans et al. 2018Hermans et al. , 2019)).Similarly, the thickness and continuity of the confining layer throughout the campus should be confirmed.If it were partly absent, the studied aquifer layer would constitute an unconfined aquifer together with the Quaternary layer, which would likely modify the hydraulic behaviour and the conclusion of the study.It can be estimated through cone penetration tests or other geophysical well logging methods such as gamma-ray or electromagnetic methods.

Conclusion
A medium-permeability aquifer with limited thickness located on Campus Sterre (Ghent, Belgium) was investigated as a possible candidate for aquifer thermal energy storage even though it is conventionally disregarded because of its limited transmissivity.Based on field experiments, it was shown that the maximum pumping and injection rate within this aquifer is only 5 m 3 /h, which is much smaller than most operating ATES systems.Nevertheless, the aquifer seems suitable as the injection rate could be sustained for a long period of time, and no clogging of the injection well could be detected, which is one of the main concerns in a medium-permeability aquifer.
For low-transmissivity aquifers, the energy demand must be covered by increasing the number of well pairs, which results in new challenges for well placement.Based on the energy demand on the campus, an ATES system using 22 well pairs should be operated.Simulation with calibrated groundwater models showed that well arrangements in clusters or in lanes, based on estimation of the hydraulic and thermal radii from a single well pair, were not adequate because they resulted in an excessive water pressure and drawdown in the injection and pumping zones, respectively.This is a direct consequence of the superposition principle, as neighbouring wells interact with each other.Instead, the wells should be arranged in a checkerboard pattern of alternately warm and cold wells with each a pumping/injection rate of 5m 3 /h, with a distance of at least 2× the estimated hydraulic radius.This shows that the hydraulic efficiency rather than the thermal recovery efficiency is the limiting factor in the design.Although such an ATES system would operate in suboptimal conditions, ATES is certainly an option to consider in comparison with BTES in a detailed design and cost analysis.
This study shows that, so far, the potential of low-transmissivity aquifer for ATES systems has likely been underestimated.With the increase in energy prices and the longterm objectives to reduce greenhouse gas emissions, the interest in ATES systems will likely increase in the future.In the absence of accessible productive aquifers, either because of their absence or because they are used for drinking water production, low-transmissivity aquifers can constitute suitable alternatives, although suboptimal and strongly dependent on local conditions.Their potential should be confirmed by more field studies targeting specifically ATES systems (long-term injection, clogging, heterogeneity).It is recommended to investigate these aquifers in two phases.In the first phase, the maximum injection and pumping rates can be estimated using one pair of wells to check the technical and economic feasibility using a numerical model.In a second phase, the interaction between pairs should be investigated with the installation of a second pair of wells to confirm that it will not hinder the efficiency of the ATES system.

Fig. 1
Fig. 1 Graphical representation of a an aquifer thermal energy storage (ATES) and b borehole thermal energy storage (BTES) system in winter (left) and summer (right) seasons (after Bloemendal 2018)

Fig. 2 a
Fig. 2 a Location of the study area in Belgium and on the suitability map for ATES systems in Flanders (after WTCB 2017).b Outline of the study area, Campus Sterre

Fig. 3
Fig. 3 Hydrostratigraphy across Campus Sterre with indication of the filter placement of the wells used for the pumping tests (after Vlaamse Overheid n.d.; Lebbe et al. 1992).mTAW = groundwater level in relation to sea level

Fig. 4
Fig. 4 Location of the well area used for the pumping tests on Campus Sterre

Fig. 5
Fig. 5 Results of the pumping tests carried out on Campus Sterre: a sites PP4 and PP6, and b sites PB JE and PP2 (the peaks can be attributed to the disruption of the measurements to carry out membrane filtering index (MFI) tests).Dates given as yyyy-mm-dd

Fig. 8 Fig. 9 Fig. 7
Fig. 8 Hydraulic head distribution after 6 months for different well arrangements, arranged by a cluster of 20-40 m, and b cluster of 40 m.The dotted line indicates Campus Sterre

Fig. 10
Fig. 10 Hydraulic head distribution after a simulation period of 20 years in a checkerboard pattern

Fig. 12 aFig. 13 a
Fig. 12 a Observation locations for the checkerboard arrangement, and (below) the temperature at the warm well area (red), the cold well area (blue), and the temperature difference (grey) after 20 years at the observation locations in the NE part of the well area b and the SW of the well area (c)

Table 2
Vandenbohede et al. (2011) used in the heat transport equation.The values are adopted fromVandenbohede et al. (2011)