Development and evaluation of a building integrated aquifer thermal storage model

(cid:1) The previous aquifer thermal energy storage (ATES) models are investigated. (cid:1) An ATES model was developed using the co-simulation between MATLAB and COMSOL. (cid:1) The developed ATES model was validated with an experimental data from an operational ATES in the Netherlands. An aquifer thermal energy storage (ATES) in combination with a heat pump is an excellent way to reduce the net energy usage of buildings. The use of ATES has been demonstrated to have the potential to provide a reduction of between 20 and 40% in the cooling and heating energy use of buildings. ATES systems are however a complex system to analyse as a number of ground conditions inﬂuence heat losses within the aquifer. ATES is also not conﬁned from the sides and is therefore vulnerable to heat losses through con- duction, advection and dispersion. The analyses of ATES system is even further complicated when the dynamic of a building is considered. When connected to a building, the temperature in the aquifer is inﬂuenced by the amount of heat exchange with the varying building load. Given the energy saving potentials of ATES systems in building operation, detailed understanding of the inﬂuence of buildings on the ATES systems and vice versa would facilitate improved operation and efﬁciency of ATES and building coupled systems. Therefore, taking into account the variations in the building and below ground con- ditions, there is the need for the development of a model that can potentially handle the dynamics on both sides. Finite element and ﬁnite volume methods are frequently used in the development of ATES models and proven as adequate tools for modelling complex ground conditions, however, most developed ATES models are often analysed independent of the building. Therefore, in this study, an ATES model that also integrates building dynamics is developed using the ﬁnite element method (FEM). The developed model was validated using data from an ATES and building in the Netherlands. The developed model was shown to have an absolute mean error of 0.17 (cid:1) C and 0.12 (cid:1) C for the cold and warm wells respectively. (cid:3) 2017 The Authors. Published by Elsevier Ltd. ThisisanopenaccessarticleundertheCCBYlicense(http://


a b s t r a c t
An aquifer thermal energy storage (ATES) in combination with a heat pump is an excellent way to reduce the net energy usage of buildings. The use of ATES has been demonstrated to have the potential to provide a reduction of between 20 and 40% in the cooling and heating energy use of buildings. ATES systems are however a complex system to analyse as a number of ground conditions influence heat losses within the aquifer. ATES is also not confined from the sides and is therefore vulnerable to heat losses through conduction, advection and dispersion. The analyses of ATES system is even further complicated when the dynamic of a building is considered. When connected to a building, the temperature in the aquifer is influenced by the amount of heat exchange with the varying building load. Given the energy saving potentials of ATES systems in building operation, detailed understanding of the influence of buildings on the ATES systems and vice versa would facilitate improved operation and efficiency of ATES and building coupled systems. Therefore, taking into account the variations in the building and below ground conditions, there is the need for the development of a model that can potentially handle the dynamics on both sides. Finite element and finite volume methods are frequently used in the development of ATES models and proven as adequate tools for modelling complex ground conditions, however, most developed ATES models are often analysed independent of the building. Therefore, in this study, an ATES model that also integrates building dynamics is developed using the finite element method (FEM). The developed model was validated using data from an ATES and building in the Netherlands. The developed model was shown to have an absolute mean error of 0.17°C and 0.12°C for the cold and warm wells respectively.

Introduction
In an attempt to decrease CO2 emissions in heating and cooling applications, the ground has been introduced as an efficient cooling/heating option for buildings. One of those applications: the 'aquifer thermal energy storage system (ATES)' uses energy stored in readily available underground water to exchange heat or cold with a building. ATES is also widely recognized as one of the most energy efficient heat/cool options [1]. Compared to traditional cooling/heating systems, ATESs can achieve 90-95% energy saving on heating/cooling when used directly, and 60-85% energy saving when coupled to a heat pump [2].
When used in buildings, the temperature of ATES varies throughout the year due to the influence of heat losses to the surroundings and the amount of injection/extraction of heat/cold to/ from the ground. Since buildings have varying heating and cooling demand patterns all year round, it influences the amount of heat/cold injection to the ground. Given that the rate and amount of heat and cold water injection and extraction to the well also influence heat loss and gain in the well, it is, therefore, imperative that this variation is considered in the evaluation of ATES connected with buildings.
A few studies [3][4][5] have been conducted to investigate the performance of ATES in connection with a building load. The authors in [3] for example investigated the performance of ATES connected to solar thermal and heat pump, while the authors in [4] analysed the influence of various operational settings in a building on ATES performance. On the other hand, the authors in [5] performed a simulation for high-temperature ATES connected to combined heat and power (CHP) system and heat pumps to assess the energy performance on the district level. These studies have shown the relationship between the building load with ATES. In these studies, the ATES model was developed using TRANSAT (ATES model in TRNSYS) [3,5,6] and C code [4] which is based on the finite difference method (FDM) that has limited capability to deal with the complex and dynamic interactions between the various parameters that influence ATES performance, in an attempt to have lightweight computational load.
Most other studies using more advanced simulation technique based on the finite element method (FEM) and finite volume method (FVM) [7][8][9][10][11][12][13][14][15][16][17] have mainly evaluated the performance of ATES independent of the connected load in the buildings. Thermal models of ATES have been discussed and analysed under various subtopics including buoyancy [7][8][9], preferential pathways [10], dispersion [11][12][13], thermal interaction [14][15][16] and natural groundwater flow [17]. These subtopics can be discussed extensively depending on the potential influence on heat recovery and the temperature distribution around the well. In general, most authors [6][7][8][9][10][11][12][13][14][15][16][17] argue that the ground is a complex medium to model since it is influenced by various parameters. For instance, buoyancy is an important factor that influences the temperature distribution for high-temperature ATES [7][8][9]. Dispersion is not negligible in heterogeneous ground layers and high-velocity water regions [18]. Thermal interaction is a concern for the places where the ATESs are densely installed and heavily affected by the operational parameters as well as certain ground conditions such as hydraulic conductivity [15]. In most of the developed FEM/FVM ATES models, the injected/extracted heat introduced is often considered constant for a certain period of time (usually days of the year), which do not dynamically change depending on building load, to analyse the influence of certain ground parameters. This results in discrepancies between the reality and the simulation since in reality operational parameters are sensitive to (halfhourly/hourly) changes in the building load or other dynamics within whole HVAC system. In addition, these models are mostly developed using strong computational fluid dynamics software such as SHEMAT, FEFLOW, MT3DMS, MODFLOW [19][20][21], which have limited capability to integrate building simulation tools (such as Energy-plus, TRNSYS, MATLAB, Modelica).
The subsequent sections of this paper are outlined as follows: Section 2 provides a literature survey on the previous models of ATES using FEM and FVM method, Section 3 provides a description of the simulation tools and the ATES model. Section 4 provides a discussion of the results and validation of the ATES model, while Section 5 provides the conclusions.

Previous modelling studies
As noted in the preceding section, models that integrate into building with ATES have been limited to those based on FDM. Most studies on ATES have often focused on below ground conditions [7][8][9][10][11][12][13][14][15][16][17] and the very few studies [3][4][5][6] that combine both above and below ground conditions are often based on the inadequate FDM. Therefore, in this section, previous modelings of ATES were collected to determine what the FEM and FVM based modelings could potentially solve in the past studies in order to address the possible challenges that a user can face with FDM method. Most of the studies presented in Table 1 dealt with various heat transfer phenomena such as dispersion, thermal interference, buoyancy and natural ground water flow to have a better prediction of the expansion of the thermal front and the temperature distribution around the well. Although the expansion of thermal front differs depending on the capabilities of the ground, the rapid increase in the utilization of the ground is expected to contribute to the scarcity of ground space so that the potential use of aquifers in the area of interest can be maximized. The proper definition of heat transfer function enables the optimal layout as well as a well-defined model for the integration into building load [22].
Previous studies shown in Table 1, depending on the applied location, the mentioned heat transfer functions were significant. Dispersion effect was often studied with heterogeneity (layering) in the ground, where the water was not uniformly distributed due to the differences in hydraulic conductivity of the layers. Dispersion results in a decrease in heat recovery as well as significant uncertainty in the thermal plume expansion [23,24]. Visser et al. [25] has predicted a maximum of 80% plume extent in comparison to homogenous porous media, which has increased the risk for thermal interaction between well groups. In the absence of natural groundwater flow, the effect of dispersion on the heat recovery is considerably low (around 0 and 3%). The presence of groundwater flow of 100 (m/year) has intensified the effect of dispersion to 0-12% since heat losses increase through advection [17].
Thermal interaction between well groups has potentially decreased the heat recovery by 20% or improve by 25% depending on the layout of the well position [14]. The risk for thermal interaction changes depending on the parameters such as the distance between well groups, injection/extraction rates, hydraulic conductivity and natural ground water flow [15,16]. In case the ATES is bounded by the sides with a concrete wall, concrete can act as good insulator and improves the heat recovery [8].
Natural ground water speed is a very significant parameter that can affect the heat recovery [3,10]. High natural velocity can make the high-temperature ATES impracticable due to the migration of high amount of heat. Chevalier and Banton [10] has predicted around 37% and a 78% decrease in heat recovery for the natural ground water velocity of 5 Â 10 7 ðm=sÞ and 10 À6 ðm=sÞ respectively. The heat recovery is influenced not only by the intensity but also the direction of the natural ground water flow [17]. Gao et al. [12] has shown the layout of the well arrangement poses different risk levels for thermal interaction.
Buoyancy has been a significant parameter for the hightemperature ATESs. Liu et al. [7] defined the critical limit temperature to take into account buoyancy effect. Nagano et al. [27] has shown that buoyancy has less influence on temperature distribution under high-velocity regions. High-temperature ATES also influences the hydraulic conductivity of the aquifer. An increase in water temperature from 16.2°C to 50°C led an increase of hydraulic conductivity from 1:97 Â 10 À6 to 2:82 Â 10 À6 ðm 2 =s PaÞ [27].
As summarized from the literature, the temperature in ATES is influenced by various transfer functions (Table 1), domain [23,24], boundary [8,15,16] and transient [27] conditions, which pose challenges for the use of FDM method. It has been shown that FEM and FVM methods were popularly used deal with the complexity in the past studies (Table 1). However, the ATES models integrated into building was limited with inadequate FDM [3][4][5][6], where certain conditions such as dispersion, buoyancy, thermal interaction and natural ground water flow, which influence ATES performances are often ignored.

COMSOL for ATES modelling
The ATES was modeled using COMSOL subsurface and heat transfer module [32]. Heat losses are influenced by two operational entities: temperature and flow in the same medium; which results in mass and heat transfer in the porous media. Convective heat transfer occurs through transportation of particles via advection due to the hydraulic head difference during the charging/discharging period. In addition to the operational parameters, there are several physical ground parameters that affect the heat transfer such as porosity, hydraulic conductivity, the thickness of aquifer and some other known parameters such as heat capacity, density and thermal conductivity.
Porosity determines the thermal capacity of a volume of the aquifer. As the porosity increases, the thermal capacity of the aquifer increases since water is able to hold more heat in comparison to the ground. It allows the ATES to keep the thermal energy close to the injection/extraction point [3]. Porosity determines the specific heat capacity from: where C p aq is the specific heat capacity of the aquifer, £ is the porosity, C pwater is the specific heat capacity of water, C psand is the specific heat capacity of sand. Hydraulic conductivity is used to determine the resistance of porous media to the movement of water. Hydraulic conductivity should be high enough to yield a certain amount of water flow to neutralize the net head increase, whereas it should be low enough to ensure low thermal migration from the injection point. High hydraulic conductivity increase the risk for thermal interaction between well groups [15]. The thickness of aquifer determines the thermal capacity of a confined ATES. Thermal energy injection should be high enough to satisfy the capacity since the well is going to experience slow changes in temperature. The authors in [3] have shown that increase in the thickness lowers the losses through conduction to the surroundings by means of low area/volume ratio.

Governing equations
Mass transfer in an aquifer is represented by forced and natural convection. Forced convection is the transport of mass in the presence of external forces (injection/extraction). Natural convection is the transport of solute through dynamic nature of fluid which can be driven by the buoyancy, diffusion or advection. Following equation satisfies the flow of water in porous media. The flow flux of water is written as a function of head gradient and coefficient of hydraulic conductivity.
whereq is Darcy flow, K is hydraulic conductivity (m/s), h is hydraulic head (m). Transient drawdown of the injection/extraction is presented with the following equation.
where S s is the storage coefficient, q is the specific heat capacity, Q s is the source/sink term. The average extraction flow rate of cold well and warm well is calculated as 40 and 18 ðm 3 =h), respectively. System does not cause significant drawdown (less than a meter) (see in Fig. 1) to take into account for low temperature ATES [3,33]. The drawdown is calculated based on the parameters in Table 2 and using analytical solution (Eq. (4)).
As a result, storage coefficient is considered as zero, which means that there is no transient drawdown of water by time; thus, the pressure distribution around the well is considered steady. Analytic translation of injection/extraction rate into hydraulic head change is shown as follows [3]: where A is the thickness, r a is actual point, r f is reference point. As it is seen from the equation increase in the thickness lower the hydraulic change and therefore, the mixing rate within the aquifer. Mass transfer in a porous media is defined by Darcy law. For where k aq is the thermal conductivity of aquifer, ðpcÞ f specific heat capacity of fluid, ðpcÞ aq specific heat capacity of the aquifer.

Meshing
The model consists of a total of 16 boundaries and 5 domains that are numerically solved by the finite element method (Fig. 2). The computational domain was quadratically meshed, and the complete mesh consists of 84,933 elements domain elements and 666 boundary elements. The minimum element quality was 0.512. The simulation time was set to 1 h. The aquifer is confined by two 5 m wide impermeable layers of clay.

Boundary conditions
Source and sink terms were introduced on the same boundary and vary time dependently. Two wells were simulated in two different domains with the same conditions. The active boundary was defined as source term on charging period, and, sink term on extraction period. Pressure set for the lateral, upper and lower sides were set to zero. The following equation (7) was used to convert the flow rate into the hydraulic head difference, on the active boundary [34]: where Q is the extraction/injection from/to aquifer ðm 3 =sÞ, B is the transmissivity ðm 2 =sÞ, r is the radius of domain ðmÞ, r 0 is radius of borehole ðmÞ. Temperature boundaries were chosen as symmetry boundary condition which means that there is no heat transfer beyond the domain. The mathematical representation of boundary conditions are shown as follows:

Assumptions
The factors that have less effect on heat recovery on ATES can be neglected for the integration into the HVAC system model. Besides, Dutch ground conditions for the case study are quite suitable to simplify the model. By referring the previous models [3,17], the factors that have less effect on heat recovery can be eliminated. The following assumptions were implemented in the model: 1. The computational domain is homogenous and isotropic. Therefore, the flow of water in each dimension was considered same.  2. ATES is confined by two impermeable layers; therefore, the vertical infiltration of water is neglected. 3. There is an only radial movement of water. (Dupuis approach) 4. Since there is no vertical infiltration, thermal interaction is neglected due to the separation of wells with a clay layer. 5. Although natural ground water flow has a considerable effect on the temperature, natural ground water flow is neglected in this study since it is very low (few meters in a year) in the Netherlands. Thus, the ATES model could be simplified into 2D Axisymmetric model [3]. 6. Dispersion and buoyancy effects are neglected, dispersion effect is very limited in cases with no natural ground water flow (0-5% decrease in heat recovery) [17].

Sizing the ATES
Sommer et al. [22] defined the maximum distance that thermal front can be away from the injection well in a homogenous ground conditions and neglecting the vertical infiltration, natural ground water flow, dispersion and thermal conduction (Eq. (9)).
where V is the volumetric injection rate, c w is the heat capacity of water, c aq is the heat capacity of the aquifer, b is the thickness of the aquifer. The Dutch society reported that the triple of the defined R th is enough to avoid thermal interaction, which can be translated into no more heat transfer beyond that radius. Since no heat flux boundary condition is defined at the end of the domain, the domain size based on the triple of (R th Þ. Therefore, First, the total injected volume of water was calculated for the projected 4 year data to determined (V storage Þ, and Eq. (9) was applied to size the thermal radius (|AB|) of the domain. Taking into account the maximum total injected volume in one season which is 176,000 ðm 3 Þ for warm well and 57,000 ðm 3 Þ for cold well (Fig. 3), the total calculated radius is 120 (m) for both of the wells. In order to confirm the size of domain, the maximum amount of volume in one season applied for both wells and it was observed that the size of the domain was sufficient to apply no  heat transfer boundary at the end of the domain. The water is injected with a fixed temperature of 16°C for warm well and 8°C for cold well (Figs. 4 and 5).

MATLAB link COMSOL
There is a developed tool for the communication between MATLAB and COMSOL. Using MATLAB Live-Link tool, COMSOL models can be controlled through MATLAB and COMSOL performs in the batch mode. All globally predefined variables, parameter, geometries, meshes and boundaries on COMSOL can be controlled on MATLAB (Fig. 6). For the time-dependent solutions, it is also possible to control the time-steps for each iteration, which allows the user to adapt various time steps depending on the level of changes in the particular period of time. Adapting time-steps is important since calculation time is heavily influenced by the amount cycles during the simulation, maybe, even more than the complexity of the model. To give an example, the simulation can solve 100 time steps of simulation in 7 s, while, it takes approximately 500 s to solve 100 steps in the 100 cycle. For instance, for this study, there are 3 different modes of operation for warm and cold well. While the cycles are applied for each time step during injection and extraction period, the cycles can be decreased for the resting mode, where there is only the heat is conducted to the surroundings. There is no difference between calculating the conduction equation, for instance, in 1 iteration for the period of 10 h or 1 hourly 10 iterations; however, the calculation time can be improved significantly.
The boundary connected to the building is activated based on the mode of operation, where the hydraulic head is defined as negative during extraction, positive during injection and zero during the resting period. There are two physics coupled in the domain which is heat transfer in porous media module and subsurface flow module. While subsurface flow module is solving pressure distribution within the domain, heat transfer module is calculating temperature distribution depending on the pressure distribution. Since it is not possible to specify the temperature directly on the boundary, heat transfer boundary is defined as heat flux depending on the pressure information.

ATES model validation
The validation study was conducted using the hourly data from an office building located in Amsterdam, the Netherlands using 4 years of data. The simulation was conducted using MATLAB in connection with COMSOL. The system is a doublet, using two wells of approximately 50 and 100 m depth for cold and warm well respectively. The ground layer information as the ground is drilled are stored on an online database [35]. From the stored information, the aquifer containing layer can be determined to be sand with a thickness of 20 m and 60 m, which are separated by the clay layer [35]. The parameter estimations for the ground layers is obtained from [32,14] ( Table 2). The temperature sensors are located on the extraction/injection point of the ATES which was filtered based on the signal coming from the active pump. In addition, as mentioned earlier (see Fig. 3), the minimal flow rate of cold extraction and warm extraction are 20 and 6 ðm 3 =hÞ, respectively; therefore, the temperature values are not taken into account for the lower values of minimal flow rate. There are moments where the extraction temperature experiences sudden peaks. These are the moments when the data is received right after the pump is activated. The sensor basically measures the temperature of the conducted heat during the resting mode, not the temperature of the ground water, which is quite visible by looking at the cold well temperature. There are many points where the water temperature go above natural ground temperature of 11.5°C, which is supposed to be maximum temperature for cold well (Fig. 7). In addition, the trend of the extraction temperature for cold and warm well is stable during the extraction period. The temperature data logger is Testo 176H2, which has an accuracy of ±0.2°K (between À20 and 70°C). The measurements are executed with a calibrated 'TA-Hydronics Scope' using pressure differential measurements over balancing valves. TA-Hydronics Scope has an accuracy of 0.1°K and the flow accuracy of ±0.1 (l=s). The sensitivity of the temperature sensor is 0.1°C. The temperature of the top layer soil has the yearly average ambient temperature, which is 10.4°C in Amsterdam The ground temperature increases by 2°C for every 100 m [31]. Therefore, ground temperature was taken as 11.5°C and 12.5°C as an initial condition for cold well and warm well domain respectively. The variable speed pump connected to the ATES was working in 6 steps ( Table 2) depending on the temperature output of the heat exchanger.
In Figs. 8 and 9, the black line shows the temperature change on the injection/extraction point. The extraction temperature is  significant to predict for the building integration; therefore, the intention was to match the extraction temperature (red 1 and blue line) with simulated temperature (black line). During the injection period, time-dependent values of flow rates and temperatures were defined on the boundary and, extraction temperature was predicted based on the extracted amount during the extraction period. The yearly mean absolute temperature deviation is as high as 0.8°C and 0.5°C for warm and cold well respectively on the first extraction temperature (Fig. 10) in the first season, the deviation for both cold and warm wells is higher due to uncertainties in the previous years. The fact is that ATES has been operational 4 years prior to that year. Therefore, the left over heat in the well is uncertain prior to the simulated dates. The measured temperature was available after the 4th year of operation which means that the measured temperature was expected to have a low-temperature gradient in comparison to the simulated temperature in the first couple of years due to the change in heat recovery by years. It is known that the heat recovery increases year by year [11,14] and becomes relatively steady in the 4th year, as a result, simulation results experience higher temperature gradient in comparison to the experimental data. The trend of temperature matches closely in the 3rd and 4th year of simulation. Specifically, looking at the final year, where the heat recovery rates close to each other and the initial discharge temperature is close to each other, The yearly mean temperature deviation in the final year is as low as 0.17°C and 0.12°C for cold well and warm well (Fig. 11), respectively. In the second year, the warm well temperature sensor on the discharge pump malfunctions at the beginning of the discharge, therefore; it deviates more at the beginning of the 2nd extraction period (Fig. 10).

Discussion
The integration of ATES into building applications has been found limited based on the review of previously studied ATES systems. On the one hand, there are various conditions of the ground, on the other hand, there are dynamics of the building that can influence the temperature distribution in ATES. Frequently, ATES models were decoupled from the building, the effect of ground parameters on the heat recovery was investigated based on the fixed boundary conditions of heat flux into the domain. Therefore, there was no model validated with the experimental data.
As it is proven in the review of the models [3][4][5][6], the heat recovery of ATES can be very site-specific. None of them [3][4][5][6] defined any different physics or domain conditions, in an attempt to make the model as simple as possible. TRNATES model which is frequently used is limited with solving the dynamics of the underground. For instance, Drenkelfort et al. [6] applied the TRNATES for high-temperature ATES, where the buoyancy effect was neglected. The changes in thermal properties such as density, hydraulic conductivity were not considered. The authors in [3,4] have taken into account the natural ground water flow; however, have not considered or evaluated the possible dispersion effect. All ATES systems were considered as individuals and completely isolated from the other ATES systems. However, in reality, ATESs can be situated in an area surrounded by a group of ATESs, where the system performance can be significantly influenced by thermal interaction. In that case, adding dispersion and thermal interaction effect to the physics of numerical solution can be necessary due to the increase of thermal plume bigger than expected. This is where a finite element method become functional to define new dynamics, domains and boundary conditions to the model. The drawback of this model is the computational load. Although the presented model is relatively simple in terms of a number of elements, the calculation period takes more time in cycles. For instance, the simulation can solve 100-time step of simulation in 7 s at once, while, it takes approximately 500 s in the cycle mode (1 cycle for each time-step).
A validation study has proven that the developed model using MATLAB-COMSOL can predict the behavior of ATES accurately and capable of handling hourly changes in the injection/extraction heat rates. More importantly, since the model is controlled from MATLAB, it can be easily introduced to the models developed in MATLAB or coupled to the co-simulation with building energy software (such as TRNSYS, Modelica, Energy plus) [36][37][38]. The studied softwares (such as SHEMAT, FEFLOW, MT3DMS, MODFLOW) [19][20][21] in previous models (presented in Table 1) have a limitation in integrating into dynamic building simulation tools.

Conclusion
This paper introduced the existing models based on the ground characteristics for low and high-temperature ATES to summarize the influence of various heat transfer function on the heat recovery. The results with these models confirmed the need for a model that is adaptable into various ground conditions and dynamics of the buildings. The co-simulation between MATLAB and COMSOL was introduced to simulate FEM model based on the varying flow rate and temperature information from the building side. In order to prove adaptability and reliability of the FEM model, a validation study performed using MATLAB in connection with COMSOL. Performed validation study has shown its potential for handling the dynamics of the building. The yearly absolute mean deviation between the simulated and the measured values for the well extraction temperature was as low as 0.17°C and 0.12°C for cold well and warm well, respectively.