Modeling a large ultra-low-temperature district heating system with a focus on an aquifer thermal energy storage

After a hesitation in the development of shallow geothermal energy systems in Belgium between 2014 and 2017, a clear revival of low-temperature district heating (LTDH) systems and thermal energy storage (TES) is now observed [1]. LTDH systems typically belong to the concept of 5th-generation district heating and cooling (5GDHC) systems [[2], [3], [4]].
Shallow geothermal energy systems have emerged as a promising solution for achieving sustainable and low-carbon heating and cooling in urban environments. There are 2 main types of geothermal TES applicable for seasonal usage: borehole TES (BTES) and aquifer TES (ATES). These systems exploit the stable subsurface temperature to store and recover thermal energy seasonally, offering a means to balance fluctuating energy demands while reducing reliance on fossil fuels. ATES has received significant attention due to its high storage capacity and potential for integration into district energy networks. ATES systems operate by injecting and extracting water at different temperatures to store excess thermal energy during periods of low demand and recover it when heating or cooling needs arise, enabling an efficient and renewable form of energy management.
BTES is mainly used for single-family houses or larger buildings where the hydrogeological conditions are unfavorable when placing an ATES. BTES consists of a borehole through which a plastic pipe is placed and a carrier fluid is pumped. Since the loop fluid is not in contact with the groundwater, the exchange of thermal energy (heat or cold) between the fluid and the surrounding soil occurs via conduction. Unlike BTES, ATES is an open-circuit system that directly uses aquifer water as a heat medium. In Flanders (the northern part of Belgium), its application could be traced to the regulations, which introduced energy performance labels and the condition to generate at least 15 kWh/m2 of renewable energy to receive governmental subsidies. This has also been promoted in Brussels’ capital region, where passive building and energy consumption of less than 15 kWh/m2 have been requirements since 2015.
The evolution of district heating technologies provides an important context for understanding the role of ATES within modern energy systems. Traditional third- and fourth-generation district heating (3GDH and 4GDH) networks are characterized by centralized heat production and relatively high supply temperatures, which limit their ability to integrate low-grade renewable sources. In contrast, 5GDHC systems operate at near-ambient temperatures and rely on decentralized, bidirectional thermal energy exchange. This paradigm shift allows for simultaneous heating and cooling, the direct use of waste and renewable heat, and stronger coupling with electrical and control systems [5]. However, it also introduces new challenges, including the need for advanced hydraulic balancing, management of temperature stratification in thermal storage, and optimization of interactions between aquifer systems, energy management systems (EMS), and end-user demands.
Given these developments, it is essential to raise several key research questions that define the motivation for this study. How can 5GDHC systems efficiently exploit aquifer-based thermal storage under spatially variable and heterogeneous subsurface conditions? What modeling strategies can capture the coupled hydraulic and thermal dynamics of multi-well configurations while remaining computationally feasible for large-scale network analysis? And how can the integration of ATES with EMS improve system flexibility, efficiency, and economic performance? Addressing these questions is critical for advancing the understanding and design of next-generation district energy systems that fully utilize shallow geothermal resources.
Previous research has primarily focused on single-well or small-scale ATES models, which provide valuable theoretical insights but often fail to capture the complexity of practical urban-scale 5GDHC networks. Real systems involve multiple interacting wells, varying geological layers, and dynamic control mechanisms governed by EMS. Consequently, there remains a clear need for comprehensive modeling tools that bridge the gap between simplified hydrothermal simulations and the operational realities of large-scale 5GDHC deployment. This study responds to that gap by developing a computational framework that links the thermal and hydraulic behavior of ATES wells with system-level energy management, offering a scalable approach for evaluating the performance and design of large 5GDHC systems.
This system-level approach captures the collective performance of multiple subsystems operating under diverse conditions. The resulting figures will present aggregated trends for the 5GDHC network, illustrating how recovery efficiency and heat exchange evolve at scale. Compared with previous analysis, this multi-well modeling approach highlights the emergent behavior of low-temperature district systems and supports practical upscaling to neighborhood or district levels.
These relationships quantify the dynamic balance between covering imbalances (e.g., via Combined Heat and Power, CHP) and storage (via EMS) components. Incorporating these metrics allows the model to link subsurface heat recovery dynamics with above-ground energy management strategies. The following literature review will support these statements.

Shi et al. [6] assess the heat losses as a main component of their paper, which can represent aquitards’ heat preservation ability. However, it doesn’t contribute to identifying promising locations for ATES systems and places of ATES wells. They also conclude the heat of aquitards is possibly collected by other technologies and methods, although its extraction is technically and financially unfeasible. In [7], the 3D hydrogeological model for an aquifer is developed. It includes 3 layers, 98 columns, 36 rows, and a grid size of 200 × 200 m, based on the geological and hydrological conditions. The difference is the model is primarily developed to track particles and their distribution to compare with operational data. Thermal properties such as longitudinal dispersion, diffusion, and sorption are simply evaluated by parametrical estimation. Schmidt et al. [8] work with the ATES concept but mostly compare it with other types. They conclude that the construction cost of the storage concepts varies significantly, although all types show a considerable effect on the economy of scale, which means the investments decrease with increasing TES volume. As an example of a cost-effective solution, they list large-scale pit TESs constructed in Denmark in 2010–2020 with specific investment costs of 20–40 EUR/m3.
The difference with [9] is a simplified model where the temperature linearly increases in warm and cold wells in an annual manner. That might be true for years when average seasonal outdoor temperature consistently increased (e.g., between 2014 and 2016). However, colder ambient temperatures would not increase the temperature of the cold well; they would either have no effect or a cooling effect. On the other hand, in summer, the temperatures increase almost every year, which validates the scenario we used in this paper. To study the possible interaction between two ATES systems, De Paoli et al. [10] also used ground-water and heat transport numerical models. Compared to the present paper, they simulate much larger ATES systems, which include 4 cold and 4 warm wells each. The third system, which they report, is even larger, with 5 cold and 5 warm wells extracting water from the deep Paleozoic layer.
Possemiers et al. [11] consider the groundwater chemistry next to 7 ATES systems in Flanders. However, the focus is on the risks exposed to the phreatic aquifers, which are less protected against contamination. The difference is that they consider the quality of the water pumped in a nearby public drinking water supply well field, especially when the well screens of the drinking water wells are situated deeper than the screens of the ATES wells. In [12], a cut out of the ATES plan maps in Delft and Amstelveen (the Netherlands) is reported. It indicates search areas for warm and cold wells and existing well locations. In addition to the case study area, another difference is the ratio between hydraulic and thermal radius, which limits the propagation of a hot plume as the square root of the thermal retardation.
Salenbien et al. [13] depict a simplified district heating network diagram showing a central heating source, possibly represented by a circle, with connections to consumers, represented by squares. It showcases a potentially cost-effective district heating network design. Combining a CHP plant and underground mine storage offers flexibility and potentially lower operational costs. We further developed this research and used the presented indicators to compare 4GDH and 5GDHC systems.
Perera et al. [14] investigated the efficiency of heat recovery in ATES systems. Their work focused on modeling the overlap area between the injected hot fluid and the extracted fluid within the aquifer. This overlap area is crucial for determining the amount of heat that can be recovered from the system. Perera et al. [14] run simulations and analytical solutions based on these equations to understand how factors like storage period, flow rate, and aquifer properties influence the overlap area and the recovery efficiency of the ATES system.
Every ATES has at least 2 wells: one for the injection/ejection of the cooled liquid and another for the warm. Wellbore placement optimization is the focus of Beernink et al.’s work [15]. They calculate the thermal resistance of the borehole itself. It considers factors like the volumetric heat capacity of water within the borehole, the volumetric heat capacity of the aquifer material, and the borehole’s geometry (volume and length). A higher thermal resistance indicates less efficient heat transfer between the borehole and the surrounding aquifer.
After that, they define the horizontal and vertical distances from the center of the wellbore. The parameters are the wellbore diameter and a reduced thermal resistance term. These reduced thermal resistances might represent the influence of the surrounding aquifer formation on heat transfer in the horizontal and vertical directions. By analyzing these distances, Beernink et al. [15] investigate the impact of wellbore placement on heat exchange efficiency within the ATES system.
Although Lanahan et al. [16] focus on the same topic of avoiding excessive heat losses; they optimize the ratio between maximum volume and surface area by changing borehole depth. The capacity of BTES may differ, ranging from 2 pipes for a single-family house up to more than 500 pipes for a residential district. This paper considers the number of boreholes within the state’s range, although the more typical number of 60 is considered. For this configuration, an acceptable level of heat loss losses was found; therefore, the same with [17], grout with a higher thermal conductivity around U-pipes may significantly improve the thermal behavior of BTES to achieve an efficiency of 40–50 %. To compare, one of the most successful BTES systems in operation is situated in Alberta, Canada. Its annual efficiency for the 3rd to 5th years of operation was 42 % [18].
In addition to the measured data, available simulation tools represent the actual performance quite well with an acceptable computational time. For instance, Puttige et al. [19] present indicators for CompPower software compared with the test data. CompPower software is applicable to simulate heat pump and TES operation. The mean average errors for a heat pump and borehole heat exchanger were 7.8 and 19.1 %, respectively, while the computational time for simulating a 4-year operation was approximately 9 h. As an example of robust multi-objective optimization, the one based on the multi-objective evolutionary algorithm is presented in [20].
Heat transfer to the 250 m-deep BTES is presented in [21]. The difference with the presented approach is the varying volumetric flow rate of water to change the amount of heat transferred to the BTES. Haq et al. [21] fix the flow rate at 1 l/s. Since BTES is a dynamic system, heat transfer intensity heavily depends on time and temperature differences. In [17], the temperature difference between the enhanced grout and the others becomes significant only after 60 h of heat injection.
In [22], the design parameters and simulation results of the BTES system in Chifeng (China) were set as a benchmark for the system parametric study. For 17 factors with fractional design (fractional coefficient equal to 9) and full factorial design on 5 levels, Wołoszyn et al. [23] carried out 109 numerical experiments. The same with [22], the following design parameters were defined as input variables: the thermal conductivity of the soil, the specific heat capacity of the soil, the thermal conductivity of the grout, the total borehole depth, and the height-to-diameter ratio of the system. In [21], the parametric study demonstrates the impact of volumetric flow rate on heat transfer to the BTES. Another example of sensitivity analysis is given in [18], where the main indicator was the BTES efficiency. Unlike [18,21,22], we did not iterate the spacing of boreholes, the number of series-connected boreholes, or the inner circulation flow rate per borehole.
For the short and medium term, Skarphagen et al. [24] suggest a simulation approach based on the ‘line source’ heat equation, first suggested by Whitehead in 1927 and later applied by Clarence Lubin and Charles Theis to the problem of groundwater heat transfer. To compute the Nusselt number in [25], the Gnielinski correlation with the friction correlation by Pethukov is applied.
To the best of the author’s knowledge, models to simulate a large 5GDHC system bringing energy engineering and geological aspects together are rare and not well-established. In contrast, such a model is a building block for understanding and managing heat transfer within a system. By incorporating them into a more comprehensive model, we aim to decrease energy usage further while maintaining comfortable conditions and ensuring the efficient operation of HVAC systems [26].