Towards Sustainable Heat Supply with Decentralized Multi-Energy Systems by Integration of Subsurface Seasonal Heat Storage

: In the energy transition, multi-energy systems are crucial to reduce the temporal, spatial and functional mismatch between sustainable energy supply and demand. Technologies as power-to-heat (PtH) allow ﬂexible and effective utilisation of available surplus green electricity when integrated with seasonal heat storage options. However, insights and methods for integration of PtH and seasonal heat storage in multi-energy systems are lacking. Therefore, in this study, we developed methods for improved integration and control of a high temperature aquifer thermal energy storage (HT-ATES) system within a decentralized multi-energy system. To this end, we expanded and integrated a multi-energy system model with a numerical hydro-thermal model to dynamically simulate the functioning of several HT-ATES system designs for a case study of a neighbourhood of 2000 houses. Results show that the integration of HT-ATES with PtH allows 100% provision of the yearly heat demand, with a maximum 25% smaller heat pump than without HT-ATES. Success of the system is partly caused by the developed mode of operation whereby the heat pump lowers the threshold temperature of the HT-ATES, as this increases HT-ATES performance and decreases the overall costs of heat production. Overall, this study shows that the integration of HT-ATES in a multi-energy system is suitable to match annual heat demand and supply, and to increase local sustainable energy use. HT-ATES feeds heat pump mode provides a efficient HT-ATES system


Introduction
To limit global warming, governments aim to reducing greenhouse gas (GHG) emissions caused by the use of fossil fuels [1][2][3][4] and transition to renewable energy sources. As a result of the transition to renewable sources, the energy system will, in part, become more decentralized with energy production by e.g., photovoltaic systems (PV) and, on a small scale, wind, brought closer to consumers [2,3]. These types of renewable energies are intermittent and thus storage facilities are required to overcome the temporal mismatch between availability of and demand for energy [5][6][7]. These mismatches can (partly) be overcome by the introduction of multi-energy systems (MES), as they provide possibilities for system integration [8][9][10][11]. In such integrated systems, production of sustainable electricity and heat, as well as storage and conversion of these commodities, are integrated with the goal to efficiently maximize the utilisation of available sustainable energy and to balance supply and demand. A MES consists of sustainable energy sources and conversion and storage facilities. Different forms of conversion can be applied in a MES, such as power-to-heat and power-to-gas (i.e., hydrogen). Consequently, traditionally separate operating sectors have become connected, such as the heat, transport and electricity sectors [12][13][14][15][16][17].
In order to overcome seasonal mismatches between supply and demand, storage of gas from power-to-gas technologies is often applied. Additionally, great potential is attributed to underground thermal energy storage techniques, like aquifer thermal energy storage (ATES) systems, as these techniques offer efficient and large storage capacities at relatively low costs [18]. Yet, the authors notice that seasonal heat storage is often overlooked in the design and evaluation of sustainable multi-energy systems [12,14,19,20], despite the large potential benefits for overall system efficiency and GHG emission reduction [13,[21][22][23]. Probably, there is a lack of attention for including seasonal storage in MES because methods for the evaluation of integration based on cost-effectiveness and robustness are lacking [12,14,19,20]. In previous work [13,24], we evaluated the potential benefit of a high temperature aquifer thermal energy storage (HT-ATES) system that is used for seasonal heat storage in a multi-energy system, and compared four MES designs [25].
Here it was shown that a scenario with HT-ATES led to the most balanced energy demand profile compared to an all-electric and hydrogen scenario, the highest amount of local electricity use, and lower costs for households (250-300 €/year) compared to an all-electric scenario [25]. To further explore the novel application of HT-ATES in a MES, a more detailed description and evaluation of different innovative designs for the combined HT-ATES and heat pump system is necessary. This was not included in earlier publications. Therefore, in this research, we elaborate on the connection between the HT-ATES and the heat pump of the MES, to provide methods for integration and foster the use of HT-ATES within multi-energy systems.

Goal and Approach
The goal of this study is to develop novel methods for integration and control of a HT-ATES system in an innovative power-to-heat energy system. Moreover, we assess how these methods affect the costs and provision of sustainable heat year-round. The methods are tested for a case-study neighbourhood to assess how conditions and the different components in the integrated HT-ATES system affect the efficiency and ability to sustainably meet the heat demand.
More specifically, three aspects are studied; firstly, the effect of the heat pump design (size, condenser temperature, modes of operation between heat pump and HT-ATES) on the heat delivery by the HT-ATES system is assessed. Secondly, the performance of the HT-ATES system is assessed to evaluate the effectiveness of this storage component. Thirdly, an analysis is done on the levelized costs of the heating system (in €/GJ) as a whole.
To do this, a detailed model of a multi-energy system is expanded and integrated with a numerical hydro-thermal model to simulate the functioning of a HT-ATES system. The integration of the multi-energy system and the HT-ATES system is applied in a case study of a neighbourhood in Nieuwegein (The Netherlands) that has been the focus of earlier work on a multi-energy and water system [13,[24][25][26][27]. In this decentralized multi-energy system, renewable electricity (PV) is converted to high temperature heat and stored in the HT-ATES system during summer. In winter, the HT-ATES system is used to fulfil the heat demand of the neighbourhood. Compared to LT-ATES, the higher storage temperatures require consideration of temperature-dependent density and viscosity to simulate the heat transport [30,31]. Therefore, the simulations for heat injection, storage and extraction are performed using the coupled groundwater flow model MODFLOW and the multi-species transport code MT3DMs in connection with SEAWAT [32][33][34]. In this study, we use Flopy, a software package to operate SEAWAT from within a Python environment [35]. The multi-energy system model is also programmed in Python code to facilitate coupling with the ATES model.

Modelling Approach
Depending on the needed maximal capacity of a specific HT-ATES system (e.g., MWth), the productivity of the aquifer and the well screen length (dependent on the aquifer thickness), multiple hot and warm wells could be needed. In practice, these wells should be placed in such a configuration that interaction between wells has a positive effect. Beernink et al. (2020) showed that positive interaction between the hot and warm Compared to LT-ATES, the higher storage temperatures require consideration of temperature-dependent density and viscosity to simulate the heat transport [30,31]. Therefore, the simulations for heat injection, storage and extraction are performed using the coupled groundwater flow model MODFLOW and the multi-species transport code MT3DMs in connection with SEAWAT [32][33][34]. In this study, we use Flopy, a software package to operate SEAWAT from within a Python environment [35]. The multi-energy system model is also programmed in Python code to facilitate coupling with the ATES model.

Modelling Approach
Depending on the needed maximal capacity of a specific HT-ATES system (e.g., MW th ), the productivity of the aquifer and the well screen length (dependent on the aquifer thickness), multiple hot and warm wells could be needed. In practice, these wells should be placed in such a configuration that interaction between wells has a positive effect. Beernink et al. (2020) showed that positive interaction between the hot and warm well(s) of a HT-ATES system causes a relatively small positive effect (<7%) on overall system performance compared to a situation without interaction [36]. Therefore, to simplify the modelling exercise, we assume that the total thermal volume is stored and recovered from one hot and one warm well which do not interact. Our simulations, compared to a 3D Energies 2021, 14, 7958 4 of 31 model, thus represents a slight underestimation of actual overall performance. Because of these simplifications, an axisymmetric hydrogeological model could be used, as was done previously [30,37]. For both the hot and the warm well, a separate axisymmetric model is initialized. From monitoring data, Lopik et al. (2016) [30] established a calibrated axisymmetric model of a high temperature (80 • C) ATES system. In this study, we therefore use a similar axisymmetric model set-up and parameter values. Several modelling aspects of the model are described below: -Spatial discretization. The spatial discretization in the horizontal direction is 0.5 m close to the well location. Because flow velocity decreases with the radial distance from the well, the cell size may increase outwards. This was done logarithmically to a maximum of 50 m at the model domain boundary. To prevent boundary conditions from influencing simulation results, several test runs were carried out. The test runs showed that outcomes did not change significantly (<1%) with an increasing grid size of 1500 m around the well. This was therefore set as the outer model boundary distance. To allow for sufficient detail/insight in the vertical (free convection) flow component, the vertical layers are also discretized at high resolution; 0.5 m thickness. Further refining of the modelling grid did not result in different results. The discretization was therefore assessed to be sufficient. -Temporal discretization. Internal time steps of the SEAWAT model are limited because of the courant number that is implemented in the advection package. We use a value of 0.8, meaning that the maximal distance that advection will be allowed in one internal transport time step is 0.8 cell length; the length of the time step is automatically adjusted accordingly. Model input and output is changed and transferred with the multi-energy system model with a daily time step. This means that the flow in/out of the wells is changed once a day and evenly divided for this timespan.  [30] and are corrected for axisymmetric flow according to Langevin (2008) [37], and given in Table 1.
The groundwater flow is solved using the Preconditioned Conjugate Gradient 2 package. The standard finite-difference method with upstream or central-in-space weighting is applied in the advection package.

Implementation of Variable Viscosity and Density
The viscosity and density of groundwater decrease with increasing temperature. This has an effect on the behaviour of the groundwater at elevated temperatures compared to the ambient groundwater. Implementation of this relationship with temperature, therefore, is important for this modelling exercise. In reality, both viscosity and density have a non-linear relationship with temperature. In this study, we use the viscosity and density dependency of water as implemented by Langevin et al. (2008) [34]. Here, the nonlinear viscosity-temperature relationship of VOSS (1984) [39] was used to determine the groundwater viscosity accurately. For density, however, a linear relationship was used to calculate the density at different temperatures. This was done because SEAWAT does not simply allow for a non-linear density function. Therefore, we used an approximation based on the injection temperature to include the possible influence of buoyancy flow according to the following linear relationship, similar to previous research [40]: where T is the water temperature in the aquifer • C and ρ is the density in kg/m 3 . Here, the density change gradient per temperature difference (δρ/δT) is changed according to the maximal injection temperature, resulting in a value of −0.22 for T inj = 50 • C and a value of −0.29 for T inj = 65 • C.

Multi-Energy System Model
To analyse different layouts for future multi-energy systems (in neighbourhoods) a simulation model was developed that matches supply and demand on an hourly basis and integrates different renewable energy sources, conversion and storage mechanisms [13,25]. In Figure 2, an overview of the model components is given, with a highlight on the aspects that are the focus of this study. Most of the modelling methodology has been explained in detail earlier, including both energy and water aspects [13,24,25]; here we focus on the coupling between the multi-energy system model and HT-ATES integration. The Python simulations of different model setups are done with an hourly time step for 10 years. The most important aspects of the model related to heat production and storage are described in the following paragraphs. explained in detail earlier, including both energy and water aspects [13,24,25]; here we focus on the coupling between the multi-energy system model and HT-ATES integration. The Python simulations of different model setups are done with an hourly time step for 10 years. The most important aspects of the model related to heat production and storage are described in the following paragraphs.

Figure 2.
Overview of the components in the multi-energy system model, with a highlight on the model parts that are the focus of this paper.

COP Heat Pump
The coefficient of performance (COP) curve of the heat pump (HP) is based on supplier information (GEA: cooling agent R717, extra heat recovery by oil cooling, 1.8 MWel installation at 65 °C and 50 °C condenser temperature) for heat pump condenser temperatures of both 65 °C and 50 °C. From literature [41], we know that the COP of a heat pump (COPHP) is mainly based on the temperature lift and less on the exact condenser temperature (heat sink). Therefore, we can use one equation for different condenser output temperatures (THP,cond) (heat sinks), as long as the temperature lift between condenser and evaporator temperature (THP,evap) is taken into account: (2)

Heat Demand Pattern
The model input for heat demand is a fixed yearly heat demand that varies per type of household. Based on the degree-day method [42], the total heat demand is distributed over the runtime with the adaptation of using degree-hours, according to: with hDH_t being the number of degree-hours in a certain hour and Tair,t being the outside air temperature, based on KNMI data from De Bilt [43]. The base temperature Tbase is set to be 14 °C to calculate the degree-hours for a certain moment in time. In the degree days method, it is important to choose the right base temperature, which is often set around 18 °C [44]. We have chosen a value of 14 °C because the houses we consider in this

COP Heat Pump
The coefficient of performance (COP) curve of the heat pump (HP) is based on supplier information (GEA: cooling agent R717, extra heat recovery by oil cooling, 1.8 MW el installation at 65 • C and 50 • C condenser temperature) for heat pump condenser temperatures of both 65 • C and 50 • C. From literature [41], we know that the COP of a heat pump (COP HP ) is mainly based on the temperature lift and less on the exact condenser temperature (heat sink). Therefore, we can use one equation for different condenser output temperatures (T HP,cond ) (heat sinks), as long as the temperature lift between condenser and evaporator temperature (T HP,evap ) is taken into account: COP HP = −0.00007· (T HP,cond − T HP,evap ) 3 + 0.0097· (T HP,cond − T HP,evap ) 2 − 0.5311 T HP,cond − T HP,evap + 14.68 (2)

Heat Demand Pattern
The model input for heat demand is a fixed yearly heat demand that varies per type of household. Based on the degree-day method [42], the total heat demand is distributed over the runtime with the adaptation of using degree-hours, according to: with h DH_t being the number of degree-hours in a certain hour and T air,t being the outside air temperature, based on KNMI data from De Bilt [43]. The base temperature T base is set to be 14 • C to calculate the degree-hours for a certain moment in time. In the degree days method, it is important to choose the right base temperature, which is often set around 18 • C [44]. We have chosen a value of 14 • C because the houses we consider in this case are well isolated and probably have a low-temperature heating system (such as floor heating) and large building inertia, and therefore space heating is not necessary above 14 • C outside air temperature. The heat degree hours are also weighted based on the season. This means  [44]. The next step is to calculate the heat demand of the neighbourhood for every hour of the run period, based on the specific fraction of the total amount of heat degree hours: Here, Q spaceheat,t is the heat demand of the neighbourhood in a specific hour (t) in MJ, ∑HDH is the sum of all heat degree hours over the total run period (in this case 10 years), and E spaceheat,total is the total heat demand for the run period in MJ.
Using this method implies that the distribution of the heat demand varies with the run period chosen. On the one hand, this is an advantage as yearly variations in outside air temperature are taken into account. On the other hand, this implies that the space heat demand of a particular hour can vary slightly for a run period of i.e., 2010-2012 vs. 2017-2019. It is therefore important to only compare scenarios with the same run period; in our case data from 2010-2019 are used for all scenarios.

Heat Supply Pattern
To ensure enough heat is stored to cover the heat demand in winter, a certain amount of heat needs to be stored in spring and summer. This should be the heat demand plus the expected loss during storage. The model applies a pattern with weekly values that represent a percentage of the total yearly heat demand that should be stored in the months from March to October, as elaborated upon in an earlier publication [13]. To compensate for the heat loss, a factor is applied to the total heat demand. The model will adapt this factor every year based on the difference in volume that is stored either in the hot or warm well. If the difference in volume is more than 15%, the factor is either increased by 0.1 (if the volume stored in the hot well is too small) or decreased by 0.15 (if the volume in the warm well is too small). The starting value of the loss factor is 1.8 for this study, meaning that the heat storage target for the HT-ATES is 1.8 times the yearly heat demand of the neighbourhood at the start of the run time. This value of 1.8 is higher than the actual expected heat loss over time, because in the first years of system operation a shell with elevated temperature around the stored volume is developed. After this period, the loss factor will start to decrease and will finally stabilize around the yearly heat loss value of the wells (which is more around 1.1-1.4).

Coupling of the HT-ATES Model and the Multi-Energy System-Model
The input of the SEAWAT model (see Section 2.1) is adjusted with a daily time step, while the multi-energy system model calculates the energy flows on an hourly basis. To bridge the hourly and the daily models, the hourly data of the multi-energy system is stored and combined into daily input data for the SEAWAT model. This can be interpreted as having an above-ground buffer tank that controls the daily fluctuations of the energy system. From the hourly multi-energy system model, both the weighted average temperature and the corresponding net flow from/to the HT-ATES are used as input for the SEAWAT model. The SEAWAT model subsequently returns the temperature of the extracted groundwater of both the warm and the hot well after simulation. These values are then used for the next 24 h of the overall energy system model calculations. In order to apply our multi-energy system with HT-ATES integration, we present a field case located in the eastern part of Nieuwegein (The Netherlands). We assume 2000 houses will be connected to a low-temperature heat grid and the houses will be suitable for low-temperature floor heating. The energy for domestic hot water production is provided with a booster heat pump. Assumptions for heat demand are given in Table 2 and are based on the Dutch energy label B, which is not the most efficient building type but is on a European level in the range of a Nearly Zero Energy Building (40-70 kWh/m 2 ) [45]. Energy demand related to domestic hot water consumption is set at 3.3 GJ/person/year [46]. In this publication, we use the case study only to determine the heat demand of the neighbourhood and the role of heat storage in the fulfilment of heat demand. For exact temporal energy flows (local use, import and export) and comparison of the MES with a HT-ATES with other possible system designs, the reader is referred to an earlier publication focussing on the same neighbourhood case [25]. The main heat source of the neighbourhood is heat from surface water (aquathermal energy) provided by the Lekkanaal, a 100 m wide channel located east of the case study area. The surface water heat is extracted in summer (April-September) and the temperature is raised by a large scale heat pump. Part of the heat is used to provide heat directly, and the remainder of the heat is stored in a HT-ATES system to be utilized during winter. A HT-ATES system is chosen as the most suitable heat storage option, because of the large storage size required for 2000 houses, combined with a relatively low amount of space available on the surface level. The large scale heat pump will use excess electricity from local sources such as a solar farm and home solar PV systems on houses, supplemented with electricity from the grid to ensure enough heat is stored to fulfil the total heat demand of the neighbourhood in winter. This system layout could assist in peak reduction on the electricity grid by PV systems in summer and prevent high peaks in electricity demand in winter if houses would be heated by i.e., an (air-sourced) heat pump. In this way, it is thus possible to decouple supply and demand, which gives more flexibility in a future energy system based on renewable energy.
Hourly weather data on outside temperature from the period 2010-2019 are used in this analysis, from the weather station De Bilt (closest to the project location) [47].

Modes of Operation for the Heat Pump and HT-ATES
The heating system has different modes of operation defined in the model, depending on the situation. In the explanation, we have chosen a heat pump condenser temperature of 50 • C and heat delivery via a district heating network (DHN) with a minimum of 40 • C. Schemes for a heat pump condenser temperature of 65 • C are given in the Supplementary Information. There are five possible modes of operation:

1.
Heat production and storage in the HT-ATES (no heat demand) In this mode of operation, heat is produced with surface water (14-26 • C) as input during the summer months (see Figure 3A). Heat is injected into the hot aquifer at 48.5 • C (assumed 1.5 • C loss over the heat exchanger).

2.
Simultaneous heat production and delivery If there is a heat demand at a time when the heat pump is producing heat (during summer), the heat demand is fulfilled directly from the heat pump (50 • C) via the hot distributor (see Figure 3B). After fulfilling the heat demand, the remainder of the heat is injected into the hot aquifer at 48.5 • C.

3.
Heat delivery from the HT-ATES above the inlet temperature of the district heating network (hot well > 43 • C) In the period from October-April (approximately), the heat pump is not able to produce heat with surface water as the surface water temperature is too low (<14 • C). Heat demand from the multi-energy system is now fulfilled by extracting heat from the hot aquifer of the HT-ATES system (see Figure 3C). Water HT-ATES is shut off. (Hot well under threshold temperature < 43 • C) In this mode of operation, the temperature of the hot aquifer has dropped below the threshold temperature (43 • C) to guarantee heat delivery at 40 • C. The HT-ATES system is now shut off and cannot be used to fulfil the heat demand of the multi-energy system. In this case, heat must be delivered from another (external) source, such as a peak boiler or a sustainable heat source. This is not included in this analysis.
For this study, we added an extra mode of operation to prolong the heat delivery from the HT-ATES system and to assess its effect on the efficiency of the heat production system and the heat production costs.

5.
HT-ATES feeds heat pump mode; heat delivery from the HT-ATES with alternative threshold temperature (hot well temperature between 30 • C and 43 • C) To prolong heat delivery from the HT-ATES system, we designed an extra mode of operation for the heat pump and HT-ATES. In this mode of operation, the heat delivery from the HT-ATES continues when the temperature of the hot well is below 43 • C but above 30 • C. This is possible because of additional heat exchange with the HT-ATES in combination with a heat lift by the heat pump (see Figure 3D). First, the return flow of the DHN is exchanged with water from the hot well (30-43 • C) and thereby the return flow is heated to 28.5-41.5 • C. This already heated return flow of the DHN is then fed to the condenser side of the heat pump and raised to 50 • C. The flow from the hot aquifer is thus first exchanged with the return flow of the DHN, which results in a temperature of 26.5 • C. Then, the same flow is led to a second heat exchanger, to extract more heat by a (closed) loop that is connected to the evaporator side of the heat pump (inlet temperature of 25 • C). The return flow from the evaporator side of the heat pump is then finally exchanged with the flow that is injected in the warm well (at 16-22.5 • C). In this way, the heat from the hot well is exchanged in two stages before being injected into the warm well, making the heat exchange more efficient. Moreover, the injection temperature in the warm well is decreased, creating a larger ∆T between the wells which could be beneficial for the HT-ATES efficiency. The electricity demand of the heat pump in this mode is retrieved directly from local RES (renewable energy systems) when possible, or from the grid when no electricity from local RES is available.

Scenarios
The goal of this study is to identify methods for integration and control of a HT-ATES within a multi-energy system, as well as assess how these methods affect the performance of the energy system. In Section 2.3, we have introduced an additional mode of operation for the heat pump/HT-ATES system. We hypothesise that this mode of operation prolongs heat delivery from the hot well and could increase the overall performance of the heat pump/HT-ATES system. Therefore, we have selected ten scenarios (Table 3) to examine this hypothesis, as well as obtain general insights on how the HT-ATES system performs within the multi-energy system.
The heat demand and heat delivery temperature are kept constant for all scenarios. The amount of heat that should be provided to the neighbourhood houses is 55.2 TJ per year at a temperature of at least 40 • C. We have chosen to vary the condenser temperature, the threshold temperature (HT-ATES feeds heat pump mode) and the size of the heat pump. Table 3 gives an overview of the ten scenarios that are investigated in this study. The description code is built up as: condenser temperature|HT-ATES threshold temperature|heat pump size.
We compare two condenser temperatures because we hope to learn more about their actual effects on heat delivery, recovery efficiency and costs. We know that a higher condenser temperature of 65 • C leads to a lower COP of the heat pump compared to 50 • C, decreasing the amount of energy that can be stored with a given heat pump size. Additionally, the heat pump will be more expensive at a higher condenser temperature and more heat loss can occur at a higher storage temperature. On the other hand, storage at 65 • C leads to a smaller storage volume for the HT-ATES system, which saves costs and space. By comparing scenarios with both condenser temperatures, we obtain more insight into the combined result of these contradictory effects on the chosen performance parameters.
To ensure heat delivery at 40 • C, the threshold value for the hot well is set at 43 • C. In Section 2.3, we introduced an HT-ATES feeds heat pump mode that allows a lower threshold temperature of the hot well of 30 • C, as the heat pump is then used for an extra temperature lift. Scenarios 65|30|2, 65|30|1.5, 50|30|2 and the 50|30|1.5 will give insight into the effects of including this HT-ATES feeds heat pump mode, and can be compared to their respective twins where the HT-ATES feeds heat pump mode is not used and the threshold temperature is set at 43 • C.
Lastly, we investigate the effect of different heat pump sizes. A smaller heat pump will save costs, but it should be large enough for the HT-ATES system to provide sufficient heat in winter. We hypothesise that the scenarios with an extra heat pump mode lead to a more efficient HT-ATES system, which could imply that the amount of heat stored in the HT-ATES could be reduced. For a good comparison, we included scenarios with a 2 MW el , 1.5 MW el and 1 MW el heat pump. We expect the 1 MW el heat pump to be too small to fulfil the heat demand of the neighbourhood in the multi-energy system (55.2 TJ/y). With an average COP of 4 and approximately five months of heat production, a maximum of 51.8 TJ/y could be stored with a 1 MW el heat pump, which is less than the heat demand of the neighbourhood (55.2 TJ/y). However, to investigate the effect of a relatively small heat pump on the HT-ATES system, we included two scenarios with a 1MW el heat pump, both with a lower (30 • C) threshold temperature.

Fulfilment of Heat Demand
The main goal of the combined heat pump and HT-ATES system is to fulfil the heat demand of the neighbourhood in the multi-energy system in a reliable and affordable way. To fulfil the heat demand at a certain moment, there are multiple options. In our system, we first check if the heat demand can be fulfilled by direct heat production, which mainly happens during summer. Next, it is checked if the heat demand can be fulfilled by the HT-ATES system, which mostly happens in the winter season. If the HT-ATES feeds heat pump mode (see Section 2.3.2) is switched on, it is possible to increase the temperature of the HT-ATES with the heat pump, and the heat demand is fulfilled by the HT-ATES system with some additional electricity use. Lastly, if the temperature of the hot aquifer is below the set threshold temperature (43 • C or 30 • C in our case), the HT-ATES system will not be used for heat production. In this case, an external source of heat is necessary. This external source is not defined in this paper, but could, for example, be a peak gas boiler, or a renewable option.
In the first few years of operation, it is likely that the HT-ATES system has not yet stored enough heat to ensure continuous heat delivery throughout the year. However, the goal would be to reach continuous operation within a few years, so no external heat source is necessary and the heat pump-HT-ATES system can function independently.
Because the heat pump can work on 100% renewable electricity with surface water as a source, there are no direct GHG emissions of the system if it operates independently from an external source. Energy use during production and recycling of system components is not taken into account here. The GHG emissions of the external heat source are not taken into account in this study as we do not define a particular heat source.

HT-ATES Performance
The performance of the HT-ATES system is assessed with two main parameters, the well recovery efficiency of stored energy and the volume balance between the wells.
Regarding the recovery efficiency, we use three different indicators that give insight into the efficiency of the heat system as part of the multi-energy system. The performance of the HT-ATES system is calculated in terms of hot well recovery efficiency, warm well recovery efficiency and system recovery efficiency. For each indicator, a separate ∆T is used, see below. The average recovery efficiency (η) over the total run period (10 years) is calculated and analysed in the results section.
T hot and T warm are the temperatures of the hot and warm well (in • C) per day, v in and v out are the flows in or out of the hot well (m 3 /h), and c w is the heat capacity of water (4180 kJ/K/m 3 ). T ambient is the ambient background temperature of the aquifer (12 • C). Depending on the ∆T applied in Formula (5), either the system recovery efficiency (6), hot well recovery efficiency (7), or warm well recovery efficiency (8) is calculated.
Secondly, the volume balance ratio (r VB ) between the hot and warm wells is calculated as described by Beernink et al. (2019) [48]: where V heat,storage is the yearly total storage volume (m 3 ) and V heat,production is the yearly total production volume (m 3 ). Here, a positive ratio means that more volume is used to store heat in the hot well. Therefore, this also means that less volume is stored in the warm wells. Oppositely, a negative ratio means that more volume is stored in the warm wells.

MES Performance
For improved integration of HT-ATES within a multi-energy system, it is relevant to know the efficiency of the MES as a whole, under varying modes of HT-ATES integration. This indicator gives insight into the final efficiency by which the heat is delivered from the MES to the houses. Thus, this parameter is the amount of total delivered energy (heat) to the neighbourhood divided by the used total energy (electricity) for heat production: With η MES as the MES efficiency, En del,10 y avg is the average amount of delivered heat to the neighbourhood (in TJ/year over the first ten years), which can thus be lower than the actual heat demand. COP 10 y avg is the average COP of the heat pump over the run time (10 years), and El dem,10 y avg is the average electricity demand of the total heat system over the run time (10 years).
It includes the efficiency of both directly produced as well as stored heat and heat loss in the district heating network. In case a certain system design is not able to provide the total heat demand, this is corrected for by only taking the delivered heat into account and not the total heat demand of the neighbourhood. The part of the heat that must be delivered from an external source is thus not taken into account by this parameter. The MES efficiency purely gives the efficiency of heat delivery of the MES including a heat pump and HT-ATES system, without other sources.

Economic Analysis/Levelized Cost of Energy
We determine the levelized costs of heat production (LCOE) in €/GJ, including the costs for heat storage, and excluding the heat delivery (district heating network), according to: in which: where α is the capital recovery factor (no unit), representing a fraction of the total CAPEX costs depending on the yearly depreciation. r is the discount rate (as a fraction of 1) and L i is the lifetime of a system component i. The CAPEX i are the total capital expenditures for a particular system component i, such as the heat pump, heat exchanger or HT-ATES wells. OPEX i are the operational expenditures and maintenance in €/year for a system component i, expressed as a percentage of the CAPEX. E cost,i are the average (in this publication, for 10 years) electricity costs for the heat system (in €/year). These costs are thus only for the modelled heat system, and not for any external source that might be necessary to fulfil the total heat demand. Q heat,delivered is the part of the heat demand delivered by the multi-energy system with HT-ATES in GJ/year. We thus only take into account the heat that can actually be delivered by the multi-energy system with HT-ATES and not by an external source. This ensures a fair LCOE comparison between different scenarios. Cost data is given in Table 4. A discount rate (r) of 6% (0.06) is applied, and electricity costs are set at 60 €/MWh e based on CBS data for non-households with an electricity demand of >70,000 MWh e in 2019 [49], and is similar to Wesselink et al. 2018 [22]. Lastly, the electricity use of the heat pump is calculated, averaged over the run period (here 10 years).

Results
In the previous section, we described the models, the coupling of the models and the scenarios that are used as input for the coupled model infrastructure. In this section, we describe the results of the scenario model runs assessed on the heat fulfilment by the (HT-ATES) system to the neighbourhood, the performance of the HT-ATES, the differences in costs (LCOE) between scenarios, and finally we combine these results in an overall analysis of the performance.

Fulfillment of Heat Demand with HT-ATES
Simulation results show that in all cases, 12-14% of the heat demand in summer is delivered to the neighbourhood directly by the heat pump for all scenarios (Figure 4) during the ten-year simulation (A) and the last five years of the simulation (B). This percentage is the same for each scenario because the net demand is the same and in all cases, the heat pump capacity is much larger than the heat demand in summer, when the heat pump produces heat that otherwise would be stored in the HT-ATES.
The simulation results show that at the heat delivered from the HT-ATES, a distinction is made between scenarios that are able to reach heat delivery by the HT-ATES system within five years, and scenarios that take longer to reach 100% delivery of demand by the heat pump and HT-ATES system. In scenarios that take longer than five years to reach 100% delivery, the HT-ATES system is shut off for >20% of the time, which negatively influences the amount of heat delivered from the HT-ATES system (40-65%). The other scenarios that reach a stable heat supply within five years have 70-80% of the heat demand delivered from the HT-ATES system. In Figure 4B, the origin of the heat delivered during the last five years of operation of the heat pump-HT-ATES system are shown. Five scenarios can operate independently of an external source after an initial start-up period.
For this case study, we show that the addition of a HT-ATES to a multi-energy system with power-to-heat results in >90% fulfilment of the heat demand during the first 10 years of operation, while after 3-5 years, the combined system can supply the full heat demand of the neighbourhood in the multi-energy system sustainably.
The HT-ATES feeds heat pump mode of operation allows a lower threshold temperature. As a result, the HT-ATES system can be operated longer during winter. Hence, with a 2 MW el heat pump, the amount of heat provided from the HT-ATES is increased by 6-7% when using the HT-ATES feeds heat pump mode (see Figure 4A). When specifically comparing the scenarios at 50 • C condenser temperature (see 50|43|2 and 50|30|2 in Figure 4A), the HT-ATES feeds heat pump mode decreases the time needed to reach 100% heat delivery from 5 to 3 years. This positive effect of the HT-ATES feeds heat pump mode is even more pronounced with a 1.5 MW el heat pump and leads to 9-11% more heat delivery from the HT-ATES system.
For the 65 • C and the 1MW el and 1.5 MW el heat pump scenarios, it is not possible to reach 100% heat delivery in 10 years. Without the HT-ATES feeds heat pump mode, an external heat source is still necessary after 10 years of operation (see 65|43|1.5 and 65|30|1.5 in Figure 4A).
At 50 • C, the HT-ATES feeds heat pump mode is able to reduce the time it takes to reach stable operation to five years, compared to nine years for a scenario without the HT-ATES feeds heat pump mode (see 50|43|1.5 and 50|30|1.5 in Figure 4A). Thus, the HT-ATES feeds heat pump mode makes it possible to create a reliable HT-ATES system with a smaller heat pump (1.5 MW el vs. 2 MW el ) at 50 • C condenser temperature.
The amount of electricity necessary for this mode comes from local RES when possible, or from the grid otherwise. In Figure 5, the electricity demand for heat production is shown over time in relation to the heat demand. This graph makes it clear that the electricity demand for the extra mode of operation is small in relation to the total electricity demand of the HT-ATES system. In all scenarios, it is never more than 2-6% of the total electricity demand of the HT-ATES system (6% for the 65|30|1.5 scenario in Figure 5). This small amount of electricity used in winter decreases the utilization of an external heat source for fulfilment of heat demand by 6-12%, depending on the scenario (12% for the 65|30|1.5 scenario). Furthermore, Figure 5 clearly shows the temporal decoupling of heat production and demand made possible by the integration of HT-ATES in the MES.
Energies 2021, 14, x FOR PEER REVIEW Figure 4. Stacked bar plots of the distribution of heat delivery from direct production of h from the HT-ATES system. In grey, the amount of heat that cannot be fulfilled by the hea and HT-ATES system is shown. The left plot (A) shows the distribution for the first ten operation. The numbers on top of the bars indicate the number of years when full heat delive the heat pump-HT-ATES system is reached, with no external source requirement. The right shows the distribution of heat delivery for the last five years of operation, with five scenar able to fulfil the heat demand. . Stacked bar plots of the distribution of heat delivery from direct production of heat and from the HT-ATES system. In grey, the amount of heat that cannot be fulfilled by the heat pump and HT-ATES system is shown. The left plot (A) shows the distribution for the first ten years of operation. The numbers on top of the bars indicate the number of years when full heat delivery from the heat pump-HT-ATES system is reached, with no external source requirement. The right plot (B) shows the distribution of heat delivery for the last five years of operation, with five scenarios fully able to fulfil the heat demand.

Well Temperature Development
Results show that the fulfilment of the heating demand by the HT-ATES system varied significantly between the 10 scenarios (Figure 4), to a large extent as a result of the varying performance of the HT-ATES system. The HT-ATES system cannot supply the entire demand if the threshold temperature is reached during heat extraction from the hot well. Figure 6 shows the well temperature over time of the relatively large 2 MW el heat pump size scenarios. The hot well temperature drops below the threshold temperature during recovery in the first years. After year 4, the hot well temperature always stays above the threshold temperature, both for the T threshold = 43 • C and the HT-ATES feeds heat pump mode (T threshold = 30 • C). The effect of the HT-ATES feeds heat pump mode is visible as a sharp drop in warm well temperature when the hot well temperature drops below the first threshold of T = 43 • C. For the 65 • C storage temperature, the HT-ATES feeds heat pump mode is no longer activated after the system has sufficiently warmed up after the first years, as shown in Figure 6A. For the 50 • C storage scenarios ( Figure 6B), the HT-ATES feeds heat pump mode is active for most of the modelled years, because the hot well temperature often drops below the threshold of 43 • C at the end of extraction season. This is also visible in the warm well temperature of this scenario in Figure 6B; the warm well temperature drops once the HT-ATES feeds heat pump mode becomes active, which results in a very irregular warm well temperature pattern. Figure 4. Stacked bar plots of the distribution of heat delivery from direct production of he from the HT-ATES system. In grey, the amount of heat that cannot be fulfilled by the heat and HT-ATES system is shown. The left plot (A) shows the distribution for the first ten y operation. The numbers on top of the bars indicate the number of years when full heat deliver the heat pump-HT-ATES system is reached, with no external source requirement. The right p shows the distribution of heat delivery for the last five years of operation, with five scenario able to fulfil the heat demand. Figure 5. Electricity demand (left axis) for heat production is shown together with the heat d (right axis) for the 65|30|1.5 scenario, which is the scenario with the highest utilization of t ATES feeds heat pump mode. The electricity demand is divided into heat production for t Figure 5. Electricity demand (left axis) for heat production is shown together with the heat demand (right axis) for the 65|30|1.5 scenario, which is the scenario with the highest utilization of the HT-ATES feeds heat pump mode. The electricity demand is divided into heat production for the HT-ATES (Elec-HT-ATES), the direct production of heat for the MES (Elec-direct heat) and the electricity for the HT-ATES feeds HP mode (Elec-HP mode).
For the smaller heat pump scenarios (1.5 and 1 MW el ), the temperature in the hot and warm wells is presented in Figure 7. Here we observe the same behaviour, but the temperature drops in both the hot and the warm wells are stronger. Hence, for these smaller heat pump capacity scenarios, less heat is stored compared to a 2 MW el heat pump, resulting in a more depleted heat storage at the end of the extraction season. As a result, the threshold temperature of 43 • C is reached for both the 65 • C scenarios and the 50 • C scenarios in almost all years. However, at the end of the simulated period, we see that higher temperatures are maintained and the threshold temperature is reached less often. This indicates that the losses of previous years have created a shell with elevated temperature around the stored volume over time. As a consequence, the heat loss decreases and the system moves towards an equilibrium situation.

HT-ATES Performance Results
The hot well, warm well and system recovery efficiency are presented in Table 5.
In the previous section, we showed that the heat that is delivered by the multi-energy system is influenced strongly by the size of the heat pump. With a large heat pump size, more energy is stored in the subsurface in summer, which results in more energy that is recovered in winter. Table 5 shows that, while this might feel counter-intuitive, the system recovery efficiency is lowest (50-55%) for the scenario where most heat is provided to the multi-energy system (largest heat pump scenarios) and vice versa. Because the amount of stored energy varies between the scenarios, the delivery of heat from the HT-ATES system cannot be directly coupled to the recovery efficiency of the HT-ATES system.

HT-ATES Performance Results
The hot well, warm well and system recovery efficiency are presented in Table 5.
In the previous section, we showed that the heat that is delivered by the multi-energy system is influenced strongly by the size of the heat pump. With a large heat pump size more energy is stored in the subsurface in summer, which results in more energy that is recovered in winter. Table 5 shows that, while this might feel counter-intuitive, the system recovery efficiency is lowest (50-55%) for the scenario where most heat is provided to the multi-energy system (largest heat pump scenarios) and vice versa. Because the amount o stored energy varies between the scenarios, the delivery of heat from the HT-ATES system cannot be directly coupled to the recovery efficiency of the HT-ATES system.
In the large heat pump scenarios, more volume is stored in summer than required (and thus recovered) in winter, resulting in low recovery efficiencies. Oppositely, for the small heat pump scenarios, more volume is extracted in winter compared to what was stored in summer. The balance between these two volume flows, the volume balance ratio In the large heat pump scenarios, more volume is stored in summer than required (and thus recovered) in winter, resulting in low recovery efficiencies. Oppositely, for the small heat pump scenarios, more volume is extracted in winter compared to what was stored in summer. The balance between these two volume flows, the volume balance ratio (r VB ) varies from −0.15 for small heat pump scenarios to 0.21 for the large heat pump scenarios (Table 5).
In Figure 8, the recovery efficiency and r VB during the last five years of simulation are shown. The r VB has an opposite effect on the hot and the warm well efficiency. For the hot well, a negative r VB means that relatively large amounts of water is extracted, compared to what was initially injected, resulting in high recovery efficiency. The opposite is true for positive r VB . In Figure 8, the recovery efficiency and rVB during the last five years of simulation are shown. The rVB has an opposite effect on the hot and the warm well efficiency. For the hot well, a negative rVB means that relatively large amounts of water is extracted, compared to what was initially injected, resulting in high recovery efficiency. The opposite is true for positive rVB. Figure 8. HT-ATES System, hot well and warm well recovery efficiency (last five years) (y-axis) vs. the volume balance ratio (rVB) during the last five years (x-axis) for the ten modelled scenarios. The average recovery efficiency during the last five years of simulation is shown to take out any effects of the first start-up years of HT-ATES operation.
The system recovery efficiency is a result of the recovery efficiency of the hot and the warm well. However, the hot well has a stronger influence, compared to the warm well, because relatively more heat is stored in this side of the system. Consequently, the highest system recovery efficiency is observed for the highest hot well recovery efficiencies at negative rVB. However, for the lowest rVB ratios, the system recovery efficiency stops increasing because the warm well recovery decrease has a stronger effect compared to the hot well recovery efficiency increase. It thus seems that maximal system recovery is obtained at slightly negative rVB.
The recovery efficiency is not only influenced by the volume balance ratio; the stor- Figure 8. HT-ATES System, hot well and warm well recovery efficiency (last five years) (y-axis) vs. the volume balance ratio (r VB ) during the last five years (x-axis) for the ten modelled scenarios. The average recovery efficiency during the last five years of simulation is shown to take out any effects of the first start-up years of HT-ATES operation.
The system recovery efficiency is a result of the recovery efficiency of the hot and the warm well. However, the hot well has a stronger influence, compared to the warm well, because relatively more heat is stored in this side of the system. Consequently, the highest system recovery efficiency is observed for the highest hot well recovery efficiencies at negative r VB . However, for the lowest r VB ratios, the system recovery efficiency stops increasing because the warm well recovery decrease has a stronger effect compared to the hot well recovery efficiency increase. It thus seems that maximal system recovery is obtained at slightly negative r VB .
The recovery efficiency is not only influenced by the volume balance ratio; the storage temperature affects the recovery efficiency as well. For higher storage temperature, more losses due to buoyancy flow occur, resulting in lower recovery efficiency. This is clearly visible for the scenario operating in volume balance (last five years) in Figure 8. Here, the warm well recovery efficiency is higher compared to the hot well recovery efficiency. Moreover, we observe two linear relationships for the hot well efficiencies, corresponding to the two varied storage temperatures, which is due to higher energy losses during storage for the 65 • C storage scenarios. Due to the limited difference in storage temperatures (65 • C vs. 50 • C) and the strong effect of volume balance, the recovery efficiency is dominated by the volume balance ratio for the simulated scenarios.
The system recovery efficiency of the HT-ATES only indirectly provides insight into the amount of heat that is delivered to the multi-energy system, because this depends on the absolute amount of energy that is actually stored in the HT-ATES system. The largest heat pump scenarios can store and deliver most heat, but they do this with a relatively small system recovery efficiency, shown in Figure 9. Consequently, this means that a relatively high amount of the stored heat is lost to the subsurface. Additionally, we observe here that for equal heat pump size, the 50 • C scenarios can provide more heat, with a lower system recovery efficiency. This is caused by the fact that more heat is stored for the 50 • C scenario (higher COP of heat pump) and that less heat is lost during storage due to the lower storage temperature and larger storage volume.
Energies 2021, 14, x FOR PEER REVIEW 22 of 33 Figure 9. The absolute amounts of delivered heat (y-axis) vs. the HT-ATES system efficiency (x-axis) for the ten modelled scenarios. An optimal system would have both a high amount of delivered heat (TJ) and a high HT-ATES system efficiency. The HT-ATES feeds heat pump mode provides a more efficient HT-ATES system for all scenarios.

LCOE of Heat Production and Storage
The levelized costs of energy (in this case heat) in €/GJ are shown in Table 6. All scenarios with a 50 °C condenser temperature have a higher LCOE (2.5-4.5 €/GJ or 11-18%) than at 65 °C condenser temperature with the same heat pump size. An important reason Figure 9. The absolute amounts of delivered heat (y-axis) vs. the HT-ATES system efficiency (x-axis) for the ten modelled scenarios. An optimal system would have both a high amount of delivered heat (TJ) and a high HT-ATES system efficiency. The HT-ATES feeds heat pump mode provides a more efficient HT-ATES system for all scenarios. To get more direct insight into the actual efficiency of the delivered heat within the MES, we have included the system delivery efficiency in the calculations (see Table 5). The system delivery efficiency of the MES is similar (1-2% lower) to the system recovery efficiency of the HT-ATES for most scenarios. This outcome indicates that heat loss in the district heating network and the direct heat production and delivery more or less outweigh each other in these scenarios. Three scenarios that take longer than eight years to reach 100% delivery and have a negative r VB (65|30|1.5, 65|30|1 and 50|30|1) have a slightly lower MES efficiency (3-6%) than HT-ATES system recovery efficiency. Still, the same trend is observed that scenarios with a small heat pump have a high MES efficiency, but a relatively low amount of heat delivered, and vice versa.
The influence of the HT-ATES feeds heat pump mode when recovery temperatures reach the threshold is observed in Figure 9. For the four scenarios for which both the normal and the extra mode of operation were simulated, we observe that both the absolute amount of energy that can be delivered to the multi-energy system and the system recovery efficiency strongly increase. The HT-ATES feeds heat pump mode thus has a positive impact on the performance of the HT-ATES, resulting in a higher recovery efficiency of the system (4-10% depending on the scenario) which increases the amount of heat delivered (in TJ) to the houses (7-15% depending on the scenario).
In Figure 9, the optimal scenario regarding heat delivery is the 50|30|2 scenario, because most heat is delivered (difference of 25 TJ), while the same system recovery efficiency is observed as the 50|30|1.5 and 65|30|2 scenarios.

LCOE of Heat Production and Storage
The levelized costs of energy (in this case heat) in €/GJ are shown in Table 6. All scenarios with a 50 • C condenser temperature have a higher LCOE (2.5-4.5 €/GJ or 11-18%) than at 65 • C condenser temperature with the same heat pump size. An important reason is the higher COP for the 50 • C scenarios, which results in more heat injection in the HT-ATES with the same heat pump size. The difference in COP is proportional to the difference in Carnot efficiency based on the extra temperature lift when the heat pump condenser temperature is raised to 65 • C instead of 50 • C. This means that it is 30% more efficient to produce heat from surface water at 50 • C (average COP = 5.5) than at 65 • C (average COP = 4). Moreover, the heat pump investment costs per kW th are higher for 65 • C (600 €/kW th ) than for 50 • C condenser temperature (400 €/kW th ). On the other hand, this cost difference is almost cancelled out as the COP is higher at 50 • C, resulting in a higher thermal heat pump capacity for the same electrical heat pump capacity. Eventually, the investment costs for the heat pump are slightly lower (around 6%) with 50 • C condenser temperature. The investments for the HT-ATES depend on the required capacity and with that the number of doublets (shown in Table 6). There is an increase of one doublet required at 50 • C vs. 65 • C condenser temperature, as the capacity delivered from the HT-ATES is quite similar over all scenarios, resulting in limited variations in the HT-ATES investment costs across the scenarios. The main decrease in costs between 50 • C and 65 • C condenser temperature is caused by the lower electricity demand (because of higher COP) at 50 • C vs. 65 • C.
Next to the condenser temperature, the HT-ATES feeds heat pump mode has an impact on the LCOE. The differences in results between 65|43|2 and 65|30|2, as well as the 50|43|2 and 50|30|2 scenarios illustrate this, as the only difference here is the extra mode of operation. In both cases with the additional mode of operation, there is an increase in costs for an additional heat exchanger, but a decrease in electricity costs while the amount of delivered heat increases. For the 65 • C condenser temperature, the electricity costs increase by 2%, while 6% more heat is delivered from the heat pump/HT-ATES system, which leads to a reduction of 1.0 €/GJ for the LCOE; while at 50 • C condenser temperature, the HT-ATES feeds heat pump mode leads to a decrease in electricity costs of 5%, and 7% more heat delivered, resulting in a total LCOE decrease of 1.4 €/GJ. Table 6. Energetic and financial performance of the heat production and storage system as part of a multi-energy system. The % of heat delivered by the HT-ATES/heat pump is the sum of the yellow and orange bars in Figure 4. The scenario name abbreviation consists of: condenser temperature|HT-ATES threshold temperature|heat pump size. The HT-ATES feeds heat pump mode thus results in extra heat delivered at limited extra costs, hence the low LCOE for these scenarios. With a smaller heat pump size of 1.5 MW el , the LCOE decreases further. For the 65 • C scenarios, the electricity demand does significantly decrease as well, although the % of heat demand delivered decreases with a smaller heat pump. Still, a reduction of 1.1 €/GJ is realized compared to a 2 MW el heat pump with an LCOE of 18.5 €/GJ. The cost reduction is further enhanced by the HT-ATES feeds heat pump mode with 12% more heat delivery resulting in an LCOE of 16.6 €/GJ. For a 50 • C condenser temperature, the same effects are visible, although in this case, the HT-ATES feeds heat pump mode does not lead to an increase in electricity demand, but raises the heat demand delivered to 90%, which is almost similar to the 50|30|2 scenario, but with a smaller heat pump (1.5 MW el vs. 2 MW el ). In this scenario (50|30|1.5), the LCOE is 13.8 €/GJ, which is 1.7 €/GJ lower than without the HT-ATES feeds heat pump mode.

Scenario
A 1 MW el heat pump with a 65 • C condenser temperature (65|30|1) does not lead to a further decrease in LCOE compared to a similar scenario with a 1.5 MW el heat pump (65|30|1.5). The savings on heat pump investment and electricity costs thus do not compensate for the decrease in heat delivery (56% vs. 80%) by the heat pump/HT-ATES system.
Overall, the 50|30|1 scenario with a 1 MW el heat pump has the lowest LCOE, with 75% of the heat demand delivered from the heat pump/HT-ATES system. The lower investment costs for the heat pump do make this scenario the most feasible. Furthermore, these results show that by installing a HT-ATES system, the heat pump size is not the determining factor for heat delivery anymore. As the HT-ATES is the main source of heat delivery during winter, the capacity of the HT-ATES system is more important than the size of the heat pump. HT-ATES systems thus help to reduce the size of the heat pump.
These results indicate a tipping point in heat pump size in relation to the amount of heat that can be delivered from the HT-ATES system to the multi-energy system. In the cases where the heat pump size is becoming a limiting factor, the positive effect of the extra mode of operation is most pronounced. The actual heat pump size tipping point will depend on the local conditions (maximum required heating power, temperature level) and thus vary per project. This analysis showed that for identifying the lowest LCOE scenario, it is worthwhile to identify this tipping point when a HT-ATES system in a multi-energy system is realized.

Overall Performance: An Integral View on Heat Delivery, System Efficiency and Costs
In this study, we considered several scenarios with varying heat pump designs and storage temperature levels and assessed how these variations impact the fulfilment of the neighbourhood heat demand, the performance of the HT-ATES system, and the costs of heat production including storage (LCOE). The overall performance of the varied scenarios is assessed here by combining all assessment criteria, as shown in Figure 10. Here, the arrows indicate the influence of the HT-ATES feeds heat pump mode compared to the normal operation mode. The most optimal system would have: 1. a large amount of delivered heat (y-axis); 2. a high HT-ATES system efficiency (x-axis); and 3. low LCOE (colour). The following observations stand out: -At a small heat pump size and a high storage temperature (65|30|1), a relatively small amount of the total heat demand can be delivered to the multi-energy system (0.56) and LCOE does not decrease compared to a 1.5 MW el heat pump scenario (both around 16.6 €/GJ). However, decreasing the storage temperature to 50 • C has a large positive effect for this small heat pump scenario (50|30|1); the LCOE strongly decreases (16.7 to 12.2 €/GJ) and the delivered amount of heat increases to 75%. -Overall, the 50 • C storage temperature ensures lower costs and higher fulfilment of heat demand to the neighbourhood. This is probably caused because the heat pump is cheaper and more effective, which results in larger amounts of stored heat. This has a stronger effect than the increase in costs due to the needed number of wells. -For all scenarios, the HT-ATES system feeds heat pump mode (43 • C vs. 30 • C shutoff temperature) has a strong positive effect on all assessment criteria: the LCOE decreases, the fulfilment of heat demand increases, and the MES efficiency increases. Hence, this mode of operation has a high potential to be integrated into future HT-ATES systems. With only a small amount of electricity use in winter (not more than 2-6% of the total electricity use), the amount of heat delivered from the HT-ATES increases by 6-12%. - The most optimal scenario for this neighbourhood, based on all assessment criteria, is the 50|30|1.5 MW scenario. This scenario can provide 100% of the heat demand yearround to the neighbourhood within the first five years, while the LCOE is 13.8 €/GJ, the second cheapest system.

The Functioning of a HT-ATES within a Multi-Energy System
We have shown that seasonal heat storage, HT-ATES, is crucial to be able to increase the amount of locally provided sustainable energy in a multi-energy system. Hence, HT-ATES is probably a better option for seasonal heat storage instead of the hot water storage tanks that are usually applied in a multi-energy system [12,14,15,19]. However, HT-ATES is not applicable everywhere, as it requires the presence of an aquifer.
The use of HT-ATES including the heat pump booster mode enables better utilisation of the stored heat, and with that, increases HT-ATES performance. We showed that this mode of operation has a high potential to be integrated into future HT-ATES systems, because the costs strongly decrease and HT-ATES performance increases, leading to an increase in local sustainable heat delivery.
The performance of the HT-ATES system is assessed via the recovery efficiency of Figure 10. The relationships between the three main aspects investigated in this study; the fulfilment of heat demand by HT-ATES + heat pump to the MES (10-year average), the 10y MES efficiency and the costs of heat production and storage (LCOE). The * added to the scenario code indicates that the HT-ATES can deliver 100% of the heat to the multi-energy system within the first five years.

The Functioning of a HT-ATES within a Multi-Energy System
We have shown that seasonal heat storage, HT-ATES, is crucial to be able to increase the amount of locally provided sustainable energy in a multi-energy system. Hence, HT-ATES is probably a better option for seasonal heat storage instead of the hot water storage tanks that are usually applied in a multi-energy system [12,14,15,19]. However, HT-ATES is not applicable everywhere, as it requires the presence of an aquifer.
The use of HT-ATES including the heat pump booster mode enables better utilisation of the stored heat, and with that, increases HT-ATES performance. We showed that this mode of operation has a high potential to be integrated into future HT-ATES systems, because the costs strongly decrease and HT-ATES performance increases, leading to an increase in local sustainable heat delivery.
The performance of the HT-ATES system is assessed via the recovery efficiency of heat from the HT-ATES system. The HT-ATES system recovery is determined by the recovery efficiency of heat from the hot and the warm well, which are strongly related to the volume (im)balance between the two wells. When the volume imbalance is negative (net extraction from the hot well), the recovery efficiency of the hot well increases and decreases at the warm well. However, the total amount of energy that can be provided to the MES can be larger than expected based on a simulation with a volume balance. For example, this could have a significant improvement for the recovery efficiency of HT-ATES wells that experience high heat losses during storage. Previous studies [31,52,53] evaluated the effect of several HT-ATES storage conditions on the recovery efficiency at V in = V out , and showed that energy losses can vary significantly. Therefore, it is recommended to analyse the combined effect of both varying storage conditions and the effect of volume (im)balance on recovery efficiency to get insight in the optimal combination of both factors for a given scenario.
In this study, we separately simulated the hot and warm wells of the HT-ATES systems. In reality, the wells will be placed next to each other, in such a way that positive interaction from the hot and warm wells is optimal. Previous studies showed that this leads to increased HT-ATES system performance [36].

The Impact of HT-ATES in a Multi-Energy System
In a multi-energy system with power-to-heat but without a HT-ATES, the heat pump size would be determined by the maximum peak in heat demand. By including a HT-ATES system in the multi-energy system, the heat demand and supply are mostly decoupled and the (peak) heat demand is delivered from the HT-ATES system during winter. The heat pump is employed in an optimal way by producing most of the heat in summer, the most likely time of abundance of (local) solar energy, but is also, to a small extent, employed in winter for the HT-ATES feeds heat pump mode when necessary. Figure 11 shows the maximum capacity delivered to the heat users and the condenser capacity of the applied heat pump. The scenarios which end up in 100% heat delivery by the HT-ATES in the multi-energy system (*) indicate that the peak heating demand is 11.8 MW th . For these scenarios, the ratio of the maximum demand and the applied condenser capacity shows how much smaller heat pump capacity can be installed to optimally utilise heat produced by power surplus. The 50|43|2 and 50|30|2 scenarios have a positive ratio, indicating the heat pump is over-dimensioned. From the 65|43|2, 65|30|2, 50|43|1.5 and the 50|30|1.5 scenario results, it is concluded that the heat pump can have a 21-25% smaller capacity than the actual peak heating demand. The other scenarios (not leading to 100% heat delivery) have ratios of 0.6-0.7, indicating that 25% lower heat pump capacity is the maximum. Besides a smaller capacity, the heat pump is employed both in summer to store the heat and in winter to provide additional heat when the HT-ATES heat production is not sufficient. This fosters year-round heat pump utilisation, although the energy use in winter is a fraction (<6%) of the total energy use of the HT-ATES system (see Figure 5). Thus, most energy demand is still during summer, with a likely oversupply of sustainable local electricity from PV.
Energies 2021, 14, x FOR PEER REVIEW 27 of 33 Figure 11. The maximum heat pump capacity vs. the heat delivery at the peak demand is shown for all scenarios. Because of the addition of the HT-ATES system, it is possible to deliver the peak demand from the HT-ATES, which explains why the maximum thermal heat pump capacity is lower than the maximum capacity delivered. The * sign shows which scenarios reach 100% heat delivery by the HT-ATES to the MES (within the first five years of operation).
Following the basic heat pump dynamics, a simplified calculation for a multi-energy system with only direct heat production will result in a total average electricity use of 12.8 TJ/year with a 3.1 MWel heat pump (assuming an average heat pump evaporator temperature of 8 °C, and 50 °C condenser temperature), with an electricity demand as shown in Figure 12. In the first 10 years of operation, a MES with HT-ATES will consume around 16-17 TJ/year of electricity and a maximum of 90% of the heat demand will be fulfilled. Yet, after the first two to four years of operation, the heat demand can be provided completely from the HT-ATES and the electricity demand of the system decreases. For example, in the case of the most optimal 50|30|1.5 scenario, the electricity demand in the 10th year of operation is 13.6 TJ. This is about 6% higher than a MES with direct heat production. The total electricity use thus slightly increases when HT-ATES is included. The total electricity use by the HT-ATES + HP system is influenced in two ways. Firstly, the source of heat (canal) is mainly harvested when its temperature is relatively high. This is favourable for the heat pump COP and results in relatively less electricity use. Secondly, more heat needs to be produced than the total yearly heat demand because a part of the produced heat will be lost during storage. Hence, the second effect is of slightly more influence on the total electricity demand in this study, for the first 10 years of operation. However, as pointed out in the previous section, the energy losses of HT-ATES systems, and thus their electricity use, may vary strongly for different HT-ATES storage conditions (e.g., storage volume, hydrogeological conditions). Consequently, under more favourable storage conditions, yearly electricity use of the MES may be lower compared to direct heat production.
More importantly, the integration of HT-ATES has a strong effect on the distribution of electricity use over the year. The electricity use is decoupled by the integration of HT-ATES, and the electricity use peak switches from winter to summer, following the availability of sustainable electricity by solar energy (Figure 12). Here, we compare the Figure 11. The maximum heat pump capacity vs. the heat delivery at the peak demand is shown for all scenarios. Because of the addition of the HT-ATES system, it is possible to deliver the peak demand from the HT-ATES, which explains why the maximum thermal heat pump capacity is lower than the maximum capacity delivered. The * sign shows which scenarios reach 100% heat delivery by the HT-ATES to the MES (within the first five years of operation).
The addition of HT-ATES in a power-to-heat MES affects the total yearly electricity use, and, the yearly distribution of electricity use, as energy demand and availability are decoupled. For a multi-energy system with power-to-heat without HT-ATES, all heat needs to be produced directly by the heat pump on demand, which can become a problem in winter because (a) the temperature of the heat source (canal) becomes too low to use, and (b) sustainable power to run the heat pump is not available. This results in the need for larger heat exchangers and the need for other sources of heat and associated costs, as well as a reduction in the use of local sustainable energy.
Following the basic heat pump dynamics, a simplified calculation for a multi-energy system with only direct heat production will result in a total average electricity use of 12.8 TJ/year with a 3.1 MW el heat pump (assuming an average heat pump evaporator temperature of 8 • C, and 50 • C condenser temperature), with an electricity demand as shown in Figure 12. In the first 10 years of operation, a MES with HT-ATES will consume around 16-17 TJ/year of electricity and a maximum of 90% of the heat demand will be fulfilled. Yet, after the first two to four years of operation, the heat demand can be provided completely from the HT-ATES and the electricity demand of the system decreases. For example, in the case of the most optimal 50|30|1.5 scenario, the electricity demand in the 10th year of operation is 13.6 TJ. This is about 6% higher than a MES with direct heat production. The total electricity use thus slightly increases when HT-ATES is included. The total electricity use by the HT-ATES + HP system is influenced in two ways. Firstly, the source of heat (canal) is mainly harvested when its temperature is relatively high. This is favourable for the heat pump COP and results in relatively less electricity use. Secondly, more heat needs to be produced than the total yearly heat demand because a part of the produced heat will be lost during storage. Hence, the second effect is of slightly more influence on the total electricity demand in this study, for the first 10 years of operation. However, as pointed out in the previous section, the energy losses of HT-ATES systems, and thus their electricity use, may vary strongly for different HT-ATES storage conditions (e.g., storage volume, hydrogeological conditions). Consequently, under more favourable storage conditions, yearly electricity use of the MES may be lower compared to direct heat production.
Energies 2021, 14, x FOR PEER REVIEW 28 of 33 electricity use of a MES with and without HT-ATES. The electricity use of scenario 50|30|1.5 is compared to a reference system that uses a HP for direct heat production only. Thus, the HT-ATES system enables the use of more sustainable and local electricity when available (in summer) and reduces the size of the heat pump compared to systems with direct heat production. The yearly electricity use slightly increases for the modelled scenarios in this study; however, this may change under more favourable HT-ATES storage conditions.

The LCOE Placed in Perspective
The LCOE for the different scenarios of heat pump size, condenser temperature and threshold temperature varies between 16.6-19.6 €/GJ for 65 °C heat pump condenser temperature, and between 12.2-17.1 €/GJ for 50 °C condenser temperature. Choosing a lower condenser temperature (50 °C instead of 65 °C ) leads to lower costs as the COP of the heat pump is higher and the costs of the heat pump are assumed to be lower. However, exact cost curves of (large scale) heat pumps are not available, so the cost difference is not known exactly and needs further investigation. Additionally, at 50 °C condenser temperature, the total stored volume in the HT-ATES is approximately twice as large as at 65 °C , and therefore requires more hydraulic pumping for the injection and abstraction.
The costs of the assessed scenarios in this paper are comparable to the costs (20.5 €/GJ) that were determined by Wesselink et al. [22], who considered a HT-ATES system in combination with geothermal energy. In an overview of different HT-ATES feasibility studies, the LCOE of heat varied between 14.7-29.3 €/GJ, with larger systems (> 350,000 m 3 ) resulting in a LCOE below 20 €/GJ [54]. The LCOE calculated in this paper thus fits within the price range known for HT-ATES systems. Heat delivered by a gas boiler costs on average 10-12 €/GJ [49]. For the scenarios considered in this paper, calculated prices are higher with a minimum difference of 0.2-2.2 €/GJ (50|30|1). This gap could be bridged by an increase in the CO2 price or subsidies on the installation or heat price of 100% renewable power-to-heat systems with HT-ATES. With a CO2-price of 55.6 €/ton CO2 and CO2-emission factor of 35.97 kg CO2/GJ of heat [55], the LCOE of a gas boiler would rise to 12-14 €/GJ, which would then fall in the same range as the 50|30|1 and 50|30|1.5 scenarios. At a CO2-price of 100 €/tonne CO2, the price range of heat from a natural gas boiler More importantly, the integration of HT-ATES has a strong effect on the distribution of electricity use over the year. The electricity use is decoupled by the integration of HT-ATES, and the electricity use peak switches from winter to summer, following the availability of sustainable electricity by solar energy (Figure 12). Here, we compare the electricity use of a MES with and without HT-ATES. The electricity use of scenario 50|30|1.5 is compared to a reference system that uses a HP for direct heat production only.
Thus, the HT-ATES system enables the use of more sustainable and local electricity when available (in summer) and reduces the size of the heat pump compared to systems with direct heat production. The yearly electricity use slightly increases for the modelled scenarios in this study; however, this may change under more favourable HT-ATES storage conditions.

The LCOE Placed in Perspective
The LCOE for the different scenarios of heat pump size, condenser temperature and threshold temperature varies between 16.6-19.6 €/GJ for 65 • C heat pump condenser temperature, and between 12.2-17.1 €/GJ for 50 • C condenser temperature. Choosing a lower condenser temperature (50 • C instead of 65 • C) leads to lower costs as the COP of the heat pump is higher and the costs of the heat pump are assumed to be lower. However, exact cost curves of (large scale) heat pumps are not available, so the cost difference is not known exactly and needs further investigation. Additionally, at 50 • C condenser temperature, the total stored volume in the HT-ATES is approximately twice as large as at 65 • C, and therefore requires more hydraulic pumping for the injection and abstraction.
The costs of the assessed scenarios in this paper are comparable to the costs (20.5 €/GJ) that were determined by Wesselink et al. [22], who considered a HT-ATES system in combination with geothermal energy. In an overview of different HT-ATES feasibility studies, the LCOE of heat varied between 14.7-29.3 €/GJ, with larger systems (>350,000 m 3 ) resulting in a LCOE below 20 €/GJ [54]. The LCOE calculated in this paper thus fits within the price range known for HT-ATES systems. Heat delivered by a gas boiler costs on average 10-12 €/GJ [49]. For the scenarios considered in this paper, calculated prices are higher with a minimum difference of 0.2-2.2 €/GJ (50|30|1). This gap could be bridged by an increase in the CO 2 price or subsidies on the installation or heat price of 100% renewable power-to-heat systems with HT-ATES. With a CO 2 -price of 55.6 €/ton CO 2 and CO 2 -emission factor of 35.97 kg CO 2 /GJ of heat [55], the LCOE of a gas boiler would rise to 12-14 €/GJ, which would then fall in the same range as the 50|30|1 and 50|30|1.5 scenarios. At a CO 2 -price of 100 €/tonne CO 2 , the price range of heat from a natural gas boiler would rise to 13.6-15.6, increasing the feasibility of a HT-ATES system with power-to-heat further.

Prolonged Lifetime of the Multi-Energy System with HT-ATES
In this study, only the first ten years of operation of the HT-ATES have been simulated. The lifetime of the multi-energy system and HT-ATES is longer than ten years, and the results show that most scenarios can reach a 100% heat delivery within the first ten years of operation (see Figures 6 and 7). After this time, the electricity demand is more representative compared to the first few years of operation as the system is still stabilizing and not able to provide the total heat demand. Hence, we performed an extra financial analysis assuming the electricity demand of the 10th year of operation would stay constant for the next 30 years of operation. This is a conservative estimation because the systems will most likely continue to improve in performance and system recovery efficiency over the years. On the other hand, we assumed that all systems would be able to provide the total heat demand of the neighbourhood after 10 years of operation, although we are not completely sure about this for i.e., scenario 65|43|1.5, 65|30|1.5, 65|30|1, and 50|30|1. The results can thus not be used to reliably compare the scenarios with each other, but give an idea of how the LCOE would develop over time. The results are shown in Table 7 and show a reduction of 1.2-5.6 €/GJ, which is a reduction in LCOE of 4-28% and 13% on average. The LCOE presented in this paper are thus a conservative estimation and will be significantly lower when considering a lifetime of 30 years, which would be further enhanced when a reduction in electricity consumption is taken into account as well.

Conclusions
In this study, it is shown that with the novel integration of a HT-ATES system to a multi-energy system, it is possible to provide >90% of the heat demand during the first ten years of operation, and supply the full heat demand after 3-5 years. The system recovery efficiency of the HT-ATES system was 50-70% depending on the system design chosen. By including a HT-ATES system in a multi-energy system, the size of the heat pump could be reduced by maximally 25% as the peak in heat demand can be delivered from the HT-ATES system, which reduces both the impact on the surface level (smaller installation) and costs. The integration of HT-ATES allows for the decoupling of energy demand and availability,