Numerical modeling and validation of a large-scale borehole thermal energy storage system in Finland

.


Introduction
The energy challenge and global warming trend have exerted pressure on reducing carbon dioxide emission.Currently, European Union countries undertakes a task of achieving climate neutrality by 2050 (European Commission, 2020).In this context, borehole thermal energy storage (BTES) system is increasingly used as a popular renewable energy technique in connection with ground source heat pump (GSHP) for heating and cooling.Modelling of BTES have been investigated for years in BTES design, in-situ ground thermal response test data analysis, and integrated building energy system simulation (Rees and He, 2013).Various analytical and numerical borehole heat exchanger models with different complexities have been developed and used for design and research.In order to verify the reliability of these models, experimental monitoring data are often used to validate the simulation results.For example, Cai et al. (2021) used measured inlet and outlet brine temperatures of deep borehole heat exchangers to validate a BTES field model containing 5 2000 m-deep boreholes in Xi'an, China.Başer and McCartney (2020) validated a COMSOL-developed model of a solar assisted BTES field consisting of 13 boreholes with 15 m depth in California by using monitored transient ground temperature at several different depths.Fadejev and Kurnitski (2015) used IDA ICE to model and analyse the performance of a large-scale GSHP system with 196 15 m-deep borehole heat exchangers of in Helsinki, Finland.In their study, the borehole heat exchanger model was validated by monitored outlet brine temperature.However, the aforementioned studies were conducted on grout filled boreholes.According to author's best knowledge, there are not many studies using monitoring data to validate BTES field model with groundwater-filled boreholes.In this paper, a large-scale BTES field consisting of 74 groundwater-filled boreholes used for an educational building in Finland was developed in IDA ICE and validated by measured 1.5-year BTES field inlet and outlet temperature data.This BTES field model will be combined in a hybrid building energy system model for hybrid GSHP system optimization in further work.

Description of the BTES system
The studied BTES field is located under an educational building in Otaniemi, Finland.The BTES field was already activated during the building construction phase in March 2018, while the building was fully open into use in June 2019.The BTES field is connected to nine centralized heat pump modules integrated in series.The total nominal heat capacity of the heat pump modules is 790 kW at the rating conditions (0/35 ℃).
The educational building is a 4/5-story building containing over 47500 m 2 of floorspace which comprises educational area, office, restaurants, gym, workshop, shopping center and a metro station.The building is equipped with heating and cooling systems.The heating is delivered through air handling unit (AHU) heating (heating is used for ventilation air), space heating and snow melting circuits.The supply water temperature for space and AHU heating is controlled within 30-45℃ according to the outdoor air temperature.The space heating is delivered through radiant ceiling panels.The space and AHU heating are mainly supplied by heat pumps in wintertime.During peak heat load, the building uses district heating for auxiliary heating.The district heating is also used for domestic hot water generation.
The cooling is provided through AHU cooling (cooling is used for ventilation air) and space cooling circuits.The space cooling is delivered through radiant panels with supply/return temperatures of 12-16 ℃/15-18 ℃.The annual cooling demand of radiant ceiling panels is mostly meet by borehole free cooling.The AHU cooling is provided mostly by heat pump evaporators operating only in peak cooling load time with average supply/return temperatures of 10 ℃/16 ℃.The BTES field consists of 74 groundwater-filled boreholes installed with single U-pipe heat exchangers.The average borehole depth is around 310 m.All 74 boreholes are drilled in homogenous granite bedrock on which 10 m soil covers.The groundwater level is at 8.2-9 m below the ground surface.The boreholes are arranged in an irregular shape shown in Fig. 1.The borehole heat exchangers are all connected in parallel.The inlet and outlet brine temperatures of the BTES field was monitored by PT100 sensors.The total brine flow rate and injected/extracted heat of the BTES field were measured by thermal energy meters (Todorov, Alanne, et al., 2021).The detailed borehole parameters are listed in Table 1.
Figure 1.BTES field layout, the black line represents the outline of the building (Todorov, Vallin, et al., 2021).

Numerical modelling of BTES system
In this study, a 3-D numerical model of the BTES field was developed in IDA ICE 4.8.IDA ICE 4.8 is a multizone simulation software for the study of thermal indoor climate and energy consumption of entire building with variable simulation time step.The IDA ICE GHX module was used for BTES field modelling.The IDA ICE GHX module uses finite difference method to calculate 2-D temperature fields.The 2-D temperature fields are then combined by means of superposition to generate a 3-dimensional ground field.The model can only be used for boreholes with equal depth.The ground, boreholes and U-pipes are divided into several layers, respectively.Temperatures of the borehole and the brine fluid are assumed as homogeneous.The model is restrained to use U-pipe borehole heat exchanger.Besides, the model does not consider groundwater flow movement.The actual ground temperature at the borehole wall is calculated by superposition.
For each borehole, the following temperature fields are calculated (Fadejev and Kurnitski, 2015): • 1-D heat transfer in downward and upward flow in Upipe with heat transfer to the borehole filling material and the ground.• 1-D heat transfer in borehole filling material with heat transfer to the brine fluid and the ground.• 2-D heat transfer in cylindrical coordinates around the borehole with heat transfer to the borehole filling material and the brine fluid.• 1-D heat transfer in an undisturbed ground field in the same way to the heat tranfer in the axial ground field around the boreholes.In multiple borehole calculation, first the thermal behavior of the ground surrounding each single borehole is calculated.Then, the effect of all boreholes on the ground temperature is summed up through superposition.
Since the borehole is filled with ground water in the target BTES field, the natural convection of water will affect the heat transfer from the working fluid to the ground.In order to consider this effect, Johnsson and Adl-Zarrabi ( 2019) proposed to use effective thermal conductivity of water in simulation.In their study, it is recommended that the effective thermal conductivity of the ground water should be within 1.2-2.0W/mK which is about 2-3 times larger than the actual water thermal conductivity to achive the effective borehole thermal resistance obtained from thermal response test.In this study, the effective thermal conductivity of the water was used in the borehole modelling to consider the natural convection effect.The heat transfer for downward and upward flows in the U-pipe are described in Eq. ( 1) and ( 2) (Fadejev and Kurnitski, 2015): (2) where ρf is the density of the brine fluid, kg/m 3 ; cp,f is the specific heat capacity of the brine fluid, J/kgK; Vf is the volume of each pipe element, m 3 ; mi is the mass flow of the brine fluid in borehole i, kg/s; Td,i,j and Tu,i,j are temperatures of the down and up flowing brine fluid respectively in node j of borehole i, ℃; Tg,d,i,j and Tg,u,i,j are temperatures of the borehole filling materials around down and up flow pipe(s) respectively at layer j of borehole i, ℃; Tw,i,j is the temperature at the borehole wall for borehole i in node j, ℃; Kfg,i is the heat transfer coefficient between the brine fluid and borehole filling material in borehole i, W/K; and Kfr,i is the heat transfer coefficient between the brine fluid and the surrounding bedrock near borehole i, W/K.The heat transfer for the borehole filling material is described as Eq. ( 3)-( 5) (Fadejev and Kurnitski, 2015).
2 , , ,  , 1 , , , , , , , , , , , , , , , , g u i j p g gg g i j g u i j fg i u i j g u i j grr w i j g u i j (5) where Mcp,g1 and Mcp,g2 are absolute heat capacities of inner and outer borehole filling materials respectively, J/K; Tg,i,j is the temperature of outer borehole filling material in node j of borehole i, ℃; Tg,u,i,j and Tg,d,i,j are borehole filling material temperatures around upflow and downflow pipe(s) respectively at node j of borehole i, ℃; Tw,i,j is the bedrock temperature at borehole i in node j, ℃; Kgg is the heat conductivity coefficient between inner the outer borehole filling material rings, W/K; Kfg,i is the heat transfer coefficient between the brine fluid and borehole filling material in borehole i, W/K; Kgrr is the heat conductivity coefficient between the outer borehole filling material ring and the bedrock, W/K.
The heat transfer for the ground ring is described as Eq. ( 6): where Mcp,r,j,k is the absolute heat capacity of the ground ring associated to nodes (j, k), J/K; DTi,j,k is the temperature change due to conditions in borehole i at node (j, k), ℃; Krr,k is the heat conductivity coefficient in radial direction between radial node k and k+1 at axial nodes j, W/K; Kzr,k is the heat conductivity coefficient in axial direction between axial nodes j and j+1 at radial nodes k, W/K.
The diagram of the BTES field model is shown in Fig. 2. The borehole model is connected to a simplified building model through a pump.The simplified building model is connected to daily BTES field injected/extracted heat power data.The injected/extracted heat profile during validation period is shown in Fig. 3, which is processed based on measured data (Todorov, Alanne, et al., 2021).In Fig. 3, the positive value refers to the injected heat to the BTES field, which is used for cooling of the building; the negative value refers to the extracted heat from the BTES field, which is used for heating of the building.The brine pump flow rate variates according to the mass flow rate data profile, which is shown in Fig. 4. The mass flow rate data profile is processed based on measured brine flow rate data (Todorov, Alanne, et al., 2021).The inlet and outlet brine temperatures of BTES follow the energy balance on the building side as following: ( ) where Tin and Tout are the inlet and outlet brine temperatures of BTES, ℃; Qfield is the BTES injected/extracted heat power, W; qm,f is the total mass flow rate of the BTES field, kg/s; cp,f is the specific heat capacity of the brine fluid, J/kgK.
In the borehole model, as all 74 boreholes are connected in parallel, the total brine flow is assumed to be evenly distributed among all boreholes.The downward brine temperature in each borehole at the surface of the ground (the depth of 0 m) is prescribed as equal to the inlet brine temperature to the BTES field: where Td,i,1 is the downward brine temperature in borehole i at the surface of the ground (depth of 0 m), ℃.At the bottom of the U-pipe, the downward and upward fluid temperatures are prescribed to be the same: , , , , where Td,i,nzHo is the bottom of the downward brine fluid in the U-pipe, ℃; Tu,i,nzHo is the bottom of the upward brine fluid in the U-pipe, ℃.
The boundary condition for the ground surface was given by the ambient air temperature which varies hourly.The ambient air temperature data was taken from data monitored by the weather observation station of the Finnish Meteorological Institute in Tapiola, Espoo, which is around 2 km away from the studied BTES field.In the BTES model, the thermal conductivity of U-pipe was set as 0.42 W/(mK), based on Reuss (2015).The thermal conductivity, specific heat capacity and density of the brine fluid were estimated at the borehole heat exchanger average temperature of 5.65 °C.The thermal conductivity, specific heat capacity and density of the brine fluid were estimated to be 0.42 W/m•K, 4243 J/kg•K and 960.9 kg/m 3 respectively by using polynomial interpolation equation of Melinder (2007).The effective thermal conductivity in groundwater-filled borehole was enhanced to 1.6 W/(mK) to consider natural convection in boreholes (Johnsson and Adl-Zarrabi, 2019).The density and specific heat capacity of the groundwater were estimated to be 999 kg/m 3 and 4200 J/kg•K, respectively.The effective thermal conductivity, specific heat capacity and density of the ground were 3.3 W/(mK) (Todorov, Vallin, et al., 2021) , 725 J/(kgK) (Janiszewski et al., 2018) and 2500 kg/m 3 (Janiszewski et al., 2018).The undisturbed ground temperature and the geothermal temperature gradient were 8.7 ℃ and 0.0119 ℃/m, obtained from on-site thermal response test.In heat extraction mode when the brine temperature is lower than the borehole wall temperature, the borehole thermal resistance was set as 0.095 mK/W, based on on-site thermal response test.In heat injection mode when the brine temperature is higher than the borehole wall temperature, the borehole thermal resistance was set slightly higher as 0.0977 mK/W, to mitigate the effect from variable brine flow rate.
The IDA ICE model used the real dimensioning of the BTES field.The borehole array arrangement was set the same to the real BTES field shown in Fig. 1.The geometries of boreholes and heat exchangers were set as the same parameters listed in Table 1.The simulation period was set from March 2018 when the BTES field was already activated to February 2021.As the target building was fully into use in June 2019, the latter 1.5 year of the simulation period (from August 2019 to February 2021) was taken as the validation period.

Results
Figs. 5 shows the IDA ICE simulated and measured inlet and outlet brine temperatures of the BTES field, respectively.The model estimated the inlet and outlet brine temperatures of the BTES field well with nonsignificant difference between heating and cooling seasons.However, the deviation in some days in August 2019 and January 2021 was relatively larger, which could be due to the real-time measurement error of the mass flow rate.Table 2 shows the average inlet and outlet brine temperature differences between IDA ICE simulation and measurement during different time periods.The differences against the measured inlet brine temperature were higher in cooling season than heating season, while the differences against measured outlet brine temperature was higher in heating season than cooling season.The average deviation of the inlet and outlet brine temperatures in the whole validation period were 0.78℃ in and 0.47℃ respectively, which shows good prediction accuracy.

Discussion
The ground water movement was not taken into account in the IDA ICE simulation of the studied BTES field.In this study, the actual groundwater level is monitored at 8.2-9 m below the ground surface.It could slightly affect the heat transfer of the top part of boreholes and lead to deviation between the simulation and measured temperature result.However, this effect is still unknown for long-term predictions.The effect may occur slowly and needs several years to discover.Due to lack of study, the ground water was supposed to have ignored effects on the simulation results.However, for other region of the world, the geology and groundwater flow rates in the bedrock might show various difference.Therefore, whether the groundwater movement effects can be neglected or not depends on the actual geology conditions.The borehole thermal resistance can be affected by the brine mass flow rate which varies over time.In this study, the borehole thermal resistance was set as two constant values for heat injection and extraction modes respectively.This approach reduces the difference of estimation accuracy between heating and cooling seasons as the result shows.However, more accurate prediction can be realized by using the real-time updated thermal resistance.This could be the next step work for the future study.
In heating-dominated regions, the building has much more heating demand than the cooling demand.There will be more heat extracted from the ground in heating seasons than the heat injected to the ground in cooling seasons.In this study, the annual extracted heat from the ground was 8.4 times of the injected heat to the ground.This underground thermal imbalance will account for a temperature drop in the ground.Before this study, a longterm ground simulation was conducted for the target BTES field by Earth Energy Designer (EED) during the design phase, which has predicted a ground temperature drop of 6-7℃ over 25 years (Nadas, 2020).Therefore, a ground thermal recharging strategy could be considered to maintain a stable GSHP system performance in the future.

Conclusion
In this study, a large-scale BTES field in Finland with 74 groundwater-filled boreholes was modelled in IDA ICE.This model was validated by 1.5-year brine temperatures at inlet and outlet of the BTES field measured by PT100 temperature sensors.The simulated temperatures matched well with the measured data.The average difference against measured inlet and outlet brine temperatures of the BTES field were within 1 ℃.The BTES field model will be used for further energy optimization in a hybrid GSHP system in the follow-up work.
government.Authors would like to thank Antti Säynäjoki from Aalto University Campus & Real Estate Ltd. for fruitful cooperation.In addition, authors would like to specially thank Oleg Todorov and Markku Virtanen for the arrangement of measurement data for this study.

Figure 3 .
Figure 3. Injected/extracted heat power during the validation period.

Figure 4 .
Figure 4. Brine mass flow rate during the validation period.

Figure 5 .
Figure 5.Comparison of simulated and measured inlet and outlet brine temperatures of BTES field.Table2shows the average inlet and outlet brine temperature differences between IDA ICE simulation and measurement during different time periods.The differences against the measured inlet brine temperature were higher in cooling season than heating season, while the differences against measured outlet brine temperature was higher in heating season than cooling season.The average deviation of the inlet and outlet brine temperatures in the whole validation period were 0.78℃ in and 0.47℃ respectively, which shows good prediction accuracy.

Table 2 :
Average brine temperature difference between IDE ICE simulation and measurement during heatingand cooling seasons and the validation period.