City-scale heating and cooling with aquifer thermal energy storage (ATES)

Sustainable and climate-friendly space heating and cooling is of great importance for the energy transition. Compared to conventional energy sources, Aquifer Thermal Energy Storage (ATES) systems can significantly reduce greenhouse gas emissions from space heating and cooling. Hence, the objective of this study is to quantify the technical potential of shallow low-temperature ATES systems in terms of reclaim-able energy in the city of Freiburg im Breisgau, Germany. Based on 3D heat transport modeling, heating and cooling power densities are determined for different ATES configurations located in an unconsolidated gravel aquifer of varying hydrogeological subsurface characteristics. High groundwater flow velocities of up to 13 m d −1 cause high storage energy loss and thus limit power densities to a maximum of 3.2 W m −2 . Nevertheless, comparison of these power densities with the existing thermal energy demands shows that ATES systems can achieve substantial heating and cooling supply rates. This is especially true for the cooling demand, for which a full supply by ATES is determined for 92% of all residential buildings in the study area. For ATES heating alone, potential greenhouse gas emission savings of up to about 70,000 tCO 2 eq a −1 are calculated, which equals about 40% of the current greenhouse gas emissions caused by space and water heating in the study areas’ residential building stock. The modeling approach proposed in this study can also be applied in other regions with similar hydrogeological conditions to obtain estimations of local ATES supply rates and support city-scale energy planning.


Introduction
One of the overarching goals formulated at the COP 26 climate change conference in 2021 is the global net zero emission of greenhouse gases by mid-century to limit climate warming to 1.5 °C compared to pre-industrial levels (COP26 2021).Besides international and national policies, climate change protection is also driven forward at the city and municipal level as many cities are developing sustainable and climate-friendly energy planning concepts (Epting et al. 2020;Lim et al. 2019, Pulselli et al. 2021).Integrated urban energy planning strategies are a suitable tool to achieve municipal climate protection plans and to reduce greenhouse gas (GHG) emissions.Since space heating and cooling in the building sector alone make up more than 30% of Germany's final energy consumption (AGEB 2021), climate-friendly heating and cooling solutions are of great importance.
One alternative to conventional space heating and cooling based on fossil fuels and cooling machines, respectively, are geothermal applications using the shallow subsurface and groundwater as a renewable source of thermal energy, such as groundwater heat pump (GWHP) systems.By seasonally reversing the pumping direction of GWHP systems, groundwater can also serve as a seasonal storage medium for heated and cooled water.This type of shallow geothermal energy is known as aquifer thermal energy storage (ATES) and allows to reduce seasonal mismatches between demand and availability of thermal energy by storing waste heat in summer and excess cooling capacities in winter (Bakr et al. 2013;Dickinson et al. 2009;Fleuchaus et al. 2018;Schüppler et al. 2019).This results in a more efficient operation of the heat pump system.For cooling, it is often possible to utilize the cooled water directly without any heat pump operation (Banks 2009;Bloemendal et al. 2018;Fleuchaus et al. 2018;Sommer et al. 2014).
A quantification approach revealing the potentially achievable heating and cooling power can facilitate the integration of ATES into city-scale energy planning (e.g.Bayer et al. 2019).Thus, this study uses the concept of power density, which relates the amount of power generated by a specific technology to the required horizontal Earth surface area to quantify the ATES potential.In recent years, the unit of power density has been increasingly used to highlight space requirements as a potential limiting factor for the transition to renewable energies, which typically have much lower power densities than fossil fuels or nuclear energy (Kammen and Sunter 2016;Smil 2015).Power density can also serve as a universal mean for comparing different electricity generation technologies (Kammen and Sunter 2016;van Zalk and Behrens 2018).For instance, according to van Zalk and Behrens (2018), nuclear power plants and coal-fired power plants have median power densities of around 241 W m −2 and 135 W m −2 , respectively.Lower median power density values of around 7 W m −2 and 2 W m −2 are stated for solar and wind power, respectively.
With regard to thermal energy, the power density concept was previously used to quantify the technical potential of shallow geothermal applications (Bayer et al. 2019;Kammen and Sunter 2016).The technical potential, as referred to in this study, relates to a specific extraction technology, such as open GWHP or closed ground source heat pump (GSHP) systems.It is, therefore, constrained by technical factors, such as space restrictions and temperature limits (Hähnlein et al. 2013;Tissen et al. 2019;2021).Bayer et al. (2019) compiled an overview of relevant studies from the literature and revealed a wide range of normalized power densities for GSHP systems with values from less than 10 W m −2 up to more than 400 W m −2 .These discrepancies result from a variety of underlying assumptions and approaches.
Other studies regarding the technical potential of shallow geothermal energy often do not explicitly use the term power density while having similar objectives.For example, Tissen et al. (2019;2021) used guideline values of achievable energy extraction rates to quantify the potential of closed geothermal systems with respect to determined space requirements on the district-and city-scale.Other studies present similar quantitative calculation approaches to estimate the thermal potential of groundwater and compare it to the energy demand as a means for subsurface thermal planning in urban settings (Epting et al. 2018;Miocic and Krecher 2022;Zhu et al. 2010).Böttcher et al. (2019) and Epting et al. (2020), for example, determined the technical geothermal potential of open GWHP systems based on 2D numerical box models considering groundwater flow conditions and different pumping rates as well as temperature changes of the extracted groundwater.However, their box models only considered hydraulic effects of GWHP systems on the aquifer, while plume propagation of thermal anomalies was not considered.For the determination of meaningful space requirements of open geothermal installations, however, modeling the thermal plume propagation is crucial to avoid thermal interactions.
In this study, we develop a novel methodology to assess the technical potential of low-temperature aquifer thermal energy storage (LT-ATES), which is commonly characterized by storage temperatures between 5 °C and 25 °C (Fleuchaus et al. 2018).To this end and for the first time, heating and cooling power densities of LT-ATES are determined considering advective heat transport and multiple adapted ATES well configurations in dependence on the local ambient groundwater flow velocity.For the power density determination, thermo-hydraulic 3D numerical box models are created using the German city of Freiburg im Breisgau (hereafter referred to as Freiburg) as an exemplary application region.These models simplify the modeling process compared to comprehensive city-scale models using a simplified geometry and representative hydrogeological and thermal underground characteristics.The 3D box models are checked against the city-scale model of the study area to evaluate the box models' representativeness.Their simple design and short simulation runtimes make the box models suitable for potential future applications in other study areas.Comparing the obtained power density values from the box models to the existing heating and cooling demand in the city of Freiburg then allows estimating heating and cooling supply rates that could be realized by ATES applications.Furthermore, this study compares GHG emissions of the potential ATES application in the city of Freiburg to those from conventional technologies, which are currently in operation.

Study area
The study area is located within the municipality of Freiburg in Southwest Germany with a population of about 230,000.The city of Freiburg has been regarded as a 'green city' role model for more than three decades due to the city's efforts to promote ecological urbanization, environmental policies and high quality of life (Fastenrath and Braun 2018;Medearis and Daseking 2012;Rohracher and Späth 2014).It covers a total area of about 153 km 2 and is located at the transition of the Upper Rhine Graben (URG) to the Black Forest mountain range (Fig. 1a).The upper graben fill sediments consisting of Pliocene and Quaternary gravel deposits from the Alps and uplifted rift flanks form productive porous aquifers, and thus provide a major portion of the regional drinking water supply as well as industrial and irrigation water demands (Geyer and Gwinner 2011;Villinger 1999).In the area of Freiburg, the upper graben fill developed as an alluvial fan of the Dreisam River (Fig. 1b) and by glacial meltwater from the Black Forest.The shallow aquifer in this area consists of two unconfined groundwater bodies, which are hydraulically connected (Fig. 1c).
The upper groundwater body in the study area is formed by the Neuenburg Formation consisting of predominantly unweathered and loosely bedded gravels with varying sand and low silt contents.The underlying Breisgau Formation consists of partially weathered sandy-silty gravels and has a lower hydraulic conductivity (Table 1).However, there is no distinct transition between these two Pleistocene formations which have a combined thickness of mostly less than 100 m (Geyer and Gwinner 2011; LUBW 2006; Villinger 1999; Wirsing and Luz 2005).In accordance with the hydraulic head contour lines in Fig. 1b, the direction of regional groundwater flow is Northwest towards the river Rhine.The contour lines are interpolated from a total of 118 groundwater monitoring wells using each well's five-year mean hydraulic head.
The two areas in Fig. 1b marked as hydraulically isolated are known as 'Lehener Bergle' (in the north) and 'Honigbuck' (in the south).They are tectonic horst structures that remained at the surface during subsidence of the surrounding rift system.These Mesozoic sedimentary rocks are hydraulically not connected to the Pliocene and Quaternary sand and gravel deposits (Villinger 1999).

City-scale model
A city-scale numerical 3D finite element method (FEM) subsurface model of the Freiburg study area is built in COMSOL.Details on the thermo-hydraulic numerical modeling in this study including the required basic equations are given in the Additional file 1 (section SD1).The Freiburg subsurface flow and heat transport model discretized by about 1.5 million tetrahedral elements covers an area of about 72 km 2 and has a vertical extent of about 290 m (Fig. 1).In this study, the city-scale model serves as a baseline benchmark for the evaluation of the box models' representativeness when determining the power density in the city of Freiburg.
Transverse dispersivity 1 m Baden-Württemberg (2009), Beims (1983) Fluid heat capacity ratio 1 -COMSOL (2020) Table 1 shows the hydraulic and thermal parameters assigned to the city-scale model including the hydraulic conductivity of the Neuenburg Formation, which is implemented as a spatially varying parameter and used for model calibration.Section SD2 in the Additional file 1 provides further information regarding the subsurface model's geometry and boundary conditions (BCs), as well as its calibration.
The characteristics of longitudinal and transverse dispersivities in transport phenomena are a field of active and current research (Di Dato et al. 2022;Park and Lee 2021;Pophillat et al. 2020b;Younes et al. 2020).For this study focusing on the development of a novel methodology for city-scale assessment of the technical ATES potential, we assume commonly used thermal dispersivities (Table 1).

Model geometry and hydrogeological subsurface data
To simplify and speed up the modeling process aimed at determining the power density of ATES systems within the study area, we utilize simplified numerical 3D finite element box models based on the complex city-scale subsurface model.Based on the spatial distribution of the hydraulic conductivity and hydraulic gradient in the calibrated city-scale subsurface flow model, the study area is divided into four homogeneous hydrogeological regions with the aim of minimizing differences of these parameters within a single region (Fig. 2, Table 2).These two parameters, the multiplication of which results in the Darcy velocity, are chosen for the delineation since the groundwater velocity highly influences thermal plume spreading (Piga et al. 2017;Pophillat et al. 2020a;2020b).
Due to very high groundwater flow velocities (29.1 m d −1 as calculated from horizontal hydraulic conductivity, hydraulic gradient and porosity) for the hydrogeological region 4 and the anticipated detrimental influence of the river Dreisam, ATES applications are assumed to be not feasible in region 4. Accordingly, no box models are created for that region.For regions 1-3, numerical box models are generated as parallelepipeds with a length of 3500 m and a width of 600 m (Fig. 3).They consist of two layers, the upper of which represents the Neuenburg Formation with a uniform thickness of 20 m.The lower layer represents the Breisgau Formation, which is implemented with a uniform thickness of 50 m.The hydraulic conductivities for both formations are set to the region-specific values given in Table 2.The slopes of the box models' surfaces and layer boundaries correspond to the respective region's representative hydraulic gradient (Table 2), while accounting for the assumption of a uniform groundwater table depth of 3 m throughout  the box models.The hydraulic gradient along a model's longitudinal extent is implicitly implemented using 1st kind constant-head BCs on both sides of the box models.As for the city-scale model, a 2nd kind no heat flux BC, i.e., a thermal insulation BC at the top of each box model (Fig. 3) leads to more conservative values for the power density (Ohmer et al. 2022).2nd kind no heat flux BC are also applied at the model bottom as well as at the upstream and downstream sides.The remaining hydraulic and thermal parameters populating the box models correspond to those from the city-scale model (Table 1).The spatial discretization of each box model comprises about 55,000 tetrahedral elements with a finer discretization around the ATES wells.

ATES configurations and model implementation
Under high groundwater flow velocities, substantial loss of stored thermal energy can occur, caused by the displacement of the injected water volume along the hydraulic gradient leading to low ATES efficiencies (Bloemendal and Hartog 2018;Bloemendal and Olsthoorn 2018).Installing two well doublets per ATES system in a line parallel to the direction of the ambient groundwater can reduce these thermal energy losses (Fig. 4).The appropriate pumping scheme then involves injecting the heated or cooled water at the upstream wells.In the following season, the stored and since displaced heated or cooled water is extracted from the corresponding downstream wells.To further improve recovery of stored thermal energy, ATES configurations with three doublets are also possible.In this case, the ATES system consists of an upstream injection well doublet, a downstream extraction doublet and a middle well doublet operating in an alternating way comparable to single-doublet systems.For 2-doublet and 3-doublet systems, the iterative adaptation of the distance between the individual ATES doublets achieves the highest possible recovery rates of stored thermal energy in each box model.Six different ATES configurations were simulated for each of the regions 1 to 3. ATES systems with one, two and three doublets were either placed in the Neuenburg Formation or in the Breisgau Formation leading to a total of 18 distinct box models.All wells are implemented with fully penetrating well screens over the entire thickness of the respective formation as proposed by Bloemendal et al. (2018).Each box model is run for 30 years according to the typical expected lifetime of ATES applications (Bloemendal et al. 2014;Sommer et al. 2015).
The distance between the warm and the cold wells of a well doublet equals two times the thermal radius R th of the individual wells which can be computed as (Doughty et al. 1982): Here c w and c aq represent the thermal capacities of water and the aquifer.V marks the volume of water that is injected during one injection period.The filter screen length of the ATES wells is represented by L. The lateral inter-well distance of two times R th ensures that no thermal interference between the warm and the cold wells of a well doublet occurs, which would lead to storage losses.
Typical seasonal ATES systems are used in heating mode in winter and in cooling mode in summer.However, short-term variations in energy demand may cause the system operation to shut down temporarily.Frequent switching between heating and cooling operation can also occur due to diurnal variations.In this study, these shortterm fluctuations are not regarded since they presumably do not affect the long term, overall characteristics of the thermal impact on the aquifer (Sommer et al. 2015).Accordingly, an ATES pumping scheme consisting of a 4 months period of heating during winter and a cooling period of the same length during summer is implemented in the box models.During the 2 months interim periods, the simulated ATES systems are not in active operation.This operation scheme corresponds to existing Dutch ATES systems (Sommer et al. 2013;2014).In the box models, the pumping scheme is implemented as time varying 2nd kind specified-flux BCs with flow rates of 600 m 3 d −1 according to typical existing GWHP systems in the city of Freiburg (Table 3). (1) The injection temperature at the warm and cold wells are defined as 1st kind BCs.Warm water is injected with a constant temperature of 18 °C, while cold water injection is set to a temperature of 6 °C.These temperatures result from assumed temperature differences during ATES operating of ± 6 K with respect to the ambient temperature of 12 °C.In Germany, this difference of ± 6 K is considered as the maximum acceptable change of groundwater temperature caused by open geothermal installations such as ATES and GWHP systems (Hähnlein et al. 2011;2013).

Calculation of thermal recovery
The efficiency of individual ATES wells or doublets of ATES wells of the same kind (i.e.warm or cold) in terms of storage loss is commonly quantified by the thermal recovery TR (Gao et al. 2017).This quantity is the ratio of extracted thermal energy and the thermal energy injected during the previous injection period, both with respect to the ambient aquifer temperature.TR accordingly describes thermal loss due to advective, conductive and dispersive heat transport as well as potential thermal interferences (Abuasbeh et al. 2021;Birhanu et al. 2015;Fleuchaus et al. 2020;Sommer et al. 2014Sommer et al. , 2013)).It can be calculated as (Abuasbeh et al. 2021): Here, E extr and E inj represent the extracted and injected thermal energy with respect to the ambient aquifer temperature T amb .T extr and T inj indicate the temperatures of the extracted and injected groundwater, respectively.V̇e xtr and V̇i nj are the extraction and injection flow rates, which are identical and constant over time according to Table 3.The values of extracted and injected thermal energy are calculated in COMSOL for both the warm and the cold ATES well(s).The corresponding energy ratio TR is then determined for each complete pumping cycle consisting of injection, passive storage and extraction.
Since the temperature of the surrounding aquifer progressively adopts the comparatively higher or lower temperature of the injected water, early storage and recovery cycles typically exhibit higher conductive storage loss and, therefore, show lower TR values (Sommer et al. 2013).After the tenth cycle, however, no further increase of TR was observed during simulations.This is in good agreement with statements from previous studies (Bakr et al. 2013;Duijff et al. 2021;Sommer et al. 2013).Accordingly, (2) the representative TR value used for further evaluation equals the average of TR for the warm and the cold well(s) for the tenth complete cycle, i.e., in the tenth year of operation.

Calculation of ATES power density
ATES power density values, which relate the amount of thermal power supplied by ATES systems to the required horizontal Earth surface area, are calculated for the city of Freiburg using the 18 box models based on the assumption that ATES systems are installed in the city as dense as possible without individual systems thermally influencing adjacent systems.For this purpose, we use the so-called thermally affected zone (TAZ) around the ATES wells in each box model after 30 years of ATES operation (Lo Russo et al. 2012).In the literature, the TAZ is commonly defined as the area where the absolute value of the temperature increase or decrease caused by ATES or GWHP systems exceeds 1 K, i.e., by the ± 1 K-plumes (Gizzi et al. 2020;Lo Russo et al. 2012;Piga et al. 2017).However, this limit for adverse thermal interferences lacks scientific justification and seems to be chosen almost arbitrarily (Pophillat et al. 2020a).Given these uncertainties, we choose the more spacious ± 0.5 K-isotherms after 30 years to delineate the TAZ in the block models.While this approach ensures even smaller thermal influences between adjacent systems, it also leads to more conservative power density values.For calculating the power density in this study, the space requirements with respect to the horizontal earth surface, i.e., the TAZ surface area, are determined as the smallest possible rectangle around the ± 0.5 K-isotherms after 30 years as shown in Fig. 5.The impact of choosing the ± 0.5 K-isotherms is briefly evaluated by also calculating the power densities based on the ± 1 K-isotherms for comparative purposes.
Based on the TAZ surface area A TAZ the power density values for heating mode PD heating and cooling mode PD cooling can be calculated as: Here, it should be noted that the calculated power densities are the annual mean power densities for heating and cooling via ATES and therefore also incorporate the periods of a year in which the system does not operate in the respective mode.This is accounted for by the availability factor f ava resulting from the 4 months' time period per year in which the system operates in heating or in cooling mode, respectively.Thus, the factor is f ava = 4/12.Calculating the power density this way allows for a meaningful comparison with the existing heating and cooling demand, which is given as the total energy demand per area and year.
The mean values of heating power Pheating and cooling power Pcooling during heating and cooling periods in eqs.( 3) and ( 4), are calculated according to: Here, the inclusion of the thermal recovery TR allows to utilize the constant injection temperatures T inj,warm and T inj,cold at the warm and the cold ATES storage, respectively, instead of the respective extraction temperatures which vary throughout the extraction phase.
The factor f HP considers the operation of a heat pump during heating mode, which adds heating power originating from the electricity grid to the thermal energy stored in the groundwater.Cooling, on the other, is assumed to be feasible without the operation of a heat pump (i.e., direct cooling).The factor f HP can be determined from the coefficient of performance COP of the heat pump according to: In this study, we assume a typical coefficient of performance COP = 3.5 (Bayer et al. 2012;Born et al. 2022;Duijff et al. 2021;Saner et al. 2010).This results in the factor f HP = 1.4.

Determination of ATES heating and cooling supply rates
This study follows the nomenclature by Tissen et al. (2019), who defined ATES heating and cooling supply rates as the shares of residential heating and cooling energy demands, respectively, that can potentially be supplied by ATES applications.
The residential heating energy demand for the city of Freiburg is available at building block level, calculated based on building characteristics, such as energetic classification, living space and building age (LUBW 2017).Besides space heating demand, the heating demand also includes thermal energy needed for residential water heating.The data set presents the heating demand referring to the original building conditions as well as lower demand values assuming building refurbishment.Since energetic building refurbishment is an important pillar of the German energy transition in the building (4) sector (Grossmann 2019), the supply rates in this study refer to the heating demand after refurbishment.These demand values are related to each building block's areal extent resulting in the heating energy demand given in MWh ha −1 a −1 .
In contrast to the heating energy demand, no such detailed data are available for the cooling energy demand in the city of Freiburg.Thus, we estimate the cooling energy demand using a ratio of 5 to 1 for heating demand to cooling demand.This demand ratio is obtained from the study by Werner (2016), who determined living space specific cooling demands of residential buildings in the European Union countries including Germany.The ratio of 5 to 1 also corresponds well to the heating and cooling energy demand in Freiburg from the European hotmaps project datasets (Mueller 2019;Mueller and Fallahnejad 2020).
The calculation of ATES heating and cooling supply rates requires the conversion of the energy demand data to the power density's physical unit W m −2 .Due to its universal nature, the power density PD determined with the ATES box models allows for its straightforward comparison with the thermal energy demand ED to calculate possible supply rates SR as follows (similar to e.g.Epting et al. (2018), Tissen et al. (2019)):

Thermal recovery
The thermal recovery TR is an important parameter describing the storage efficiency of ATES systems and quantifying thermal energy loss in the subsurface.In this study, it is determined using numerical 3D box models, each of which represents a distinct hydrogeological region of the city-scale subsurface model of Freiburg.For each ATES configuration the thermal recoveries are calculated according to Eq. ( 2).As shown in Fig. 6 thermal recoveries of ATES systems in the deeper aquifer, i.e., the Breisgau Formation, are consistently higher for all ATES configurations and all regions compared to systems in the upper aquifer (Neuenburg Formation), which is characterized by higher groundwater flow velocities (Table 4).To mitigate the detrimental effects of the groundwater flow on the energy storage efficiency, different ATES configurations with one, two or three well doublets are modeled as described above.For the Breisgau Formation, 3-doublet ATES configurations show the highest thermal recoveries with up to TR = 59% in hydrogeological region 1 (Fig. 6).In contrast, 1-doublet systems recover the lowest share of thermal energy.This system configuration also shows the lowest recovery values for ATES in the Neuenburg Formation, while the highest thermal recoveries for the Neuenburg Formation can be observed for systems with 2 well doublets with up to 42%.These results demonstrate the strong influence of the ambient groundwater flow velocity on the thermal recovery and the suitable ATES design.High advective heat transport rates caused by high groundwater flow velocities of up to 12.7 m d −1 entail significant subsurface energy loss (Fig. 6, Table 4).Thus, thermal recovery values for 1-doublet systems in the Neuenburg Formation are very low with the maximum being TR = 4% in hydrogeological region 1.By adding a second downstream extraction well doublet, thermal recoveries substantially increase up to TR = 42%.However, the thermal recovery values for 2-doublet ATES configurations are still significantly lower than typical recovery values reported in the literature and discussed below.This is especially true for regions 2 and 3 with TR = 27% and TR = 13%, respectively, indicating that at high groundwater flow velocities the 2-doublet configuration type can mitigate thermal losses only to a limited extent.The use of three well doublets in the upper aquifer does not improve thermal recovery but instead leads to recovery values which are slightly lower than for 2-doublet systems.Compared with the groundwater flow velocity, variations of other parameters, such as the ambient groundwater temperature and the solid thermal conductivity, have been shown to have a negligible effect on the thermal recovery.Further details on this are provided in the Additional file 1: Section SD3.
Previous studies of LT-ATES systems revealed thermal recoveries of 47% up to 90% for ambient groundwater flow velocities ranging from 0.01 m d −1 to 1.6 m d −1 further demonstrating this parameter's influence on the storage efficiency (Abuasbeh et al. 2021;Bloemendal and Olsthoorn 2018;Kangas and Lund 1994;Sommer et al. 2014).Figure 7 shows the relation between thermal recoveries of ATES systems from the literature as well as from this study and the corresponding ambient groundwater flow velocities.The recovery values from this study correspond to the maximum recoveries of each hydrogeological region as shown in Fig. 6.Sommer et al. (2014) analyzed an LT-ATES in Utrecht, the Netherlands, where the flow velocity is a mere 0.01 m d −1 and, therefore, significantly lower than in the Freiburg study area (Table 4).A thorough monitoring showed mean thermal recoveries during a seven-year operation period of 68 and 82% for the warm storage area and the cold storage area, respectively.Drijver et al. (2012) state values between 70 and 90% as typical range of thermal recoveries for LT-ATES systems in aquifers with low flow velocities.Similar recovery values are reported in Bakr et al. (2013) for multiple densely placed LT-ATES systems in the Dutch city of The Hague using thermohydraulic modeling.The minimum thermal recovery amongst all ATES systems was 68% in the first year.Over the course of the ten years modeling period, the recovery values showed an increasing trend towards steady state values with the maximum thermal recovery being 87% in the tenth year.The above-mentioned studies all refer to Dutch systems operating under conditions of very low ambient groundwater flow velocities.
Other publications studied ATES systems located in aquifers with higher ambient groundwater flow velocities.Numerically computed thermal recoveries of mostly between 73 and 80% are stated for 2-doublet ATES systems and an ambient groundwater flow velocity of 0.1 m d −1 in Bloemendal and Olsthoorn (2018).For a flow velocity of 0.3 m d −1 , which entails higher subsurface heat loss, the same study gives recovery values ranging mostly from 65 to 75%.Flow velocities of about 0.1 m d −1 and 0.2 m d −1 are stated for a Swedish LT aquifer storage system, which result in mean thermal recovery values of 47% and 60% for the warm and the cold storage areas, respectively (Abuasbeh et al. 2021).These numbers are in line with this study's thermal recoveries for 3-doublet systems in the Breisgau Formation ranging from TR = 49% to TR = 59% corresponding to groundwater flow velocities between  4).In the literature, it was also demonstrated that LT-ATES systems can be possible with ambient flow velocities of up to 1.6 m d −1 when using multi-doublet configurations (Kangas and Lund 1994).

Power density of ATES
The heating and cooling power densities achievable with ATES systems in the city of Freiburg are calculated according to eqs. ( 3) and ( 4), respectively, for the ATES designs with the highest thermal recoveries (Fig. 6).For the Neuenburg and the Breisgau Formation, these are the 2-doublet and 3-doublet configurations, respectively.
The highest power densities for ATES systems placed in the Neuenburg Formation are calculated for hydrogeological region 3 with PD heating = 1.6 W m −2 and PD cool- ing = 1.2 W m −2 , while aquifer storage systems in region 1 lead to the lowest power densities of PD heating = 1.3 W m −2 and PD cooling = 0.9 W m −2 (Fig. 8).This is in contrast to the thermal recovery, which is highest for region 1 and lowest for region 3 (Fig. 6) and shows the high influence of the thermally affected zone's area on the power density according to eqs. ( 3) and ( 4).Due to the Neuenburg Formation's high ambient groundwater flow velocities in region 2 and 3 relative to region 1 (Table 4), the injected thermal plume undergoes a more pronounced dispersive spread in regions 2 and 3 leading to shorter ± 0.5 K-isotherms.In region 1, on the other hand, the comparatively low flow velocity results in more stable and thus much more elongated ± 0.5 K-isotherms.Accordingly, the TAZ for region 1 is about 58% larger than for region 3 and about 30% larger than for region 2 resulting in the lowest power density for ATES systems in region 1.
These effects are similar for the utilization of the deeper Breisgau Formation.There, however, the ambient groundwater flow velocity and accordingly the extent of the TAZ vary much less between regions 1 to 3. The highest and lowest heating power density of ATES systems placed in the Breisgau Formation are calculated to PD heating = 3.2 W m −2 Fig. 8 ATES power densities for Neuenburg and Breisgau Formations.ATES systems in the Neuenburg Formation are simulated as 2-doublet configurations.ATES systems in the Breisgau Formation use a 3-doublet configuration and PD heating = 2.7 W m −2 for hydrogeological regions 2 and 1, respectively (Fig. 8).As for the systems in the Neuenburg Formation, due to free cooling without the use of a heat pump, the power density values for cooling mode are uniformly smaller by the factor f HP = 1.4.
Previous studies on the technical potential of shallow geothermal applications were compared regarding their power density in Bayer et al. (2019).However, they only covered studies on closed geothermal systems such as GSHP systems.The compiled power density values range from about 7 W m −2 up to a 460 W m −2 reflecting a large variety of underlying assumptions and methodological approaches (Bayer et al. (2019) and references therein).Power densities of GSHP systems in an urban quarter ranging from 14 W m −2 to 93 W m −2 can be inferred from Tissen et al. (2019).GSHP systems typically induce much smaller thermal anomalies in the subsurface explaining these consistently higher power densities compared to the values of the open systems in this study (Perego et al. 2022).
Power density values for open GWHP systems in an urban quarter described in the literature were calculated similar to the approach described in this study, i.e., using the areal extents of the thermal plumes (Tissen et al. 2019).However, these were determined analytically, and may therefore deviate significantly from numerically simulated thermal plumes (Pophillat et al. 2020a).Nevertheless, with values of about 3 W m −2 , the power densities for two GWHP scenarios are similar to the power density values from this study (Fig. 8).It should be noted, however, that Tissen et al. ( 2019) used the ± 1 K-isotherms to delineate the TAZ compared to our more conservative approach of using the ± 0.5 K-isotherms.For further comparisons, Fig. SD4 in the Additional file 1 presents the less conservative power density values calculated from the box models using the ± 1 K-isotherms after 30 years.On average, these power densities are about 2.1 times as high as the values shown in Fig. 8.The strong sensitivity of the thermal recoveries of open systems to the groundwater flow velocity discussed above also contributes to the differences in reported power densities.
Open GWHP systems were also evaluated regarding their thermal potential in Epting et al. (2020) in the Swiss city of Basel.In an exemplary city quarter, power densities between 12 W m −2 and 720 W m −2 can be inferred for various well distances and temperature differences of up to 8 K.These power densities are significantly higher than in this study and in Tissen et al. (2019).This is because Epting et al. (2020) only considered hydraulic effects resulting from groundwater extraction and injection rather than accounting for thermal plume propagation.The power densities' reference to surface area therefore also refers to the much smaller inter-well distances of between 10 and 50 m leading to higher power densities.
The heating and cooling power densities calculated in this study are annual mean values resulting from the used ATES pumping scheme, i.e., a 4-month period each for the heating and cooling modes and two 2-month passive periods in between.This is in contrast to the cited studies, where the power density only refers to the periods when the system is actually in operation.In the mentioned studies of closed systems, this period is a fixed operating time per year of 2400 h a −1 or 1700 h a −1 (Bayer et al. 2019;Tissen et al. 2019), while the power density for the open GWHP systems studied in Tissen et al. (2019) relates to a year-round operation.
Our study's results refer to a specified technology used in well-defined hydrogeological conditions.Model-implemented characteristics, such as specified pumping rates or hydraulic gradients, lead to more realistic power density values constrained by technical restrictions.This is in contrast to many previous studies on the topic of shallow geothermal power densities, several of which are based on more general assumptions or ignore the influence of groundwater flow (Bayer et al. 2019).The present study also for the first time distinguishes between power densities for heating and cooling modes.In this respect, it is also important to stress the influence of the 2nd kind no heat flux BC applied to the top side of the numerical models.As shown by Ohmer et al. (2022) this can significantly increase the lateral thermal plume propagation as it impedes any dissipation or input of thermal energy to or from the atmosphere.

Heating and cooling supply rates with ATES
The spatial distribution of ATES power densities is shown in Fig. 9 for ATES systems placed in the Breisgau Formation, which have higher thermal recoveries (Fig. 6) and power densities (Fig. 8), indicating that ATES systems in this formation are more feasible regarding heating and cooling supply rates.The figure also presents the supply rates calculated according to Eq. ( 8) for each block of residential buildings.The bar charts in Fig. 9 show the percentages of residential buildings in the study area for which ATES systems placed in the Breisgau Formation could supply a certain share of their heating or cooling energy demand.For about 28% of the residential buildings, the heating supply rate is above 80%.Heating demand supply rates of 60% or more could be achieved for about 50% of the residential buildings.Especially for residential buildings located in suburban or commercial and industrial districts, such as the commercial and industrial area in the northern part of the study area, ATES could supply the entire heating demand.In contrast, ATES systems in the Breisgau Formation can completely supply about 92% of all residential buildings with cooling energy (Fig. 9).
Similar to this study, Schiel et al. (2016) showed in a previous study on urban space heating with GSHP systems that the highest supply rates were located in suburban areas rather than the city's center.For Southwest Germany, calculations of the technical potential of GSHP systems showed that the heating demand of 65% to 93% of all buildings could be completely supplied by such systems with the highest supply rates again concentrating on rural and suburban areas (Miocic and Krecher 2022).This is in line with a study of western Switzerland revealing that a complete supply of heating demand by GSHP systems is possible in many rural and suburban areas as opposed to densely populated cities for which a supply rate deficit is to be expected (Walch et al. 2021).These findings from previous publications match the spatial distribution of the ATES heating supply rates in Freiburg, which are lower in more densely populated areas (Fig. 9).
The above-mentioned studies refer to heating energy supplied by GSHP systems.Besides GSHP systems, Tissen et al. (2019;2021) also calculated heating supply rates of GWHP systems for urban settings in two different cities.For both cities, the supply rates of GSHP systems are significantly higher than of GWHP systems.GSHP systems could potentially supply the total heating demand, i.e., 100%, after refurbishment in more than half of all districts in one of the cities.In contrast, the highest supply rate of GWHP systems was determined to be 83% among all districts.The numbers for GWHP systems are in line with this study's heating supply rates of ATES systems, which are also open shallow geothermal systems.

Greenhouse gas emission savings with ATES
In the year 2020, the final energy demand for residential space and water heating in the city of Freiburg was 1000 GWh a −1 (GEF Ingenieur AG 2021).The energy was mostly supplied by a mix of natural gas, district heating and heating oil (further details in Additional file 1: section SD5).To evaluate possible GHG emission savings achievable with ATES, two different scenarios are considered using ATES systems in the Breisgau Formation.These calculations are based on the emission factors of the Freiburg heating energy mix as well as on an ATES emission factor of 0.083 tCO 2 eq MWh −1 adopted from Stemmle et al. (2021).
Scenario 1 assumes an ATES supply limited to the residential building blocks for which a heating supply rate of 100% or more was calculated in the previous chapter.This is true for about 15% of all residential buildings (Fig. 9).In scenario 2, ATES systems supply all residential buildings according to the building block specific supply rates.Buildings located in the study area's hydrogeological region 4 are excluded.
Figure 10 shows the resulting annual GHG emission savings.Scenario 1 yields savings of about 17,023 tCO 2 eq a −1 , which equals about 10% of the estimated current annual GHG emissions caused by space and water heating of all considered residential buildings in the study area.Higher GHG emission savings result for scenario 2. They amount to about 70,398 tCO 2 eq a −1 or 40% of the total residential space and water heating related GHG emissions in the study area.Installing individual ATES systems for buildings with small supply rates is unlikely to be implemented in practice.However, integration of ATES in the existing district heating network in Freiburg poses a promising option to make use of the full technical potential.
For space cooling, lack of information on emission factors and the current structure of supply in Freiburg prevents calculation of GHG emission savings.The high cooling supply rates (Fig. 9), however, indicate a large unused potential for sustainable space cooling and high GHG emission savings compared to conventional cooling technologies, such as compression chillers.
Like ATES, air source heat pumps (ASHP) are systems that can provide both, heating and cooling energy.In contrast to ATES, ASHP systems typically have lower coefficients of performance (COPs) due to much stronger temperature variations of the heat source, i.e., the outside air, throughout the year (Esen et al. 2007;Gao et al. 2021;Sarbu and Sebarchievici 2014).Thus, ATES systems typically require less electricity for the same amount of thermal output than ASHP systems.These energy savings are another advantage of ATES besides the possibility of using otherwise unused waste thermal energy.

Limitations of the box model approach
The box model approach used in this study prevents any thermal interferences between ATES systems, which impedes detrimental effects on the thermal recoveries.Various studies showed, however, that a holistic planning framework based on the coordinated placement of ATES systems can maximize subsurface utilization in contrast to separate planning of individual systems.This is of great importance in dense urban areas with a high ATES adoption rate and other subsurface infrastructure, such as sewage and traffic tunnels.A central component of such energy planning framework is the identification of an optimal trade-off between the combined energetic benefits across all ATES systems and the efficiency of individual systems.Furthermore, deliberate thermal interferences between wells of the same type, i.e., wells of either the warm or the cold storage areas, can decrease thermal loss at the storage volume boundaries due to a better ratio of storage surface area to storage volume (Bloemendal et al. 2018;Duijff et al. 2021;Pellegrini et al. 2019).Expanding on this study's modeling approach aimed at ATES potential determination, a city specific optimization of ATES placement in Freiburg could account for such positive thermal interferences, while also including existing shallow GWHP systems.
Utilizing a uniform mean ambient groundwater temperature of 12 °C across all box models, i.e., in all hydrogeological regions, does not allow to study the influence of different thermal regimes on thermal plume propagation.Due to the ambient groundwater flow velocity's strong influence on the plume propagation, using a mean value for this parameter seems reasonable in the context of a city-level potential assessment.This assumption is further supported by Fig. SD3 in the Additional file 1, which shows a very small influence of ambient groundwater temperatures ranging between 10 °C and 13 °C on the thermal recovery.Nevertheless, future studies using a similar approach of representative box models could potentially achieve more accurate potential estimations by including information on ambient groundwater temperature variations.Reasons for spatial variations of the ambient groundwater temperature in urban areas include heat fluxes caused by anthropogenic heat sources and sinks, such as basements, subway infrastructure, sewage systems and district heating grids (Benz et al. 2015).The presented model approach does not account for these heat fluxes due to the box models not being spatially localized within the study area besides their allocation to one of the hydrogeological regions (Fig. 2).
Figure 11 exemplarily shows the temperature distribution after 30 years in the cityscale model with ATES systems implemented in the Neuenburg Formation.The distances between individual ATES systems along the groundwater flow direction correspond to the lengths of the TAZ in the box models and thus depend on the hydrogeological region (Fig. 2).Along the Dreisam River the ATES placement is less dense, and hydrogeological region 4 is again excluded.
Most of the modeled ± 0.5 K-isotherms shown in Fig. 11 are shorter than the length of the TAZ in the respective box model, while some of them are longer.This is to be expected since the box models are only approximations for each hydrogeological region.The differences of the thermal plume lengths between the city-scale subsurface model and the box models could be reduced by increasing the number of hydrogeological regions or the number of variable parameters in the box models, e.g., the site-specific formation thickness of the aquifers.
The use of representative box models has a number of advantages compared to using the city-scale subsurface model.Provided that a reasonable delineation of sufficiently homogeneous hydrogeological subsurface regions is possible, the approach allows a fast approximation of power density values for ATES operation in different cities or regions.In this study, they took an average of about 55 min to simulate the ATES operation over 30 years.In contrast, the Freiburg city-scale model took about 50 h to complete the task on the same computer (8 CPU cores with a base clock of 3.6 GHz and 128 GB of RAM).This also means that adaptation and evaluation of different ATES system configurations are less time consuming.

Conclusions
Using representative numerical 3D thermo-hydraulic models of a simple geometry, this study assesses the technical potential of low-temperature aquifer thermal energy storage (LT-ATES) applications in the city of Freiburg.For this purpose, the power density, which relates the amount of power generated by a specific technology to the required surface area, of space heating and cooling energy supply using ATES systems is quantified.
Simulating various ATES configurations with multi-doublet designs reduces energy storage loss.Using two or three well doublets consistently increases thermal recoveries compared to single-doublet systems.Nonetheless, relatively high groundwater flow velocities of up to about 13 m d −1 considerably reduce the recovery of stored thermal energy.This is especially true for ATES in the upper aquifer, the so-called Neuenburg Formation, where the maximum thermal recovery is 42%.For ATES operation in the deeper Breisgau Formation with lower groundwater flow velocities of less than 1 m d −1 , thermal recovery values of up to 59% are obtained.
Accordingly, ATES in the Breisgau Formation leads to higher power densities of up to 3.2 W m −2 , while the highest power density for the Neuenburg Formation is only 1.6 W m −2 .For the Breisgau Formation, this also enables considerably higher ATES supply rates with respect to the existing residential heating and cooling demand in Freiburg.While heating energy supply rates of larger than 60% are determined for about 50% of all residential buildings in the study area, the cooling energy demand could be supplied entirely by ATES systems for 92% of the buildings.
Based on the calculated supply rates and today's final energy mix for space and water heating in Freiburg, potential GHG emission savings of up to about 70,000 tCO 2 eq a −1 for ATES heating alone can be estimated.This equals about 40% of the current overall GHG emissions caused by space and water heating in the study area's residential buildings.While the extensive utilization of all the available subsurface space using ATES systems is not realistic, these numbers still show promising opportunities for ATES applications in the city of Freiburg.
In the future, this modeling approach could be expanded upon with the aim of integrating ATES into more specific and practice-oriented urban energy planning.This way, the technology could face an increasing use and could help to achieve climate protection goals at the municipal level and beyond.A brief comparison of our study results with power densities and supply rates from the literature reveals that closed GSHP systems should also be considered in urban energy planning scenarios since this type of shallow geothermal systems can potentially lead to higher power densities and supply rates.Numerical models as proposed in this study could help to identify the suitable type of shallow geothermal system for regions with similar hydrogeology.

Fig. 1
Fig. 1 a, b: Location of the city of Freiburg im Breisgau with the highlighted study area in Southwest Germany.c Profile section A-A*.The Riegel Horizon (RH) shown in the profile section is not regarded in the numerical models.Data from GDI-BW (2015), Geofabrik (2022), USGS (2017).Hydraulic head data from the Environmental Protection Authority Freiburg and the Baden-Württemberg State Institute for the Environment, Survey and Nature Conservation (LUBW).Profile section modified from Wirsing and Luz (2005)

Fig. 2
Fig. 2 Delineated representative hydrogeological regions within the study area of Freiburg.Region delineation is based on the horizontal hydraulic conductivity and the hydraulic gradient

aFig. 3
Fig. 3 Exemplary box model of a 2-doublet ATES system in the Neuenburg Formation in hydrogeological region 1.The injection wells are implemented via 1st kind BCs (temperature) and 2nd kind BCs (mass flow rate).The extraction wells are implemented using 2nd kind BCs (mass flow rate).The black lines mark the ± 1 K-isotherms and ± 0.5 K-isotherms, respectively

Fig. 4
Fig. 4 Schematic ATES configurations showing top view illustrations of ATES systems with one, two and three well doublets

Fig. 5
Fig. 5 Exemplary box models of two 2-doublet ATES systems in the Neuenburg Formation in hydrogeological region 1 (left) and hydrogeological region 2 (right).The TAZ around the ± 0.5 K-isotherms after 30 years of operation are highlighted.The smaller ± 1 K-isotherms are also shown (black solid lines)

Fig. 6
Fig. 6 Thermal recovery values for various ATES configurations in different groundwater flow regimes.Ambient groundwater flow velocities v in the Neuenburg and the Breisgau Formations in hydrogeological regions 1-3 are also shown

Fig. 7
Fig. 7 Thermal recoveries of ATES systems from this study and previous publications plotted against the corresponding ambient groundwater flow velocities.Data points from Abuasbeh et al. (2021),Bloemendal  and Olsthoorn (2018),Sommer et al. (2014)

Fig. 9
Fig. 9 ATES supply rates in the city of Freiburg for heating (left) and cooling (right) using the Breisgau Formation.The bar charts illustrate the share of the total residential buildings within each individual supply rate class

Fig. 10
Fig. 10 Possible GHG emission savings from ATES heating compared to the current heating energy mix in Freiburg.Calculated for two scenarios based on data from GEF Ingenieur AG (2021)

Fig. 11
Fig. 11 Top view of the city-scale subsurface model of Freiburg with ± 0.5 K-isotherms after 30 years of ATES operation in the Neuenburg Formation.The smaller ± 1 K-isotherms are also shown

Table 1
Subsurface parameters and their corresponding values used in the numerical Freiburg cityscale model and box models

Table 2
Representative hydraulic conductivities and hydraulic gradients of the four defined hydrogeological regions in Freiburg

Table 3
ATES design parameters and respective values used in the box models