# Assessment of short-term aquifer thermal energy storage for demand-side management perspectives: Experimental and numerical developments

Guillaume de Schepper, Claire Paulus, Pierre-Yves Bolly, Thomas M Hermans, N. Lesparre, Tanguy Robert

<figure><img src="image_1.png" /></figure>

## ► To cite this version:

Guillaume de Schepper, Claire Paulus, Pierre-Yves Bolly, Thomas M Hermans, N. Lesparre, et al.. Assessment of short-term aquifer thermal energy storage for demand-side management perspectives: Experimental and numerical developments. Applied Energy, 2019, 242, pp.534-546. 10.1016/j.apenergy.2019.03.103. hal-02444112

HAL Id: hal-02444112

https://hal.science/hal-02444112v1

Submitted on 27 Sep 2025

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Assessment of short-term aquifer thermal energy storage
for demand-side management perspectives: Experimental
and numerical developments

Guillaume De Schepper¹*, Claire Paulus²,³, Pierre-Yves Bolly¹,³, Thomas Hermans⁴,
Nolwenn Lesparre⁵, Tanguy Robert⁶,⁷

* Corresponding author: g.deschepper@aquale.com

1. AQUALE SPRL, R&D Department, Rue Ernest Montellier 22, 5380 Noville-les-Bois, Belgium.

2. Raco BVBA, Meylandtlaan 39, 3550 Heusden-Zolder, Belgium.

3. UC Louvain, Louvain School of Engineering, Rue Archimède 1, 1348 Louvain-la-Neuve, Belgium.

4. Ghent University, Department of Geology, Campus Sterre, Krijgslaan 281, 9000 Gent, Belgium.

5. University of Strasbourg, Laboratoire d'hydrologie et de géochimie de Strasbourg (LHyGeS), CNRS UMR 7517, 1 Rue Blessig, 67084 Strasbourg, France.

6. University of Liège, Urban & Environmental Engineering, Quartier Polytech 1, Allée de la découverte 9, B52 Building (Sart-Tilman), 4000 Liège, Belgium.

7. F.R.S. - FNRS (Fonds de la Recherche Scientifique), Rue d'Égmont 5, 1000 Brussels, Belgium.

The final version of this manuscript has been published in *Applied Energy*, **242**, 534-546.
https://doi.org/10.1016/j.apenergy.2019.03.103
## Highlights

* Shallow alluvial aquifers are suitable to perform short-term thermal energy storage.

* It has a high development potential for demand-side management applications.

* Energy recovery rates are high for typical demand-side management frequencies.

* Preheating shallow alluvial aquifers for demand-side management is feasible.

## Key words

Aquifer thermal energy storage; Demand-side management; Energy efficiency; Numerical modelling; Geothermics; Geophysics

## Abstract

In the context of demand-side management and geothermal energy production, our proposal is to store thermal energy in shallow alluvial aquifers at shorter frequencies than classical seasonal aquifer thermal energy storage. We first conducted a one-week experiment in a shallow alluvial aquifer, which is characterized by a slow ambient groundwater flow, to assess its potential for thermal energy storage and recovery. This experiment has shown that up to 90 % of the stored thermal energy can be recovered and would therefore suggest that aquifer thermal energy storage could be considered for demand-side management applications. We then conceptualized, developed, and calibrated a deterministic 3D groundwater flow and heat transport numerical model representing our study site, and we simulated 77 different scenarios to further assess this potential. This has allowed us to demonstrate that low-temperature aquifer thermal energy storage (temperature differences of -4 K for precooling and 3, 6, and 11 K for preheating) is efficient with energy recovery rates ranging from 78 to 87 %, in a single aquifer thermal energy storage cycle. High-temperature aquifer thermal energy storage (temperature differences between 35 and 65 K) presents lower energy recovery rates, from 53 to 71 %, with
all other parameters remaining equals. Energy recovery rates decrease with increasing storage duration and this decrease is faster for higher temperatures. Retrieving directly useful heat (without upgrading with a groundwater heat pump) using only a single storage and recovery cycle appears to be complicated. Nevertheless, there is room for aquifer thermal energy storage optimization in space and time with regard to improving both the energy recovery rates and the recovered absolute temperatures.
# Nomenclature

## Abbreviations

<table><tr><td>ATES</td><td>Aquifer Thermal Energy Storage</td></tr><tr><td>DHW</td><td>Domestic Hot Water</td></tr><tr><td>DSM</td><td>Demand-Side Management</td></tr><tr><td>ERT</td><td>Electrical Resistivity Tomography</td></tr><tr><td>GWHP</td><td>Ground Water Heat Pump</td></tr><tr><td>HT-ATES</td><td>High-Temperature Aquifer Thermal Energy Storage</td></tr><tr><td>LT-ATES</td><td>Low-Temperature Aquifer Thermal Energy Storage</td></tr></table>

## Parameters

<table><tr><td>c<sub>s</sub></td><td>Volumetric heat capacity of the porous medium [J m<sup>-3</sup> K<sup>-1</sup>]</td></tr><tr><td>K<sub>XY</sub></td><td>Horizontal hydraulic conductivity [m s<sup>-1</sup>]</td></tr><tr><td>K<sub>Z</sub></td><td>Vertical hydraulic conductivity [m s<sup>-1</sup>]</td></tr><tr><td>n</td><td>Pore-size distribution of the van Genuchten model [-]</td></tr><tr><td>n<sub>e</sub></td><td>Effective porosity of the porous medium [-]</td></tr><tr><td>n<sub>t</sub></td><td>Total porosity of the porous medium [-]</td></tr><tr><td>S<sub>r</sub></td><td>Residual saturation [-]</td></tr><tr><td>α</td><td>Inverse of the air entry suction of the van Genuchten model [m<sup>-1</sup>]</td></tr><tr><td>α<sub>L</sub></td><td>Longitudinal dispersivity [m]</td></tr><tr><td>α<sub>T</sub></td><td>Transverse dispersivity [m]</td></tr><tr><td>δ</td><td>Error criterion value [-]</td></tr><tr><td>λ<sub>s</sub></td><td>Thermal conductivity of the porous medium [J m<sup>-1</sup> s<sup>-1</sup> K<sup>-1</sup>]</td></tr></table>

## Variables

<table><tr><td>i</td><td>Hydraulic gradient [-]</td></tr><tr><td>Δt</td><td>Aquifer storage duration [h]</td></tr><tr><td>ΔT</td><td>Temperature difference between the initial and stored water [K]</td></tr><tr><td>T</td><td>Absolute temperature of groundwater [°C]</td></tr><tr><td>η</td><td>Energy recovery rate [-]</td></tr></table>
# 1 Introduction

Interest in the diversification of energy sources has become a driving force for energy-transition political decision-making. Recent researches in the Netherlands [1], India [2], and Italy [3] demonstrated that the aging of the population and their increasing wealth clearly offset improvements in the energy efficiency of the building stock. As a consequence, countries need to promote renewable and sustainable energies in addition to the latter, as the European Union did through the Energy Efficiency Directive 2012/27/EU [4]. In this context, demand-side management (DSM) has played a crucial role in the development and integration of new (often renewable, decentralized, and intermittent) energy production techniques [5].

With the increasing use of electrically-driven heating, ventilation, and air-conditioning systems in buildings, smaller consumers are seen as key players for DSM in the low-voltage grid [6]. Space heating and cooling as well as domestic hot water (DHW) production with heat pumps in particular can be combined with various thermal energy storage strategies for flexibility purposes [7]. It is possible to take advantage of the ability of those systems to decouple electricity consumption and heat/cold demand in buildings, by partly or totally rescheduling the electrical consumption (thermostatically-controlled load-shifting) from peak hours to off-peak hours. Office and residential buildings (and by extension, the tertiary sector) can therefore offer a potential for flexibility by interacting with the electrical grid to improve the system-level efficiency [8] or the system reliability [9]. A detailed review of interactions between heat pumps and smart grids can be found in Fischer and Madani [7].

In this context, and provided that a productive aquifer is present, groundwater heat pumps (GWHP) offer a significant potential for flexibility when we consider their thermal (hundreds to thousands of kW) and electrical power (hundreds of kW), together with the thermal inertia of buildings [10]. Until now, thermostatically-controlled load-shifting has been achieved using the thermal envelope of the building or using water tanks [11]. Whereas water tanks contain
finite volumes of water, a GWHP system is connected to an aquifer which possesses an almost
infinite (and insulated) storage volume inside the porous media [12]. With this study, we not
only propose to further consider strategies centered on aquifers for thermal energy storage to
provide a new potential flexibility tool but also to improve the overall energy efficiency of
GWHP systems.

In a GWHP system, groundwater is pumped from a production well and delivered to a heat pump (or heat exchanger). Groundwater acts as a heat source for space heating and DHW production, while it is a chilled source for space (free) cooling. Often, groundwater is directly reinjected in the aquifer (with a lower or higher temperature depending on whether the system is used for heating or cooling) using an injection well. Such systems are called open-loop geothermal systems. Most of those systems are designed with well doublets, namely an injection well distinct from the production well [13], working in cyclic or continuous mode [14].

Aquifer Thermal Energy Storage (ATES) systems further also use GWHP to store thermal
energy in the subsurface with heat being stored during space cooling and cold being stored
during space heating for example [15]. The interested reader is referred to Hesaraki et al. [16]
for a comprehensive review of subsurface thermal energy storage, to Bayer et al. [17] for a
review of the environmental impact of geothermal energy production, and to Haehnlein et al.
for its international legal status [18] and its sustainability [19]. At present, ATES are mainly
used on a seasonal timescale [16] due to the basic functioning of GWHP. Moreover, ATES is
neither optimized in space, nor in time [20]. Some authors have considered other types of
thermal energy source for ATES, including waste heat [21], solar heat [22], and power-to-heat
converted electricity [23] to improve the efficiency of such systems. In this study, we do not
focus on the energy source itself but on the ability of the subsurface to efficiently store heat and
cold.
The efficiency of thermal energy recovery depends on various processes, namely thermal conduction and dispersion, regional groundwater flow, or even density-driven flow [24]. Regional groundwater flow is also related to the regional hydraulic gradient (driving force) and the hydraulic conductivities of the area of interest [25]. Bloemendal and Olsthoorn [26] demonstrated for example that recovering stored thermal energy in high ambient groundwater flow velocity aquifers is manageable when the recovery wells are located over the heat transport path. The heterogeneity of the porous media, in terms of hydraulic conductivity distribution, significantly impacts thermal energy storage [27]. Short-term storage appears to provide better results in high than in low permeability media, with higher recovery rates, while long-term storage has poorer recovery rates in high permeability media compared to low ones [28]. In addition, Ferguson [29] also stated that, in a model, the heterogeneity of the hydraulic parameters has a stronger effect than the heterogeneity of the thermal parameters.

The main idea behind our work is to further consider short-term ATES for their potential for flexibility, more specifically with the implementation of two different strategies. The first one consists in slightly preheating the aquifer ($T < 30°C$) during off-peak periods and recovering the stored thermal energy during peak periods. Note that precooling the aquifer for space cooling also fits in this tactic. This strategy is called low-temperature ATES (LT-ATES) and takes profit of improving the coefficient of performance (COP) of the GWHP for space heating and of the heat exchanger for space cooling. The second strategy consists in storing thermal energy with higher temperatures in order to retrieve heat ($T > 50°C$) that can be directly used for space heating without the need for upgrading (no GWHPs needed) [12]. This strategy is called high-temperature ATES (HT-ATES) and enables the use of various higher temperature heat sources which are present near cities. Among them, we can cite fatal heat from industry, waste heat (from incineration or inside landfills), solar heat, power-to-heat, and, to a certain extent, sewers.
In the following sections, we first introduce our study site, which is located in Belgium, and
justify our choice of this type of aquifer (section 2). We then present our ATES experiment run
over an entire injection, storage and recovery cycle (section 3). We also show how to
conceptualize, build and calibrate a deterministic predictive numerical model for groundwater
flow coupled to heat transport using the data gathered during this experiment (section 4). We
then use this predictive model to forecast 77 different short-term ATES scenarios (consisting
of single injection, storage and recovery cycles) related to typical DSM frequencies (real-time,
intraday, and interday). These scenarios are grouped to reflect the three above-mentioned
categories: (1) precooling of the aquifer for space cooling, (2) preheating of the aquifer (LT-
ATES) for space heating and DHW production with the help of a GWHP, and (3) direct use of
the recovered heat (HT-ATES) for space heating (section 5). We discuss the energy recovery
rates for all scenarios and more particularly discuss the recovered temperatures for the third
group to meet a no-upgrading objective (section 6). At the end of this paper, we draw our
conclusions and propose perspectives to improve energy recovery rates for both LT- and HT-
ATES, and obtain higher recovered temperatures in HT-ATES cases (section 7).
## 2 Study site

In Belgium, most of the major cities are built on shallow alluvial aquifers [30] composed of sand and gravel. Such aquifers are characterized by their high productivity, even though some clay lenses may be locally responsible for the low productivity of certain areas [31]. These aquifers are seen as ideal targets for low-enthalpy geothermal systems, as pointed out by Allen and Milenic [32]. These authors demonstrated that a pumping rate of 72 m³/h (reachable with one or two production wells in such alluvial aquifers) and a temperature reduction of 8 K in the GWHP can generate a heating power of 672 kW.

For these reasons, we conducted our study in the alluvial plain of the Sambre River (Error! Reference source not found.a). The topography of the study site is roughly flat. The area of interest is located in sparse woods, which extend from east to west and is bordered by a dirt track (north) and a grass field (south) (Figure 1b). The local geological setting is typical of an alluvial plain: a 2 to 3 m thick layer of Quaternary clayey loam soils, partially mixed with backfill soils, is underlain by ~5 m thick Quaternary sandy gravel alluvial deposits. A clay layer, several meters thick and originating from the weathering of the subjacent Carboniferous coal shales from the regional scale Hainaut coalfield [33] is found deeper (Error! Reference source not found.).
<figure><img src="image_3.png" /><figcaption>**Figure 1.** Location of the study site (a) in Belgium (Wallonia), and location of the largest cities built on alluvial deposits of the Sambre and Meuse Rivers. The experimental setup (b) shows the implementation of the pumping and injection wells (W1 and W2), together with observation piezometers PzA, PzB and PzC. The 4D ERT setup is composed of 6 ERT profile lines and 189 electrodes in total.</figcaption></figure>
From a hydrogeological perspective, the shallow clayey loam layer confines the subjacent
fully saturated alluvial aquifer. Hydraulic conductivities in the aquifer were estimated by
pumping tests to be 10⁻³ m/s on average [10]. The deeper clay layer is almost impervious and
is considered as the bottom of the targeted aquifer. The local hydraulic gradient (driving force
of the groundwater flow) in the aquifer is extremely low (i = 0.0001 in average) with a very
slow groundwater flow from north to south. Groundwater temperature in the vicinity of our
experiment is 10 °C [10].
### 3 Aquifer thermal energy storage experiment

In May 2015, an ATES experiment designed to test the efficiency of a well doublet was conducted during 4 days in the study area, focusing on short-term heat storage in the alluvial aquifer (Figure 2). Groundwater was pumped from the pumping well at a rate of 2.55 m³/h (W1, Figure 1b), at a temperature of 10°C (aquifer ambient temperature). The water was heated at 40 °C (ΔT = 30 K) and injected in the aquifer through the injection well (W2, Figure 1b) at 2.55 m³/h during 5 h 20 min and stored for 72 h.

Following the storage phase, the stored energy was recovered by pumping water from the injection well during 4 h 30 min, at an average pumping rate of 7.6 m³/h. The very low local hydraulic gradient allowed the energy to be stored efficiently, avoiding the occurrence of any convection phenomenon. Yet, conduction and thermal dispersion led to the formation of a plume that was monitored by 4D ERT (see Appendix). All wells and piezometers referred to in the study site had dataloggers; they recorded water pressure and temperature during the entire experiment on a minute-to-minute basis.
<figure><img src="image_4.png" /><figcaption>**Figure 2.** ATES well doublet experiment. Water is extracted from the aquifer at 10°C in the pumping well (step 1), it is heated with the water heater (step 2), injected into the aquifer through the injection well (step 3), stored (step 4) and recovered from the injection well (step 5).</figcaption></figure>

The observed peak temperature did not exceed 36.5 °C during the injection phase, due to heat loss along the injection pipe and mixing of injected water and groundwater (Figure 3A). During the storage phase, a significant temperature decrease was measured (13 K loss) during the first 24 h and a gentler decrease afterwards. On day 3, the pump was activated to test the submerged pump, causing the quick temperature rise and decline seen on Figure 3. At the beginning of the pumping phase, a very brief decrease in temperature occurred, followed by a 4 K rise before declining again. This temperature rise through pumping is believed to be linked with groundwater flowing towards the well. Since the heat plume globally ascends in the aquifer, it is assumed that the temperature recorded by the datalogger in W2 (placed at a fixed depth of 5 m) is not representative of the aquifer during the storage phase. Radial flow towards the
injection well brings warmer water into the well. Then, the temperature decreases following the energy recovery process.

During the recovery phase (Figure 3B), the temperature drops due to recovery from 20.8 to 14.8 °C. After 4 h 30 min of recovery, the pump was stopped and started again a bit later during 30 min and definitely stopped. During the pump interruption a slight temperature rise was observed, followed by a gentle decrease over time during 10 days until it reached the constant simulated temperature.

<figure><img src="image_5.png" /><figcaption>**Figure 3.** Observed and simulated temperature evolution in the injection well (W2) during the ATES experiment (A), with a zoom on the recovery phase (B) highlighting the observed and simulated heat recovery efficiency.</figcaption></figure>

The energy recovery efficiency was estimated based on the injected energy plotted on this figure. The energy recovery efficiency of the observed trend is 61.7 % of the injected energy. The volume of injected water was recovered after 1 h 51 min, with 30 % of energy recovered
at that time (Table 1). At the end of the recovery process, 280 % of the injected water volume was pumped.

<table><caption>**Table 1.** Evolution of recovered water volume and energy (observed and simulated), and temperature (observed and simulated) during the recovery phase.</caption><thead><tr><th rowspan="2">Recovery duration</th><th rowspan="2">Injected water volume recovered (%)</th><th colspan="2">Energy recovered (%)</th><th colspan="2">Temperature (°C)</th></tr><tr><th>Obs.</th><th>Sim.</th><th>Obs.</th><th>Sim.</th></tr></thead><tbody><tr><td>0h</td><td>0</td><td>0</td><td>0</td><td>20.8</td><td>20.9</td></tr><tr><td>1h 51'</td><td>100</td><td>29.6</td><td>30.3</td><td>16.8</td><td>13.7</td></tr><tr><td>4h 30'</td><td>280</td><td>61.7</td><td>58.2</td><td>14.8</td><td>12.1</td></tr></tbody></table>

No significant temperature variations were observed and simulated in well W1 and piezometers PzA, PzB and PzC.
# 4 Numerical model

A hydrogeological model was created to reproduce the ATES experiment by means of the subsurface flow and transport FEFLOW model code. In this study, coupled subsurface water flow under variably saturated conditions and heat transport were simulated, taking variations of density and viscosity into account. The simulated thermal plume was then compared to observation data collected through 4D ERT monitoring experiments (see Appendix).

## 4.1 Subsurface flow

Constant hydraulic head boundary conditions were prescribed along the northern and southern lateral boundaries of the model to represent the local hydraulic gradient (i = 0.0001), with values of -1.13 m and -1.33 m respectively. The eastern and western lateral boundaries of the model are no-flow boundaries. The same applies to the bottom boundary of the model.

The ATES experiment extraction and injection rates were assigned to the pumping well (W1) and the injection well (W2) respectively. Both wells were given a 0.24 m diameter and were screened along the entire aquifer thickness, assuming radial convergent or divergent flow during the pumping or injection phases.

The three hydrogeological units embodied in the model were given initial hydraulic parameters based on field measurements and values found in the literature (Table 2). The clayey loam and shale units were given average values for their hydraulic conductivity, porosity and residual saturation. The hydraulic conductivity values of the aquifer were determined based on data from a step-drawdown pumping test performed prior to the ATES experiment. Porosity and residual saturation values were taken from Epting et al. [34], who developed a heat transfer model in an aquifer similar to our alluvial aquifer. The modified Mualem-van Genuchten formulation was used to account for unsaturated hydraulic properties [35]. The hydraulic parameters were set with initial values, and some were calibrated against observed data (Table 2).
<table><caption>**Table 2.** Hydraulic parameters of the numerical model for the 3 hydrogeological units.</caption><thead><tr><th>Hydrogeological unit</th><th>Parameter</th><th>Initial value</th><th>Calibrated value</th><th>Source</th></tr></thead><tbody><tr><td rowspan="5">Aquitard<br>(clayey loam)</td><td>Horizontal hydraulic conductivity, K<sub>XY</sub> (m s<sup>-1</sup>)</td><td>1.0·10<sup>-6</sup></td><td>-</td><td>-</td></tr><tr><td>Vertical hydraulic conductivity, K<sub>Z</sub> (m s<sup>-1</sup>)</td><td>1.0·10<sup>-7</sup></td><td>-</td><td>-</td></tr><tr><td>Porosity, n<sub>t</sub> (-)</td><td>0.45</td><td>-</td><td>[36]</td></tr><tr><td>Residual saturation, S<sub>r</sub> (-)</td><td>0.08</td><td>-</td><td>[36]</td></tr><tr><td>van Genuchten parameter α (m<sup>-1</sup>)</td><td>1.58</td><td>-</td><td>[36]</td></tr><tr><td rowspan="5">Aquifer<br>(sandy gravel)</td><td>van Genuchten parameter n (-)</td><td>1.41</td><td>-</td><td>[36]</td></tr><tr><td>Horizontal hydraulic conductivity, K<sub>XY</sub> (m s<sup>-1</sup>)</td><td>1.0·10<sup>-3</sup></td><td>5.0·10<sup>-5</sup> - 5.0·10<sup>-3</sup></td><td>[10]</td></tr><tr><td>Vertical hydraulic conductivity, K<sub>Z</sub> (m s<sup>-1</sup>)</td><td>1.0·10<sup>-4</sup></td><td>2.0·10<sup>-5</sup> - 2.0·10<sup>-4</sup></td><td>[10]</td></tr><tr><td>Porosity, n<sub>t</sub> (-)</td><td>0.38</td><td>0.38</td><td>[34]</td></tr><tr><td>Residual saturation, S<sub>r</sub> (-)</td><td>0.12</td><td>-</td><td>[34]</td></tr><tr><td rowspan="5">Aquiclude<br>(weathering clay)</td><td>van Genuchten parameter α (m<sup>-1</sup>)</td><td>3.99</td><td>-</td><td>[34]</td></tr><tr><td>van Genuchten parameter n (-)</td><td>3.50</td><td>-</td><td>[34]</td></tr><tr><td>Horizontal hydraulic conductivity, K<sub>XY</sub> (m s<sup>-1</sup>)</td><td>1.0·10<sup>-8</sup></td><td>-</td><td>-</td></tr><tr><td>Vertical hydraulic conductivity, K<sub>Z</sub> (m s<sup>-1</sup>)</td><td>1.0·10<sup>-9</sup></td><td>-</td><td>-</td></tr><tr><td>Porosity, n<sub>t</sub> (-)</td><td>0.42</td><td>-</td><td>[36]</td></tr><tr><td></td><td>Residual saturation, S<sub>r</sub> (-)</td><td>0.10</td><td>-</td><td>[36]</td></tr></tbody></table>

## 4.2 Heat transport

With regard to simulated heat transport, meteorological parameters input such as air temperature or wind speed are applied with recharge [37]. A time-varying heat-transport boundary condition (Cauchy type), based on the air temperature recorded during the ATES experiment, was applied to the nodes of the top model surface.

The eastern, southern, western, and bottom boundaries were given no heat-flux boundary conditions. A 10 °C fixed temperature boundary condition (Dirichlet type) was applied to the northern boundary of the aquifer and aquiclude units, since water is flowing through the model domain from the northern boundary towards the southern one.

A heat nodal source was set in the injection well (W2) 5 m below ground surface, as implemented in the field. A power of 7.69·10⁹ J/d was applied to that node during the injection phase, which corresponds to an injection of water with a ΔT of 30 K at 2.55 m³/h during 5 h 20 min. Before and after the injection phase, no boundary condition was applied to that node.
Heat transport parameters are listed in Table 3; most of the parameters were selected from similar soils properties found in the literature. No specific thermal parameters measurements were taken in this case. The properties found in Epting et al. [36] were chosen to define initial parameter values in our model since the hydrogeological context of the aquifer they investigated is similar to the one presented here. The initial longitudinal and transverse dispersivity values were evenly assigned to all three units, with respective values of 5 m and 0.5 m. As with the hydraulic parameters, the heat transport parameters were calibrated (Table 3). The standard FEFLOW water thermal parameters were not modified.

The convective form of the thermal transport equation was chosen [38], with fluid viscosity dependency on temperature based on an empirical equation [39]. Since the ΔT of the injected water was significant but did not exceed 40 °C, nonlinear fluid density dependency on temperature was considered, as illustrated by Molson et al. [37].

<table><caption>**Table 3.** Heat transport parameters for the 3 hydrogeological units.</caption><thead><tr><th>Hydrogeological unit</th><th>Parameter</th><th>Initial value</th><th>Calibrated value</th><th>Source</th></tr></thead><tbody><tr><td rowspan="5">Aquitard (clayey loam)</td><td>Volumetric heat capacity of solid, c<sub>s</sub> (J m<sup>-3</sup> K<sup>-1</sup>)</td><td>3.0·10<sup>6</sup></td><td>-</td><td>[40]</td></tr><tr><td>Thermal conductivity of solid, λ<sub>s</sub> (J m<sup>-1</sup> s<sup>-1</sup> K<sup>-1</sup>)</td><td>1.9</td><td>-</td><td>[40]</td></tr><tr><td>Effective porosity, n<sub>e</sub> (-)</td><td>0.06</td><td>-</td><td>[41]</td></tr><tr><td>Longitudinal dispersivity, α<sub>L</sub> (m)</td><td>5</td><td>-</td><td>[42]</td></tr><tr><td>Transverse dispersivity, α<sub>T</sub> (m)</td><td>0.5</td><td>-</td><td>[42]</td></tr><tr><td rowspan="5">Aquifer (sandy gravel)</td><td>Volumetric heat capacity of solid, c<sub>s</sub> (J m<sup>-3</sup> K<sup>-1</sup>)</td><td>2.87·10<sup>6</sup></td><td>2.0·10<sup>6</sup></td><td>[10]</td></tr><tr><td>Thermal conductivity of solid, λ<sub>s</sub> (J m<sup>-1</sup> s<sup>-1</sup> K<sup>-1</sup>)</td><td>2.7</td><td>3.94</td><td>[43]</td></tr><tr><td>Effective porosity, n<sub>e</sub> (-)</td><td>0.12</td><td>0.10</td><td>[34]</td></tr><tr><td>Longitudinal dispersivity, α<sub>L</sub> (m)</td><td>20</td><td>5</td><td>[34]</td></tr><tr><td>Transverse dispersivity, α<sub>T</sub> (m)</td><td>2</td><td>0.5</td><td>[34]</td></tr><tr><td rowspan="5">Aquiclude (weathering clay)</td><td>Volumetric heat capacity of solid, c<sub>s</sub> (J m<sup>-3</sup> K<sup>-1</sup>)</td><td>2.3·10<sup>6</sup></td><td>-</td><td>[40]</td></tr><tr><td>Thermal conductivity of solid, λ<sub>s</sub> (J m<sup>-1</sup> s<sup>-1</sup> K<sup>-1</sup>)</td><td>1.5</td><td>-</td><td>[40]</td></tr><tr><td>Effective porosity, n<sub>e</sub> (-)</td><td>0.05</td><td>-</td><td>[44]</td></tr><tr><td>Longitudinal dispersivity, α<sub>L</sub> (m)</td><td>5</td><td>-</td><td>[42]</td></tr><tr><td>Transverse dispersivity, α<sub>T</sub> (m)</td><td>0.5</td><td>-</td><td>[42]</td></tr></tbody></table>
## 4.3 Initial conditions

When performing transient simulations, setting proper initial conditions are of prime importance to prevent simulations results from being distorted [45]. The model was first run for a period of 2 years to set the initial conditions, by applying boundary conditions of the very first set of time-varying data available, starting from a fully saturated situation. The resulting calculated hydraulic head distribution was then set as initial conditions for water flow. Regarding heat transport, an initial temperature of 10 °C was set in the entire model in accordance with the temperature measured in the wells and piezometers prior to the ATES experiment.

## 4.4 Temporal discretization

Simulations were run on a 64 Go RAM computer with 2 multi-core processors at 2.20 GHz (20 physical cores, 40 logical processors). The model was run over a 15-days period, starting with a warm-up day, followed by the 4-days ATES field experiment and a 10-days additional period. The automatic time-step control option implemented in FEFLOW was chosen, with which the model defines the time-steps size by controlling the numerical calculations with a Euclidian $L_2$ integral root mean square error criterion [46]. Convergence criteria are met when an error criterion value of $\delta = 5 \cdot 10^{-4}$ is reached; the model then initiates the next time-step calculation. Conservative time-step constraints are to be set to avoid numerical oscillations and instability [47]. An initial time-step size of 0.0001 day was chosen here, with a growth factor between subsequent time-steps of 1.2, and the maximum time-step size was limited to 0.1 day.
## 4.5 Model calibration

The model was automatically calibrated using FEPEST, the FEFLOW integrated version of the model-independent parameter estimation code PEST [48], by means of the state-of-the art pilot points inverse modeling method [49]. The FEPEST workflow is documented by Goretzki et al. [41]. The calibration process consisted of two steps. First, the groundwater flow model was calibrated against the hydraulic head data recorded in wells and piezometers by means of 200 pilot points. Then, the heat transport model was calibrated against the temperature measured in the same wells and piezometers.

The groundwater flow calibration process focused on the adjustment of the aquifer unit horizontal hydraulic conductivity ($K_{xy}$). The vertical hydraulic conductivity ($K_z$) was calculated in accordance with $K_{xy}$ by applying a 0.1 anisotropy factor. Calibrated values of $K_{xy}$ range from $5.0 \times 10^{-3}$ to $5.0 \times 10^{-5}$ m/s (Figure 4). The other hydraulic parameters found in Table 2 were insensitive and therefore not included in the calibration process. In fact, Hermans et al. [27] showed that in similar heat storage and recovery simulations, the most sensitive parameters (in a global sensitivity analysis) were the hydraulic conductivities of their model, confirming our findings.
<figure><img src="image_6.png" /><figcaption>**Figure 4.** Model structure, with the horizontal (a) and vertical (b) discretization exposed. A zoom on the experimental area shows the mesh refinement along ERT profile lines, and around wells and piezometers. Vertically, the 16 model layers are displayed. The aquitard, aquifer and aquiclude are respectively made of 5, 8 and 3 layers of various thicknesses.</figcaption></figure>

During the heat transport calibration phase, 4 parameters of the aquifer unit were calibrated: the heat capacity of the porous medium matrix, its thermal conductivity, the porosity and the longitudinal dispersivity. The transverse dispersivity value was constantly linked to the longitudinal dispersivity value with a 0.1 anisotropy factor. Since the heat plume did not reach PzA or PzB during the injection phase, no other observation data than the injected warm water in W2 were exploitable. The use of the pilot points method to calibrate the heat transport parameters mentioned previously was therefore unnecessary and we considered the thermal
properties of the aquifer as homogeneous. The calibrated values of these parameters are listed
in Table 3.

The 4D ERT monitoring measurements of the thermal affected zone were used for validating
the heat transport model as described in the Appendix.

The simulated temperature curve is shown in Figure 3. The model was calibrated to peak at
the same temperature. The storage phase observed and simulated temperature curves are
significantly different because the model was calibrated to fit the observed recovery peak
temperature at 20.8 °C, as seen on Figure 3B, having in mind that that temperature is
representative of the actual groundwater temperature. As for the recovery phase (Figure 3B),
both curves are somewhat different. The slight temperature rise observed 5 h after starting the
recovery phase could not be reproduced in the model. Hence the temperature drop due to
recovery (from 20.8 to 12.3 °C) was of higher amplitude than what was observed (from 20.8 to
14.8 °C). The temperature rise after 4h30 of recovery cannot be reproduced by the model. After
the final stop, the simulated temperature was a bit lower but constant. Snapshots of the thermal
plume vertical extension at the end of the injection, storage and recovery phases are exhibited
in Figure 5. The heated water injection spot along W2 is clearly seen at the end of the injection
phase. The plume extension snapshots were taken at temperatures higher than 11 °C. It
continues to extend even after the injection phase due to heat conduction along the aquifer full-
thickness. Heat convection is the main transport process that occurs during injection and
recovery, while heat conduction is the dominant process during the storage phase [50].
Conduction also leads to partially warming the confining aquitard and underlying aquiclude.
The influence of the aquifer hydraulic conductivity heterogeneities, as displayed for K<sub>XY</sub> on
Figure 4, is seen on Figure 5 with an asymmetric thermal plume extending in the X direction
on either side of well W2. The vertical extension of the plume is limited by the aquitard and
aquiclude units low hydraulic conductivities.
No significant temperature variations were observed and simulated in well W1 and piezometers PzA, PzB, and PzC.

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Simulated heat plume in 2D vertical view along the Y = 12.5 m profile, which includes the injection well (W2). Three snapshots are shown at the end of the injection, storage and recovery phases.</figcaption></figure>
# 5 Predictive simulations set-up

With the help of our calibrated numerical model, we ran 77 different predictive simulation scenarios to assess short-term ATES, the scope being DSM applications. We constrained our simulations within the same area, duration (less than one week), and season of the year (May 2015) as our experiment to gain advantage of more robust predictions.

All scenarios follow the same sequence. First, groundwater is pumped from well W1 at a fixed pumping rate of 15 m³/h. The initial groundwater temperature is 10 °C, as it was measured in-situ. This pumped groundwater is then reinjected in the aquifer in well W2 at an equivalent flow rate of 15 m³/h. The only difference between the simulations lies in the injected water temperature since we impose a specific temperature difference ΔT (-4, 3, 6, 11 K for LT-ATES, and from 35 to 65 K with a 5 K increment for HT-ATES), which consequently allows us to estimate the total amount of energy stored in the aquifer. After a specific storage period Δt (0.25, 1, 6, 12, 24, 48, and 72 h, corresponding to real-time, intraday, and interday DSM frequencies) varying from one scenario to another, the stored warm water is recovered from well W2 at a constant flow rate of 15 m³/h and for a constant duration of 48 h. Energy recovery rates are then calculated.

The pumping and injection flow rates were set at 15 m³/h for all scenarios because it corresponds to the average critical flow rate of the wells installed on the study site. Long-term pumping rates higher than 15 m³/h could not be sustained in the local alluvial aquifer because of a low aquifer recharge; the aim here being to simulate scenarios that are consistent with real-world applications. Nevertheless, pumping and injection flow rates of 15 m³/h, associated with a temperature reduction of 3, 6, and 11 K, could potentially generate heating powers of respectively 52, 104, and 192 kW that is typical of operational LT-ATES systems. Note that the imposed 48 h duration of the recovery phase is arbitrary but long enough to recover most of the energy that was previously stored.
We simulated 3 groups of ATES systems:

1. LT-ATES for space heating or DHW production with the help of a GWHP in an aquifer preheated at 3, 6, or 11 K, corresponding to regular operational ΔT values found in the literature, and mostly related to state regulations [18].

2. LT-ATES for air conditioning with cold storage at ΔT = -4 K, also based on operational references [18].

3. HT-ATES for space heating or DHW production without upgrading (no use of GWHP) by trying to recover the highest absolute temperature (ideally, 45°C or higher) and handling ΔT values of 35, 40, 45, 50, 55, 60, and 65 K.

Note that ATES is not optimized in our work since we focus only on a single heat (or cold) storage and recovery cycle.
# 6 Results and discussion

The energy recovery ($\eta$) results of our 77 simulations (7 different $\Delta T$, and 11 different storage durations $\Delta t$) are compiled in Figure 6 ($\eta$ vs. $\Delta T$), Figure ($\eta$ vs. $\Delta t$), Figure ($T_{Recovered}$ vs. $T_{Injected}$), and Table 4 (summary).

In terms of energy recovery rates, values vary between 53 % and 87 % (Table 4). The highest value of 87 % is obtained for scenarios simulating LT-ATES ($\Delta T$ = -4, 3, 6, and 11 K) and real-time to intraday frequencies ($\Delta t$ = 0.25, 1, 6, and 12 h). Energy recovery remains constant in these 16 scenarios. The scenario with the highest injected temperature ($\Delta T$ = 65 K, $T$ = 75 °C) and an interday frequency ($\Delta t$ = 72 h) result in the lowest value of 53 %.

As seen in Figure 6, energy recovery rates first remain constant at low temperatures ($\Delta T$ = -4, 3, 6, and 11 K) and then decrease with increasing temperature ($\Delta T$ = 35 to 65 K). This behavior remains similar throughout the 7 different storage periods considered (plots of Figure 6 are parallel).

<table><caption>**Table 4.** Energy recovery rate values ($\eta$, in %) for the 77 DSM simulated scenarios (7 storage periods $\Delta t$, 11 temperature differences induced $\Delta T$).</caption><thead><tr><th rowspan="2">η (%)</th><th rowspan="2"></th><th colspan="11">ΔT (K)</th></tr><tr><th>-4</th><th>3</th><th>6</th><th>11</th><th>35</th><th>40</th><th>45</th><th>50</th><th>55</th><th>60</th><th>65</th></tr></thead><tbody><tr><td rowspan="7">Δt (h)</td><td>0.25</td><td>87</td><td>87</td><td>87</td><td>87</td><td>71</td><td>71</td><td>70</td><td>68</td><td>67</td><td>66</td><td>65</td></tr><tr><td>1</td><td>87</td><td>87</td><td>87</td><td>87</td><td>71</td><td>71</td><td>69</td><td>68</td><td>67</td><td>66</td><td>65</td></tr><tr><td>6</td><td>87</td><td>87</td><td>87</td><td>87</td><td>71</td><td>71</td><td>70</td><td>68</td><td>67</td><td>66</td><td>65</td></tr><tr><td>12</td><td>86</td><td>87</td><td>87</td><td>87</td><td>71</td><td>70</td><td>69</td><td>68</td><td>66</td><td>65</td><td>64</td></tr><tr><td>24</td><td>84</td><td>84</td><td>84</td><td>83</td><td>67</td><td>67</td><td>65</td><td>64</td><td>63</td><td>61</td><td>61</td></tr><tr><td>48</td><td>84</td><td>83</td><td>83</td><td>84</td><td>67</td><td>66</td><td>65</td><td>63</td><td>62</td><td>60</td><td>59</td></tr><tr><td>72</td><td>78</td><td>78</td><td>78</td><td>78</td><td>62</td><td>60</td><td>58</td><td>57</td><td>55</td><td>54</td><td>53</td></tr></tbody></table>

The energy recovery rates first remain constant for real-time to intraday frequencies ($\Delta t$ = 0.25, 1, 6, and 12 h) and then decrease when the storage duration is increased ($\Delta t$ = 24, 48, and 72 h, representing interday frequencies) (Figure 6). This behavior is similar for the 11 different injected temperature considered (plots of Figure 6 are parallel).
<figure><img src="image_8.png" /><figcaption>**Figure 6.** Comparison between energy recovery rate values and temperature differences.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Comparison between energy recovery rate values and storage periods.</figcaption></figure>

In terms of absolute temperatures, even after a very short storage duration (i.e. 0.25 h), the recovered temperature already decreases by a few K. The temperature drop (between the injected and recovered temperatures) increases in parallel to the injected temperature but also when storage duration is increased (Figure ).
<figure><img src="image_10.png" /><figcaption>**Figure 8.** Comparison between injected and right before recovery absolute temperatures.</figcaption></figure>

The analysis of the results shows that:

- Energy recovery rates essentially decrease when storage duration increases, and energy recovery rates are therefore higher in real-time than at intraday and higher at intraday than at interday frequencies.

- Energy recovery rates decrease when stored temperature increases and LT-ATES is more efficient in terms of energy recovery rates than HT-ATES (for a single storage and recovery cycle and without considering exergy).

- Recovered temperatures decrease rapidly when stored temperature and storage duration increase.

In quantitative terms, LT-ATES systems (ΔT of -4, 3, 6, and 11 K in our simulations) operated at real-time (15 min) and intraday (1, 6, and 12 h) frequencies present an almost constant energy recovery rate value between 86 and 87 %. The latter value drops to 83 % (after 48 h) and 77 % (after 72 h) for interday frequencies. HT-ATES systems (ΔT of 35 to 65 K with an increment of 5 K in our simulations) present lower energy recovery rate values, with a decrease from 70 to 63 % for real-time to interday frequencies (0.25, 1, 6, and 12 h), from 67
to 59 % after 24 and 48 h, and from 62 to 53 % after 72 h (interday frequency). It should be
noted that, since energy recovery was underestimated for the real case experiment, the
simulated energy recovery of our 77 scenarios are likely to have been underestimated too.
# 7 Conclusions and perspectives

Short-term aquifer thermal energy storage should be further investigated through experimental and numerical developments for flexibility purposes by preheating the aquifer (T < 30 °C) to improve the performance coefficient of groundwater heat pumps, and by directly storing potentially useful heat (at higher temperatures for water recovery at T > 50 °C) for space heating or domestic hot water production (no upgrading with groundwater heat pumps). Our study, which is representative of productive shallow alluvial aquifers with slow ambient groundwater flow (a few meters per year), demonstrates that warm or cold water can be stored during off-peak periods and recovered during peak periods (at real-time, intraday and interday frequencies) with energy recovery rates up to 90 %. For a single aquifer thermal energy storage and recovery cycle, low-temperature storage presents higher energy recovery rates (from 78 to 87 %, according to our predictive scenarios) than high-temperature storage (from 53 to 71 %, according to our predictive simulations). In addition, energy recovery rates decrease with increasing storage durations and with increasing working temperatures.

Before aiming to apply this research to real-world, further studies should be undertaken to investigate the behavior of other types of aquifer (e.g. shallow alluvial aquifers with faster ambient groundwater flow, fractured rock reservoirs, lower porosity aquifers, etc.). Future developments aiming for the optimization of thermal energy storage and recovery in space (system sizing), in time (cyclicality), and in terms of absolute temperatures should also be investigated. Automatic tools to control those systems should also receive attention. The interested professional community would highly benefit from more research being carried out to address the technical and economical constraints, as well as from the implementation of a legislative context, if needed, and life-cycle assessment of short-term ATES systems.
Acknowledgments

The present work was funded by the BEWARE Fellowships Industry and BEWARE Fellowships Academia programs of the Wallonia-Brussels Federation (Belgium), as part of the SMARTMODEL and SUITE4D projects under contracts No. 1610056 and No. 1510466, respectively. Both projects were also co-financed by the FP7 Marie Skłodowska-Curie COFUND program of the European Union (Project No. 600397). The authors also wish to thank the journal editors and three anonymous reviewers for their insightful comments and suggestions that helped improving this paper. The data used to produce the results of this paper can be obtained from the corresponding author upon request.
# 8 Appendix

## 8.1 Numerical model spatial discretization

The model was built by creating a squared 2D mesh that is 2000 m x 2000 m and centered on the study site. As the lateral boundaries were distant from the study site, they were expected not to influence the simulations results. The model was discretized with the FEFLOW integrated Triangle mesh generator [51], which ensures mesh quality by respecting the 2D Delaunay in-circle criterion [52]. The mesh was refined along the ERT profile lines, with a maximal inter-nodal distance of 2 m which was the distance between two successive electrodes along the profile lines, in order to be able to compare thermal plume simulated results and observed ERT data. The mesh was also refined around the wells and piezometers of the study site. Based on the drilling diameter (0.24 m) of the piezometers and wells logs, an inter-nodal distance of 0.10 m was chosen to avoid numerical dispersion errors in the calculations of heat transport [37]; because heat convection calculations are highly sensitive to spatial discretization [47], especially since the thermal radius is known to be smaller than the hydraulic radius around heat injection wells in ATES systems [25]. Overall, the triangular elements edges range from 0.1 m (close to the wells and piezometers) to 75 m (along the model lateral boundaries).

The resulting 2D triangular mesh contained 7,267 nodes and 14,375 triangular elements. This 2D mesh was replicated in the third dimension with elevations ranging from 0 m to -10 m, and a total of 16 layers created (Figure 4). The 3D model was consequently composed of 123,539 nodes and 230,000 elements. The top clayey loam horizon was discretized with 5 layers of variable thicknesses. The alluvial aquifer was made of 8 layers while the 3 remaining layers constituted the deeper weathering clay aquiclude layer. A 0.1 m thick layer was set at the top of the model domain to avoid water infiltration calculations instabilities due to time-varying boundary conditions applied to the model surface. Thinner layers (0.1 m thick) were accounted for on both sides of lithological interfaces (at -2 m and -7.5 m depths) to avoid coarse
discretization, and consequently thermal calculation instabilities [47], at sharp hydraulic conductivity changeover fronts.

## 8.2 4D geophysical monitoring for model validation

Monitoring the injection and extraction of warm water in an ATES system is easily feasible through injection, production, and observation wells. However, tracking the extension of a heat/cold plume in an aquifer through boreholes only would be representative of the vicinity of the boreholes at a specific moment in time [53]. Recent research demonstrated the ability of geophysical techniques (analogous to medical imaging techniques but designed to auscultate the subsurface), and especially electrical resistivity tomography (ERT), to monitor physical processes occurring in the subsurface [54], including heat storage and recovery or heat convection [55]. Because bulk electrical resistivity (imaged with ERT) depends on temperature [56] among other parameters, ERT monitoring combined with classical hydrogeological measurements in wells were demonstrated by Hermans et al. [57] to adequately track a heat plume evolution in space (heterogeneity) and in time (quantitatively). The interested reader is referred to Hermans et al. [58] and Lesparre et al. [59] for detailed information on ERT monitoring of heat transport in the subsurface.

Bearing this in mind, the full ATES experiment was monitored by surface 4D ERT (4D stands for the 3 dimensions in space and the temporal dimension), in addition to the dataloggers in boreholes, as such devices are able to image the thermal affected zone [59]. The ERT system consisted of 9 profile lines with an equally spaced inter-distance of 3 m, each line being 61 m long. Each line had 21 electrodes with a 2-m inter-spacing, except for the two last pairs of electrodes at both ends that had a 4-m inter-spacing (Figure 1b). This experimental setup covered an area of 225 m². Electrical resistivity data were acquired with a combination of dipole-dipole and multiple gradient electrode arrays following the work of Van Hoorde et al.
[60], which were optimized for multichannel acquisition. Two background images were acquired to assess the background noise in ERT images, before the heat injection process. A noise level of 0.5 % was estimated after reading Robert et al. [61] and subsequently used to filter ERT images. Above this threshold value, electrical resistivity variations were estimated to express temperature variations accurately [57].

Several ERT images were recorded during the storage phase; 6 of these are displayed in Figure A, 25, 27 and 47 hours after the injection phase with 2 temperature distribution plots in the horizontal plane at -4 m and -5.5 m depth (dashed contours + solid line at 11 °C). The simulated temperature is also displayed for comparison purposes (color scale). ERT measurements provide bulk conductivity data that are dependent on the electrode spacing and therefore representative of a given volume of aquifer (matrix and pore space). Since the local hydraulic gradient is very low, the monitored and simulated thermal plumes are centered on the injection well W2; no convection process is involved. Yet, an offset of -0.5 m in the ordinate's direction is visible at the highest temperature of the monitored plume when compared to that of the simulated plume.

The highest simulated temperature at the injection point was reached at the end of the storage phase (36.5 °C). On the other hand, the monitored ERT thermal plume peaked at 12 °C, between 27 to 47 h after injection. In addition, the maximum extension of the monitored plume was reached between 27 and 47 h after injection (maximum temperature: 12 °C), while the simulated plume continued to expending gently during the storage phase (Figure A).

When looking closely at the 11 °C contours in Figure A, the monitored and simulated thermal plumes appear similar, regardless of the little offset and extension variations in the monitored data. This suggests that the initial simplifications proposed for the geophysical data inversion and the model conceptualization and calibration were correctly hypothesized. The calibrated groundwater flow and heat transport model was validated with the help of these 4D
ERT monitoring data, following the example of Hermans et al. [58]. A linear relation allows the establishment of a link between electrical resistivity variations and temperature variations in the 10 – 40 °C interval [56]. In addition, coefficients coupling electrical conductivity to temperature variations were defined in previous studies [55]. Temperature variations can be derived from ERT measurements using these relations. Yet, the smoothing of geophysical measurements in their inversion process tends to overestimate the plume extension. As observed in Figure A, the monitored plume was slightly larger than the simulated one. In addition, temperature data derived from 3D ERT images can be misinterpreted due to background local electrical conductivity values reported to be higher than natural. Using ERT data as hard data for our model calibration by tracking the heat plume extension in time-lapse was extremely difficult since the error caused by the high groundwater electrical conductivity values was unknown. The data quality however was good enough to validate our calibration procedure.
<figure><img src="image_11.png" /><figcaption>**Figure A.** Heat plume ERT monitoring during the storage phase at - 4 m, - 5.5 m and - 7 m depth. The simulated heat plume is displayed in the background.</figcaption></figure>
References

[1] Brounen D, Kok N, Quigley JM. Residential energy use and conservation: Economics and demographics. European Economic Review 2012;56:931-45. doi:10.1016/j.euroecorev.2012.02.007.

[2] Filippini M, Pachauri S. Elasticities of electricity demand in urban Indian households. Energy Policy 2004;32:429-36. doi:10.1016/S0301-4215(02)00314-2.

[3] Besagni G, Borgarello M. The determinants of residential energy expenditure in Italy. Energy 2018;165:369-86. doi:10.1016/j.energy.2018.09.108.

[4] European Council. Directive 2012/27/EU of the European Parliament and of the Council of 25 October 2012 on energy efficiency, amending Directives 2009/125/EC and 2010/30/EU and repealing Directives 2004/8/EC and 2006/32/EC (Text with EEA relevance). Official Journal of the European Union 2012:1-56.

[5] Palensky P, Dietrich D. Demand Side Management: Demand Response, Intelligent Energy Systems, and Smart Loads. IEEE Transactions on Industrial Informatics 2011;7:381-8. doi:10.1109/TII.2011.2158841.

[6] Arteconi A, Hewitt NJ, Polonara F. State of the art of thermal storage for demand-side management. Applied Energy 2012;93:371-89. doi:10.1016/j.apenergy.2011.12.045.

[7] Fischer D, Madani H. On heat pumps in smart grids: A review. Renewable and Sustainable Energy Reviews 2017;70:342-57. doi:10.1016/j.rser.2016.11.182.

[8] Arteconi A, Patteeuw D, Bruninx K, Delarue E, D'haeseleer W, Helsen L. Active demand response with electric heating systems: Impact of market penetration. Applied Energy 2016;177:636-48. doi:10.1016/j.apenergy.2016.05.146.

[9] Baetens R, De Coninck R, Van Roy J, Verbruggen B, Driesen J, Helsen L, et al. Assessing electrical bottlenecks at feeder level for residential net zero-energy buildings by integrated system simulation. Applied Energy 2012;96:74-83. doi:10.1016/j.apenergy.2011.12.098.

[10] Paulus C, Bolly P-Y, Hermans T, Koo Seen Lin E, Robert T. Stockage de chaleur en aquifère et flexibilité de la demande électrique : quelles possibilités ? Proceedings of the 34èmes Rencontres de l'AUGC, University of Liège, Belgium 2016:18.

[11] Lund JW, Boyd TL. Direct Utilization of Geothermal Energy 2015 Worldwide Review. Proceedings of the World Geothermal Congress 2015, Melbourne, Australia: 2015, p. 1-31.

[12] Wesselink M, Liu W, Koornneef J, van den Broek M. Conceptual market potential framework of high temperature aquifer thermal energy storage - A case study in the Netherlands. Energy 2018;147:477-89. doi:10.1016/j.energy.2018.01.072.

[13] Vanhoudt D, Desmedt J, Van Bael J, Robeyn N, Hoes H. An aquifer thermal storage system in a Belgian hospital: Long-term experimental evaluation of energy and cost savings. Energy and Buildings 2011;43:3657-65. doi:10.1016/j.enbuild.2011.09.040.

[14] Kim J, Lee Y, Yoon WS, Jeon JS, Koo M-H, Keehm Y. Numerical modeling of aquifer thermal energy storage system. Energy 2010;35:4955-65. doi:10.1016/j.energy.2010.08.029.

[15] Lo Russo S, Civita MV. Open-loop groundwater heat pumps development for large buildings: A case study. Geothermics 2009;38:335-45. doi:10.1016/j.geothermics.2008.12.009.

[16] Hesaraki A, Holmberg S, Haghighat F. Seasonal thermal energy storage with heat pumps and low temperatures in building projects—A comparative review. Renewable and Sustainable Energy Reviews 2015;43:1199-213. doi:10.1016/j.rser.2014.12.002.
[17] Bayer P, Rybach L, Blum P, Brauchler R. Review on life cycle environmental effects of geothermal power generation. Renewable and Sustainable Energy Reviews 2013;26:446-63. doi:10.1016/j.rser.2013.05.039.

[18] Haehnlein S, Bayer P, Blum P. International legal status of the use of shallow geothermal energy. Renewable and Sustainable Energy Reviews 2010;14:2611-25. doi:10.1016/j.rser.2010.07.069.

[19] Haehnlein S, Bayer P, Ferguson G, Blum P. Sustainability and policy for the thermal use of shallow geothermal energy. Energy Policy 2013;59:914-25. doi:10.1016/j.enpol.2013.04.040.

[20] Bloemendal M, Olsthoorn T, Boons F. How to achieve optimal and sustainable use of the subsurface for Aquifer Thermal Energy Storage. Energy Policy 2014;66:104-14. doi:10.1016/j.enpol.2013.11.034.

[21] Rybach L, Kohl T. Waste heat problems and solutions in geothermal energy. Geological Society, London, Special Publications 2004;236:369-80. doi:10.1144/GSL.SP.2004.236.01.21.

[22] Gao L, Zhao J, Tang Z. A Review on Borehole Seasonal Solar Thermal Energy Storage. Energy Procedia 2015;70:209-18. doi:10.1016/j.egypro.2015.02.117.

[23] Chen H, Cong TN, Yang W, Tan C, Li Y, Ding Y. Progress in electrical energy storage system: A critical review. Progress in Natural Science 2009;19:291-312. doi:10.1016/j.pnsc.2008.07.014.

[24] Doughty C, Hellström G, Tsang CF, Claesson J. A dimensionless parameter approach to the thermal behavior of an aquifer thermal energy storage system. Water Resour Res 1982;18:571-87. doi:10.1029/WR018i003p00571.

[25] Bloemendal M, Hartog N. Analysis of the impact of storage conditions on the thermal recovery efficiency of low-temperature ATES systems. Geothermics 2018;71:306-19. doi:10.1016/j.geothermics.2017.10.009.

[26] Bloemendal M, Olsthoorn T. ATES systems in aquifers with high ambient groundwater flow velocity. Geothermics 2018;75:81-92. doi:10.1016/j.geothermics.2018.04.005.

[27] Hermans T, Lesparre N, De Schepper G, Robert T. Bayesian Evidential Learning: A field validation using push-pull tests. Hydrogeology Journal 2019;In Press.

[28] Bridger DW, Allen DM. Influence of geologic layering on heat transport and storage in an aquifer thermal energy storage system. Hydrogeology Journal 2014;22:233-50. doi:10.1007/s10040-013-1049-1.

[29] Ferguson G. Heterogeneity and Thermal Modeling of Ground Water. Ground Water 2007;45:485-90. doi:10.1111/j.1745-6584.2007.00323.x.

[30] Robert T, Hermans T, Lesparre N, De Schepper G, Nguyen F, Defourny A, et al. Towards a subsurface predictive-model environment to simulate aquifer thermal energy storage for demand-side management applications. Proceedings of the 10th International Conference on System Simulation in Buildings 2018:P017.

[31] Fossoul F, Orban P, Dassargues A. Numerical simulation of heat transfer associated with low enthalpy geothermal pumping in an alluvial aquifer. Geologica Belgica 2011;14.

[32] Allen A, Milenic D. Low-enthalpy geothermal energy resources from groundwater in fluvioglacial gravels of buried valleys. Applied Energy 2003;74:9-19.

[33] Ruthy I, Dassargues A. Carte hydrogéologique de Wallonie, Fleurus - Spy 47/1-2: Notice explicative [in French] 2014.

[34] Epting J, Händel F, Huggenberger P. Thermal management of an unconsolidated shallow urban groundwater body. Hydrology and Earth System Sciences 2013;17:1851-69. doi:10.5194/hess-17-1851-2013.
[35] Schaap MG, van Genuchten MT. A Modified Mualem-van Genuchten Formulation for Improved Description of the Hydraulic Conductivity Near Saturation. Vadose Zone Journal 2006;5:27-34. doi:10.2136/vzj2005.0005.

[36] Schaap MG, Leij FJ, van Genuchten MT. ROSETTA: a computer program for estimating soil hydraulic parameters with hierarchical pedotransfer functions. Journal of Hydrology 2001;251:163-76. doi:10.1016/S0022-1694(01)00466-8.

[37] Molson JW, Frind EO, Palmer CD. Thermal energy storage in an unconfined aquifer. 2. Model development, validation, and application. Water Resour Res 1992;28:2857-67. doi:10.1029/92WR01472.

[38] Diersch H. Feflow: finite element modeling of flow, mass and heat transport in porous and fractured media. New York, NY: Springer Berlin Heidelberg; 2013.

[39] Hassanizadeh SM. Modeling species transport by concentrated brine in aggregated porous media. Transport in Porous Media 1988;3. doi:10.1007/BF00235333.

[40] Somogyvári M, Bayer P, Brauchler R. Travel-time-based thermal tracer tomography. Hydrology and Earth System Sciences 2016;20:1885-901. doi:10.5194/hess-20-1885-2016.

[41] Goretzki N, Inbar N, Kühn M, Möller P, Rosenthal E, Schneider M, et al. Inverse Problem to Constrain Hydraulic and Thermal Parameters Inducing Anomalous Heat Flow in the Lower Yarmouk Gorge. Energy Procedia 2016;97:419-26. doi:10.1016/j.egypro.2016.10.038.

[42] Dassargues A. Hydrogeology: groundwater science and engineering. Taylor & Francis CRC Press. Boca Raton, USA: 2018.

[43] Markle JM, Schincariol RA, Sass JH, Molson JW. Characterizing the Two-Dimensional Thermal Conductivity Distribution in a Sand and Gravel Aquifer. Soil Science Society of America Journal 2006;70:1281. doi:10.2136/sssaj2005.0293.

[44] McWhorter DB, Sunada DK. Ground-water hydrology and hydraulics. Water Resources Publications, LLC. Highlands Ranch, Colorado, USA: Water Resources Publications, LLC; 1977.

[45] Refsgaard JC. Parameterisation, calibration and validation of distributed hydrological models. Journal of Hydrology 1997;198:69-97. doi:10.1016/S0022-1694(96)03329-X.

[46] Diersch H-JG. Error norms used in FEFLOW. White Pap II WASY GmbH, Institute for Water Resource Planning and System Research, Berlin 2009:109-15.

[47] Kaiser BO, Cacace M, Scheck-Wenderoth M. 3D coupled fluid and heat transport simulations of the Northeast German Basin and their sensitivity to the spatial discretization: different sensitivities for different mechanisms of heat transport. Environmental Earth Sciences 2013;70:3643-59. doi:10.1007/s12665-013-2249-7.

[48] Doherty J. PEST: Model-independent parameter estimation, user manual. 6th ed. Watermark Numerical Computing. Brisbane: 2016.

[49] Doherty J. Calibration and Uncertainty Analysis for Complex Environmental Models. Watermark Numerical Computing. Brisbane: 2015.

[50] Miotliński K, Dillon PJ. Relative Recovery of Thermal Energy and Fresh Water in Aquifer Storage and Recovery Systems. Groundwater 2015;53:877-84. doi:10.1111/gwat.12286.

[51] Shewchuk JR. Triangle: Engineering a 2D quality mesh generator and Delaunay triangulator. In: Lin MC, Manocha D, editors. Applied Computational Geometry Towards Geometric Engineering, Springer Berlin Heidelberg; 1996, p. 203-22.

[52] Gargallo-Peiró A, Roca X, Peraire J, Sarrate J. Optimization of a regularized distortion measure to generate curved high-order unstructured tetrahedral meshes. International Journal for Numerical Methods in Engineering 2015;103:342-63. doi:10.1002/nme.4888.
[53] Rau GC, Andersen MS, McCallum AM, Roshan H, Acworth RI. Heat as a tracer to quantify water flow in near-surface sediments. Earth-Science Reviews 2014;129:40–58. doi:10.1016/j.earscirev.2013.10.015.

[54] Binley A, Hubbard SS, Huisman JA, Revil A, Robinson DA, Singha K, et al. The emergence of hydrogeophysics for improved understanding of subsurface processes over multiple scales: The Emergence of Hydrogeophysics. Water Resources Research 2015;51:3837–66. doi:10.1002/2015WR017016.

[55] Hermans T, Nguyen F, Robert T, Revil A. Geophysical Methods for Monitoring Temperature Changes in Shallow Low Enthalpy Geothermal Systems. Energies 2014;7:5083–118. doi:10.3390/en7085083.

[56] Hayley K, Bentley LR, Gharibi M, Nightingale M. Low temperature dependence of electrical resistivity: Implications for near surface geophysical monitoring. Geophysical Research Letters 2007;34. doi:10.1029/2007GL031124.

[57] Hermans T, Wildemeersch S, Jamin P, Orban P, Brouyère S, Dassargues A, et al. Quantitative temperature monitoring of a heat tracing experiment using cross-borehole ERT. Geothermics 2015;53:14–26. doi:10.1016/j.geothermics.2014.03.013.

[58] Hermans T, Vandenbohede A, Lebbe L, Nguyen F. A shallow geothermal experiment in a sandy aquifer monitored using electric resistivity tomography. Geophysics 2012;77:B11–21.

[59] Lesparre N, Robert T, Nguyen F, Boyle A, Hermans T. 4D electrical resistivity tomography (ERT) for aquifer thermal energy storage monitoring. Geothermics 2019;77:368–82. doi:10.1016/j.geothermics.2018.10.011.

[60] Van Hoorde M, Hermans T, Dumont G, Nguyen F. 3D electrical resistivity tomography of karstified formations using cross-line measurements. Engineering Geology 2017;220:123–32. doi:10.1016/j.enggeo.2017.01.028.

[61] Robert T, Caterina D, Deceuster J, Kaufmann O, Nguyen F. A salt tracer test monitored with surface ERT to detect preferential flow and transport paths in fractured/karstified limestones. Geophysics 2012;77:B55–67. doi:10.1190/geo2011-0313.1.