Effects on the Unit Commitment of a District Heating System Due to Seasonal Aquifer Thermal Energy Storage and Solar Thermal Integration

Abstract

The ongoing transformation of district heating systems (DHSs) aims to reduce emissions and increase renewable energy sources. The objective of this work is to integrate solar thermal (ST) and seasonal aquifer thermal energy storage (ATES) in various scenarios applied to a large DHS. Mixed-integer linear programming (MILP) is used to develop a comprehensive model that minimizes operating costs, including heat pumps (HPs), combined heat and power (CHP) units, electric heat boilers (EHBs), heat-only boilers (HOBs), short-term thermal energy storage (TES), and ATES. Different ATES scenarios are compared to a reference without seasonal TES (potential of 15.3 GWh of ST). An ATES system with an injection well temperature of about 55 °C has an overall efficiency of 49.8% (58.6% with additional HPs) and increases the integrable amount of ST by 178% (42.5 GWh). For the scenario with an injection well temperature of 20 °C and HPs, the efficiency is 86.6% and ST is increased by 276% (57.5 GWh). The HOB heat supply is reduced by 8.9% up to 36.6%. However, the integration of an ATES is not always economically or environmentally beneficial. There is a high dependency on the configurations, prices, or emissions allocated to electricity procurement. Further research is of interest to investigate the sensitivity of the correlations and to apply a multi-objective MILP optimization.

1. Introduction

The use of renewable heat from geothermal, solar, and biomass sources in urban district heating systems (DHSs) has gained importance [1]. As more (renewable) heat sources are integrated, the importance of thermal energy storage (TES) increases [2]. Short-term TES in DHS is mainly used to shape daily or weekly heat load variations [3,4]. Mathematical modeling and optimization of DHSs including short-term TES usually show the economic benefits of combined heat and power (CHP) plants when producing at times of high wholesale electricity prices [5,6]. Furthermore, the use of short-term TES in combination with heat pumps (HPs) or electric heat boilers (EHBs) offers the potential for decarbonizing the heat supply in DHS and for more efficient control of power systems with high shares of renewable energy [7,8]. In contrast, long-term TES is particularly suitable for a seasonal mismatch between the availability of renewable energy sources, e.g., solar thermal (ST), and the heat demand of the buildings connected to a DHS [1]. In general, there are four types of long-term TES: tank thermal energy storage (TTES), underground or pit thermal energy storage (PTES), borehole thermal energy storage (BTES), and aquifer thermal energy storage (ATES) [9]. For DHS applications, TTES and PTES are state of the art, e.g., storing thermal energy in large basins has been implemented commercially in Danish DHSs [10,11]. The integration of other seasonal TES is not yet widespread and is instead more scientifically studied with different levels of detail regarding specific technologies [12] or focusing on, e.g., city-scale applications [13].

1.1. Literature Review

Usually, BTES is researched at a smaller scale, e.g., single-family buildings or residential districts with less than 100 buildings [14]. However, there are some interesting methodological approaches related to heat supply concepts including BTES. Miglani et al. [15] developed a multi-objective genetic algorithm considering costs and carbon emissions for concept design as well as mixed-integer linear programming (MILP) optimization of typical demand days estimating the operation of ground source HPs and ST for a single house. Based on this bi-level optimization framework, they further optimized the energy system for the same use case also considering solar photovoltaic, EHBs, and small-scale TES [16]. Renaldi and Friedrich [17] conducted a multi-year parametric study for the techno-economic design of ST and BTES without HPs on a district level including 52 buildings and an annual heat demand of approximately 611 MWh. Reed et al. [18] performed a financial model calculation with about 10,000 scenarios resulting in discounted savings per building for a 50-house development with DHS including BTES and ST compared to individual heating with natural gas-fired heat-only boilers (HOBs).

More suitable for large-scale use is ATES, which is a geological system of rocks used as a storage medium containing groundwater, which is used as a heat transfer fluid and is extracted from the aquifer through multiple groundwater wells and reinjected after being heated or cooled [19]. However, particular DHS concepts including ATES depend heavily on local conditions, which are determined by hydrogeological requirements of the underground rock as well as individual DHS temperature levels, often leading to design concepts with HPs that boost the ATES temperature for utilization in high-temperature DHS [9]. Pellegrini et al. [14] analyzed the barriers to ATES implementation in Europe and derived possible solutions to transform the identified barriers into opportunities for the development of ATES. According to Guelpa and Verda [9], the characteristics of ATES are low construction costs with relatively good thermal capacity compared to other seasonal TES.

However, research often focuses on the thermal performance of ATES design using process engineering models or monitoring real-world systems, concluding that ATES can balance heating and cooling demands efficiently [20,21,22,23]. Hereby, e.g., Dickinson et al. [24] and Godinaud et al. [25] also considered HPs, but the ATES application was only on a small scale for a few buildings not including DHS. Moreover, detailed simulation studies are often performed in research. Gonzalez-Ayala et al. [26] performed finite element method (FEM) simulations of an ATES using accurate thermal maps of the subsoil in order to determine storage efficiencies. Li et al. [27] computed multilateral-well ATES without other energy system components in a multi-year approach and used FEM simulations to improve the utilization of underground space. In contrast, Beernink et al. [28] considered several ATES design improvements on a larger city scale including DHS but also with an isolated view of the ATES itself. This is also the situation with other current literature. Stemmle et al. [29] developed a numerical three-dimensional thermo-hydraulic model of a simple low-temperature ATES geometry in order to determine the technical potential for applications in Freiburg, Germany. Collignon et al. [30] developed a detailed numerical model estimating thermal losses and capacities of high-temperature ATES for DHS in Geneva, Switzerland, but neglecting fluctuating renewable heat supply from, e.g., ST. Hirvonen et al. [31] performed an energy system design study by combining ST and seasonal TES for different use cases and DHS sizes but did not use an optimization approach.

1.2. Resulting Problem and Contribution of This Paper

In summary, most current works outline energy system analyses with seasonal TES in different depths of detail. However, there is a lack of novel seasonal TES applied to DHS because either renewable heat from, e.g., ST, is only partially included or innovative methodological approaches using MILP optimization do not include ATES, even though ATES is more suitable for large-scale seasonal TES in urban DHSs with limited space than BTES and sometimes even PTES. In general, the market penetration of large-scale ST and ATES in Germany’s DHS is very low [19]. However, seasonal TES is expected to be a central cornerstone of future DHSs as well as cross-sectoral smart energy systems [1,2]. Thus, novel design approaches considering existing large-scale DHSs should also consider ATES. However, the majority of research does not analyze their operation potentials in DHSs with an interacting view of the current heat supply and future options regarding renewable supply using MILP optimization. Thus, the questions of economic and ecological potential remain unanswered so far.

Consequently, there is a research gap in the design and modeling of existing heat supply in DHSs toward a decarbonized energy system including ATES and other renewables. This paper presents the implementation of an MILP optimization algorithm to simulate the annual operation of a DHS, especially including ATES. The intention of this paper is to answer the question of how large a TES in combination with ST is economically and ecologically to be favored in urban DHS. This work reveals the potential of durable cheap heat from renewable heat sources, which is a significant contribution to the related topics. More precisely, the model developed is applied and analyzed with real-world data from a German energy utility. The conclusions derived shall help to identify design and operation parameters that support a sustainable and secure energy supply in urban areas with DHS.

Aiming to address this research issue, a novel modeling approach using MILP optimization is developed in the Section 2. The approach is applied to a real-world use case for a DHS located in Berlin, Germany, and the results are presented in the Section 3. The drivers of the results are identified in the Section 4. Furthermore, relevant system boundaries and model limitations are examined. Finally, the most significant outcomes are outlined in the Section 5.

2. Materials and Methods

First, a model of DHS operation is developed for the integration of ST energy through ATES. Based on MILP, a mathematical optimization model is computed in Python using the Gurobi 11.0.0 solver. Due to the computational time, the time horizon of such complex optimization problems is usually only a few weeks. Therefore, the next step is to cluster the use case data, e.g., heat load profiles, in order to simulate the DHS operation with an acceptable duration considering different heat suppliers and TES. Additional information about the use case and further assumptions are given in the last part of this section.

2.1. Linear Optimization Approach

The field of mathematical optimization is divided into several distinct areas, with MILP representing one of them. Similar to continuous linear optimization, it concerns the optimization of linear objective functions over a set constrained by linear equations and inequalities. The distinction lies in the fact that some of the variables are integer-valued and others are continuous. The method is applied to a multitude of disparate NP-hard complex problems, such as combinatorial problems, numerous optimization problems in the energy industry like unit commitment problems, transportation logistics such as the vehicle routing problem, or machine scheduling in production [32,33].

In this work, MILP is used to describe a cost optimization problem for the scheduling of production plants in DHS to achieve beneficial operation for both heat generation and the electricity grid, similar to the unit commitment problem [7]. The objective function shown in Equation (1) minimizes the aggregated cost of the operation ($C_{op}$) for all intervals of the simulation horizon ($i=1,\dots ,I$). Capital costs are not considered in optimized unit commitment problems.

$$\min C_{op}=\sum _{i=1}^I\left(C_{a,i}+C_{d,i}+\sum _{b=1}^B{C}_{b,i}+\sum _{h=1}^H{C}_{h,i}+\sum _{l=1}^L{C}_{l,i}+\sum _{k=1}^K{C}_{k,i}\right)$$

(1)

These total operation costs contain the costs for the different DHS units: the costs for the ST plant $C_{a,i}$, the costs for the ATES $C_{d,i}$, and the costs for the various controllable plants, EHB $C_{b,i}$, HOB $C_{h,i}$, large HP $C_{l,i}$, and CHP units $C_{k,i}$. The sum of all cost components over the entire time horizon $I$ is minimized.

Hereby, the costs for the ST plant are calculated from (2) through the ST heat load and the heat-related maintenance costs $c_{mh,a}$. The ST heat load depends on the length of a simulation interval $\mathsf{\Delta }I$, the collector area $A$ of the ST, and the specific collector yield for each simulation step ${qa}_{a,i}$. The size of the collector field is a decision variable.

$$C_{a,i}=\mathsf{\Delta }I\times A\times {qa}_{a,i}\times c_{mh,a} \forall a,i$$

(2)

The specific collector yield is determined by given parameters (solar irradiation, outside air temperature, wind speeds, flow temperatures, etc.) and is therefore not a decision variable. Solar yields must always be fed into the DHS in order to avoid stagnation of the ST system.

Equation (3) represents the operating costs for the ATES system $d$, based on the decision variable, the number of drillings $D$, the interval length $\mathsf{\Delta }I$, the specific heat supply to the DHS per drilling ${qd}_{d,i}$, the heat-related maintenance costs for the ATES $c_{mh,d}$, the specific electric demand of the ATES ${pd}_{d,i}$, and the electricity-related costs (electricity price $c_{el,i}$ and costs for the purchase of power from the grid $c_{pp}$). Here, one ‘drilling’ is used as a summarized expression for the pair of an injection well (cold side) plus a production well (warm side) of the ATES. For the charging period of the ATES, the specific heat supply is negative (treated like an additional heating demand).

$$C_{d,i}=\mathsf{\Delta }I\times D\left({qd}_{d,i}{c}_{mh,d}+{pd}_{d,i}\left(c_{el,i}+c_{pp}\right)\right) \forall d,i$$

(3)

Both the specific heat supply (or demand) ${qd}_{d,i}$ and the specific electricity demand ${pd}_{d,i}$ are calculated from the assumed magnitudes of the ATES. The heat supply (or demand) in Equation (4) equals the specific heat source of the storage ${hd}_{d,i}$ in addition to the potential electricity demand of the ATES-HP ${pd}_{hp,d,i}$, while $N_{hp,d}$ is a binary parameter for considering an ATES system with or without additional HPs to boost the DHS supply temperature.

$${qd}_{d,i}={hd}_{d,i}+N_{hp,d}{pd}_{hp,d,i} \forall d,i$$

(4)

If an HP is part of the ATES system configuration, the power demand for the HP ${pd}_{hp,d,i}$ is calculated with the HP’s coefficient of performance ${COP}_{d,i}$ and the available heat source ${hd}_{d,i}$, while the heat source has a positive value, as shown in (5). This means that the ATES system is in a discharge state.

$${hd}_{d,i}\ge 0⟹{pd}_{hp,d,i}={\displaystyle \frac{{qd}_{d,i}}{{COP}_{d,t}}}={\displaystyle \frac{{hd}_{d,i}}{{COP}_{d,t}-1}} \forall d,t,i$$

(5)

The ${COP}_{d,t}$ (6) for the ATES HP is determined by the exergy efficiency ${\eta }_{ex}$ and the Carnot efficiency ${COP}_{carnot,t}$, while the Carnot COP is based on the hot condensation temperature $T_H$ and the cold vaporization temperature $T_{C,t}$ of the HP cycle [34].

$${COP}_{d,t}={\eta }_{ex}{COP}_{carnot,t}={\eta }_{ex}\left({\displaystyle \frac{T_H}{T_H-T_{C,t}}} \right) \forall d,t$$

(6)

The total electricity demand of the ATES ${pd}_{d,i}$ in (7) is then the power supply of the potential HP ${pd}_{hp,d,i}$ plus the constant power demand for the deep groundwater well pumping ${pd}_{dd,d,i}$.

$${pd}_{d,i}=N_{hp,d}{pd}_{hp,d,i}+{pd}_{dd,d,i} \forall d,i$$

(7)

Similar to the ST heat supply, the ATES heat source, which is the basis for the ATES heat supply, is calculated from assumed or given parameters according to (8). The volumetric flow rate ${\dot{V}}_{d,i}$, the density of water $\rho$, the specific heat capacity of water $c_p$, and the temperature difference between the upper and lower temperature levels ($T_{ul,d,i}$ and $T_{ll,d,i}$) of the ATES side are needed.

$${hd}_{d,i}={\dot{V}}_{d,i}\times \rho \times c_p\left(T_{ul,d,i}-T_{ll,d,i}\right) \forall i$$

(8)

The model distinguishes between boilers powered by electricity and boilers using fuel combustion. First, in (9), the costs for a respective EHB unit $b$ are calculated from the associated heat supply ${\dot{Q}}_{b,i}$, the thermal efficiency ${\eta }_{th,b}$, the electricity price $c_{el,i}$ and costs for power purchase $c_{pp}$, and the heat-related $c_{mh,b}$ and hourly $m_b$ maintenance costs. Taking into consideration the interval length $\mathsf{\Delta }I$ and the binary operating variable $X_{b,i}$, which allows the shutdown and starting of the plants, the thermal efficiency (10) is determined by the nominal heat load ${\dot{Q}}_{max,b}$ and nominal input load, which is the input power $P_{max,b}$ for EHB units.

$$C_{b,i}=\mathsf{\Delta }I\times X_{b,i}\left({\dot{Q}}_{b,i}\left({\displaystyle \frac{c_{el,i}+c_{pp}}{{\eta }_{th,b}}}+c_{mh,b}\right)+m_b\right) \forall b,i$$

(9)

$${\eta }_{th,b}={\displaystyle \frac{{\dot{Q}}_{max,b}}{P_{max,b}}} \forall b$$

(10)

Equation (11) shows, in parallel, the hourly costs $C_{h,i}$ for a HOB unit $h$. Different from the EHB, the driving costs are calculated by the fuel-related costs instead of electricity. The specific fuel cost $c_{f,h,i}$, and the costs for total emissions, based on the specific emissions for the fuel ${\epsilon }_h$ and the specific emission costs $e_{h,i}$, have to be considered. In addition, also the interval length, the binary operating variable $X_{h,i}$, the heat supply ${\dot{Q}}_{h,i}$, the thermal efficiency ${\eta }_{th,h}$, the heat-related maintenance cost $c_{mh,h}$, and hourly maintenance cost $m_h$ form the variables and parameters. As in (10) the thermal efficiency for the HOB (12) depends on the nominal heat load ${\dot{Q}}_{max,h}$ and the input load, as well as the fuel consumption ${\dot{Q}}_{f,h}$.

$$C_{h,i}=\mathsf{\Delta }I\times X_{h,i}\left({\dot{Q}}_{h,\mathrm{i}}\left({\displaystyle \frac{c_{f,h,i}+e_{h,i}{\epsilon }_h}{{\eta }_{th,h}}}+c_{mh,h}\right)+m_h\right) \forall h,i$$

(11)

$${\eta }_{th,h}={\displaystyle \frac{{\dot{Q}}_{max,h}}{{\dot{Q}}_{f,h}}} \forall h$$

(12)

The costs for HP $C_{l,i}$ are calculated according to (13) using the same variables and parameters as mentioned in detail before, adapted to the particular HP $l$. Furthermore, the coefficient of performance ${COP}_{l,i}$ and potential electricity-related maintenance costs $c_{me,l}$ are considered, while the ${COP}_{l,i}$ is calculated, simplified by the nominal heat load ${\dot{Q}}_{max,l}$ and the nominal power consumption $P_{max,l}$ shown in (14).

$$C_{l,i}=\mathsf{\Delta }I\times X_{l,i}\left({\dot{Q}}_{l,i}\left({\displaystyle \frac{c_{el,i}+c_{pp}+c_{me,l}}{{COP}_{l,i}}}+c_{mh,l}\right)+m_l\right) \forall l,i$$

(13)

$${COP}_{l,i}={\displaystyle \frac{{\dot{Q}}_{max,l}}{P_{max,l}}} \forall l$$

(14)

In (15), the thermal efficiency of the CHP ${\eta }_{th,k}$, calculated by the nominal load ${\dot{Q}}_{max,k}$ and the fuel consumption ${\dot{Q}}_{f,k}$ in (16), and the electric efficiency ${\eta }_{el,k}$, calculated by the power output $P_{max,k}$ of the CHP in (17), are needed for the determination of the CHP costs $C_{k,i}$.

$$C_{k,i}=\mathsf{\Delta }I\times X_{k,i}\left({\dot{Q}}_{k,i}\left({\displaystyle \frac{c_{th,k,i}+{\eta }_{el,k}{c}_{el,k,i}}{{\eta }_{th,k}}}+c_{mh,k,i}\right)+m_{k,t}\right) \forall k,i$$

(15)

$${\eta }_{th,k}={\displaystyle \frac{{\dot{Q}}_{max,k}}{{\dot{Q}}_{f,k}}} \forall k$$

(16)

$${\eta }_{th,k}={\displaystyle \frac{P_{max,k}}{{\dot{Q}}_{f,k}}} \forall k$$

(17)

Furthermore, all fuel- or heat-related costs are summarized in (18) and all electricity-related costs are calculated from (19). For high electricity prices, the term of the electricity-related costs can be negative, associated with the revenue from feeding electricity into the public power grid.

$$c_{th,k,i}=c_{f,k,i}+e_{k,i}{\epsilon }_k \forall k,i$$

(18)

$$c_{el,k,i}=c_{me,k,i}-c_{el,i} \forall k,i$$

(19)

In addition to the objective function, constraints are required to define the entire MILP problem so that the optimizer can appropriately set the decision variables. One of the most important constraints is the energy balance (20). On the basis of each simulation interval, the balance of all heat supply and heat demand flows has to be maintained. The left-hand side of the equation shows the heat demand parameter ${\dot{Q}}_{de,i}$ of the DHS. In each simulation step $i$, this heat load must be provided by the ST heat, $A·{qa}_i$, the charging (negative) or discharging (positive) of the ATES, $D·{qd}_i$, the heat supply of controllable plants, EHB, HOB, HP, and CHP, and the charging (negative) or discharging (positive) of the short-term TES, ${\dot{Q}}_{s,i}$.

$${\dot{Q}}_{de,i}=A\times {qa}_i+D\times {qd}_i+\sum _{b=1}^B{X}_{b,i}{\dot{Q}}_{b,i}+\sum _{h=1}^H{X}_{h,i}{\dot{Q}}_{h,i}+\sum _{l=1}^L{X}_{l,i}{\dot{Q}}_{l,i}+\sum _{k=1}^K{X}_{k,i}{\dot{Q}}_{k,i}+\sum _{s=1}^S{\dot{Q}}_{s,i} \forall i$$

(20)

All controllable plants, EHB, HOB, HP, and CHP, have individual upper and lower load limits defined in (21)– (24). The heat supplied by the units cannot exceed the nominal load ${\dot{Q}}_{max}$. At the same time, many systems rarely operate below a limited load ${\dot{Q}}_{min}$ for a very short time (when ramping up or down). Therefore, the main operation is restricted by this value.

$${\dot{Q}}_{min,b}\le {\dot{Q}}_{b,i}\le {\dot{Q}}_{max,b} \forall b,i$$

(21)

$${\dot{Q}}_{min,h}\le {\dot{Q}}_{h,i}\le {\dot{Q}}_{max,h} \forall h,i$$

(22)

$${\dot{Q}}_{min,l}\le {\dot{Q}}_{l,i}\le {\dot{Q}}_{max,l} \forall l,i$$

(23)

$${\dot{Q}}_{min,k}\le {\dot{Q}}_{k,i}\le {\dot{Q}}_{max,k} \forall k,i$$

(24)

Further important restrictions are made for short-term TES. First, the storage level of TES ${QS}_{s,i}$ has to be below the maximal storage capacity ${QSC}_s$. As also shown in (25), the storage level is calculated by the storage level of the previous interval ${QS}_{s,(i-1)}$, the interval length $\mathsf{\Delta }I$, the charging (negative) or discharging (positive) of TES ${\dot{Q}}_{s,i}$, and the thermal losses ${\dot{Q}}_{lo,s,(i-1)}$. The thermal losses are computed on the previous interval’s storage level, its temperature, and other given parameters of TES, such as the geometry, location, ambient temperature, insulation material, or insulation layer thickness.

$${QS}_{s,i}={QS}_{s,(i-1)}-\mathsf{\Delta }I\left({\dot{Q}}_{s,i}+{\dot{Q}}_{lo,s,(i-1)}\right)\le {QSC}_s \forall s,i$$

(25)

Similar to the DHS heat supplies of the controllable plants, the charging (negative) and discharging (positive) of TES ${\dot{Q}}_{s,i}$ is limited in (26) by the maximal charging load ${\dot{Q}}_{cc,max,s,i}$ and the maximal discharging load ${\dot{Q}}_{cd,max,s,i}$.

$${\dot{Q}}_{cc,max,s,i}\le {\dot{Q}}_{s,i}\le {\dot{Q}}_{cd,max,s,i} \forall s,i$$

(26)

Equation (27) represents the law of energy conversion. The TES storage level at the end of the simulation horizon ${QS}_{s,I}$ has to equal the starting storage level ${QS}_{s,0}$.

$${QS}_{s,I}={QS}_{s,0} \forall s$$

(27)

The last step for the MILP problem definition is the summary of all variables used in the mathematical model: the DHS heat supply by all controllable plants, shown in (28), their binary operating variables in (29), and TES charging or discharging, the TES storage level, and the ST and ATES dimensioning in (30). All continuous and discrete variables of the problem are summarized in

Table A1. Continuous and discrete variables of the problem.

Continuous	Discrete
$A\in {\mathbb{R}}^{+}$: collector area of ST plant	$X_{b,i}\in \left(\mathrm{0,1}\right)$: binary operating variable of electric heat boilers
$C_{a,i}\in {\mathbb{R}}^{+}$: operating costs of ST	$X_{h,i}\in \left(\mathrm{0,1}\right)$: binary operating variable of heat only boilers
$C_{b,i}\in {\mathbb{R}}^{+}$: operation costs of EHB	$X_{k,i}\in \left(\mathrm{0,1}\right)$: binary operating variable of combined heat and power plants
$C_{d,i}\in {\mathbb{R}}^{+}$: operation costs of ATES	$X_{l,i}\in \left(\mathrm{0,1}\right)$: binary operating variable of heat pumps
$C_{h,i}\in {\mathbb{R}}^{+}$: operation costs of HOB
$C_{k,i}\in {\mathbb{R}}^{+}$: operation costs of CHP
$C_{l,i}\in {\mathbb{R}}^{+}$: operation costs of HP
$C_{op}\in {\mathbb{R}}^{+}$: total operation costs
${\dot{Q}}_{b,i}\in {\mathbb{R}}^{+}$: heat supply of EHB
${\dot{Q}}_{h,i}\in {\mathbb{R}}^{+}$: heat supply of HOB
${\dot{Q}}_{k,i}\in {\mathbb{R}}^{+}$: heat supply of CHP
${\dot{Q}}_{l,i}\in {\mathbb{R}}^{+}$: heat supply of HP
${\dot{Q}}_{s,i}\in {\mathbb{R}}^{+,-}$: charging (negative) or discharging (positive) of the short-term TES
${QS}_{s,t}\in {\mathbb{R}}^{+}$: storage level of the TES

(Appendix A).

$${\dot{Q}}_{b,i}\in {\mathbb{R}}^{+}; {\dot{Q}}_{h,i}\in {\mathbb{R}}^{+}; {\dot{Q}}_{l,i}\in {\mathbb{R}}^{+}; {\dot{Q}}_{k,i}\in {\mathbb{R}}^{+} \forall i,b,h,l,k$$

(28)

$$X_{b,i}\in \left(\mathrm{0,1}\right); X_{h,i}\in \left(\mathrm{0,1}\right); X_{l,i}\in \left(\mathrm{0,1}\right); X_{k,i}\in \left(\mathrm{0,1}\right) \forall i,b,h,l,k$$

(29)

$${\dot{Q}}_{s,i}\in {\mathbb{R}}^{+,-}; {QS}_{s,i}\in {\mathbb{R}}^{+}; A\in {\mathbb{R}}^{+}; D\in {\mathbb{N}}_0 \forall i,s$$

(30)

Based on this mathematical scheduling problem description, a simulation model was developed in Python and solved with the Gurobi 11.0.0 optimizer for different frameworks, e.g., optimized or defined number of ATES drillings.

In order to compare the simulation results economically and ecologically, the specific heating costs $c_{op}$ (31) based on the total operating costs $C_{op}$ and the total heat demand $Q_{de}$ and the specific emissions of the energy supply ${\epsilon }_{op}$ in (32) are used. To determine the specific emissions, the total emissions ${\mathrm{E}}_{op}$ are calculated by the emissions of each interval step ${\mathrm{E}}_{op,i}$.

$$c_{op}={\displaystyle \frac{C_{op}}{Q_{de}}}$$

(31)

$${\epsilon }_{op}={\displaystyle \frac{{\mathrm{E}}_{op}}{Q_{de}}}={\displaystyle \frac{\sum _{i=1}^I{\mathrm{E}}_{op,i}}{Q_{de}}}$$

(32)

(33) shows the emissions for each simulation interval ${{\rm E}}_{op,i}$, while the emissions for each technology have to be calculated based on their power (ATES, EHB, and HP) or fuel (HOB and CHP) consumption and the specific electricity emissions ${\epsilon }_{el,i}$ or the specific fuel emissions ${\epsilon }_f$.

$${{\rm E}}_{op,i}=\mathsf{\Delta }I\left({\epsilon }_{el,i}\left(D\times {pd}_i+\sum _{b=1}^B{P}_{b,i}+\sum _{l=1}^L{P}_{l,i}\right)+{\epsilon }_f\left(\sum _{h=1}^H{\dot{Q}}_{f,h,i}+\sum _{k=1}^K{\dot{Q}}_{f,k,i}\right)\right) \forall i,b,l,h,k$$

(33)

2.2. Data Hierarchical Clustering Approach

Regarding the various framework conditions for the simulations, the computational times of the optimization problem should not exceed a reasonable duration. By considering seasonal TES through the ATES investigations, time horizons of at least one year have to be solved. Usually, short-term scheduling or unit commitment problems based on operating objective functions are applied for several hours or days [33].

To meet this challenge, available hourly data will be clustered according to the hierarchical method described below. The method combines arithmetic average determination with descending data sorting. The goal is to obtain an averaged data distribution without heavy flattening of positive or negative peaks.

First, the original data ${data\_o}_i$ are sorted in a descending order ${data\_sort}_i$. Both are the basis of the clustering. The arithmetic average of the interval steps $i$ is now calculated from these data, and the size of the dataset is reduced by the clustering factor $fc$. Equation (34) determines the arithmetic average of the original data, ${data\_o\_ar}_i$, and reduces the number of intervals within the simulation horizon ($i=1,\dots ,I$).

$${data\_o\_ar}_i={\displaystyle \frac{1}{fc}}\sum _{n=1+\left(i-1\right)·fc}^{i·fc}{data\_o}_n \forall i$$

(34)

For an original dataset of length $I_{data\_o}=8760$ (when hourly data are available for the simulation) and an interval length of $fc=8$, the new averaged dataset has a length of $I_{data\_o\_ar}=1095$, because for $i=1095$, the final value of the sum sign reaches the maximal value for the original dataset ($i\times fc=1095\times 8=8760$).

The average sorted values ${data\_sort\_ar}_i$ are calculated similarly in (35).

$${data\_sort\_ar}_i={\displaystyle \frac{1}{fc}}\sum _{n=1+(i-1)\times fc}^{i\times fc}{data\_sort}_n \forall i$$

(35)

Both resulting new datasets, ${data\_o\_ar}_i$ and ${data\_sort\_ar}_i$, have (as requested) a lower number of containing values than the original sets, and the total average value is not influenced. The averaged original data represent the general chronological sequence of the data, while the averaged sorted original data have a better representation of the individual values (higher maximal and lower minimal values). In order to distribute the better-fitting individual values based on the averaged arrangement, the averaged original data must also be sorted, ${data\_o\_ar\_sort}_i$. Next, the assignment is made according to (36).

$${data\_o\_ar}_i={data\_o\_ar\_sort}_j⟹{data\_fc}_i:={data\_sort\_ar}_j \forall i,j$$

(36)

The final clustered data ${data\_fc}_i$ of a specific interval $i$ is equal to the specific value $j$ of the averaged sorted data ${data\_sort\_ar}_j$. The position $j$ of the value is identified via the search of the averaged original data ${data\_o\_ar}_i$ within the descending sorted list of these values ${data\_o\_ar\_sort}_j$.

2.3. Application Case and Assumptions

The described clustering procedure is applied to the data of the use case for the selected reference year 2023. The DHS of a German utility company Blockheizkraftwerks-, Träger- und Betreibergesellschaft Berlin (BTB) represents the use case. All associated data of the use case and further assumptions are presented in this section.

Table 1. Total nominal loads (electricity $P_{max}$, DHS heat supply ${\dot{Q}}_{max}$, and fuel consumption ${\dot{Q}}_f$), energy source, and units of the controllable DHS plants.

Index	Supply Plant			Total Nominal Loads (MW)
	Type	Units	Source	${\mathit{P}}_{\mathit{m}\mathit{a}\mathit{x}}$	${\dot{\mathit{Q}}}_{\mathit{m}\mathit{a}\mathit{x}}$	${\dot{\mathit{Q}}}_{\mathit{f}}$
CHP 1	TPS	1	Biomass	20	30	55
CHP 2	Engine	9	Natural gas	2.93	3.08	6.59
CHP 3	Engine	1	Natural gas	0.42	0.51	1.07
CHP 4	Engine	1	Natural gas	0.8	0.86	1.89
HOB 1	HOB	2	Natural gas	-	18.5	20.7
HOB 2	HOB	3	Natural gas	-	33	36.4
HOB 3	HOB	2	Natural gas	-	20	20.62
HOB 4	HOB	1	Natural gas	-	10	10.42
HOB 5	HOB	2	Natural gas	-	2	2.2
HOB 6	HOB	10	Natural gas	-	2.76	3.06
HOB 7	HOB	1	Natural gas	-	0.96	1.42
EHB 1	EHB	2	Electricity	3.3	3.3	-
EHB 2	EHB	1	Electricity	0.44	0.44	-
HP 1	HP	2	River water	1.72	4.3	-
HP 2	HP	1	Process heat	0.19	1.03	-

shows the energy source and total nominal loads for the controllable DHS supply plants of the application case.

The biomass combusted in CHP 1 is waste wood. The power-to-heat (PTH) plants are supplied with electricity from the public power grid. The heat source of HP 1 is the Spree, a river in Berlin whose temperature must be at least 8 °C to operate the HP, and the heat source of HP 2 is the process heat from CHP 3. Furthermore, the hourly maintenance costs for the CHP plants are between 12 EUR/h and 30 EUR/h, and 15 EUR/h for all HPs. The thermal-related maintenance costs are assumed to be 0–10 EUR/MWh for the HOB and EHB and 1 EUR/MWh for the HP. Electricity-related maintenance costs are 13 EUR/MWh for CHP 1. The DHS system also includes short-term TES with a total capacity of 121 MWh at 2.8 MW of maximal charging or discharging loads.

In the following, hourly use case data are presented. To shorten the calculation duration, a clustering factor $fc$ of 8 is applied in the MILP optimization, taking Equations (34)–(36) into account, so that the final simulation intervals have a length of 8 h.

The hourly and clustered thermal heat load of BTB’s DHS is shown in Figure 1. The average heat load in 2023 is equal to 67.5 MW, while the peak load is 181 MW and the average base load in summer is about 24 MW. The total heat demand in 2023 is 591 GWh.

Figure 2 illustrates the ambient air temperatures in Berlin and the river water temperatures of the Spree for the reference year 2023. The air temperature data are provided by the German weather service DWD [35], and the river water temperatures are recorded by measuring stations on the Spree [36]. The clustered air temperature time series (red line) moderately attenuates the maximum and minimum peaks observed in the original time series (black line) while maintaining a high degree of fidelity. In contrast, simple averaging over 8-h intervals would result in a greater flattening of the curve. River water temperature, on the other hand, shows relatively low temporal variability with minimal pronounced peaks. As a result, the clustered time series (blue line) closely follows the original data (light grey line), with negligible divergence.

The average air temperature in 2023 is 11.5 °C and the average river temperature is 13.1 °C. The river temperature follows the air temperature except slightly shifted and with flattened peaks. On the basis of these data, the losses of short-term TES, the COP of the river water HP, or the solar yield of the ST plant are affected.

The day-ahead electricity price of the German power grid and the related carbon emissions are shown in Figure 3 [37,38]. In 2023, the average electricity price is 95.18 EUR/MWh [37], while the average emissions related to the electricity mix are about 395 kg/MWh [38].

Additionally, PTH plants (EHB and HP) that use electricity from the public power grid have to pay taxes, surcharges, concession fees, and location-dependent grid fees. These specific procurement costs amounted to 65.83 EUR/MWh for Berlin in 2023.

Further assumptions are made regarding the wholesale prices and emission factors of biomass and natural gas [39,40,41], summarized in Table 2. Smaller plants have to pay for their emissions according to the German national emission trading system (ETS) nEHS based on the fuel emission trading act BEHG, and larger plants are assigned to the European ETS. It is assumed that these systems will converge in the long term and, therefore, only a single price for the carbon dioxide (CO₂) emissions of 50 EUR/MWh is assumed within the simulations [42,43].

The solar yield of an exemplary ST collector, illustrated in Figure 4, is calculated on the basis of weather data from the DWD (which also includes the solar irradiation on the horizontal plane or the wind speeds). The observed total solar irradiation in the area of Berlin was 1139 kWh/m² in the year 2023 [35].

To simplify the model and speed up the calculation duration, this calculation is performed outside the MILP optimization using the Solar-Keymark-Output-Calculator (ScenoCalc) of the research institute Solites [44,45]. Within ScenoCalc, a high-temperature vacuum tube collector, model Vitosol 200-T SPX-S [46] from the company Viessmann Group GmbH & Co. KG, Allendorf (Eder), Germany, was selected and simulated for Berlin, the weather conditions of 2023, and with a collector inclination of 20°. As determined by ScenoCalc, the solar irradiation on the inclined collector surface is equal to 1206 kWh/m² and the resulting solar yield per collector area is 433.5 kWh/m². The selected collector results in a system efficiency of 38%.

The results of the ScenoCalc determination and further assumptions for the ST and ATES systems are summarized in Table 2 [42,45,46,47].

All assumptions made for the ATES system are based on Hebold et al. [47], a feasibility study executed by the German company Geothermie Neubrandenburg GmbH, Neubrandenburg, Germany, and the research institute GFZ Helmholtz Centre for Geosciences, Potsdam, Germany. This study includes the development of potential technical implementation concepts for the integration of a seasonal ATES facility into the DHS operated by BTB at the Adlershof side, a part of BTB’s DHS. First, the study assesses the hydro-geological conditions and then develops suitable drilling concepts. One of the concepts assumes a higher injection temperature (HIT) of about 55 °C, so that the underground water at the injection well also heats up during the lifecycle of the ATES, and another concept is based on a lower injection temperature (LIT). The LIT is close to the average underground temperature of 20 °C and can only be achieved with an HP. In both cases, the warm well of the ATES is charged with a temperature of 95 °C. The total power demand for the ATES deep groundwater well pumping is determined as 1660 MWh/a, which results in an electricity load for the pumps of 18.95 kW [47].

Assuming a heat exchanger temperature difference of 3 K, exergy efficiency of 50% [34,48], an ATES high-production temperature of 95 °C, and a low injection temperature of 20 °C for the LIT-ATES or 55 °C for the HIT-ATES, the HIT-COP is calculated by Equation (6) to be 2.09, while it is 3.4 for the LIT-COP when the targeted HP outlet temperature is 105 °C, which is close to the DHS supply temperature in high-demand seasons.

Based on the assumptions made, four main portfolio scenarios were selected for the simulations:

• S: ST (solar yield of 433.5 kWh/m²/a) without ATES.
• A: ST with HIT-ATES (injection well temperature of 55 °C).
• B: ST with HIT-ATES and HP (DHS supply temperature of 105 °C, COP of 3.4).
• C: ST with LIT-ATES (injection well temperature of 20 °C) and HP (COP of 2.09).

While the controllable plant operation, short-term TES, and collector area of the ST system are optimized within the MILP model, the dimensioning of the ATES can be both optimized or set. The background is that optimization via the operating cost result does not necessarily lead to the use of an ATES, but the effect of different ATES capacities on the optimal dimensioning of the ST should still be investigated. The ATES number $D$ can therefore either be a decision variable or a fixed parameter. The dimensioning of the ATES system is added to the symbol of the respective scenario in the results shown in the Section 3. Hereby, the cost optimum is marked by an asterisk.

Figure 5 shows the volumetric flow rate and the upper and lower temperature levels of the ATES scenarios, assumed according to the aforementioned feasibility study.

In general, the upper temperature is kept constant during the charging period in the middle of the year, while the lower temperature level is approximately constant during the discharging period from 15 October to 14 April. The available supply temperature during the discharging period declines from 95 °C to 65 °C due to heat losses to the surrounding ground for both the HIT- and LIT-ATES scenarios. The lower temperature declines from 55 °C to 40 °C during the charging period of the HIT-ATES. While the LIT-ATES’ cold well is always close to the surrounding ground temperature, there is no decline in the temperature. The limiting factor for the actual ATES flow rate, which depends on the permeability of the water-bearing layer, is the maximum filtration capacity of the built-in filter and gravel pack. It is equal to 140 m³/h [47] for one ATES system consisting of one hot well and one cold well.

Since the volume flow during the charging period is defined as negative (the flow direction from injection to the production well of the ATES), the resulting available heat source of the ATES, shown in Figure 6, is also negative. For scenarios A and B, the heat source load of the ATES is identical (consequently, the solid blue line overlaps with the solid black line), and for scenario C (solid red line), its absolute value is higher (for both charging and discharging) due to the higher temperature difference between the injection and production wells. For A and B, the heat load of the thermal source decreases from 6.5 to 1.6 MW while discharging and the absolute heat demand for the charging arises from 6.5 to 8.9 MW. The discharging heat source load for scenario C decreases from 12.2 to 7.3 MW, while the resulting heat load during charging is constant at 12.2 MW. Normally, the resulting total available or needed heat load of the ATES system is equal to the heat source load and the power demand is equal to the deep groundwater pumping electricity demand, as observed in Scenario A. The electricity demand for groundwater pumping is only slightly above zero (black dashed line) and the resulting available heat load (black dotted line) corresponds exactly to the heat source. Consequently, the black dotted line is perfectly aligned with both the black and blue solid lines. However, for scenarios B and C there, is an additional HP. The power load increases in the discharging period and so does the resulting available total heat load, whereby in scenario B, it is up to 9.2 MW, and in C, its peak load is about 23.3 MW.

The energy flow chart in Figure 7 illustrates the interactions between the components of the overall system. All heat generators contribute to meeting the heat demand of the district heating system (DHS), represented by the red flows. The controllable heat generators (HP, EHB, CHP, and HOB) are powered by their respective energy sources. The HP and EHB units use electricity from the grid (green flows), the HOB units use only natural gas (yellow flows), and the CHP units use natural gas or biomass (brown flow) depending on their configuration. In addition to generating heat, CHP units also generate electricity, which is fed back into the electricity grid.

ST is the only non-controllable component that supplies heat to the DHS. Each scenario also includes a short-term TES, which can function both as a heat generator and a consumer, provided its storage capacity allows it. The optimization model determines when thermal energy is stored or released by the TES.

A unique feature of the system is the incorporation of seasonal ATES, which can be implemented in two primary configurations, highlighted in light blue in the diagram. In scenario A, the ATES operates at a high injection temperature, eliminating the need for an additional HP to lower the return temperature during the discharge phase. Conversely, the configuration used in scenarios B and C includes an additional HP to assist in discharging the ATES and providing heat to the DHS. This configuration allows a higher flow temperature to the heating network while maintaining a lower injection temperature to the ATES, optimizing system performance.

Table 2 summarizes all the assumptions for the parameters of the optimization problem for a better overview and completes the basic assumptions for the calculations, such as the solar irradiation or the temperature difference of the heat exchangers to determine the evaporation and condensation temperature of the ATES-HP.

Table 2. Assumptions for parameters within the optimization problem.

Index	Parameter	Value	Unit
$B$	Number of EHB units	2	-
$c_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}$	Specific fuel cost, wholesale price biomass [40]	20	EUR/MWh
$c_{\mathrm{g}\mathrm{a}\mathrm{s}}$	Specific fuel cost, wholesale price natural gas [39]	37.5	EUR/MWh
$c_{me,k}$	Electricity-related maintenance costs (CHP)	0–13	EUR/MWh
$c_{mh}$	Heat-related maintenance costs	0–10	EUR/MWh
$c_p$	Heat capacity of water (around 80 °C)	4193.7	J/(kgK)
$c_{pp}$	Costs for the purchase from the power grid	65.83	EUR/MWh
${COP}_{d,HIT}$	COP of HP of HIT-ATES system (B)	3.4	-
${COP}_{d,LIT}$	COP of HP of LIT-ATES system (C)	2.09	-
$D$	Number of drillings of the ATES system (injection and production)	0–7	-
$e$	Specific emission costs	50 (150)	EUR/MWh
${\epsilon }_{biomass}$	Specific emissions for biomass [41]	27	kg CO2/MWh
${\epsilon }_{gas}$	Specific emissions for natural gas [41]	201	kg CO2/MWh
$H$	Number of HOB units	7	-
$I$	Simulation horizon	1095	-
$\mathsf{\Delta }I$$, fc$	Simulation interval, clustering interval	8	h
$K$	Number of CHP units	4	-
$L$	Number of heat pumps	2	-
$m_k$	Hourly maintenance costs (CHP)	12–30	EUR/h
$m_l$	Hourly maintenance costs (HP)	15	EUR/h
${\eta }_{ex}$	Exergy efficiency of ATES-HP [34,48]	50	%
${pd}_{dd,d,i}$	Power demand for ATES deep pumping	18.95	kW
-	Total solar irradiation, Berlin 2023 [35]	1139	kWh/m2/a
-	Total solar irradiation on ST collector area	1206	kWh/m2/a
$\sum {qa}_{a,i}\mathsf{\Delta }I$	Solar yield per ST collector area	433.5	kWh/m2/a
${\dot{Q}}_{cc,max}$	Maximal charging load	2.8	MW
${\dot{Q}}_{cd,max}$	Maximal discharging load	2.8	MW
$Q_{de}$	Heat demand of the DHS	591,000	MWh
${\displaystyle \frac{{\dot{Q}}_{min}}{{\dot{Q}}_{max}}}$	Minimal partial load	0–13	%
${QSC}_s$	Maximal storage capacity	121	MWh
$\rho$	Density of water (around 80 °C)	972.8	kg/m3
-	Temperature difference of heat exchangers	3	K
$T_H$	Hot condensation temperature of ATES-HP	108	°C
$T_{C,HIT}$	Cold evaporation temperature HIT-ATES-HP	52	°C
$T_{C,LIT}$	Cold evaporation temperature LIT-ATES-HP	17	°C

Table 2. Assumptions for parameters within the optimization problem.

Index	Parameter	Value	Unit
$B$	Number of EHB units	2	-
$c_{\mathrm{b}\mathrm{i}\mathrm{o}\mathrm{m}\mathrm{a}\mathrm{s}\mathrm{s}}$	Specific fuel cost, wholesale price biomass [40]	20	EUR/MWh
$c_{\mathrm{g}\mathrm{a}\mathrm{s}}$	Specific fuel cost, wholesale price natural gas [39]	37.5	EUR/MWh
$c_{me,k}$	Electricity-related maintenance costs (CHP)	0–13	EUR/MWh
$c_{mh}$	Heat-related maintenance costs	0–10	EUR/MWh
$c_p$	Heat capacity of water (around 80 °C)	4193.7	J/(kgK)
$c_{pp}$	Costs for the purchase from the power grid	65.83	EUR/MWh
${COP}_{d,HIT}$	COP of HP of HIT-ATES system (B)	3.4	-
${COP}_{d,LIT}$	COP of HP of LIT-ATES system (C)	2.09	-
$D$	Number of drillings of the ATES system (injection and production)	0–7	-
$e$	Specific emission costs	50 (150)	EUR/MWh
${\epsilon }_{biomass}$	Specific emissions for biomass [41]	27	kg CO2/MWh
${\epsilon }_{gas}$	Specific emissions for natural gas [41]	201	kg CO2/MWh
$H$	Number of HOB units	7	-
$I$	Simulation horizon	1095	-
$\mathsf{\Delta }I$$, fc$	Simulation interval, clustering interval	8	h
$K$	Number of CHP units	4	-
$L$	Number of heat pumps	2	-
$m_k$	Hourly maintenance costs (CHP)	12–30	EUR/h
$m_l$	Hourly maintenance costs (HP)	15	EUR/h
${\eta }_{ex}$	Exergy efficiency of ATES-HP [34,48]	50	%
${pd}_{dd,d,i}$	Power demand for ATES deep pumping	18.95	kW
-	Total solar irradiation, Berlin 2023 [35]	1139	kWh/m2/a
-	Total solar irradiation on ST collector area	1206	kWh/m2/a
$\sum {qa}_{a,i}\mathsf{\Delta }I$	Solar yield per ST collector area	433.5	kWh/m2/a
${\dot{Q}}_{cc,max}$	Maximal charging load	2.8	MW
${\dot{Q}}_{cd,max}$	Maximal discharging load	2.8	MW
$Q_{de}$	Heat demand of the DHS	591,000	MWh
${\displaystyle \frac{{\dot{Q}}_{min}}{{\dot{Q}}_{max}}}$	Minimal partial load	0–13	%
${QSC}_s$	Maximal storage capacity	121	MWh
$\rho$	Density of water (around 80 °C)	972.8	kg/m3
-	Temperature difference of heat exchangers	3	K
$T_H$	Hot condensation temperature of ATES-HP	108	°C
$T_{C,HIT}$	Cold evaporation temperature HIT-ATES-HP	52	°C
$T_{C,LIT}$	Cold evaporation temperature LIT-ATES-HP	17	°C

3. Results

Using the developed MILP optimization model, the previously defined scenarios were solved for the minimum total cost of plant operation. First, the scenario of optimizing the ST supply to the DHS without an ATES (S) was performed as a reference for the following results. Then, different types of ATES configurations (A, B, and C) were studied, focusing on the effects of the unit commitment of controllable plants with different ATES sizing. Specifically, for scenario A, a single well (A1) and up to seven wells (A7) were investigated because this ATES has the lowest heat load per drilling. For the ATES configurations with additional HP, the cost of electricity leads to the economic deterioration of the operation, although only the effects of a single drilling (two wells) of the ATES system are presented. The goal of the optimization is not the number of wells or the dimensioning of the ATES system, but rather the investigation of the optimal operation of the existing DHS portfolio under the new conditions (integrated ST and ATES).

3.1. Scenario S: ST Without ATES

The reference for the following results is scenario S or S0 without considering any seasonal TES. Figure 8 shows the cost-optimized heat generation of all controllable suppliers and the non-controllable ST plant with the framework conditions defined in preceding Section 2. The heat demand of the DHS (red line) must be met accurately on an hourly basis, but the actual heat supply may differ from the heat demand. Any resulting discrepancy must be compensated by charging or discharging the TES. The portfolio of heat generators is represented by different colors: light green for ST, dark green for biomass CHP, various shades of orange for natural gas CHP, various shades of grey for HOB, violet for EHB and blue for HP. The absence of certain colors in the diagram indicates that these plants, although included in the portfolio, are not part of the optimal solution. This result is consistent with the intended analysis.

The scenario results in optimal dimensioning for ST of 35,246.2 m², producing 15.3 GWh heat per year, which is a share of 2.6%. The share of renewables in total is due to the biomass CHP 1 with 40.8% for an existing DHS, which is relatively high. EHB and HP only account for about 1.1%, while conventional heat producers account for the highest share. Moreover, 23.6% of the heat is produced by the other natural gas-fired CHP and 35.7% by HOB. The total annual emissions for scenario A are about 115,500 t CO₂ and the optimal costs are equal to 4.9 million EUR. These values result in specific emissions of 195.4 kg/MWh CO₂ and specific heating costs of 8.3 EUR/MWh for this specific DHS scenario. All of the results are summarized in

Table A2. DHS heat loads, CO₂ emissions, and operation costs resulting for all scenarios (A, B, C).

	S0	A1	A2	A5*	B1	B2	C1	C2
DHS heat supply in GWh/a
ST	15.29	42.53	69.98	124.98	42.49	69.98	57.73	91.81
CHP 1	219.36	224.94	227	235.86	224.96	226.67	224.88	222.73
CHP 2–4	139.2	140.75	144.04	162.25	140.74	143.34	140.4	135.22
HOB	210.7	191.99	174.94	141.11	184.64	161.69	133.51	80.41
EHB	0.11	0.11	0.11	0.11	0.11	0.11	0.11	0.11
HP	6.4	6.86	7.22	7.36	6.87	6.78	6.52	4.95
ATES	0	−16.12	−32.24	−80.61	−8.76	−17.51	27.91	55.82
Surplus 1	+0.03	+16.15	+32.28	+80.64	+8.79	+17.54	−27.88	−55.79
Heat production shares in %
ST	2.59	7.19	11.84	21.15	7.19	11.84	9.77	15.53
CHP 1	37.11	38.06	38.41	39.91	38.06	38.35	38.05	37.68
CHP 2–4	23.55	23.81	24.37	27.45	23.81	24.25	23.75	22.88
HOB	35.65	32.48	29.6	23.87	31.24	27.36	22.59	13.6
EHB	0.02	0.02	0.02	0.02	0.02	0.02	0.02	0.02
HP	1.08	1.16	1.22	1.25	1.16	1.15	1.1	0.84
Renewable	40.8	46.43	51.49	62.33	46.43	51.36	48.94	54.07
CO2 emissions in thousand t CO2
CHP 1	10.86	11.13	11.24	11.68	11.14	11.22	11.13	11.03
CHP 2–4	59.9	60.57	61.98	69.82	60.56	61.68	60.41	58.18
HOB	44.24	40.23	36.57	29.33	38.65	33.74	27.79	16.67
EHB	0.02	0.02	0.02	0.02	0.02	0.02	0.02	0.02
HP	0.48	0.53	0.55	0.57	0.53	0.52	0.5	0.38
ATES	0	0.66	1.31	3.28	3.8	7.59	17.27	34.53
Total	115.5	113.13	111.67	114.68	114.7	114.78	117.11	120.81
Specific CO2 emissions in kg/MWh CO2
CHP 1	49.51	49.48	49.52	49.52	49.52	49.5	49.49	49.52
CHP 2–4	430.32	430.34	430.3	430.32	430.3	430.31	430.27	430.26
HOB	209.97	209.54	209.04	207.85	209.33	208.67	208.15	207.31
EHB	181.82	181.82	181.82	181.82	181.82	181.82	181.82	181.82
HP	75	77.26	76.18	77.45	77.15	76.7	76.69	76.77
Total	195.42	191.41	188.94	194.03	194.07	194.2	198.15	204.41
Operation costs in million EUR
ST	0.05	0.13	0.21	0.37	0.13	0.21	0.17	0.28
CHP 1	−4.44	−4.43	−4.44	−4.56	−4.43	−4.45	−4.5	−4.73
CHP 2–4	−1.38	−1.37	−1.36	−1.23	−1.37	−1.39	−1.43	−1.57
HOB	10.47	9.52	8.66	6.94	9.15	7.99	6.58	3.95
EHB	−0.01	−0.01	−0.01	−0.01	−0.01	−0.01	−0.01	−0.01
HP	0.24	0.25	0.27	0.27	0.25	0.25	0.23	0.18
ATES	0	0.27	0.53	1.34	1.48	2.95	6.71	13.43
Total	4.93	4.35	3.86	3.12	5.19	5.54	7.77	11.52
Specific operation costs in EUR/MWh
ST	3	3	3	3	3	3	3	3
CHP 1	−20.22	−19.7	−19.56	−19.33	−19.7	−19.65	−19.99	−21.24
CHP 2–4	−9.9	−9.77	−9.46	−7.59	−9.77	−9.7	−10.17	−11.59
HOB	49.7	49.6	49.48	49.2	49.56	49.4	49.27	49.06
EHB	−85	−82.89	−85	−85	−82.65	−85	−85	−85
HP	36.76	36.92	37.02	36.07	36.83	36.58	35.96	35.35
Total	8.34	7.37	6.52	5.28	8.78	9.38	13.14	19.48

¹ The DHS heat production surplus compared to the DHS demand of 591 GWh without considering ATES heat supply or demand. * This is the optimal ATES dimensioning scenario in terms of operating costs.

in Appendix A.

3.2. Scenario A: ST with HIT-ATES

In comparison to the reference, scenario A is the first to consider seasonal TES: an ATES with a temperature level of approximately 55 °C for the injection well and a maximal temperature of 95 °C for the production well. For this configuration, the optimal ST plant has a collector area of almost 288,300 m² with five ATES drillings (A5*). The total ATES heat demand during the charging period is therefore equal to 169.1 GWh and the total DHS supply from the ATES during the discharging period is 88.5 GWh. As the system has no extra HP to raise the DHS supply temperature, the ATES power demand is equal to the demand for the deep pumping of 18.95 kW, the ATES demand varies between 32.4 MW and 44.6 MW, the maximum ATES heat supply to the DHS is 32.4 MW at the start of the discharging period, and therefore the total heat load (DHS demand minus ATES heat load) is increased much more in the charging period than the demand is decreased during discharging, as seen in Figure 9. * indicates that this is the optimal ATES dimensioning scenario in terms of operating costs.

This shift in heat demand from DHS and ATES can also be transferred to the representation of the heat production divided among the various supply plants throughout the year as shown in Figure 10. The total ST heat supply is 125 GWh, which is a share of 21.2% of the total heat production. However, at the same time, the production surplus due to storage losses (mainly ATES, not short-term TES) is equal to 80.6 GWh, which is about 13.6% of the DHS heat demand. Nevertheless, the cost-optimal solution for ATES dimensioning leads to an increase in ST heat of 717%. Of course, the production surplus has to be taken into account. If the ATES losses are assumed to be mainly served by ST, an increase in DHS supply by ST of 210% is still achieved.

In the following section, comprehensive investigations into the effects of different ATES dimensioning show that neither the increase nor decrease in a specific producer’s DHS supply is linear with a linear enlargement of the ATES system. Also, the cost-wise optimum does not represent the best solution regarding carbon emissions.

Figure 11 shows a thermal energy split between the different producers and the ATES delta (heat demand and supply). The ST heat increases the most as the ATES size increases. However, the rate of increase per extra drilling reduces at the same time. The increase in ST heat decreases in absolute terms from simulation A5* onwards. CHP, EHP, and HP also tend to deliver more heat to the DHS. The EHP and HP only generate a really small amount of heat due to their installed design and the actual price situation. However, the potential for an increase in HP supply through an ATES is evident. In scenario A1, the relative heat supply of the HP is increased by 7.3% and in A7 by 38.9%. CHP systems also experience an increase of 3.7% (A1) to 38% (A7). Particularly remarkable is the rise in heat production of natural gas-fired CHP systems, which significantly increased from scenario A5* onward (a 6.3% increase from A4 to A5*, while the increase before is about 1–3% per extra drilling). Only the heat yield from the HOB is continuously reduced by 8.9% (A1) over 33% (A5*) to 34.2% (A7). Due to the constant total efficiency of the ATES system, the storage losses (resulting in a production surplus) increase evenly.

The carbon emissions associated with the DHS heat supply are shown in Figure 12. The minimal total emissions result in scenario A3 with about 110,900 t of CO₂. This is a decrease of 4% compared to the reference. The cost-optimal result A5* has a reduction of 0.7%. Other scenarios (A6 or A7) even show an increase in total occurring emissions.

The most notable effect of the ATES sizing variation is summarized in Figure 13. The total DHS operational costs, the ATES-related costs, and the costs of ST heat generation are illustrated. Through seasonal TES, the total operation costs are reduced. As mentioned previously, the optimum is scenario A5* with five ATES systems and total costs of 3.12 million EUR. Compared to the reference, the operation costs are reduced by 36.7%. The allocation of costs to the components is conspicuous: at some point, the ATES costs dominate the costs of the heat producers. For scenario A5*, the ATES costs are almost equal to the costs of the controllable producers (except ST) and the ATES account for almost 50% of the overall operation costs in scenario A7.

3.3. Scenario B: ST with HIT-ATES and HP

After discussing the results of scenario A in detail, the other ATES configurations are also briefly presented, before summarizing the most important variants at the end of the section.

Within scenario B, the previous scenario A is only supplemented by an HP, which raises the DHS supply temperature to 105 °C. Figure 14 shows the results for one drilling (B1). Due to rising ATES costs, the cost-related optimum would be the solution without an ATES system. While the discharging available heat per drilling is 17.7 GWh for scenario A, it is enlarged through the extra HP to 25.1 GWh. For both scenarios, the heat demand during the charging period is equal to 33.8 GWh per drilling. At the same time, the electricity demand with HP is 9 GWh and, therefore, significantly higher than in scenario A. Regardless of the discharging period, which causes the only difference between A and B, the ST heat supply is increased to about 42.5 GWh. In particular, this is about 178% more than the reference and represents a share of 7.2%.

3.4. Scenario C: ST with LIT-ATES and HP

Due to the larger temperature spread in scenario C, the total heat demand and the heat supply per drilling are significantly higher than in the other two scenarios. As shown in Figure 15, the heat demand shift during the charging and discharging period through the ATES is significantly larger than in Figure 14. While the power demand for the HP is almost 40.5 GWh, the heat supply during the ATES discharging period is 81.3 GWh and thus higher than the heating demand during charging of 53.4 GWh. As a result, this ATES system acts as an overall heat supplier. The optimal ST supply is 57.7 GWh, which is slightly more than for one ATES system in scenarios A or B. Other controllable producers are only marginally affected, but there is a notable impact on HOB operation. HOB heat production is reduced by 36.6% for one LIT-ATES system or reduced by 61.8% for two ATES drillings. In comparison, there is only an 8.9% (A1) to 23.3% (B2) reduction in HOB heat production for the other scenarios.

3.5. Comparison of the Scenarios

Finally, in this part of the Section 3, the main optimization outcomes are compared regarding the DHS heat supply, emissions, and costs.

Figure 16 shows the heat generation by ST, CHP 1 (biomass CHP), other CHPs (2–3), HOB, EHB, HP, and ATES. As mentioned before, scenarios A and B generate a heat production surplus due to the equalization of storage losses. In scenario C, thermal losses are balanced within the ATES system by the ATES-HP heat production. The heat supply of the controllable producers is reduced in this case. Under the selected framework conditions, scenario A has its cost optimum with five ATES systems (A5*).

Scenarios B and C do not have an optimum with ATES (S0 would be better). The absolute CO₂ emissions and specific CO₂ emissions are shown in Figure 17. Due to the electricity consumption, which is particularly relevant for the additional HP (scenarios B and C), emissions can also be allocated to the ATES system. Renewable producers (CHP 1, EHB, HP, and ST) have a high share in the heat supply of the DHS for all scenarios, but as expected, they hardly contribute to the overall carbon emissions. Although an ATES generally leads to more ST, fewer HOBs, and a higher share of renewables, the total (and specific) carbon emissions could be higher than the reference, e.g., 117,100 t CO₂ for scenario C1 compared to 115,500 t in S0. Figure 18 shows the cost allocation and total specific heating costs of the scenarios. The contribution of the ATES to the total operating costs was already described for scenario A, but for B and C, it is even clearer.

The forced use of the seasonal ATES system means that less heat is provided by the other controllable plants. This particularly affects the HOBs, as they provide the peak load at particularly high heat generation costs. CHP, EHBs, and HPs are specifically used at the optimum times depending on the electricity market and are therefore (for many hours, if they were operating in the reference) more cost-effective than the HOB units, which are independent of the electricity market.

All details of the results are shown in Table A2 (Appendix A).

4. Discussion

This study has successfully developed an MILP optimization model that not only simulates the optimal operation of DHS supply plants but also optimizes the dimensions of ST and ATES systems from an operational cost perspective. The model provides a comprehensive evaluation of the thermal energy supply for DHS, carbon emissions, operation costs, renewable energy share, and other critical parameters such as specific emissions, costs, and heat losses. In response to the research question, the results confirm that the integration of an ATES system can indeed contribute to higher potential ST heat amounts, reduced emissions, and lower operating costs, although these benefits are not always achieved simultaneously.

However, a larger ATES size does not necessarily lead to lower operation costs or lower emissions. This is mainly because of the complex interactions between the ATES and other suppliers in the DHS, taking into account the fluctuating electricity prices and associated emissions. For example, the HOB displaced by the ATES operates more economically and has lower emissions at high electricity prices, higher electricity-related emissions, and relatively low fuel and emission costs, and therefore can represent a more sensible option from a cost and emission perspective than the ATES system with a continuously operated HP and high electricity demand. The ATES also slightly increases the operation of the CHP instead of the HOB, which may be economical but can still lead to higher emissions (if it is not biomass CHP 1). The simulations with an ATES system optimize the reduced operation of all controllable plants. Since HOB units are less flexible and often have the highest heat generation costs due to their independence from the electricity market, they are most affected by ATES and ST integration.

As noted above, ATES-HP’s high demand for electricity from the public grid is a major cause of the negative impacts. In particular, the requirement for continuous operation of the HP due to the continuous source heat load from the ATES during the discharge period means that it is not possible to optimally match the electricity demand with the situation in the electricity market. The decoupling of the ATES-HP from the electricity grid price through additional storage capacity (batteries or ATES source heat TES) could be a very useful extension of the portfolio, which could integrate the ATES more efficiently into the system. It can be concluded that the downstream heat pump must be dimensioned larger if the thermal energy of the ATES production is stored temporarily in a short-term TES. However, it is not required to run continuously, as was assumed in the present work, and can be operated at the most economical electricity prices of the day or week, resulting in lower operating costs of the ATES-HP but higher investment costs for the whole system. A battery storage system also decouples the electricity demand of the ATES-HP from highly fluctuating electricity prices. In this case, the HP can continue to operate at all times. Thanks to the battery, more electricity can be purchased in times of low prices, and in times of high prices, the storage is discharged and the electricity demand is reduced. This is then not only linked to the ATES-HP but could also influence the operation of other systems.

Moreover, the ATES-HP has a relatively low COP caused by its high temperature differences. Although the increased temperature spread in the ATES system offers the advantage of greater source heat, it is associated with a lower COP for the ATES-HP. The less auxiliary power is required, the more efficient the overall ATES system. The reduction in the DHS supply temperature provided by the ATES-HP and the acceptance of a lower ATES heat load due to higher ATES injection temperatures must be balanced for an optimal result. As renewable power capacity continues to expand, electricity prices and the carbon emissions allocated to the electricity mix will decrease, benefiting the ATES system with HPs.

The results of scenario A must also be critically evaluated. Here, only an increase in the DHS supply temperature without additional HPs was assumed, but it was not determined whether this increase could be compensated at any time by the controllable heat generators and brought to the required level.

The analysis suggests that differences in the generation portfolio could have a significant impact on the flexibility of the system to better respond to high electricity prices. A more balanced portfolio with a stronger presence of EHP, HP, and CHP could be a more effective strategy to fully exploit and complement the benefits of the ATES system without increasing the negative effects of oversizing.

Small variations, such as an ATES-HP outlet temperature lower than 105 °C, can lead to significant improvements in the COP of the HP. This is due to the smaller temperature difference between the evaporation and condensation temperatures of the HP, a factor that was evident in the comparison between scenarios B and C. The higher COP in scenario B resulted in this scenario becoming more cost-effective at a faster rate. The results indicate that external factors such as framework conditions can significantly influence the optimal size of ST systems, suggesting that while ATES is a viable solution, other impacts may also effectively increase the integration of ST into DHS, e.g., changing the CO₂ emission price distinctly changes the optimal outcome: scenario S0 without ATES reaches a higher ST supply and the optimal solution for all the ATES scenarios shifts. Figure 19 shows the heat supply for a price of 150 EUR/MWh CO₂ and scenario S0. The ST amount is 53.4 GWh, which is a 240% increase compared to only 15.7 GWh in the original scenario with a 50 EUR/MWh emission price.

With a higher carbon emission price of 150 EUR/t CO₂ instead of 50 EUR/t, scenario A would require only four ATES drillings, while the optimal solution for scenario B would be two, as opposed to none under previous conditions. Despite these variations, the overall trends remain consistent with those observed under different emission prices, as seen in Figure 20. * indicates that this is the optimal ATES dimensioning scenario in terms of operating costs.

Therefore, an important limitation is that the results are highly dependent on the assumptions made for all framework conditions, e.g., as shown above, the carbon price. This study also has several other limitations that may affect the generalizability of the results. First, investment or capital costs were not considered, which may affect the economic feasibility of the proposed solutions. Second, ATES and ST systems were not directly implemented in the MILP model; instead, precalculated time series for ST yield and ATES injection and production temperatures based on literature data were used. Finally, the model does not currently support multi-objective optimization, which would have allowed a more balanced assessment of different objectives such as cost, emissions, and renewable energy share.

Future research should focus on conducting a sensitivity analysis to better understand the influence of various assumptions and parameters on the model’s results. This could include applying the MILP model to different use cases or exploring future scenarios that include changes in electricity prices, another reference weather year, the introduction of new heat supply portfolios (such as more heat pumps), and varying fuel and emission prices, especially since these assumptions are also strongly influenced by each other. In addition, further development of the MILP model is needed to address the limitations identified in this study, such as incorporating investment costs and enabling multi-objective optimization. These improvements will provide a more robust and versatile tool for evaluating the integration of ATES and ST systems into DHS. To adapt the model to broader applications, the incorporation of diverse energy sources (e.g., industrial heat, wind, or geothermal), regional climatic frameworks and demand profiles, and economic and policy contexts must be considered. Due to the modular structure of the model, it can already be scaled relatively well to different sizes of small- and large-scale DHS, but the representation of the individual components should be improved, e.g., with regard to resource availability or partial load behavior. At this stage, the model is not easily transferable to large national scales as the spatial resolution of the heat supply is a key challenge that is not included in the developed model.

Most importantly, it should be reiterated that the studies are exclusively concerned with assessing the optimal operation of controllable DHS units under the influence of the integration of various ATES and ST systems. Since the design of ATES and ST is not the focus of the optimization, investment costs can be neglected; therefore, this work cannot be the basis for investment decisions or future strategies of the asset portfolio.

5. Conclusions

Using the developed MILP model for DHS scheduling, the effects of integrating different ATES configurations on a potential ST plant sizing could be studied, taking into account optimal operating costs. First, the MILP model was computed in Python and solved with the Gurobi 11.0.0 optimizer. In order to subsequently apply the model to the use case of the DHS of BTB Berlin, the available hourly time series (e.g., for heat load, air and river water temperature, solar irradiation, ATES temperatures, electricity prices, and associated emissions) were averaged using a hierarchical clustering method.

Four main scenarios (S, A, B, and C) were then computed. With 15.3 GWh, scenario S (ST optimization without an ATES system) has an ST share of 2.6% of the total heat production. PTH in EHB and HP account for 1.1%, CHP 60.7%, and HOB 35.6%. The specific cost of heat production is 8.34 EUR/MWh and the specific emissions are 195.42 kg CO₂ per MWh. Scenarios A and B are supplemented by an ATES at a high injection-well temperature of about 55 °C and reach 33.8 GWh of extra heat demand per drilling (one injection and one production well) during the summer charging period of the ATES. Scenario A has no additional HP and preheats the DHS return temperature with 17.7 GWh while the power demand from the deep groundwater pumping is 1.7 GWh, so the total efficiency is 49.8%. Adding an HP, which reaches a DHS supply temperature of 105 °C in scenario B, the power demand increases to 9 GWh at 25.1 discharging heat supply to the DHS (total efficiency of 58.6%). Scenario C also considers an HP and an ATES injection temperature of 20 °C and has the highest total efficiency of 86.6%. The charging heat demand increases to 53.4 GWh, the power annual power demand is 40.5 GWh, and the discharging heat supply is 81.3 GWh. Scenario A (HIT-ATES without HP) achieves an increase in ST by a factor of 8.2 (125 GWh) while reducing HOB heat production by 33% with an optimal dimensioning of 5 ATES drillings (equal to 10 wells), resulting in 5.28 EUR/MWh for the costs and 194.03 kg/MWh for the emissions. For scenario B (HIT-ATES with HP) and scenario C (LIT-ATES with HP), no ATES would be the cost-optimal solution under the given framework conditions. However, if one ATES is considered, it may not result in a cost reduction (8.78 EUR/MWh for B and 13.14 EUR/MWh for C), but rather an increase in ST: 42.5 GWh or factor 2.8 for scenario B and 57.7 GWh or factor 3.8 increase for scenario C. At the same time, the amount of HOB heat is reduced by 12.4% (B) or 36.6% (C). CO₂ emissions decline for scenario B to 194.07 kg/MWh and increase for scenario C to 198.15 kg/MWh. CHP and PTH were less affected. Important factors in reducing costs or emissions are the higher proportion of very cheap and emission-free ST heat, the displacement of the HOB heat by the ATES or ATES-HP, and the greater flexibility between CHP and PTH in the ATES discharge time with reduced DHS heat demand due to seasonal storage.

As a result, ATES integration does not automatically lead to lower operation costs or carbon emissions, but always to more renewable heat supply, e.g., from ST. Furthermore, the study shows the extreme influence of the framework conditions, which strongly change the optimal result with only small changes (such as an increase in the CO₂ emissions price). The high power demand of the HP, which raises the ATES temperature for the DHS supply, has a major impact on the economic and environmental efficiency of the system. Increasing the COP by lowering the DHS supply temperature or increasing the ATES injection temperature, as well as decoupling the ATES-HP from the public power grid (through groundwater heat TES or batteries), may be useful extensions of the research. Whether a return temperature increase as in scenario A is even feasible without an additional HP depends on the additional heat producers and is beyond the scope of this study. In addition to extending the model to include such currently limiting factors, investment costs or multi-criteria optimization could be part of future work.