Aquifer Thermal Energy Storage (ATES) for District Heating and Cooling: A Novel Modeling Approach Applied in a Case Study of a Finnish Urban District

.


Introduction
According to Eurostat, in 2018, the share of renewable energy sources (RES) used for heating and cooling in EU was 21% and several countries, like Sweden (65%), Latvia (56%), Finland (55%) and Estonia (54%), covered more than half of their heating and cooling consumption with renewable sources [1].The variability of renewable generation between heating and cooling seasons, as well as the low coincidence between supply and demand are important challenges for RES penetration, therefore short-and long-term energy storage is needed for maximizing the usage of RES.Aquifer thermal energy storage (ATES) is an attractive technological option suitable for large buildings and utilities as well as capable to enable important storage capacities [2,3].Moreover, the utilization of GSHP operating within the urban subsurface space, is an efficient and resilient alternative for sustainable generation of heating and cooling energy in a district level [4].
The potential of ATES integration as a part of sustainable heating and cooling in combination with a ground-source heat pump (GSHP) for energy recovery from the subsurface has been acknowledged worldwide.Fleuchaus et al. [3] presented a complete overview of global ATES development and Energies 2020, 13, 2478 2 of 19 application: nowadays some 3000 ATES systems are operated worldwide.The Netherlands with 85% of all ATES realizations, followed by Sweden, Denmark and Belgium, are the undisputed frontrunners.Schmidt et al. [5] revealed that there are some 100 large-scale utility ATES systems utilized in district heating (DH) and cooling (DC) networks.
Normally, the long-term impact of ATES utilization is a combination of thermal, hydrological, microbiological and chemical impact on the affected aquifer and should be thoroughly investigated [6].The regulation of shallow geothermal plants (depth below 400 m) varies significantly among countries [7].Countries like Denmark, the Netherlands and Austria limit the lower and higher storage temperatures, whereas France and Switzerland establish a maximum fluctuation of groundwater temperature.Finland has no explicit legislative references to groundwater utilization for thermal storage, thus the findings of the present work can contribute for developing a specific normative framework in the future.
In the same line, the ground-source heat pump (GSHP) is a key technology for decarbonization of existing heating and cooling, which are nowadays mostly based on the use of fossil fuels [8][9][10].The work of Paiho et al. [8] revealed the importance of large-scale heat pumps for increasing the flexibility of Finnish energy systems.Within the same research, different examples are presented for heat pump integration in Finnish DH-DC networks-including the Kakola plant in Turku utilizing heat from sewage wastewater, and the Katri Vala plant in Helsinki generating heating and cooling in a single operation.
Fleuchaus et al. [11] evaluated the performance of ATES based on different criteria and concluded that ATES integration into heating and cooling systems was rarely addressed.In order to fill this gap, the integration of GSHP in tandem with ATES within the existing DH-DC networks of a Finnish urban district is presented and developed in the current case study.The main objective of this work is to propose a mathematical modeling of the whole ATES-GSHP-DH-DC energy chain in order to improve the system's energy management, as well as to study its technical and economic feasibility and the long-term environmental impact.Finnish public data sources are available, like the Finnish Environmental Institute (SYKE) regarding the hydrological resources, the Geological Survey of Finland (GTK) on hydrogeological conditions, and the National Land Survey of Finland (NLSF) for geographical data.The present research also introduces a methodology for fetching data from the aforementioned sources in order to calibrate and validate a groundwater model of the studied area, which in turn is an indispensable tool for studying the ATES-GSHP impact in the long-term.

Materials and Methods
The modeling procedure of the combined ATES-GSHP-DH-DC system, depicted in Figure 1, is based on the following steps, namely, (i) input data of the target DH-DC networks and the nearby groundwater areas, (ii) perform mathematical modeling of combined ATES-GSHP operation, (iii) undertake technoeconomic and sensitivity analysis, and iv) study the impact of ATES operation on aquifer areas, by developing and calibrating a specific groundwater model.A groundwater model based on the finite difference method code MODFLOW (Harbaugh et al. [12]) has been adopted and developed in the present case study.The model is calibrated against long-term data (hydraulic heads of the observation wells).The particular case study is introduced in Section 2.1, while a detailed explanation and demonstration of the modeling procedure is presented in Sections 2.2 and 2.4.

Figure 1.
General modeling procedure of the system based on aquifer thermal energy storage (ATES), ground-source heat pumps (GSHP), district heating (DH), and district cooling (DC).

Input Data of the DH and DC Networks
The target district heating and cooling networks are located in the central district of Kupittaa in the town of Turku, located in the southwest part of Finland.The available data is hourly-based and the most relevant parameters of both DH and DC networks are summarized in Table 1 and Figure 2.   The target district heating and cooling networks are located in the central district of Kupittaa in the town of Turku, located in the southwest part of Finland.The available data is hourly-based and the most relevant parameters of both DH and DC networks are summarized in Table 1 and Figure 2.

Input Data of the Groundwater Areas
Available open data from Finnish public sources was retrieved for characterizing the target groundwater areas.In this research, data from Finnish Environment Institute [13] was used relative to groundwater areas, monitoring stations, and observation wells.Available information was utilized for 15 observation wells located in HK-Ruokatalo area of Kupittaa district and eight wells in the Kaarninko area, and their long-term statistical data for average head were used for steady-state model calibration.In the HK-Ruokatalo area, the average measured head was 16.7 m above sea level with standard deviation 0.7 m, while for Kaarninko area these values were 19.3 m and 0.5 m respectively.The average aquifer thickness was estimated as 10 m (Joronen [14]) and the maximum allowed daily pumping rate was 2500 m 3 /day.The undisturbed aquifer temperature in HK-Ruokatalo area was around 10 ℃, quite high due to the subsurface heat island effect observed in cities (Bayer et al. [15]).

Geographical Data
Open data from the National Land Survey of Finland [16] was used, particularly its "10 m elevation model".The elevation model was retrieved as Geo-TIFF raster file and transformed to Surfer Grid file (GRD) using QGIS software [17].The Aura River, a half-kilometer to the northwest, as well as the Baltic Sea, located several kilometers to southwest establish the hydrogeological boundaries of the groundwater model.

ATES-GSHP Integration for District Heating and Cooling
A ground-source heat pump (GSHP), operating with an abstraction and injection well (well doublet) was considered.The condenser side of the heat pump is connected to DH network while the evaporator side is connected to aquifer pumping stream.
In the base (first) scenario, the ATES pumping flow path encounters two serial exchangers-HP evaporator and cooling for DC network.In the second scenario before the HP evaporator, a precooling exchanger is added, providing a first stage cooling to the DC network.As will be shown in the result section, with this configuration the DC demand can be more efficiently covered and GSHP efficiency (COP) can be improved since heat pump inlet temperature increases several degrees after a precooling exchanger.ATES-GSHP integration within the existing DH-DC networks is

Input Data of the Groundwater Areas
Available open data from Finnish public sources was retrieved for characterizing the target groundwater areas.In this research, data from Finnish Environment Institute [13] was used relative to groundwater areas, monitoring stations, and observation wells.Available information was utilized for 15 observation wells located in HK-Ruokatalo area of Kupittaa district and eight wells in the Kaarninko area, and their long-term statistical data for average head were used for steady-state model calibration.In the HK-Ruokatalo area, the average measured head was 16.7 m above sea level with standard deviation 0.7 m, while for Kaarninko area these values were 19.3 m and 0.5 m respectively.The average aquifer thickness was estimated as 10 m (Joronen [14]) and the maximum allowed daily pumping rate was 2500 m 3 /day.The undisturbed aquifer temperature in HK-Ruokatalo area was around 10 • C, quite high due to the subsurface heat island effect observed in cities (Bayer et al. [15]).

Geographical Data
Open data from the National Land Survey of Finland [16] was used, particularly its "10 m elevation model".The elevation model was retrieved as Geo-TIFF raster file and transformed to Surfer Grid file (GRD) using QGIS software [17].The Aura River, a half-kilometer to the northwest, as well as the Baltic Sea, located several kilometers to southwest establish the hydrogeological boundaries of the groundwater model.

ATES-GSHP Integration for District Heating and Cooling
A ground-source heat pump (GSHP), operating with an abstraction and injection well (well doublet) was considered.The condenser side of the heat pump is connected to DH network while the evaporator side is connected to aquifer pumping stream.
In the base (first) scenario, the ATES pumping flow path encounters two serial exchangers-HP evaporator and cooling for DC network.In the second scenario before the HP evaporator, a precooling exchanger is added, providing a first stage cooling to the DC network.As will be shown in the result section, with this configuration the DC demand can be more efficiently covered and GSHP efficiency (COP) can be improved since heat pump inlet temperature increases several degrees after a precooling exchanger.ATES-GSHP integration within the existing DH-DC networks is depicted in the general scheme presented in Figure 3, where temperature values illustrate the second scenario setup.

GSHP Utilization for District Heating
Generally speaking, GSHP is utilized to recover and upgrade all excess heat proceeding from the DC network and inject it in DH network.In this context ATES is utilized for balancing the energy system and mitigating the variability and no-coincidence of the simultaneously dispatched heating and cooling loads.For this purpose, heat pump supply temperature is calculated, based on the demanded power fraction k (the ratio between heat supplied by the heat pump and total heat demanded in the DH branch).The flow fraction recirculated through HP condenser can be calculated as: k p , where 0 ≤ p ≤ 1 is additional exponent parameter, thus the ratio between heat pump condenser THPC and total TDH of DH network is k 1-p .Therefore, for each hour n, given that TDH,R,n and TDH,S,n are DH return and supply temperatures respectively, the heat pump supply temperature THPC,S,n can be calculated as follows: ,, =  ,, + ( ,, −  ,, ) 1− ⇒ ∆ , = ∆ ,  1− (1) The resulting supply temperature TDH,S2,n after mixing can be calculated as: In the present case study, the exponential parameter p was chosen to equal 0.6.It can be observed that there is a significant advantage in partial load operation, since there is no need to increase HP supply temperature as high as DH supply, thus boosting the COP.For example, for power fraction k = 0.4 and TDH = 40 ℃, the GSHP should elevate DH return temperature by roughly 28 ℃ instead of 40 ℃.After mixing with the supply DH flow, the maximum temperature drop of the flow TDHS is ≈ 7 ℃.

COPH Estimation Model
For industrial and large-scale processes, multiple HP units in serial connection increase overall system efficiency, and therefore the Lorentz COP [9,18] would describe more accurately the behavior of the HP configuration, since it takes into account the logarithmic mean temperature of the sink and source, as well as both inlet and outlet temperatures of the condenser and evaporator.According to

GSHP Utilization for District Heating
Generally speaking, GSHP is utilized to recover and upgrade all excess heat proceeding from the DC network and inject it in DH network.In this context ATES is utilized for balancing the energy system and mitigating the variability and no-coincidence of the simultaneously dispatched heating and cooling loads.For this purpose, heat pump supply temperature is calculated, based on the demanded power fraction k (the ratio between heat supplied by the heat pump and total heat demanded in the DH branch).The flow fraction recirculated through HP condenser can be calculated as: k p , where 0 ≤ p ≤ 1 is additional exponent parameter, thus the ratio between heat pump condenser ∆T HPC and total ∆T DH of DH network is k 1-p .Therefore, for each hour n, given that T DH,R,n and T DH,S,n are DH return and supply temperatures respectively, the heat pump supply temperature T HPC,S,n can be calculated as follows: The resulting supply temperature T DH,S2,n after mixing can be calculated as: In the present case study, the exponential parameter p was chosen to equal 0.6.It can be observed that there is a significant advantage in partial load operation, since there is no need to increase HP supply temperature as high as DH supply, thus boosting the COP.For example, for power fraction k = 0.

COP H Estimation Model
For industrial and large-scale processes, multiple HP units in serial connection increase overall system efficiency, and therefore the Lorentz COP [9,18] would describe more accurately the behavior of the HP configuration, since it takes into account the logarithmic mean temperature of the sink and source, as well as both inlet and outlet temperatures of the condenser and evaporator.According to Energies 2020, 13, 2478 6 of 19 Reinholdt et al. [18], the maximum theoretical COP of a heat pump can be estimated by calculating Lorentz COP, defined as follows: , where T lm,H = T HPC,S − T HPC,R ln T HPC,S T HPC,R ; In Equation (3), T lm,H and T lm,L are, respectively, the logarithmic mean temperature of the sink and source, where notations HPC and HPE stand for heat pump's condenser and evaporator temperatures, while notations I/O stand for inlet/outlet temperatures of the evaporator and S/R stand for supply/return temperatures of the condenser (all values expressed in Kelvin).Based on the best industrial refrigeration systems, Reinholdt et al. [18] suggested values for Lorentz efficiency between 50% and 60% of the maximum Lorenz COP calculated with Equation (3).In our case study, a more conservative value of 45% was adopted.

GSHP Utilization for District Cooling
As mentioned previously, part of DC demand can be produced by free cooling in a first stage cooling exchanger located at the beginning of ATES pumping flow.After that, GSHP is utilized in the second place for simultaneously cooling the ATES flow in the evaporator as well as supplying heat to DH network in the condenser (see Figure 3).Finally, second stage cooling is applied, and groundwater is injected into the aquifer.
For each hour of operation, it is crucial to determine the exact aquifer pumping flow rate Q [m 3 /s] since there is constraint for daily pumping of 2500 m 3 /day.Due to this limitation, the maximum heat output of the GSHP condenser is limited to 1.4 and 1.6 MW in scenario 1 and 2 respectively, and pumping flow rate is calculated according to the iterative algorithm developed below.

Computation of ATES Hourly Pumping Rate
Since there are several exchangers (two and three, respectively, for scenario 1 and 2) in the ATES flow path, the minimum needed pumping flow rate is proposed to be estimated iteratively.If Φ heat,n and Φ cool,n are, respectively, heating and cooling demand to be covered in hour n, as the first estimation of the pumping flow can be taken the maximum flow needed either for heating or cooling (notations according to Figure 3): Once the first estimation for Q n is known, it is possible to calculate separately all exchangers within the ATES flow path, in both scenarios 1 and 2, as follows.Scenario 1: Recalculation of temperature after HP evaporator: Step 2 : Scenario 2: Recalculation of first and second stage cooling demands: Step The ATES flow is recalculated again in Step 1, and if the new value deviates more than a predefined threshold from the previous one (in this case a 5% threshold is adopted), then the whole loop (Step 1/Step 2) is repeated.

Calculation of ATES Pumping Power Demand
The required pumping power [kW] for ATES operation can be calculated on an hourly basis, assuming overall pressure drop in the line ∆p = 600 kPa and standard pumping efficiency η = 0.55 [19], as follows:

Calculation of Pumping Power Demand to DH-DC Network
Similarly, pumping power [kW] to provide DH-DC through the GSHP condenser/evaporator respectively can be calculated hourly, assuming overall pressure drop between supply and return lines ∆p DH = ∆p DC = 250 kPa [20] and standard pumping efficiency η = 0.55 [21], as follows: where The volumetric heat capacity of water S VC,wat used was 4.19 and 4.1 MJ/m 3 K, respectively, for cooling and heating operation.

Numerical Model and Its Calibration for Steady State
The groundwater model is set up utilizing the finite difference code MODFLOW [12] with ModelMuse environment [22].In ModelMuse, the aquifer is discretized with a 100 x 100 m square cell grid, covering a physical extension of about 20 km 2 , delimited between the Aura River to the northwest and the Baltic Sea to the southwest.Southeast and northeast borders are assumed as no-flow boundaries (see Figure 4).
Groundwater model calibration for steady state was carried out taking into account the long-term statistical data for 15 observation wells in the Kupittaa area and eight observation wells in Kaarninko.Calibration was done according to the procedure developed by Todorov et al. [23], by using root mean squared error (RMSE) [24] and mean absolute error (MAE) [25] for the close field (Kupittaa) and far field (Kaarninko).As seen in Figure 5, the results of Kaarninko (far-field area) were more dispersed (calculated RMSE = 1.32 m/MAE = 1.07 m), since our model is intended to present better correlation between measured and simulated values (RMSE = 0.54 m/MAE = 0.29 m) within the close-field calibration.In most of Kupittaa's observation wells, this difference was within the margins of the measured long-term standard deviation.A typical horizontal hydraulic conductivity for sand/gravel aquifer was selected: K = 5 × 10 −5 m/s (Luoma [26]), and during model calibration was adjusted to 5 × 10 −4 m/s for the area containing the observation wells (small black rhombs in Figure 4 delimited by circles).The value of vertical hydraulic conductivity was chosen as K z = 0.1K.Typical values were also utilized for storativity (S = 1 × 10 −5 ), porosity (n = 0.25) and recharge rate of R = 1.3 × 10 −8 m/s [26].
adjusted to 5 × 10 -4 m/s for the area containing the observation wells (small black rhombs in Figure 4 delimited by circles).The value of vertical hydraulic conductivity was chosen as Kz = 0.1K.Typical values were also utilized for storativity (S = 1 × 10 -5 ), porosity (n = 0.25) and recharge rate of R = 1.3 × 10 -8 m/s [26].In Figure 4, the calibrated groundwater model and its steady state solution are depicted, where iso-lines represent groundwater head in meters while the color grid represent elevation between blue (min.elevation, 0.3 m above sea level) and red (max.elevation 59.6 m).As seen in Figure 4, the natural groundwater flow moves from the Kaarninko area (highest hydraulic heads 19-20 m) towards the Kupittaa area (hydraulic heads 16-17 m), and finally reaches the lowest hydraulic boundaries represented by the Aura River and Baltic Sea.In Figure 4, the calibrated groundwater model and its steady state solution are depicted, where iso-lines represent groundwater head in meters while the color grid represent elevation between blue (min.elevation, 0.3 m above sea level) and red (max.elevation 59.6 m).As seen in Figure 4, the natural groundwater flow moves from the Kaarninko area (highest hydraulic heads 19-20 m) towards the Kupittaa area (hydraulic heads 16-17 m), and finally reaches the lowest hydraulic boundaries represented by the Aura River and Baltic Sea.

Technoeconomic Evaluation of GSHP-ATES
Based on hourly calculations, different technical variables are computed, like the annual energy generation for heating/cooling, the electricity consumption and the average daily ATES pumping rate.A cost database regarding various energy generation technologies was used (after Nielsen et al. [27,28]), as well as prices for ATES well drilling, heat exchangers, and piping (Drenkelfort et al. [29]) for estimating the investment cost.Based on the annuity method, the energy generation cost is

Technoeconomic Evaluation of GSHP-ATES
Based on hourly calculations, different technical variables are computed, like the annual energy generation for heating/cooling, the electricity consumption and the average daily ATES pumping rate.A cost database regarding various energy generation technologies was used (after Nielsen et al. [27,28]), as well as prices for ATES well drilling, heat exchangers, and piping (Drenkelfort et al. [29]) for estimating the investment cost.Based on the annuity method, the energy generation cost is calculated, assigning annual investment payments (annuity) and assuming 5% interest rate as well as the investment lifetime of 20 years (Nielsen et al. [27]).The operation and maintenance (O&M) costs (1% of investment) are also included within the overall annual cost, as well as the electricity cost for GSHP and pumping (given electricity price of 100 €/MWh, including taxes, transfer and distribution fees [30]).The economic evaluation was developed according to Todorov et al. [23], including the calculation of the following variables listed in Table 3.

Technoeconomic Analysis
The main technical parameters of ATES operation for both studied scenarios are shown in Table 4.It can be acknowledged that even with 5%-6% of peak heat power for scenario 1 and 2, the GSHP coverage ratio is 18%-20% of the annual heating demand.Moreover, an important advantage of scenario 2 is shown when comparing a cooling demand covered by ATES.The scheme with two cooling exchangers in scenario 2 allows 78% coverage of DC demand annually (compared to 67% in scenario 1), from which the first stage cooling represents roughly one sixth.The estimation of economic feasibility parameters and the production cost of thermal energy are shown in Tables 5 and 6 respectively.The resulting thermal energy production cost in scenario 2 is slightly below 30 €/MWh.Overall investment cost is around 2.3 million €: 26% corresponds to GSHP/exchangers and 73% is related to the underground components (connection pipes and wells), figures close to similar ATES realization in Germany (Schüppler et al. [31]).The specific investment cost per installed heat pump capacity is 1.6/1.4€/W for scenario 1 and 2 respectively, values comparable to the 1.8 €/W reported for a similar ATES system in a Belgian hospital (Vanhoudt et al. [32]).Additionally, scenario 2 is investigated with more details, as follows.GSHP COP is 3.2 on average, slightly improving to 3.3 during the winter due to lower GSHP supply temperature (64 • C on average), while, during the summer, GSHP covers a higher heat fraction and the average supply temperature increases to 69 • C (see Figure 6).
ATES operation is based on energy conversion using electricity to cogenerate heating and cooling in a single operation.GSHP is the principal electricity consumer accounting for 90% of the annual demand, followed by ATES pumping (6%) as well as pumping needed to inject HP supply energy to DH-DC networks-respectively 1% and 3%.This is important to acknowledge since total electricity demand (4.8 GWh/a) has a significant impact on the annual cost, and, consequently, on the specific cost of generated heating and cooling energy, as seen in Table 6.The ATES system is well balanced, as seen from the average injection and abstraction temperatures that are both equal to the aquifer's undisturbed temperature of 10 • C.Moreover, the system is balanced in terms of energy, as shown in Table 4, since the annual heat demand covered is equal to cooling demand covered plus GSHP power demand (13.9 GWh). Figure 7 depicts the annual variation of all temperatures along the ATES flow path: abstraction, after first stage cooling, after GSHP evaporator, and finally injection.ATES operation is based on energy conversion using electricity to cogenerate heating and cooling in a single operation.GSHP is the principal electricity consumer accounting for 90% of the annual demand, followed by ATES pumping (6%) as well as pumping needed to inject HP supply energy to DH-DC networks-respectively 1% and 3%.This is important to acknowledge since total electricity demand (4.8 GWh/a) has a significant impact on the annual cost, and, consequently, on the specific cost of generated heating and cooling energy, as seen in Table 6.The ATES system is well balanced, as seen from the average injection and abstraction temperatures that are both equal to the aquifer's undisturbed temperature of 10 ℃.Moreover, the system is balanced in terms of energy, as shown in Table 4, since the annual heat demand covered is equal to cooling demand covered plus GSHP power demand (13.9 GWh). Figure 7 depicts the annual variation of all temperatures along the ATES flow path: abstraction, after first stage cooling, after GSHP evaporator, and finally injection.ATES operation is based on energy conversion using electricity to cogenerate heating and cooling in a single operation.GSHP is the principal electricity consumer accounting for 90% of the annual demand, followed by ATES pumping (6%) as well as pumping needed to inject HP supply energy to DH-DC networks-respectively 1% and 3%.This is important to acknowledge since total electricity demand (4.8 GWh/a) has a significant impact on the annual cost, and, consequently, on the specific cost of generated heating and cooling energy, as seen in Table 6.The ATES system is well balanced, as seen from the average injection and abstraction temperatures that are both equal to the aquifer's undisturbed temperature of 10 ℃.Moreover, the system is balanced in terms of energy, as shown in Table 4, since the annual heat demand covered is equal to cooling demand covered plus GSHP power demand (13.9 GWh). Figure 7 depicts the annual variation of all temperatures along the ATES flow path: abstraction, after first stage cooling, after GSHP evaporator, and finally injection.

Sensitivity Analysis of the System's Operation
As shown in Table 6, about 70% of energy production cost is related to electricity consumption, of which the GSHP accounts for around 90%.The heat pump's COP is an important variable to consider in order to boost the system's efficiency and decrease cost.That is why, in this section, a sensitivity analysis will be performed regarding COP and energy production cost, and how they depend on the exponent parameter p.The effect of varying p within the interval [0;1] is that, e.g., for p = 1, flow recirculated through GSHP condenser is directly proportional to power fraction k, and thus heat pump supply temperature should be equal to DH supply temperature.Figure 8 plots a ∆T HPC /∆T DH fraction of GSHP condenser calculated with Equation (1) and a temperature drop fraction after heat pump junction ∆T DH,S /∆T DH calculated with Equation ( 2) for different values of p = 0.2, 0.4, 0.6, 0.8.The comparative thermal effect for ∆T DH = 40 • C is illustrated in the secondary vertical axis.From Figure 8, it can be seen that, for lower values of p, the GSHP has higher efficiency when working at power fractions (e.g., during the winter period) since HP supply temperature is not so high as when working at lower power fractions (e.g., during the winter period) since HP supply temperature is not so high as DH supply.The drawback is that, after HP junction, DH supply temperature TDH,S2 can also present an important temperature drop ΔTDH,S (e.g., red dashed curve, for p = 0.2).It is also interesting to explore what is the maximum ΔTDH,S for each p within the interval [0;1].Let's define the following function f(k), as the ratio between ΔTDH,S and ΔTDH, according to Equation ( 2): As seen from Figure 8, in the interval [0;1], f(k) has one maximum, which can be found where the function's first derivative is zero: Similarly, it is possible to calculate the average value of f(k) within [0;1] as: As seen from Figure 8, in the interval [0;1], f(k) has one maximum, which can be found where the function's first derivative is zero: Similarly, it is possible to calculate the average value of f(k) within [0;1] as: The results for f max and f avg calculated respectively with Equations ( 8) and ( 9) are presented in Figure 9.
The results for fmax and favg calculated respectively with Equations ( 8) and ( 9) are presented in Figure 9.As previously noticed, and also presented in Figures 8 and 9, for low values of p, the temperature drop after HP junction can increase significantly.Moreover, when p = 0.2, for the maximum and average temperature drop, the fraction is as high as 0.53 and 0.33 of ΔTDH respectively, while for p = 0.8 the possible gains in low power fraction are quite limited.Therefore, a sensitivity analysis for scenario 2 was performed, for four different values of p: All cases are simulated on an hourly basis, with the same constraint for ATES average daily pumping flow, and the annual results are listed in Table 7.By decreasing p, the average GSHP supply temperature also decreases, while the average COP increases, and vice versa.From Table 7, it also can be seen that the higher the COP, the lower GSHP heat capacity needed, since the fraction supplied by the HP compressor is lower.The latter also explains why the annual heat demand covered is lower for high COP, while the annual cooling demand covered remains stable among the four cases.As previously noticed, and also presented in Figures 8 and 9, for low values of p, the temperature drop after HP junction can increase significantly.Moreover, when p = 0.2, for the maximum and average temperature drop, the fraction is as high as 0.53 and 0.33 of ∆T DH respectively, while for p = 0.8 the possible gains in low power fraction are quite limited.Therefore, a sensitivity analysis for scenario 2 was performed, for four different values of p: All cases are simulated on an hourly basis, with the same constraint for ATES average daily pumping flow, and the annual results are listed in Table 7.By decreasing p, the average GSHP supply temperature also decreases, while the average COP increases, and vice versa.From Table 7, it also can be seen that the higher the COP, the lower GSHP heat capacity needed, since the fraction supplied by the HP compressor is lower.The latter also explains why the annual heat demand covered is lower for high COP, while the annual cooling demand covered remains stable among the four cases.
The percentage variations compared to the base case 3 are plotted in Figure 10.It is important to notice how the energy production cost decreases as COP increases.On the other hand, an important drawback for cases 1 and 2 is the high temperature drop in DH supply temperature, 19.1 • C and 11.6 • C respectively, which is a confirmation that case 3 is a reasonable trade-off between the system's efficiency, economic feasibility, and technical constraints.The percentage variations compared to the base case 3 are plotted in Figure 10.It is important to notice how the energy production cost decreases as COP increases.On the other hand, an important drawback for cases 1 and 2 is the high temperature drop in DH supply temperature, 19.1 ℃ and 11.6 ℃ respectively, which is a confirmation that case 3 is a reasonable trade-off between the system's efficiency, economic feasibility, and technical constraints.

Impact on Groundwater Areas
Although the undisturbed aquifer temperature is as high as 10 ℃, first stage cooling can be used in 8736 out of 8760 h, and annually it represents 17% of cooling demand covered by ATES (about one sixth of 9.6 GWh).This configuration also increases the temperature before GSHP evaporator by 1.5 ℃ on average, which improves the COP and enhances the heat pump's capacity in the evaporator as well.The average injection temperature lies in a narrow range of roughly 10 ± 1 ℃ (see Table 4), which justifies one-way ATES operation and, consequently, the thermal impact on the aquifer remains very limited.
The long-term hydraulic impact was simulated in MODFLOW by taking a weekly-based average for ATES pumping rate and defining 52 × 20 = 1040 stress periods.The result after 20 years of one-way operation is presented in Figure 11, where hydraulic head is represented by iso-lines with resolution of 0.25 m.In order to mitigate the hydraulic impact of pumping, the injection well is placed downstream while the abstraction well is located upstream (Figure 11).The maximum simulated drawdown is 1.28 and 1.17 m for summer and winter operation respectively, which corresponds to 5.0 and 4.7 m inside the pumping well.The overall impact of ATES pumping vanishes at about 500 m from each well, thus it does not affect the surrounding groundwater areas in a significant way.

Impact on Groundwater Areas
Although the undisturbed aquifer temperature is as high as 10 • C, first stage cooling can be used in 8736 out of 8760 h, and annually it represents 17% of cooling demand covered by ATES (about one sixth of 9.6 GWh).This configuration also increases the temperature before GSHP evaporator by 1.5 • C on average, which improves the COP and enhances the heat pump's capacity in the evaporator as well.The average injection temperature lies in a narrow range of roughly 10 ± 1 • C (see Table 4), which justifies one-way ATES operation and, consequently, the thermal impact on the aquifer remains very limited.
The long-term hydraulic impact was simulated in MODFLOW by taking a weekly-based average for ATES pumping rate and defining 52 × 20 = 1040 stress periods.The result after 20 years of one-way operation is presented in Figure 11, where hydraulic head is represented by iso-lines with resolution of 0.25 m.In order to mitigate the hydraulic impact of pumping, the injection well is placed downstream while the abstraction well is located upstream (Figure 11).The maximum simulated drawdown is 1.28 and 1.17 m for summer and winter operation respectively, which corresponds to 5.0 and 4.7 m inside the pumping well.The overall impact of ATES pumping vanishes at about 500 m from each well, thus it does not affect the surrounding groundwater areas in a significant way.

Conclusion
The presented case study was successful in demonstrating and developing a mathematical model for system's management: calculation of HP recirculation flow, estimation of heat pump COP, as well as an algorithm for computation of ATES pumping flow rate based on the capacity to cover heating and cooling demand in a single operation.Additionally, the system's technoeconomic feasibility, efficiency, and the impact of GSHP-ATES operation on the nearby aquifer were evaluated.The groundwater model was developed and calibrated, utilizing different available data sources like the National Land Survey of Finland (NLSF), Finnish Environment Institute (SYKE), Geological Survey of Finland (GTK) and different tools such as MS Excel, QGIS and MODFLOW.
The dispatch of combined heating and cooling loads using annual data for an existing urban district in Finland was developed in tandem with the GSHP-ATES model.It presented an attractive economic outcome-competitive energy production cost around 30 €/MWh, far below 76.7 €/MWh, which was the average Finnish DH price in 2017 [33], as well as very limited long-term environmental impact on the nearby aquifer.The maximum drawdown within the pumping well was estimated as 5 m after 20 years of operation, and the overall hydraulic impact is limited to 500 m around the wells.Injection temperature deviates from undisturbed aquifer temperature at a level of about 1 ℃ on average, which is within the limits of Swiss (3 ℃) and French (11 ℃) legislation.Additional sensitivity analysis revealed that, by varying HP recirculation flow parameter (exponent parameter p), it is possible to influence the heat pump's COP and the energy production cost.However, the only constraint to be considered is the temperature drop after the HP junction in the DH supply line, which in the base case resulted in 6.4 ℃ on average.The future transition to low district heating networks (Guzzini et al. [34]) by the introduction of heat pumps can eventually benefit from the proposed mathematical methodology due to its capability to find a trade-off between the energy production cost and the maximum allowed temperature drop introduced by the heat pump in the DH supply line.
The environmental assessment presented in this work was carried out using available hydrogeological data of the investigated groundwater area.The lack of reliable hydrogeological information at the present prefeasibility stage was a significant constraint for simplifying the groundwater model-the aquifer was assumed as uniform, isotropic and confined to the considered domain, and model calibration was performed only for steady state.Additional pumping tests and more detailed geological exploration would be needed as future steps.

Conclusions
The presented case study was successful in demonstrating and developing a mathematical model for system's management: calculation of HP recirculation flow, estimation of heat pump COP, as well as an algorithm for computation of ATES pumping flow rate based on the capacity to cover heating and cooling demand in a single operation.Additionally, the system's technoeconomic feasibility, efficiency, and the impact of GSHP-ATES operation on the nearby aquifer were evaluated.The groundwater model was developed and calibrated, utilizing different available data sources like the National Land Survey of Finland (NLSF), Finnish Environment Institute (SYKE), Geological Survey of Finland (GTK) and different tools such as MS Excel, QGIS and MODFLOW.
The dispatch of combined heating and cooling loads using annual data for an existing urban district in Finland was developed in tandem with the GSHP-ATES model.It presented an attractive economic outcome-competitive energy production cost around 30 €/MWh, far below 76.7 €/MWh, which was the average Finnish DH price in 2017 [33], as well as very limited long-term environmental impact on the nearby aquifer.The maximum drawdown within the pumping well was estimated as 5 m after 20 years of operation, and the overall hydraulic impact is limited to 500 m around the wells.Injection temperature deviates from undisturbed aquifer temperature at a level of about 1 • C on average, which is within the limits of Swiss (3 • C) and French (11 • C) legislation.Additional sensitivity analysis revealed that, by varying HP recirculation flow parameter (exponent parameter p), it is possible to influence the heat pump's COP and the energy production cost.However, the only constraint to be considered is the temperature drop after the HP junction in the DH supply line, which in the base case resulted in 6.4 • C on average.The future transition to low district heating networks (Guzzini et al. [34]) by the introduction of heat pumps can eventually benefit from the proposed mathematical methodology due to its capability to find a trade-off between the energy production cost and the maximum allowed temperature drop introduced by the heat pump in the DH supply line.
The environmental assessment presented in this work was carried out using available hydrogeological data of the investigated groundwater area.The lack of reliable hydrogeological information at the present prefeasibility stage was a significant constraint for simplifying the groundwater model-the aquifer was assumed as uniform, isotropic and confined to the considered domain, and model calibration was performed only for steady state.Additional pumping tests and more detailed geological exploration would be needed as future steps.

Figure 1 .
Figure 1.General modeling procedure of the system based on aquifer thermal energy storage (ATES), ground-source heat pumps (GSHP), district heating (DH), and district cooling (DC).

Figure 3 .
Figure 3. General scheme of ATES integration in DH-DC networks.

Figure 3 .
Figure 3. General scheme of ATES integration in DH-DC networks.
4 and ∆T DH = 40 • C, the GSHP should elevate DH return temperature by roughly 28 • C instead of 40 • C.After mixing with the supply DH flow, the maximum temperature drop of the flow ∆T DHS is ≈ 7 • C.
The drawback is that, after HP junction, DH supply temperature T DH,S2 can also present an important temperature drop ∆T DH,S (e.g., red dashed curve, for p = 0.2).It is also interesting to explore what is the maximum ∆T DH,S for each p within the interval [0;1].Let's define the following function f(k), as the ratio between ∆T DH,S and ∆T DH , according to Equation (2):

Figure 9 .
Figure 9. Maximum and average values for f function.

Figure 9 .
Figure 9. Maximum and average values for f function.

Figure 11 .
Figure 11.Hydraulic impact after 20 years of ATES operation in scenario 2.

Figure 11 .
Figure 11.Hydraulic impact after 20 years of ATES operation in scenario 2.

Table 1 .
Relevant DH-DC network parameters of Kupittaa district in Turku.

Table 1 .
Relevant DH-DC network parameters of Kupittaa district in Turku.

Table 2 .
Technical variables of ATES

Table 2 lists
the relevant ATES technical variables.

Table 2 .
Technical variables of ATES.

Table 3 .
Variables for economic evaluation.

Table 4 .
ATES system technical parameters.

Table 6 .
Energy production cost.

Table 7 .
Sensitivity analysis based on four cases (case 3 is the base case).

Table 7 .
Sensitivity analysis based on four cases (case 3 is the base case).