Investigations into the First Operational Aquifer Thermal Energy Storage System in Wallonia (Belgium): What Can Potentially Be Expected?

Abstract

In the context of energy transition, new and renovated buildings often include heating and/or air conditioning energy-saving technologies based on sustainable energy sources, such as groundwater heat pumps with aquifer thermal energy storage. A new aquifer thermal energy storage system was designed and is under construction in the city of Liège, Belgium, along the Meuse River. This system will be the very first to operate in Wallonia (southern Belgium) and should serve as a reference for future shallow geothermal developments in the region. The targeted alluvial aquifer reservoir was thoroughly characterized using geophysics, pumping tests, and dye and heat tracer tests. A 3D groundwater flow heterogeneous numerical model coupled to heat transport was then developed, automatically calibrated with the state-of-the-art pilot points method, and used for simulating and assessing the future system efficiency. A transient simulation was run over a 25 year-period. The potential thermal impact on the aquifer, based on thermal needs from the future building, was simulated at its full capacity in continuous mode and quantified. While the results show some thermal feedback within the wells of the aquifer thermal energy storage system and heat loss to the aquifer, the thermal affected zone in the aquifer extends up to 980 m downstream of the building and the system efficiency seems suitable for long-term thermal energy production.

1. Introduction

Aquifer thermal energy storage (ATES) is a well known technique used to temporarily store warm and cold water in the subsurface in groundwater wells [1]. In most cases, this technique uses the stored warm or cold water, respectively, for space heating or cooling through seasonal scale storage schemes. Such systems, together with other underground thermal energy storage techniques [2,3], target the partial energy independence of buildings in the global context of the reduction of greenhouse gas emissions when producing energy through the use of renewable energies [4]. Thermal energy storage offers a reliable solution to reduce energy production and provides substantial energy savings [5]. ATES systems generally involve higher financial investments compared to other heat pump systems (due to drilling and associated underground studies) but their average payback times from 2 to 10 years [6,7,8,9] as their total lifetime is expected to range from 25 to 50 years [10,11] and makes them potentially promising investments. ATES systems are found for various types of facilities in the tertiary sector (office, housing, and sport buildings, stores, airports, universities, greenhouses, data centers, etc.) including, for example, Stockholm-Arlanda Airport [12], Carleton University in Ottawa [13], the Klina Hospital in Antwerp [7], and the Reichstag in Berlin [14]. Along with warm or cold water, stored and recovered energy may originate from power-to-heat electricity [15], solar heat [16] or waste heat [17] sources. ATES systems are typically bi-directional systems, with a mono-well or a well doublet allowing for the injection and recovery of warm or cold water, able to run in heating and cooling mode [1,3,18]. The interested reader is referred to Hesaraki et al. [2] and Xu et al. [19] for detailed reviews of subsurface thermal energy storage.

The interest aroused by thermal energy storage systems has increased in recent decades [20]. The first ever thermal energy storage project was implemented in the 1960s in Shanghai, China [3,21]. Since then, ATES applications have extended worldwide, with more than 2800 running systems at the present time [3]. Most of these systems were installed in Northern Europe: about 95% of all systems being in the Netherlands (85%), Sweden, Belgium, and Denmark. While the Netherlands are considered to be in the maturity market phase, Sweden, Belgium, and Denmark are still in the growth market phase and all the other countries are only in the emerging market phase. When taking a closer look at the regional level in Belgium, a striking observation is that all the running ATES systems are located in Flanders and the Brussels-Capital Region [3,22,23]; none of these were present in Wallonia before the Province of Liège decided to undergo the whole authorization process which was recently granted [24]. The development of ATES applications is known to be influenced by local geological and climatic conditions, together with economical, technical, and regulatory conditions, among others [3,25]. In Belgium, no specific regulation regarding the exploitation of geothermal energy has been documented yet, either on a national level or on a regional level, even though it is known that only 5% of shallow geothermal energy in Belgium is used through ATES systems [26]. Local decisions are taken based on the Decree on Environmental Permits from 1985. In Flanders, a decree on deep geothermal use was approved in 2016 and new regulations are expected to be implemented soon [27]. In Wallonia, efforts to develop ATES systems have been made [28,29], with several studies performed on shallow aquifer thermal injection and recovery highlighting some potential ATES target aquifers [30,31,32,33,34,35,36].

Reproducing heat transport originating from ATES systems with numerical models has been successfully achieved by various authors [37]. By using the finite element FEFLOW modeling code [38,39], Bridger and Allen [40] simulated 3D seasonal ATES in a heterogeneous unconfined aquifer, highlighting thermal short-circuit issues, which were reported in several other studies [41,42]. Milnes and Perrochet [43] set up an experiment to assess a well doublet sustainability based on thermal feedback and recycling on an annual timescale. Sheldon et al. [44] developed a catchment scale 3D model to study the impact of a groundwater supercomputer cooling system in an urban aquifer with seasonal storage cycles. Applications of ATES are found in unconfined as well as confined aquifers. In the latter case, a 3D model was recently developed at high temperatures considering the geothermal gradient [45]. Other thermal transport models were developed with other modeling codes, such as HydroGeoSphere [30,34], SEAWAT [46,47], HST3D/HSTWin [48], or METRA [49]. To the authors’ knowledge, most of the ATES numerical applications span over an annual or seasonal timescale with very few available data.

In the present study, extensive datasets were acquired to conceptualize and calibrate a coupled groundwater flow and heat transport model. The ATES system presented here is the very first completed project of its kind in Wallonia. It should hopefully be considered as a cornerstone for future shallow open geothermal systems developments in southern Belgium, since the market is currently in its growth phase [3]. Most of the major Walloon cities are located on top of highly productive alluvial aquifers similar to the one investigated here. The main objectives were (i) to study the feasibility of a future ATES system that will be implemented in the city of Liège, Wallonia (Belgium), (ii) and to understand the targeted aquifer low-temperature geothermal behavior by creating a reliable 3D groundwater flow and heat transport model working at a daily timescale.

2. Materials and Methods

2.1. Study Site

The study site is located in the 80 ha Outremeuse district of the city of Liège, Belgium (Figure 1). A wasteland plot (3.8 ha) was chosen to build a new Knowledge and Resource Centre on behalf of the administrative authority of the Province of Liège. Heating and air conditioning provided by an open-loop well doublet ATES system is planned for the future building. The Outremeuse district is fairly flat, with absolute elevations varying from 63 to 64 m ASL, and is bounded by the Meuse River to the north and the Meuse diversion to the south. Both are channeled and flow from south-west to north-east. In terms of land use, the district is highly urbanized (buildings, roads, sidewalks, etc.), except for the study site, which is currently a wasteland plot but will be urbanized once the ATES system is up and running.

The local hydrogeological setting is typical of an alluvial plain [50]. The top 1.5–5.0 m of the subsurface are composed of backfill soil. Directly below that layer, a Quaternary silty aquitard (2.0 to 5.0 m) confines the sandy gravel Meuse alluvial aquifer (semiconfined) that is about 5.0 to 8.0 m thick. The aquifer is limited underneath by a Carboniferous compact dark shale bedrock. This bedrock is usually considered as possessing aquiclude to aquifer properties since local sandstone layers or lenses and fractured zones might be found.

2.2. Aquifer Characterization through Field Experiments

2.2.1. Pumping Tests

The study site was characterized using geophysical profiles (electrical resistivity and seismic refraction tomographies), well logs, pumping tests, and tracer tests. The piezometers O1 to O5 and the ATES well doublet (W1, W2) were drilled and screened in the aquifer. W1 was screened across the entire aquifer thickness, while W2 was only screened across the lower half of the aquifer. Pumping tests were performed in W2 and monitored with multi-parameter dataloggers in O1 to O5. A step-drawdown pumping test was carried out, producing rates varying from 20 to 65 m³/h. In addition, a constant-rate pumping test set at 50 m³/h in W2 was carried out during 3 days, producing a 0.34 m drawdown in W2 which was reached a few hours after the start of the pumping test. This drawdown remained stable during the 3 days of pumping (Figure 2B). The same trend was observed at the piezometers, reflecting the highly productive nature of the alluvial aquifer, with hydraulic conductivities ranging from 1.5 × 10⁻² to 6.5 × 10⁻³ m/s, which were calculated with the Theis solution under transient-flow regime [51,52]. Those values are in line with previous studies undertaken in the alluvial aquifer of the Meuse River [32,53]. Additionally, the estimated local hydraulic gradient is very low (i = 0.00017), with groundwater gently flowing from west to east, corresponding to an average water flux of 0.09 to 0.21 m/d. Based on these tests, the hydraulic properties of the aquifer seemed suitable for its exploitation through an ATES system [33,54,55,56].

2.2.2. Tracer Tests

Heat and dye tracer tests were carried out to gain insight into the transport parameters of the aquifer. In both cases, the tracers were monitored in W2 where a constant pumping rate of 50 m³/h had been implemented. First, the heat tracer test was carried out by injecting warm water (heated to 48 °C) at a constant injection rate of 2 m³/h in O1 during 14 h. The background temperature was 14 °C, resulting in a ΔT of +34 K. The dye tracer test was initiated by pouring 1 kg of sodium naphtionate diluted in 50 L (concentration = 20 g/L) in O1, and monitored in W2 with a fluorimeter [57]. The recovery rate in W2 was 98%, with a concentration peak of about 1200 ppb. Both breakthrough curves are displayed in Figure 3, assuming t = 0 h as the initial injection time for both tracers. The tracer tests parameters are set out in

Table 1. Heat and dye tracer tests characteristics.

	Heat Tracer	Dye Tracer
Volume/Concentration	28 m3 (+34 K)	50 L (20 g/L)
First arrival	9.2 h	1.7 h
Peak time	31.2 h	3.0 h
Maximum flow velocity	1.8 m/h	9.4 m/h
Peak flow velocity	0.5 m/h	5.3 m/h

.

Due to the 50 m³/h pumping rate implemented in W2, advection was the dominant transport process in the dye tracer test. High groundwater velocities (maximum and peak), together with a short residence time, meant that no adsorption process was involved. However, tailing can be observed on the breakthrough curve (Figure 3B). The fact that the alluvial aquifer is highly heterogeneous suggests that local hydrodynamic dispersion or heterogeneous advection is happening [58,59]. In addition, the first peak is followed by a second lower one (9 h after injection). As for the heat tracer test, very low temperature variations were recorded, as seen on the breakthrough curve in Figure 3A. With a 14.65 °C maximum peak temperature, the monitored ΔT in W2 is only of +0.25 K. This is due to dilution processes and to the difference between the injection flow rate in O1 (2 m³/h) and the pumping flow rate in W2 (50 m³/h); the volume of water pumped from the aquifer was much greater than the one injected. A second lower peak appears around 51 h after injection, which could confirm the hypothesis that hydrodynamic dispersion or heterogeneous advection is happening. In this context, heat conduction in the aquifer seems unlikely to take place as long as W2 is in action; heat convection is likely the dominant process here.

Even if the tested aquifer was the same and the distance between O1 and W2 is short, both tracer tests resulted in significantly different breakthrough curves. The dye tracer arrived at the end in a very short time (1.7 h), much earlier than the heat tracer (9.2 h). Accordingly, the dye tracer peak was recorded before the heat tracer peak, respectively, at 3 h and 32 h. The dye tracer had a clear advective behavior, in which transport was dictated by water flowing through pores, while the heat tracer had a hybrid convective and conductive behavior in which transport was influenced by pores and the matrix of the porous medium. The breakthrough curves were used as references for model calibration, as seen with the dashed lines in Figure 3.

2.3. Numerical Model Parametrization

The hydrogeological model presented here was created and run with FEFLOW [38,39]. A coupled subsurface fully saturated model was set up, with coupled mass transport and heat transport processes. The model is made of 388,703 elements and 236,664 nodes, with inter-nodal distances varying from 30 m (outside the study site perimeter) to 0.05 m (around piezometers and wells of the study site).

The model area is limited by the Meuse River and its diversion (Figure 1B). Both are channeled but the deeper part of their streambed is connected to the bottom of the aquifer (Figure 4). Concrete foundation walls with a 2-m thickness are present and were added to the model along its lateral boundaries. Five layers were created in the model. The first one represents the backfill soil layer and is limited at the top by the ground surface. Layers 2 and 3 form the silty aquitard. The sandy gravel aquifer is divided in two layers (4 and 5).

The parameters presented below were chosen and set based on research in the existing literature and field data collected in situ. The three hydrogeological units seen in Figure 4 were given initial hydraulic conductivity values, as summarized in

Table 2. Hydraulic and transport parameters of the aquifer, with their initial and calibrated values.

Parameter	Initial	Calibrated	Reference
Horizontal hydraulic conductivity [m·s−1]	5.3 × 10−3	2.1 × 10−2 to 3.0 × 10−4	Pumping tests
Vertical hydraulic conductivity [m·s−1]	5.3 × 10−4	2.1 × 10−3 to 3.0 × 10−5	Pumping tests
Specific storage [m−1]	1.2 × 10−3	-	Pumping tests
Effective porosity [-]	0.12	0.16 a	[60]
Matrix heat capacity [MJ·m−3 K−1]	2.87	2.27 b	[61]
Matrix thermal conductivity [J·m−1·s−1·K−1]	2.7	2.7 b	[61]
Longitudinal dispersivity [m]	20	10	[60]
Transverse dispersivity [m]	2	1	[60]

^a Calibration based on the dye tracer test. ^b Calibration based on the heat tracer test.

and

Table 3. Hydraulic and transport parameters of the upper backfill soil and aquitard.

Parameter	Loamy Backfill Soil	Reference	Silty Aquitard	Reference
Horizontal hydraulic conductivity [m·s−1]	1.0 × 10−4	[32]	2.0 × 10−6	[32]
Vertical hydraulic conductivity [m·s−1]	1.0 × 10−5	[32]	2.0 × 10−7	[32]
Specific storage [m−1]	1.0 × 10−4	[32]	1.0 × 10−4	[32]
Effective porosity [-]	0.05	[32]	0.02	[32]
Matrix heat capacity [MJ·m−3 K−1]	3.0	[57]	2.22	[31]
Matrix thermal conductivity [J·m−1·s−1·K−1]	1.9	[57]	1.37	[31]
Longitudinal dispersivity [m]	5	[32]	5	[32]
Transverse dispersivity [m]	0.5	[32]	0.5	[32]

. The alluvial aquifer hydraulic conductivity was applied in line with the pumping test results. For these hydrogeological units, a 1/10 vertical anisotropy factor was fixed for the hydraulic conductivity values. In addition, the concrete foundation walls were given a 5 × 10⁻⁷ m/s hydraulic conductivity value. They were not assumed to be totally impermeable since leakage phenomena may occur between the Meuse and its diversion to the aquifer in both directions.

Transport parameters of the alluvial aquifer were initially set based on parameters estimated from the heat and dye tracer tests. In the loamy aquitard and backfill soil horizons, transport parameters were defined based on research in the existing literature (Table 2 and Table 3).

Reclamation drainage channels are found at the interface between the loamy backfill soil and silty horizons. They run along the lateral boundaries of the model, close to the concrete foundation walls, in order to keep the water table below the district buildings’ foundations. These channels were included in the model and were given free drainage boundary conditions, with Cauchy (third-type) boundary conditions set to take into consideration the elevation of the bottom of the pipes.

The top boundary of the model was given a 30 mm/year mean recharge rate. It was calculated by multiplying the average precipitation rate in Liège (900 mm/year) by a 1/30 factor typical of urbanized sealed surfaces, as proposed by Epting [60] in a similar urban context. In terms of thermal transfer, a Cauchy (third-type) boundary condition was applied at the upper boundary, with daily air temperature data available from Liège airport (50°39′ N, 5°27′ E). Since the unsaturated zone was not included in our model, this upper thermal boundary was applied directly to the water table. Such an assumption may lead to a significant overestimation of the influence of air temperature on groundwater temperature; the thermal transfer coefficient in the first layer of the model was given a value of 1 J/m²/s/K [62].

Even if the bedrock has low hydraulic conductivity values, vertical water flow from the fractured shale bedrock is known to recharge the aquifer in the study area, with a mean value of 100 mm/year [50]. This value was applied to the bottom boundary of the model. No data are available on the thermal flux or the temperature of the vertical water flowing in the aquifer, but a 0.06 W/m² heat flux corresponding to the mean geothermal gradient was accounted for.

Along the lateral boundaries of the model, on the lower part of the aquifer that is connected to the Meuse and its diversion streambeds, Cauchy (third-type) boundary conditions were set to take both streams into consideration. Their water levels were set based on data collected 3 km upstream of the study site (Angleur, 50°36′ N, 5°36′ E) and locally calculated by respecting the river gradient (i = 0.00017).

In terms of heat transport, third-type temperature boundary conditions were applied to the Meuse and its diversion. The data used here were from the Flémalle station (50°35′ N, 5°26′ E). The water temperature shows wide variations that are intrinsically dependent on the seasons, ranging from a few degrees in winter to more than 25 °C in the summer, 28 °C being the maximum authorized temperature for the Meuse River. A nuclear power plant is located about 30 km upstream of Liège that uses water pumped from the river for cooling processes during the energy production. When all units of the power plant produce energy, 3.5 million m³ of surface water are collected and discharged every day in the Meuse. The released water temperature is strictly controlled, with a maximum heating limit fixed to 5 K [63].

As no groundwater abstraction takes place in the Outremeuse district, no pumping well boundary condition was initially implemented in the model. Another point that should be highlighted is that no significant underground heat source is known or documented in the aquifer.

2.4. Model Calibration

2.4.1. Steady State Groundwater Flow Calibration

The model was automatically calibrated in two phases using FEPEST, the FEFLOW integrated version of the model-independent parameter estimation code PEST [64,65]. First, it was calibrated in steady state against hydrogeological data to adjust the hydraulic parameters of the model. The pilot points method was used [66,67]; the resulting horizontal hydraulic conductivity heterogeneous distribution in the lower part of the aquifer is shown in Figure 5.

The observed and simulated hydraulic heads (O1 to O5), together with their calibration residuals, are listed in

Table 4. Observed vs. calibrated hydraulic head data in the observation wells.

	Observed (m)	Simulated (m)	Residual (m)
O1	59.21	59.21	0.00
O2	59.21	59.22	0.01
O3	59.19	59.18	-0.01
O4	59.18	59.16	-0.02
O5	59.23	59.23	-0.01

. These residuals may seem very low (less than 3 cm), but one should bear in mind that the hydraulic gradient in the aquifer is very low. Moreover, it should be pointed out that the dataloggers used for hydraulic head monitoring had a measurement accuracy of ± 0.005 m.

2.4.2. Transient Mass and Heat Transport Calibration

For the second step of the calibration procedure, both mass and heat transport were considered [31]. The tracer tests carried out in the field were reproduced, with a 50 m³/h well boundary condition implemented at W1 in the model. Homogeneous transport parameters of the aquifer were then manually calibrated, taking into consideration zones of piecewise consistency. The initial transport parameters found in the literature were used to set the loamy backfill horizon and silty aquitard. The initial and calibrated transport parameters are listed in Table 3.

Breakthrough curves, both simulated and monitored in situ, are plotted in Figure 3. The simulation of the heat breakthrough curve was deemed adequate, even if the monitored temperature variation was very low, and enabled the calibration of the volumetric heat capacity and thermal conductivity of the aquifer porous matrix (Table 2). The dye tracer peak concentration was appropriately simulated by calibrating the effective porosity of the aquifer. The simulated first arrival, however, came in earlier than the one that was monitored. The tailing effect was correctly reproduced by the model, unlike the second peak, which is clearly visible on the graph. Even if only heat transport was considered in the long-term simulation, calibrating breakthrough curves of both tracers provide complementary information on the transport processes in the aquifer [68].

2.5. Transient Model Validation

The calibrated hydraulic parameters were validated under transient a flow regime by applying a 56 m³/h pumping rate in W2. Observed and simulated hydraulic head data were compared for O3, O5 and W1, as displayed in Figure 6. Drawdowns induced by the active pumping well showed a significant decrease during the first 15 min and then decreased gently over time. This was fitly simulated by the model, but a slight offset between the observed and simulated curves was noticed after 8 hours, the maximum being 0.05 m in W1 and the minimum 0.02 in O5. While the monitored hydraulic head decrease in O5 and W1 was smooth and regular, O3 revealed measurement fluctuations during the experiment. As similar devices were used to monitor the hydraulic head at the three locations, this trend could simply be explained by a technical defect. Since O3 is indeed close to the Meuse diversion, one would expect a stable decrease due to the proximity of such a hydrogeological boundary. No data were available for validating our model with mass and heat transport.

With regard to groundwater flow, the calibration and validation phases were performed against hydraulic head data only. Even if the model area is limited by two streams, no discharge data of the Meuse, its diversion or the reclamation channels were available.

2.6. Simulation Run

The calibrated model was set to run a predictive simulation over a 25-year period based on existing data collected over the last decade, from 1 January 2009 to 31 December 2018, and replicated to cover the 25-year period. The model was initialized with a 2-month warm-up period. Then, the 25-year simulation was run based on time-varying boundary conditions applied along the lateral boundaries of the model, representing the water level and temperature variations in the Meuse River and its diversion. Such data are available online free of charge on the Wallonia Public Service waterways (http://voies-hydrauliques.wallonie.be) and Aquapol project (http://aquapol.environnement.wallonie.be) websites. The initial conditions were set up with streams water level and temperature data from 1 January 2009. In addition, the hydraulic head distribution in the model was given the same value as the streams water level, and the initial groundwater temperature was 14 °C.

In terms of computing performances, the simulation ran for 2 days, 17 h and 53 min on a 64 Go RAM computer with two multi-core processors at 2.20 GHz (20 physical cores, 40 logical processors).

2.6.1. ATES System Operational Characteristics

The future running ATES system was implemented into the model, with W1 being the warm well and W2 the cold well (Figure 1). This system is designed for a maximal operational temperature difference (ΔT) of 5 K between W1 and W2. Daily pumping rate data estimated for the future building needs were applied to both wells, following the trend seen in Figure 7A. Pumping rates during the heating period do not exceed 7.2 m³/h, while they vary up to 25 m³/h during the cooling period, with the highest groundwater extraction rates planned from June to August. The related temperature difference (ΔT) need is displayed in Figure 7B. The air conditioning period during which water is abstracted from the cold well runs from the beginning of March to mid-November, while the space heating period during which water is abstracted from the warm well runs from mid-November to the end of February. As the significant pumping rates differences show, air conditioning needs are much higher than the space heating needs, due to the high efficiency thermal insulation of the future building.

2.6.2. ATES System Thermal Radii

Based on the mean monthly pumping rate and relative temperature data, monthly theoretical thermal radii (R_th) were calculated for both wells when active in their injection phase. The thermal radius is defined as the thermal footprint of an ATES well on its surrounding aquifer in a cylinder representing the volume of water injected [55]. It is calculated based on the concept introduced by Doughty et al. [69], with the following equation:

$$R_{th}=\sqrt{\frac{c_w{V}_{in}}{c_{aq}\pi L}}$$

(1)

with the volume of injected water (V_in), the well screen length (L), and the volumetric heat capacities of the water and the saturated aquifer (c_w and c_aq), the latter being calculated with the porous matrix and the water filling its pores taken into consideration:

$$c_{aq}=n\cdot c_w+(1-n)\cdot c_{matrix}$$

(2)

The resulting mean monthly thermal radii are compiled in

Table 5. Observed vs. calibrated hydraulic head data in the observation wells.

	Warm Well (W1)			Cold Well (W2)
	Vin [m3]	Rth [m]	L/Rth [-]	Vin [m3]	Rth [m]	L/Rth [-]
January	-	-	-	1890	14.5	0.24
February	-	-	-	1430	12.6	0.28
March	585	5.9	1.10	-	-	-
April	3552	14.6	0.44	-	-	-
May	6170	19.3	0.34	-	-	-
June	7588	21.4	0.30	-	-	-
July	11,261	26.0	0.25	-	-	-
August	10,213	24.8	0.26	-	-	-
September	6411	19.6	0.33	-	-	-
October	3819	15.2	0.43	-	-	-
November	1315	8.9	0.73	-	-	-
December	-	-	-	59.23	14.5	0.24

, together with the injected volume of water and the ratio between screen length and thermal radius (L/R_th).

The minimum thermal radius of the warm well is slightly wider than that of the cold well, while the difference in the maximum thermal radius is significant. This is due to the total injected volume of warm water that is considerably larger than that of cold water, as well as to the wells screen length. In W2, the well was screened in the lower part of the aquifer (L = 3 m), where gravel is the coarsest, whereas in W1, the aquifer thickness was fully screened (L = 6.5 m). Doughty et al. [69] stated that the optimal L/R_th ratio is 1.5 to avoid heat loss in the aquifer caused by conduction and dispersion processes. This optimum is, however, encompassed in a 1 to 4 minimum to maximum L/R_th ratio range that is acknowledged to limit conduction and dispersion losses [55]. The low ratios of W1 (0.25 to 1.10) and W2 (0.24 to 0.28) indicate a strong likelihood that loss processes take place in the present configuration. The maximum and minimum thermal radius values for both wells are displayed in Figure 8. Even if groundwater flows roughly in the eastern direction, the warm thermal plume that will develop around W1 could be expected to reach W2 at some point due to the larger volume of warm water injected. As groundwater flow velocities are higher than 25 m/y (i.e., 150 m/y), the recovery efficiency should be mostly affected by ambient groundwater flow advection [54]. In addition, the high injection rates applied in W1 during the heating period should boost the advective flow. Furthermore, the wider thermal radius of W1, linked to its shorter screen length, could lead to heat loss through conduction and dispersion processes in the aquifer, as demonstrated by Bloemendal and Hartog [54].

3. Results

The time span mentioned here is indicative and the related boundary conditions were set based on the data available for the 2009 to 2018 decade, which were replicated to cover the 25-year simulation period. It is important to point out that the ATES system presented here was not in operation during the afore-mentioned period.

3.1. Groundwater Flow

The hydraulic head in the aquifer is dependent on the water level of the streams imposed along the model boundaries. The data prescribed as boundary conditions show seasonal water level variations of 1.5 to 2.0 m on average, with a maximum level of 62.9 m recorded in January of the years 5 and 17 (Figure 9A). The seasonal variation trends are echoed in the ATES system production wells W1 and W2 (Figure 9B) and in the observations wells (Figure 9C), but the corresponding hydraulic head shows variations of a lower intensity.

When comparing both ATES wells, hydraulic head variations were simulated with a difference of up to 0.2 m between W1 and W2 during cooling periods (at the highest pumping rates in July and August). During the heating periods, the simulated hydraulic head in W1 and W2 shows similar values and trends due to low pumping rate values applied, unlike what happens in the cooling phases (Figure 7A).

As for the observation wells, seasonal hydraulic head variation trend and intensity are similar for all of them, regardless of their location and their proximity to the boundary of the Meuse diversion. Differences between the simulated hydraulic head curves of O1 to O5 are hardly distinguishable, even though one would expect the high injection rates applied during the cooling period in W2 to have an impact on O1 because of their proximity. Water is, however, injected in a high conductivity area in W2, unlike in O1, and it flows downgradient to a preferential flow zone.

The eastwards groundwater flow direction is respected in the model throughout the entire simulation, regardless of which air conditioning period or applied pumping and injection rates are considered (Figure 10). Yet, a gradient inversion occurs from an air conditioning period to another with groundwater flowing from W1 to W2 during the cooling period and vice versa during the heating period.

3.2. Groundwater Temperature

The Meuse River is prone to high seasonal temperature variations, as was observed in the course of the simulation. Inter-seasonal variations range up to 25 K within a year, with a maximum recorded temperature of 28 °C summer and a minimum of 2 °C in winter (Figure 11A).

The daily temperature differences in Figure 7B are consistent throughout the 25 years simulations, with ATES warm and cold wells (W1 and W2) abstraction periods clearly visible (Figure 11B). In contrast to the simulated hydraulic head values (Figure 9B), the simulated temperatures in both wells are significantly different. Overall, the temperature rises up to 20.2 °C in W1, with an inter-seasonal variation of a few K. In W2, temperature fluctuates between 14 and 15 °C during the cooling period, and it varies from slightly below 14 to 18 °C during the heating period. This is due to the applied temperature difference between both wells that show high variations during the heating period, while the applied temperature difference is generally steady throughout the cooling period. The groundwater temperature increases during the first three operational years and then reaches a plateau where the rise in groundwater temperature is limited to 20.2 °C. The first two months were considered as a model warm-up period. Thermal feedback occurs here, as seen in Figure 10: a thermal plume is generated in the aquifer and keeps on extending as the ATES system is running. It is intercepted by the cold well (W2) as soon as the first cooling period is run. Indeed, this well has not recorded any temperature lower than the initial background aquifer temperature: after the warm-up year, the simulated temperatures in W2 were mostly close to and higher than 14 °C.

When considering the observation wells (Figure 11C), O1 shows distinct temperature variations following the variations in W1 located nearby. No temperature variations are noted in O2 and O5, demonstrating that the temperature variations in the Meuse and its diversion do not influence the overall aquifer groundwater temperature. This is supported by the temperature distribution maps in Figure 10. The ATES system should not, therefore, be affected by the temperature variations of the stream water. From year 3 onwards, O4 was showing a constant temperature rise until it reached a plateau in year 6 because it is located downstream of the ATES system and the thermal plume reached O4. The calculated temperature in O3 rose by about 0.5 K during the first years and followed the seasonal trends.

The thermal plume was defined based on the thermal affected zone showing a temperature rise higher than 0.5 K (Figure 10), as investigated by Attard et al. [70]. After the warm-up period and the first ATES cycle, the thermal plume had already reached W2; the thermal feedback effect was already on. The plume expanded in the east-northeast direction during the simulation, reaching a 980 m longitudinal extension and a 175 m transverse extension at the end of year 25. Its expansion is mainly due to heat convection processes induced by the ATES system that is continuously in action and by the high ambient groundwater flow. In year 5, the thermal plume reaches the Meuse diversion boundary downstream of the ATES system, as displayed for years 15 and 25 in Figure 10. Groundwater exfiltration to the stream, together with heat transfer processes, occurred during the simulation.

4. Discussion

The model presented here was built based on field data collected through conventional hydrogeological tests (pumping tests, heat and dye tracer tests). An inverse calibration was considered for the hydraulic conductivity in the model, which is a procedure known to produce non-unique parameters distributions [71,72]. The resulting hydraulic conductivity distribution found in the aquifer simulated in this study leads to specific calculations of the predicted hydraulic head and temperature, representing one solution among others that bears inherent uncertainties, especially since data scarcity limits the accuracy of the calibration procedure. Data scarcity could be an obstacle to proper predictions, as is the hydraulic head measurement error of ±0.005 m inherent to the dataloggers used in this study. In a context with a very low ambient hydraulic gradient, such an error related to data used for calibration purposes could mislead the well informed modeler and produce uncertain calibrated parameters distributions. While the aim of this paper was to provide a thoroughly documented ATES case study, no uncertainty analysis was performed. It is believed that the hydraulic conductivity of the aquifer and its variance are the most sensitive parameters of the model, as the sensitivity analysis results of a study undertaken in the same aquifer (but at another location) have shown. This study was carried out with similar parameters and boundary conditions, by performing and simulating a heat push-pull test [34]. Inverse modeling has been widely applied in previous research projects [31,73] and recent developments of Bayesian evidential learning techniques push forward the allowance for direct predictions to be performed based on prior models [34]. This technique, or even inverse stochastic modeling with particle tracking [74], could shed light on the model weaknesses, and help identify its uncertainties and substantially enhance its predictive strengths, allowing to make fully informed decisions.

The calibration was performed automatically for groundwater flow and manually for the heat transport parameters. It is generally good practice to calibrate groundwater flow and heat transport parameters together in order to reduce uncertainties related to ill-posed models, especially with models developed in heterogeneous hydrogeological environments [64,75,76]. Due to its highly time-consuming aspect, performing a transient groundwater flow and heat transfer model calibration was hardly feasible; but not doing so bears inherent uncertainties [64]. Nevertheless, we chose to proceed because the calibrated parameters would have been representative of the porous medium located between O1 and W2 only, which represents an area covering a distance of 16 m only, i.e., the inter-well distance. The dye and heat tracer tests presented here were performed between O1 and W2, and no other transport data were available elsewhere in the aquifer investigated. Data scarcity is obviously tied to the budget allocated to the present study. A combined flow and transport calibration would have provided a more robust model but with a level of detail that would not likely be significant compared to the scale of the study area. Applying homogeneous transport parameters to the model domain based on local scale calibrated parameters between O1 and W2 imposes the unique character of these parameters and details on heterogeneous materials are lost [71,76,77]. Extensive datasets collected in all the wells by carrying out several heat and solute tracer tests would have helped to better characterize the aquifer [68].

The simulation was run before the ATES system was up and running. Results could not be stressed based on actual data storage and recovery data, which would be of particular interest once several seasonal ATES cycles are performed. They provide, instead, insight into a worst-case scenario of potential groundwater temperature distributions. This worst-case scenario was chosen as the ATES system simulated here was run for 25 years in a row in continuous mode at its full capacity, while such systems are usually run in cyclic mode [3]. A significant imbalance between stored warm and cold water is therefore observed, namely 5100 m³ of cold water stored during the heating period and 50,000 m³ of warm water stored during the cooling period. It should, however, be pointed out that the building for which the system was designed is an office building. According to the building needs, the system is expected to be in operation mostly during office hours. Late hours and week-ends should be storage phases with no groundwater abstraction, and no heat or cold recovery, accordingly. The thermal feedback effect highlighted here is driven by thermal interference between both wells and could be avoided by respecting a well-to-well distance and a thermal radius ratio of 2/1 [55] or 3/1 [46]. It is expected to be of lower amplitude and should not have such an impact on the aquifer. If this were to be the case after all, since both ATES wells are equipped with pumps and used for injection or abstraction, reversing the ATES system could be an option to limit the increase in groundwater temperature. Another solution would be to use the stored energy to heat other buildings of the Outremeuse district, or to store cold water from an external source so that warm and cold stored water volumes would be balanced [55]. According to Bloemendal and Olsthoorn [78], the use of recirculation systems in aquifers with high groundwater flow velocities with well-doublets of the same kind (warm or cold) helps in preventing heat loss to the aquifer. In the present case, the use of at least two cold wells and two warm wells should be considered. Such options would minimize the impact of the ATES system on groundwater while at the same time improve the energy efficiency of potential new systems implemented in the aquifer.

Even if the Meuse River and its diversion are channeled, their streambeds are connected to the aquifer. No data or indication regarding the losing or gaining parts of the streams could be found. Regarding the Meuse diversion, the piezometric data clearly show its draining effect on the aquifer in the vicinity of the ATES system, east of the study site. The western boundary of the aquifer is therefore assumed to be a losing portion of the Meuse River, a part of it recharging the aquifer. Third-type (Cauchy) boundary conditions were set along the aquifer lateral boundaries by applying daily water level variations of the Meuse River and its diversion. The transfer rate associated with these boundary conditions was not calibrated since no stream discharge data were available; it was estimated instead. Being able to calibrate this parameter would help to better understand the interactions between the alluvial aquifer and the streams, because they have a strong influence on the simulated ATES system production temperature in terms of groundwater flow.

Civil works and structures were previously shown to affect the Meuse river and the alluvial aquifer [79]. Such issues are not specific to the city of Liège and were assessed in previous studies [80,81,82]. Quay walls and tunnels are present along the lateral boundaries of the model but their exact location and depth are unknown and very little is known on their impact on local hydrogeological conditions. Having detailed information on low-conductive structures should help to better constrain the model inflow and outflow at the boundaries.

The upwards recharge from the bedrock is known to occur along vertical fractures in low-conductive matrix [50]. The exact location of these fractures is not known in the model area. Groundwater flowing from the bedrock to the aquifer potentially has a lower temperature than the alluvial aquifer (14 °C). It could locally decrease the background aquifer temperature, interfere with the running ATES system and affect its energy recovery efficiency. It should be noted here that no specific heat-flux boundary condition was applied at the bottom boundary of the model; the geogenic geothermal gradient was not accounted for either.

In alluvial aquifers, the groundwater temperature may show significant seasonal variations due to surface water–groundwater interactions [83] or to the seasonal thermal fluctuation zone of the subsurface profile [84]. The initial temperature set in our model (14 °C) was based on temperatures measured at the end of summer in the aquifer. Since our simulations started in winter, an artificial temperature offset may have been induced in our model that may have influenced the final results, even if a warm-up ATES heating/cooling cycle was considered. An experimental site implemented in the same hydrogeological unit has been extensively investigated, inter alia, with heat tracer tests [31,32,35,36,84]. It is located in Hermalle-sous-Argenteau, 13 km downstream of the studied site along the Meuse, in a suburban area. The average groundwater temperature in the underlying aquifer is 13 °C, with measured seasonal variations ranging from 10.5 to 13.5 °C. Land cover was demonstrated to influence groundwater temperature, inducing a temperature rise of 2.0 ± 0.7 K under artificial surfaces, as found in urban environment [85], while the temperature rise linked to natural terrestrial vegetation was limited to 0.2 ± 0.8 K.

The initial high temperature in Liège is probably due to the urban heat island effect that was assessed in many cities worldwide [86,87,88,89,90], especially since our model area has one of the highest building densities in the city. Arola et al. [91] reported elevated groundwater temperatures due to urbanization, combined with shallow seasonal thermal fluctuation zone. Underground car parks, multi-story basements, road tunnels or district heating pipes could be potential heat sources that affect the aquifer temperature [92]. Several short-distance tunnels, partially burrowed in the aquifer, are found along the lateral boundaries of the model area. Their short distance (<20 m) implies that the temperature in these tunnels is likely in equilibrium with air temperature. No data are available about the depth of car parks or buildings basements found in the Outremeuse district, and their potential penetration depth in the aquifer is not known either. Improving the model presented here with heat flux estimations from buildings based on other studies could help [70,93,94]. In unconfined aquifers, heat flux through the unsaturated zone to the saturated zone originates from anthropogenic sources and structures resulting in a significant heat transfer processes to the aquifer [83,92]. Groundwater flow and heat transport simulations were previously performed in the semi-confined alluvial aquifer in Liège by neglecting the effect of the temperature sources of urban structures [30].

The temperature of the two surface water bodies was prescribed as a third-type boundary condition with data collected at a shallow depth in the Meuse River. The temperature was assumed to be constant along the streams vertical profile. Yet, one could argue that their shallow part is more likely to be affected by seasonal temperature variations than their deeper part [95]. The simulations run here are considered worst-case scenarios since only the bottom part of the river is connected to the aquifer through its streambed, with temperature data collected close to the river surface affected by variations in seasonal temperature.

5. Conclusions

We presented a new aquifer thermal energy storage system that will be implemented in the Meuse alluvial aquifer in Wallonia (southern Belgium). A transient simulation was run over a 25-year period, and potential thermal impact on the aquifer was assessed. In its present state, the ATES system shows thermal feedback between the pumping and production wells due to cooling periods which are longer than heating periods, inducing larger storage volumes of warm water than cold water. The groundwater temperature increases during the first operational years and is then stabilized, limiting the groundwater temperature rise to 20.2 °C, and it is also higher than the initial aquifer temperature (14 °C). Water level fluctuations in the Meuse River were demonstrated to influence the simulated temperature at both ATES wells, but they were not affected by the variations in stream temperature.

The future operational system will likely be designed with the aim to prevent the thermal feedback effect, by considering the installation of recirculation systems with cold and warm wells pairs, or by injecting cold water from an external source during the heating period so that the stored volumes of warm and cold water are at equilibrium over a year. Reaching a balance between the stored volumes of cold and warm water is crucial to ensure long-term sustainable thermal energy production. Sharing the groundwater thermal resources for space heating or air conditioning is also considered to be a reliable solution to limit the impact of the temperature rise on local groundwater resources.

This case study opens opportunities for further shallow open geothermal developments in Wallonia, as the promising market of the region is still in its emerging phase. Most of the major Walloon cities are located on top of highly productive alluvial aquifers similar to the one investigated here. These aquifers are likely to be suitable reservoirs for ATES systems implementation and new systems are expected to spring up in the next decades.