Sustainable Cities and Society A method and analysis of aquifer thermal energy storage (ATES) system for district heating and cooling: A case study in Finland

Aquifer thermal energy storage (ATES) systems with groundwater heat pumps (GWHP) provide a promising and e ﬀ ective technology to match the renewable energy supply and demand between seasons. This paper analyses the integration of an ATES and GWHP system in both district heating (DH) and district cooling (DC) networks in terms of system's e ﬃ ciency, techno-economic feasibility and impact on the surrounding groundwater areas. To that end, a novel method of holistic integration of groundwater modeling is proposed and demonstrated for retrieving and analyzing data from a variety of open Finnish public data sources. A case study is presented, where the ATES integration is examined within an existing district heating network in Southern Finland. It is concluded that combining heating and cooling, with seasonally reversible ATES operation and balanced pumping volumes during summer and winter periods, had low impact on the aquifer area and is economically feasible. Finally, the study concludes that even with limited data, obtained from open public data sources, it is possible to assess the ATES integration with an acceptable accuracy.


Introduction
Buildings in Europe account for 41 % of the final energy consumption, followed by transport (32 %), and industry (25 %) (EU, 2013). In line with the Paris Agreement, carbon neutrality should be obtained by 2050 (UN, 2015). Hence, the integration of renewable energy technologies in heating and cooling of buildings and communities is a necessity. Since the variations of the availability of renewable energy between heating and cooling seasons are significant, a seasonal storage is needed to maximize the usage of renewable energy. One of the most promising technological options is Aquifer Thermal Energy Storage (ATES) due to its relative affordability and ability to enable large storage capacities (Pellegrini et al., 2019). Furthermore, by utilizing the available subsurface space in cities, ATES systems can potentially provide sustainable heating and cooling energy for different building typologies, thus the ATES integration in a district/urban level can be more efficient and resilient compared to a conventional separate heating and cooling generation (Hooimeijer & Maring, 2018). ATES heat storage is based on mature technologies, with groundwater heat pump (GWHP) as a core of an ATES system. The research interest in shallow (< 300 m of depth) geothermal energy systems is more recent, dating back to the 1970's (Sanner, Systems, Systems, & Pumps, 2001).
The potential of ATES technology to be integrated as a sub-system of sustainable space heating and cooling in tandem with the recovery of heat in the subsurface has been acknowledged worldwide. According to Fleuchaus, Godschalk, Stober, and Blum (2018), there were around 3000 ATES systems in the world in 2017, mostly concentrated in Europe. The Nordic countries and particularly the Netherlands (2500), Sweden (220) and Denmark (55) are the frontrunners of ATES application. According to Schmidt, Pauschinger, Sørensen, Snijders, and Thornton (2018), from these 3000 ATES projects, there are 100 largescale systems integrated in district heating and cooling (DH/DC) networks. The research of Paiho et al. (2018) revealed that large-scale heat pumps can increase the flexibility of the district heating system, especially when utilizing heat pumps for first-stage preheating before the CHP-plant finally adjusts the supply temperature. In the same study are given several examples for heat pump utilization in DH/DC networks in Finlandthe Kakola plant in Turku utilizing heat from waste water and the Katri Vala plant in Helsinki producing DH and DC in a single process. Bloemendal, Olsthoorn, and van de Ven (2015) developed a method for identifying the available world ATES potential combining climatic and hydro-geological data, as well as elaborated a world map for ATES suitability. The study concluded that some 50 % of world urban areas have medium potential for ATES (remaining stable among the present century), while 15 % have high potential -a figure which will decrease to 5 % in the second half of 21 st century due to climate change. Lu, Tian, and He (2019) presented similar approach for evaluating world ATES potential based on socio-economic, geo-hydrological, climatic and groundwater factors, as well as concluded that ATES potential is very good, good and moderate in 7 %, 20 % and 34 % of the zones respectively. Arola and Korkka-Niemi (2014) have concluded that especially for high density urban areas, the undisturbed groundwater temperature could be even 3−4°C higher than the average air temperature due to the heat island effect. Bayer, Attard, Blum, and Menberg (2019) coined the concept of subsurface urban heat islands (SUHI) and concluded that large-scale urban subsurface temperature could be 2−6°C higher than in the countryside. The research of Bonafoni, Baldinelli, and Verducci (2017) suggested the utilization of satellite remote sensing analysis of albedo and surface temperature changes in order to control the increase of SUHI effect in cities. Consequently, cities, at least in Nordic conditions where heat extraction from ground is dominant, are the perfect candidates for ATES development and integration. Multiple research works so far have been focused on ATES planning/monitoring in high density urban areas in terms of optimization of available subsurface space, flow/thermal interference and ATES overall efficiency (Bakr, van Oostrom, & Sommer, 2013;Bloemendal, Olsthoorn, & Boons, 2014;Bloemendal, Jaxa-Rozen, & Olsthoorn, 2018;Bozkaya, Li, Labeodan, Kramer, & Zeiler, 2017;Caljé, 2010;Fleuchaus, Schüppler, Godschalk, Bakema, & Blum, 2020;Hoving et al., 2014;Sommer, Valstar, Leusbrock, Grotenhuis, & Rijnaarts, 2015).
It has been concluded that realizing a sustainable ATES system requires long-term balancing of charging and discharging the storage. Additionally, the system's performance can be optimized by storing extra heat or cooling capacity from other (sustainable) sources e.g. by using solar thermal harvesting (Ghaebi, Bahadori, & Saidi, 2014;Kastner et al., 2017;Paksoy, Andersson, Abaci, Evliya, & Turgut, 2000). Optimization can be achieved also by controlling the energy exchange between warm or cold wells in individual ATES systems taking into account its impacts on the neighboring ATES systems or by modifying the energy demand itself (demand side management). Here, groundwater heat pumps (GWHP) take advantage of stable subsurface temperatures and the ability of groundwater areas to store the excess energy in order to use it when necessary (seasonal shifting of energy demand). Taking also into account that groundwater temperature is annually stable and, in Northern climate, several degrees above the average air temperature, heating mode operation of GWHP is much more efficient compared to conventional air-to-air heat pumps or fossil fuel boilers. Lund and Boyd (2016) revealed that geothermal energy penetration in Finland is high. However, it is still mostly driven by domestic and small-scale installations and practically inexistent in a large-scale/district level. Surprisingly, despite the availability of important aquifer areas and highly developed DH networks, Finland has still not implemented ATES systems in a large scale. Moreover, the availability of numerous groundwater resources in Finland is an important additional advantage for GWHP introduction and deployment, taking into account the variety of ATES projects implemented during the last decades in other Nordic countries like Sweden. For the above reasons, Finland has been chosen as the target country of the present study.
ATES operation results in a combined hydrological, thermal, chemical and microbiological impact on the affected groundwater areas and should be carefully evaluated (Bonte, Stuyfzand, Hulsmann, & van Beelen, 2011). The legislation of shallow geothermal installations (depth less than 400 m) is diverse among countries (Haehnlein, Bayer, & Blum, 2010). Regulations for installations of wells concern the use of hazardous materials and proper backfilling of the drilling hole to avoid hydraulic short circuiting between aquifers. Other legislation concerns protection of groundwater areas for drinking water supply. Some countries adopt limits for minimum and maximum storage temperatures, like Austria (5-20°C), Denmark (2-25°C) and the Netherlands (5-25°C) -while others adopt a maximum change in groundwater temperature, for example Switzerland (3°C) and France (11°C). In Finnish legislation, there is no explicit reference to groundwater use for energy generation and storage, the only generally related laws are the Water Act (1961) and the Environmental Protection Act (2000). Arola and Korkka-Niemi (2014) presented a theoretical study for ATES utilization in Finnish conditions, combining simulated energy demands for different building types with groundwater modeling, as well as performing separate simulations for heating/cooling loads with assumed constant heat pump COP, where finally a long-term environmental impact was assessed for each case. In Finland, there are publicly available data sources regarding the hydrological resources (Finnish Environmental Institute), geological conditions (Geological Survey of Finland) and geographical data (National Land Survey of Finland). The present work introduces a novel method for retrieving data from the aforementioned Finnish public sources and modeling the ATES system  using combined heating/cooling loads and variable COP-model for ATES integration in DH/DC networks. The present paper highlights also the importance of some calculated ATES parameters (thermal radius, heat recovery efficiency), which can be useful during pre-feasibility evaluation or in combination with more sophisticated modeling tools. The novelty of this paper is to investigate holistically the ATES application for heating and cooling, and although using limited and uncertain data, be able to prepare a sufficiently accurate feasibility study from three different perspectivestechnical rationale, economic feasibility and environmental impact. Furthermore, a simplified but holistic groundwater and techno-economic modeling approach of the energy system is introduced and developed. The present research will contribute to develop a valuable knowledge on ATES systems analysis, modeling and implementation. This work can potentially benefit engineers, energy companies and final users as well as organizations and agencies responsible for environmental legislation and energy efficiency.

Materials and methods
The modeling scheme of an ATES system is developed through the following general steps, namely -i) pre-processing of input data, ii) data processing using different tools & methods, and iii) post-processing with results presentation and analysis ( Fig. 1). In order to adequately model an ATES system, it is fundamental to develop and calibrate a specific groundwater model. Basically, both ATES and groundwater models should be connected and linked together by using different tools for data interaction. Finite difference method-based groundwater models (MODFLOW, Harbaugh, 2005) and solute/heat transport (MT3DMS, Zheng & Wang, 1999) are adopted and developed in the present work. Models are calibrated against statistical data (observation wells) of the studied groundwater areas. In Sections 2.2-2.4, the development of the modeling methodology is explained with details and demonstrated referring to the information and data collected from the Finnish case study (introduced in Section 2.1).

Characterization of the district heating network
In the design process, the first step is to characterize the target district heating network, which is in this case study located in the village of Pukkila, a Finnish municipality located in the Uusimaa region in the southern part of Finland, and its heat is mainly produced by wooden chips in a 1.5 MW nominal boiler. In addition to the base load chips boiler, the heating plant has two peak load oil boilers (2 MW), which have been used sporadically. The annual district heat generation in Pukkila was 4407 MW h in 2017 and a daily-based power demand is presented in Fig. 2 (Hynynen, 2018a, b). The efficiency of the base load boiler (chips) is estimated to be about 80 % (Hynynen, 2018a).

Characterization of groundwater areas
In the second step, available data is retrieved from open Finnish public sources in order to characterize the groundwater areas. Here, data from Finnish Environment Institute (Suomen Ympäristökeskus, SYKE) website is used (https://www.syke.fi/en), and particularly Hertta 5.7 application regarding groundwater areas, as well as monitoring stations and observation wells. Pukkila's groundwater area is composed of three different aquifer zones, of which two (Vanhalanmäki-161602 and Pukkilan kk-161601) are close to Pukkila village and its district heating plant. Porvoonjoki River is a natural border of the south-eastern part of the village and separates area 161,601 in two parts, being also a specified head boundary for the studied area (see Fig. 3).
There is available information for 5 observation wells in area 161,602 (# 105, 205, 305, 405 and 505) and 5 observation wells in area 161,601 (# 605, 705, 805, 905 and 1005). The water level variation of each well has been recorded during the last 10 years and the average values were used for the steady state model calibration.

Geographical data (elevation model)
The third step of the design process is to assess and model appropriately the groundwater flow as well as to establish correctly the boundary conditions, starting with reliable terrain model. The open data of the National Land Survey of Finland (https://tiedostopalvelu. maanmittauslaitos.fi/tp/kartta?lang=en), particularly its "10 m elevation model" is used for this purpose. The elevation model was downloaded as Geo-TIFF raster file in two separate cadastral sheets, L4111 and L4112, converted to Surfer Grid file (GRD) using QGIS (QGIS, 2019). The extracted area of roughly 5.2 km in length and 1.3 km in width contains a dynamic terrain with elevations between roughly 35 and 95 m a.s.l., as seen in Fig. 4. The Porvoonjoki River and its effluents Virenoja and Kuutinoja determine the hydrogeological contours and will play a key role as natural boundaries of the groundwater model.

ATES integration for DH and DC
Groundwater source heat pump (GWHP), operating with warm and cold well (well doublet) is considered. For this purpose, both #605 (as warm) and #705 (as cold) wells are adopted, since they are located close to the actual district heat production plant (see Fig. 5). According to the steady state head distribution and gradient (supported by measurements of observation wells and simulations), the dominant groundwater flow is from the north (higher head values) to the southeast reaching the Porvoonjoki River as lowest head boundary.  Hourly-based simulation results for annual heating and cooling demand of office building are used in order to introduce a dynamic variable cooling load (Tuominen, Holopainen, Eskola, Jokisalo, & Airaksinen, 2014). The data used is for a typical office building of 2695 m 2 net area and type D1 (according to Finnish building code, part C3-2010). The data was converted into an average daily-based heating and cooling loads and up-scaled appropriately (by a factor of 27) in order to match the existing Pukkila's DH network profile. GWHP is used to satisfy partially the heating demand (0.35 MW base load), using the existing chips boiler for peak loads. If boiler is needed, GWHP would be used first to increase DH network temperature from 40°C (assumed return temperature) to some intermediate value, and after that, the final   Todorov, et al. Sustainable Cities and Society 53 (2020) 101977 DH supply temperature would be reached by the boiler. This would improve HP efficiency in partial load mode, since COP H is better with lower HP supply temperature. Different and reversible operation during summer and winter periods is assumed, creating an ATES well doublet -warm well (#605) and cold well (#705). During the summer operation a primary ATES circuit starts from the cold abstraction well, providing district cooling. After the cooling exchanger, water at up to 14°C is utilized in GWHP evaporator, and finally injected into the warm well. During the winter period the process is reversed; water is taken from the warm well, conducted if needed through the district cooling network exchanger, used with GWHP, and finally injected into the cold well. ATES operation is simulated, calculating all relevant parameters on a daily-basis. The average pumping flow rate for each day is calculated as a maximum value between the flow needed for heating and the flow needed for cooling, using Eqs. (A.3) and (A.5). The injection temperature (after GWHP) is calculated according to Eq. (A.3) and solving for temperature drop ΔT. Finally, the electric energy demand for pumping is estimated from Eq. (A.4), assuming pumping pressure p = 600kPa and overall efficiency η estimated between 0.5 and 0.6 depending on pump's size, based on producers' data for submersible pumps Grundfos SP (Grundfos, 2019). ATES reversible operation during summer and winter period is presented in Fig. 6. It should be noted that main streams are not mixed within the central flow commutation during the winter operation.

COP H estimation model
In order to evaluate how both source temperature and heat pump O. Todorov, et al. Sustainable Cities and Society 53 (2020) 101977 (HP) supply temperature affect GWHP efficiency, a linear regression model based on GWHP producer's data presented in the work of Pero (2016) and Hynynen (2018b), is implemented. Available data is used for four source temperatures T 1 (0°, 10°, 20°and 30°) and for five HP supply temperatures T 2 (40°, 50°, 60°, 70°and 80°). In addition, data was linearly fitted in order to obtain analytical equations for each match line (see linear regression model A.10). This linear regression model is used in all simulations, interpolating for any value of source temperatures T 1 in order to estimate heat pump COP H as a function of T 1 and T 2 .

Groundwater flow model
Groundwater model is a simplified quantitative tool developed in order to synthesize and represent the actual hydrological processes, as well as to be able to describe and predict their future development and trends. In order to study both the steady state and transient behavior (due to additional stresses produced e.g. from pumping), the groundwater model was developed using the finite difference code MODFLOW-2005 (Harbaugh, 2005) under ModelMuse (ModelMuse, 2019) as graphical user interface. Alternatively, a simple steady state model was created with Microsoft Excel spreadsheet in order to validate the simulated results.

Boundary conditions
As mentioned previously, the studied aquifer area has natural boundaries with "specified head boundary" (Dirichlet condition). Porvoonjoki River limits the area from south-east, Virenoja stream from north and Kuutinoja stream crosses the area in the middle. Aquifer north-east and south-west limits would be considered as "no flow" boundaries, special cases of "specified flow boundary" (Neumann condition) applied with zero flow. Abstraction and injection wells would be represented as point sources and modeled with "specified flow boundary". Similarly, the recharge rate is represented by a "specified flow boundary" distributed evenly over the total model area.

Numerical models and steady state calibration
In MODFLOW, the aquifer area is discretized using 100 × 100 m square cell and grid of 40 columns by 13 rows, covering a physical area  O. Todorov, et al. Sustainable Cities and Society 53 (2020) 101977 of roughly 3 km 2 , comprised between the aquifer north-west border and the natural boundary, Porvoonjoki River, from the east (Fig. 7). The MODFLOW model is divided in 2 layers. Terrain's topography is introduced as "point average interpolation" for "Model_Top" parameter of the upper layer. The lower layer is a confined aquifer between elevations 20 and 40 m (see Fig. 8). A standard value for horizontal hydraulic conductivity is chosen (sand/gravel aquifer) K x = 1.10 −4 m/ s (Arola, Okkonen, & Jokisalo, 2016). Vertical hydraulic conductivity is set as K z = 0.1K x . Standard values are also used for porosity (n = 0.25) and storativity (S = 1.10 -5 ). Initially, a typical value for recharge rate of R = 5.10 -9 m/s is used (Arola et al., 2016), and after model's calibration, adjusted to R = 6.10 -9 m/s. Specified head boundaries are designated for Porvoonjoki River between 43.5 m (north) and 40.5 m (south), Virenoja stream between 48 m (east) and 48.8 m (west) and Kuutinoja from 46−47 m (east part) to 53 m (west part).
A simplified steady state model is implemented and calibrated using Microsoft Excel spreadsheet. For confined isotropic aquifer with recharge R and transmissivity T a two dimensional steady state groundwater model for head value h of cell (i,j) can be discretized according to Eq. (A.7). It is clear that steady state solution depends on R/T ratio, therefore in order to calibrate the model according to observation wells measured values, R/T ratio would be an important sensitive parameter to vary.
In Fig. 9, there is presented a steady state solution of a 35 × 9 cell Excel model, with heads in meters above sea level. Red border cells are specified head cells, whereas yellow cells represent a "no flow" boundary (outside model´s domain). All intermediate cells are set according to Eq. (A.7). Black border cells are the location of the observation wells (from left to right #105, 205 … 905).
Groundwater models need to be calibrated against the measured head values of the observation wells. For steady state model calibration would be used the estimated average heads between 2006 and 2016 (Anderson, Woessner, & Hunt, 2015). Trial-and-error matching was applied varying the sensitive parameter R/T between 2.10 −6 and 3.5.10 −6 using both Excel and MODFLOW groundwater models. Nearfield target heads values (observation wells 605, 705, 805 and 905 in area 161,601) would have higher matching ranking than far-field heads values (Anderson et al., 2015). Root mean squared error (Eq. (A.9)) would be used as more robust indicator than a simple average and has been computed for all residual well head values, but also separately for near-field (#161,601) and far-field (#161,602) areas. The results for the best near-field match (R = 6.10 -9 m/s / T = 2.10 -3 m 2 /s, values finally adopted in the model) are summarized in Table 1.

Solute and heat transport modeling
The solute transport model MT3DMS (Bedekar et al., 2016) is used to simulate the heat transport in shallow confined aquifers due to similarities between the mathematical formulation of solute and heat transport equations (Hecht-Méndez, Molina-Giraldo, Blum, & Bayer, 2010; Hecht-Méndez, Molina-Giraldo, Blum, & Bayer, 2010). The main parameters and correlation coefficients used in MT3DMS simulation are summarized in Table 2.

Model implementation in MODFLOW/MT3DMS
Both winter and summer periods are simulated in MODFLOW, introducing weekly-based stress periods (injection/abstraction flows and injection temperatures of warm/cold wells are weekly averaged) as well as making local grid refinement (LGR) near the wells to improve model accuracy. Initially, abstraction temperatures from warm and cold wells are assumed as 7.5°C and 6°C respectively, but since these values are arbitrarily chosen, a different and novel approach is developed for estimating iteratively a convergent solution.
LGR is adopted in the MODFLOW/MT3DMS model, where nearby areas to warm and cold wells are discretized with 50 × 50 m cell size (Fig. 10).
In Fig. 10, the green rectangle represents a 50 × 50 m LGR near the wells, red and blue squares represent warm and cold wells areas affected by a calculated thermal radius of roughly 75 m (Eq. (A.1)). They also are used to estimate the average abstraction temperature (during the abstraction periods). Observation point (pink) is located some 280 m downstream the cold well.

Techno-economic analysis of ATES operation
In principle, it would be beneficial to limit DH network supply temperature to a maximum value of 80°C (2018b) and reduce it according to the annual heating demand in order to limit the network losses and improve the overall efficiency. The minimum temperature of domestic hot water is recommended to be in the range 55−60°C (in order to prevent the risk for Legionella formation, according to Banks, 2012), and therefore, the minimum recommended DH supply temperature could be established in 70°C (Finnish Energy, 2013).
Daily-based Excel spreadsheet calculations were performed and based on them, it is possible to calculate the annual energy demands for heating, cooling and electricity as well as the average daily pumping flow rate. All relevant technical parameters are summarized in Table 3.
Cost database for different HP generation technologies is used, according to Nielsen and Möller (2013), as well as prices for piping, heat exchangers and ATES well drilling (Drenkelfort, Kieseler, Pasemann, & Behrendt, 2015) in order to determine the investment cost. The cost of produced energy is calculated based on the annuity method, assigning annual payments for the investment and taking into account the investment's lifetime of 20 years (Nielsen & Möller, 2013) as well as adopting a value of 5 % for interest rate. The annuity also includes O&M cost (assumed as 1 % of investment) and electricity cost for GWHP and ATES pumping (assumed as 100€/MWh, including taxes, transfer/distribution fees, Nordpool, 2019). The economic analysis includes the estimation of the following parameters summarized in Table 4.

Initial ATES model settings
The annual variation of all relevant parameters is presented in Fig. 11. The annual combined demand are presented on the left axis respectively as red curve (heating) and cyan curve (cooling). Heat pump COP varies from 2.6 to 4.2 (3.4 on the average). The average injection temperature is 2.8°C/7.9°C (cold/warm well) and the average pumping flow rate is 0.015 m 3 /s. Summer period lasts from week 18 to week 36 (both inclusive), abstracting from cold well (#705) and injecting into the warm well (#607) respectively. Inversely, winter period comprises weeks 37-52 and 1-17, and the system operates abstracting from warm well and injecting into cold well respectively.

ATES model iterations
Additionally, an 8-year period simulation was performed with the previously developed MODFLOW model in order to study the abstraction temperature variation of the warm and the cold well, as well as to estimate the charged and discharged thermal energy. Warm/cold abstraction temperatures are calculated as an average of warm/cold well area (defined according to Fig. 10), only during the well's abstraction period (summer for the cold well and winter for the warm well). Since heating demand is dominant, charge and discharge annual cycle calculations are performed for the cold well and results are summarized in Table 5. It can be noted that after roughly 3-4 years of operation, the abstraction temperatures and heat recovery factor (HRF) converge.

Table 1
Steady state model calibration with the ratio of R/T is 3.10 −6 .
In the last column are highlighted the best RMSE values for MODFLOW far-field, near-field and overall results respectively.      Todorov, et al. Sustainable Cities and Society 53 (2020) 101977 From the weekly-based calculations it has also been observed that ATES summer operation should be reversed already in week 34 instead of 37. This is implemented in the next model iterations. Finally, after four iterations, warm/cold abstraction temperatures converge with good accuracy to 7.2-7.3°C/3.9-4.2°C, as shown in Table 6.

Techno-economic analysis
The fourth ATES model iteration is used for the analysis (warm/cold abstraction temperatures 7.3/3.9°C). The main technical parameters of the system are summarized in Table 7. It should be noticed that even with 11 % of peak heat power, the GWHP coverage ratio is roughly 38 % of the annual heating demand. The average COP is 3.4 slightly increased to 3.6 during the winter period due to the lower GWHP supply temperature. Since heating is mostly generated during the winter period roughly 2/3 of GWHP power consumption is in winter. As the amount of pumped water is roughly balanced between winter and summer operation, the calculated thermal radius of the warm/cold wells are approximately equal (75 m). On the other hand, the excessive maximum calculated drawdown of the cold well (over 14 m) due to more intensive pumping (in only 16 weeks) during the summer period should be noted.
The economic feasibility estimation for investment costs as well as the energy production cost are presented in Tables 8 and 9 respectively. Investment costs do not include the necessary DC network implementation. The resulted overall energy production cost is roughly 41.5 €/MWh. The total investment cost is roughly 1.06 million € from which 22 % correspond to HP/exchangers and 75 % is related to the underground part (wells, pipes). These figures are in line with the research of Schüppler, Fleuchaus, and Blum (2019), where similar ATES system in Germany resulted in total investment cost of roughly 1.28 million € from which 23 % correspond to HP/exchangers and 60 % is related to wells/piping.

Impact on groundwater areas
ATES model iteration #4 is adopted in order to study the long-term effect (20 years) of warm and cold wells´thermal interaction. The thermal front simulation of ATES operation for year 1, 2, 4, 8 and 20 is presented in the following Fig. 12, for the week when the cold and warm plumes achieve their maximum expansion.
In Fig. 12, both wells are represented by a 50 × 50 m pink cell. Left images depict the maximum annual cold well plume expansion (end of the winter period, after week 11), while right images are the maximum annual warm well plume expansion (end of the summer period, after week 31).
It can be observed that the thermal plume of the warm well maintains more or less within its thermal radius of roughly 75 m, since heat injected in the aquifer is less compared to the heat abstracted from it. Moreover, the heat plume around the warm well almost vanishes at the end of the winter period (left images), while the plume around the cold well increases slowly over the years. After 20 years of ATES operation, the thermal plume of the cold well expands several hundreds of meters to the south-east, following the dominant groundwater flow direction. All in all, it can be concluded that the locations of the wells (cold well located downstream) and the separation between them is favorable for a correct long-term ATES operation.
The temperature evolution of the warm and cold well areas over the 20-year period of ATES operation (including the observation point temperature), as well as the pumping flow rate are presented in Fig. 13. It can be acknowledged that after the third year the ATES system converges with cyclically varying temperatures: warm well 279−281 K, cold well 276−279 K. Weekly averaged pumping rates are varying between roughly 0.01 m 3 /s (winter period) and −0.025 m 3 /s (summer period) with some peaks when heating or cooling loads are exceptionally high. Cold plume reaches the downstream observation point after roughly 8 years of operation and after 20 years its effect is slowly attenuated to roughly −0.8°C compared to aquifer's undisturbed temperature.   Todorov, et al. Sustainable Cities and Society 53 (2020) 101977 4. Discussion All analyses presented in this work were carried out using quite limited and uncertain information regarding the hydrogeology of the studied case area. The fundamental assumptions -normally simplified during the prefeasibility phase -were that the aquifer layer is uniform, confined and isotropic in the considered area. For more accurate aquifer parameters estimation, additional geological survey and tests for non-equilibrium (transient) flow conditions should be conducted as next steps for model calibration and validation.
The presented case study for ATES integration within DH/DC networks was economically feasible and had limited environmental impact even within a 20-year horizon of ATES operation. Combined heating and cooling demands, with seasonally reversible ATES operation and roughly balanced pumping volumes during summer and winter periods, had low impact on the aquifer area and seemed to be environmentally equilibrated.
The introduction of combined heating and cooling demand (in addition to the existing DH demand) using simulated data for typical Finnish office building was coupled with the ATES model. It showed promising economic results (energy production cost 41.5 €/MWh, compared to the average DH price in Finland of 76.7 €/MWh according to DH, 2017) and low long-term environmental impact on the surrounding groundwater areas, less than 1°C after 20 years of operation in a radius of 280 m from the pumping wells, which is below the limits set by Swiss (3°C) and French (11°C) legislation. The average cold and warm storage temperatures, computed for a 20-year period of ATES operation, are 4.1°C and 7.1°C respectively, fulfilling the Danish limits of 2−25°C. Thus, it can be concluded that the presented ATES integration for combined heating and cooling proves to be a feasible solution as well as environmentally equilibrated, fulfilling most European existing regulations.
Another important parameter is the imbalance ratio (IR), defined by Bozkaya et al. as a ratio of the difference between the heat injected and extracted from the ground and the maximum of these values. In the aforementioned research was highlighted the influence of thermal imbalance in the ground temperature change presenting variations between −10°C and +10°C for heating and cooling dominated systems respectively (Bozkaya, Li, & Zeiler, 2018). The current work is a case of heating dominated system with IR = −27 % and ground temperature change of roughly −1°C after 20 years of ATES operation.
The research of this work could be continued and the result accuracy could be improved for example by the following steps: • additional geological survey and slug & pumping tests in order to improve groundwater model quality • study of some additional energy sources integration, such as thermosolar and industrial waste heat • more detailed study on how ATES system is affected by energy prices • more efficient ATES optimization and control strategies

Conclusions
The present work was successful in demonstrating that using limited and uncertain data (mostly obtained from public open-data sources), it is possible to holistically integrate groundwater models, energy demand and ATES system as well as to study system's performance efficiency, techno-economic feasibility and the impact of ATES operation on the surrounding groundwater areas. Groundwater flow and thermal models were developed and calibrated, using a variety of available data sources (National Land Survey of Finland, Finnish Environment Institute) and tools (MS Excel, QGIS, MODFLOW, MT3DMS). Heat pump COP estimation analytical model was also implemented and coupled with the groundwater models. The purpose was to develop a case study for ATES integration within the existing Pukkila's district heating network, as well as to assess the long-term environmental flow and thermal impact generated to aquifer groundwater areas.
All in all, ATES systems prove to be an efficient and a sustainable alternative for traditional fossil fuel boilers, due to their capacity to annually store and recover cooling & heating energy from the subsurface. Significant technical and economic improvement could be achieved when simultaneous or seasonable cooling and heating loads are dispatched, within integrated district energy networks.

Declaration of Competing Interest
None.

Acknowledgments
This work is part of GESATES project, an innovative initiative for groundwater heat pump application using aquifer thermal energy storage (ATES), funded by Business Finland and promoted by different participating companies and organizations, such as HELEN, Turku Energy, Nastola Foundation, Nivos Energy. It also included different research locations in Finland -Helsinki, Turku, Nastola and Pukkila-Mäntsälä region.  • Heat recovery factor (HRF) for heating/cooling stored energy (reversible ATES operation) is defined as a ratio between the annually discharged and charged energy (Kranz and Bartels, 2010). , Bloemendal and Hartog (2018) also utilize a term of recovery efficiency, referred to a warm or cold stored energy over the whole charge/discharge cycle (normally one year): where indices out/in stand for discharge/charge cycles, warm/cold stand for warm/cold wells and E is the stored/recovered aquifer energy.
• Pumping flow rate, where Φ i is the heat demand for a day i (condenser side), ΔT is water's temperature drop in heat pump evaporator side and COP H is heat pump COP: • Electric pumping energy consumption P i : • Pumping flow rate, where Φ is the heat demand (evaporator side) and ΔT is water's temperature drop in heat pump evaporator side: • Annuity factor AF, where r is the interest rate and k is the number of years • Real head in cell (i,j) with simulated head h i,j containing the pumping well (with radius r w ) can be computed additionally, applying Thiem equation ( O. Todorov, et al. Sustainable Cities and Society 53 (2020)