# Long-Term Investigation of a Borehole Heat Exchanger (BHE) Array with the Influence of Same Inlet Temperature (SIT) and Same Heat Load (SHL) Boundaries

Authors: Yuping Zhang; Xiao Yang; Jun Liu; Boyang Liu; Fujiao Tang; Hossein Nowamooz
Journal: Advances in Civil Engineering
 (2023)

---

## Abstract

Borehole heat exchanger (BHE) array is a key element in the ground source heat pump system (GSHPS). The BHE array is traditionally simplified by assuming the same heat exchange rate for each BHE, which normally is not correspondent to the real applications. In this investigation, a TRT was firstly conducted for a site located in Shanxi Province (China) to get the local ground thermal properties. A numerical simulation framework was next verified by the TRT results. Subsequently, the ground thermal properties obtained by the TRT were used in the numerical simulation model to investigate the influence of two operation strategies on the ground temperature development of the surrounding soil and the BHE performance. The results showed that both ground temperature and BHE array performance decreased during the 8-year operation period of the BHE array in the heating mode. The ground temperature around the center BHE was lower in the same heat load (SHL) scenario (same heat load for the BHEs in the array) than in the same inlet temperature (SIT) scenario (same carrying fluid inlet temperature for the BHEs in the array), and this difference can be 2 °C in the 7th year. The annual average outlet temperature in the SHL scenario was also lower than that in the SIT scenario, especially in the long term. The investigation shows that SIT working strategy superiors the SHL strategy, providing a reference for the upcoming investigation on the design aspect of the GSHPS.


---

## Body

## 1. Introduction

Ground source heat pump (GSHP) is currently a successful technology for heating and cooling due to its stability, environmental friendliness, and wide availability [1], which possesses a considerable potential in reducing the carbon emission [2, 3]. The GSHP can offer higher energy efficiency for air conditioning than conventional systems due to the higher temperature for heating and lower temperature for cooling the underground environment provided [4, 5]. Energy demand is normally provided by multiple borehole heat exchanger (BHE) since a single BHE traditionally cannot provide adequate energy for heating and cooling [6–9].During the service period of a BHE array, heat would be accumulated in the vicinity of boreholes in the cooling season, and cold would be gathered around the boreholes in the heating season [10]. If the heat and cold around the BHE could not be compensated during the recovery season, ground temperature imbalance occurs, which is a major threat deteriorating the performance of a GSHP system [11, 12]. The temperature imbalance varies along the BHE [13] and can lead to mechanical problems [14]. The imbalance increases with the continuous operation of the BHE array, and changing parameters such as the intermittent ratio and separation distance can relieve the ground imbalance [15].Heat load influences the ground temperature imbalance, especially in the heating or cooling only applications [16, 17]. The average soil temperature could decrease by 0.5–1.0 °C per year in the long-term observation of heating only scenario [18], depending on the borehole designs, GSHP performance, load profiles, and soil conditions [19]. The heating performance of the GSHP system could lower down by 2–4% with a 1.0 °C carrying fluid temperature decrease [20, 21]. Generally, the ground temperature imbalance problem increases by year when the seasonal heat load was kept constant, while the increase rate decreases [10].Thermal interaction is another major factor influencing the ground temperature imbalance, which in the single BHE scenario could be neglected. The larger the spacing between the BHEs, the smaller the effect of heat interaction. To avoid substantial interference between BHEs, a distance of 7-8 m is recommended by Signorelli et al. [22]. However, the BHEs installation area is normally limited, and the optimal spacing should be adopted during the design stage. In practical applications, 5 m of distance is in favorable [23]. Gultekin et al. [24] concluded that a 4.5 m of spacing is enough to keep the total performance losses within 10% for a period of 2400 h nonstop operation. The BHE layout is another important parameter influencing the thermal interaction of BHEs. Chen et al. [11] found that the subarray with a larger number of installed BHEs is shifting its thermal load towards the area with less BHEs installed, and the heat load on the central BHEs is gradually shifted towards those located at the edge. Giordano and Raymond [25] concluded that square-shaped BHEs have a lower heat loss compared to a conventional circular design. Zhang et al. [26] obtained that the performance of a BHE array increased by 14.82% and 12.10% in cooling mode, and 20.62% and 16.98% in heating mode, compared with staggered and square arrangements. There are also investigations of groundwater flow on the influence BHE array thermal interactions. Adopting higher thermal load for the outer and downstream BHEs with the existence of groundwater flow is helpful to balance the ground temperature [27].Considerable efforts have been made to study the ground temperature imbalance of BHE array. However, the investigations usually simply use the constant carrying fluid inlet temperature or average heat load as the boundary condition of each BHE. In real working conditions, the fluid inlet temperatures and flow rates are usually assumed to be identical for all boreholes. The heat pump outlet temperature is equal to the BHE-carrying fluid inlet temperature and is distributed to all boreholes at nearly the same temperature. Their outlet temperatures are different due to ground temperature variation surrounding the BHEs, leading to the difference in the heat exchange rate for each borehole [28, 29]. In practice, the heating/cooling load of the BHE array can be controlled, indicating the variation of carrying fluid inlet temperature with the operation time. The ground temperature change can also lead to the variation of carrying fluid temperature. Therefore, simply using the constant carrying fluid inlet temperature or average heat load as the boundary condition of each BHE is not able to reflect the real working condition of a BHE array.In this investigation, a TRT is firstly conducted in-situ to get the local geotechnical parameters. A numerical simulation model is furthermore verified by the TRT. Afterwards, the working mode of the BHEs with the same carrying fluid inlet temperature is going to be studied (the fluid inlet temperature is going to change with time according to the heat load and ground temperature) and compared with the simplified method by considering the same heat load for each BHE. The investigation is able to optimize the design of a BHE array, especially for a long-term operation of a GSHP system.


## 2. Methodology and Numerical Simulation Verification

Numerical simulation is an approach to studying effectively and economically the performance of the BHE array [30], and a robust numerical simulation model should be verified with the in-situ experiment.


### 2.1. Methodology

This investigation is based on the numerical simulation, which mainly contains the heat transfer in the surrounding soil, in the grout, in the pipe, and in the circulating fluid.The heat transfer in the surrounding soil, in the grout, and in the pipe is dominated by heat conduction, while heat transfer in the carrying fluid is mainly determined by heat convection. It should be noted that the pipe is normally with the thickness of several millimeters, and heat storage in the pipe material is often neglected. The heat conduction is given by(1)ρC∂T∂t=∇⋅k∇T+Q,where ρ is the medium density (kg·m−3), C is the medium-specific heat capacity (J·kg−1·K−1), T is the medium temperature (K), k is the medium thermal conductivity (W·m−1·K−1), and Q is the medium heat source (W·m−3).The heat transfer in the carrying fluid is mainly dominated by heat convection, and the governing equation is given by(2)AρfCp_f∂Tf∂t+AρfCp_fuf⋅∇Tf=∇⋅Akf∇Tf+Qwall,where A is the pipe inner cross-sectional area (m2), ρf is the fluid density (kg·m−3), Cp-f is the fluid specific heat capacity (J·kg−1·K−1), Tf is the fluid temperature (°C), uf is the fluid flowing velocity (m·s−1), kf is the fluid thermal conductivity (W·m−1·K−1), and Qwall is the energy from the surrounding media (W·m−1).A TRT was conducted to get the initial temperature and the soil thermal properties (thermal conductivity and volumetric heat capacity). This section firstly introduces the in-situ TRT, and the TRT is further used to verify the capacity of the numerical simulation framework. The TRT results are also used in the next section.


### 2.2. In-Situ TRT

The in-situ TRT was conducted in Xi’an, Shaanxi, China (Figure1). The borehole diameter was 0.09 m, with the depth of 270 m. The grout material is with the thermal conductivity of 3 W·m−1·K−1, density of 2500 kg·m−3, and specific heat capacity of 800 J·Kg−1·K−1. The test pipe is a single high-density polyethylene (HDPE) U-pipe (De40), with the thermal conductivity of 0.4 W·m−1·K−1, and water is the carrying fluid in the pipe, with the flow rate of 0.6 m·s−1. The average temperature difference between the inlet and the outlet of the pipes was kept 2.28 °C.Figure 1
Test site.The test has lasted for more than 3 days. The ambient, inlet, and outlet temperatures during the TRT are shown in Figure2.Figure 2
Recorded ambient, inlet, and outlet temperatures.According to the fluid inlet temperature, outlet temperature, and time, the average fluid temperature with time is shown in Figure3. And soil thermal conductivity can therefore be obtained as 2 W·m−1·K−1.Figure 3
The variation of carrying fluid temperature with time to get the ground thermal properties.


### 2.3. Numerical Simulation Verification

The above in-situ measurement was further adopted to verify the numerical simulation model used in the following investigation. The difference between the measurements (Tin-situ) and the predictions (Tnumierical) is presented by ε, given by(3)ε=Tnumerical−Tin−situTnumerical×100%.The comparison between the in-situ measurement and the numerical simulation is shown in Figure4. The comparison shows that the relative error between the measurement and the numerical simulation is less than 5%, proving the accuracy of the numerical simulation model.Figure 4
Comparison between the numerical simulation and in-situ measurement.


## 2.1. Methodology

This investigation is based on the numerical simulation, which mainly contains the heat transfer in the surrounding soil, in the grout, in the pipe, and in the circulating fluid.The heat transfer in the surrounding soil, in the grout, and in the pipe is dominated by heat conduction, while heat transfer in the carrying fluid is mainly determined by heat convection. It should be noted that the pipe is normally with the thickness of several millimeters, and heat storage in the pipe material is often neglected. The heat conduction is given by(1)ρC∂T∂t=∇⋅k∇T+Q,where ρ is the medium density (kg·m−3), C is the medium-specific heat capacity (J·kg−1·K−1), T is the medium temperature (K), k is the medium thermal conductivity (W·m−1·K−1), and Q is the medium heat source (W·m−3).The heat transfer in the carrying fluid is mainly dominated by heat convection, and the governing equation is given by(2)AρfCp_f∂Tf∂t+AρfCp_fuf⋅∇Tf=∇⋅Akf∇Tf+Qwall,where A is the pipe inner cross-sectional area (m2), ρf is the fluid density (kg·m−3), Cp-f is the fluid specific heat capacity (J·kg−1·K−1), Tf is the fluid temperature (°C), uf is the fluid flowing velocity (m·s−1), kf is the fluid thermal conductivity (W·m−1·K−1), and Qwall is the energy from the surrounding media (W·m−1).A TRT was conducted to get the initial temperature and the soil thermal properties (thermal conductivity and volumetric heat capacity). This section firstly introduces the in-situ TRT, and the TRT is further used to verify the capacity of the numerical simulation framework. The TRT results are also used in the next section.


## 2.2. In-Situ TRT

The in-situ TRT was conducted in Xi’an, Shaanxi, China (Figure1). The borehole diameter was 0.09 m, with the depth of 270 m. The grout material is with the thermal conductivity of 3 W·m−1·K−1, density of 2500 kg·m−3, and specific heat capacity of 800 J·Kg−1·K−1. The test pipe is a single high-density polyethylene (HDPE) U-pipe (De40), with the thermal conductivity of 0.4 W·m−1·K−1, and water is the carrying fluid in the pipe, with the flow rate of 0.6 m·s−1. The average temperature difference between the inlet and the outlet of the pipes was kept 2.28 °C.Figure 1
Test site.The test has lasted for more than 3 days. The ambient, inlet, and outlet temperatures during the TRT are shown in Figure2.Figure 2
Recorded ambient, inlet, and outlet temperatures.According to the fluid inlet temperature, outlet temperature, and time, the average fluid temperature with time is shown in Figure3. And soil thermal conductivity can therefore be obtained as 2 W·m−1·K−1.Figure 3
The variation of carrying fluid temperature with time to get the ground thermal properties.


## 2.3. Numerical Simulation Verification

The above in-situ measurement was further adopted to verify the numerical simulation model used in the following investigation. The difference between the measurements (Tin-situ) and the predictions (Tnumierical) is presented by ε, given by(3)ε=Tnumerical−Tin−situTnumerical×100%.The comparison between the in-situ measurement and the numerical simulation is shown in Figure4. The comparison shows that the relative error between the measurement and the numerical simulation is less than 5%, proving the accuracy of the numerical simulation model.Figure 4
Comparison between the numerical simulation and in-situ measurement.


## 3. Necessary Parameters for the Numerical Simulation Model

In the above section, the initial soil thermal properties were obtained by TRT. The same geotechnical and climate conditions were considered in the numerical simulation model as mentioned in the previous TRT. The detailed local temperature and other necessary parameters would be introduced in this section.


### 3.1. Local Temperature

A sinusoidal equation can be used to approximate the measured local ambient temperature (Figure5), with the highest daily average temperature of 30 °C and the lowest daily temperature of −5 °C. The sinusoidal equation is further used as the surface boundary condition for the numerical simulation model. The site is with the ground temperature gradient of 0.06 K·m−1.Figure 5
Approximated and measured ambient temperature.


### 3.2. Geometry and Mesh

There are in total 36 pipes included in the BHE array, and they are symmetrically installed, with a spacing of 5 m (Figure6). In order to save the calculation time, ¼ of the pipes are investigated. It should be noted that the selected geometry should be changed when there is groundwater flow according to the seasonal groundwater flow direction. The numerical simulation model geometry has both a width and a length of 100 m. The borehole is with the depth of 80 m.Figure 6
The whole and the selected geometry for the numerical simulation model.There are in total 188705 tetrahedra and 552000 prisms in the numerical simulation model. The mesh size could guarantee the accuracy of the numerical simulation model.


### 3.3. Other Parameters

A 25% of glycol mixing with water was considered for the carrying fluid [20], with the circulating velocity of 0.5 m·s−1. The pipe is with the thickness of 0.15 cm and the inner diameter of 4 cm.The BHEs are used for heating during the winter (120 days of operation during 1 year). And the heat load for the 9 pipes is shown in Figure7.Figure 7
Total heat load for the 9 pipes during 5 years.


## 3.1. Local Temperature

A sinusoidal equation can be used to approximate the measured local ambient temperature (Figure5), with the highest daily average temperature of 30 °C and the lowest daily temperature of −5 °C. The sinusoidal equation is further used as the surface boundary condition for the numerical simulation model. The site is with the ground temperature gradient of 0.06 K·m−1.Figure 5
Approximated and measured ambient temperature.


## 3.2. Geometry and Mesh

There are in total 36 pipes included in the BHE array, and they are symmetrically installed, with a spacing of 5 m (Figure6). In order to save the calculation time, ¼ of the pipes are investigated. It should be noted that the selected geometry should be changed when there is groundwater flow according to the seasonal groundwater flow direction. The numerical simulation model geometry has both a width and a length of 100 m. The borehole is with the depth of 80 m.Figure 6
The whole and the selected geometry for the numerical simulation model.There are in total 188705 tetrahedra and 552000 prisms in the numerical simulation model. The mesh size could guarantee the accuracy of the numerical simulation model.


## 3.3. Other Parameters

A 25% of glycol mixing with water was considered for the carrying fluid [20], with the circulating velocity of 0.5 m·s−1. The pipe is with the thickness of 0.15 cm and the inner diameter of 4 cm.The BHEs are used for heating during the winter (120 days of operation during 1 year). And the heat load for the 9 pipes is shown in Figure7.Figure 7
Total heat load for the 9 pipes during 5 years.


## 4. Numerical Simulation Results

Two operation strategies were investigated in the numerical simulation model (Figure8). The same heat load (SHL) mode assumes the circulating fluid flows through a main pipe to the branch pipes, and the outlet temperatures from the branch pipes would be collected into the same main pipe. Therefore, the BHEs have the same inlet fluid temperature, while the temperature varies when the circulating fluid flows out from the BHEs. The same inlet temperature (SIT) mode assumes that the BHEs work separately, having different inlet and outlet carrying fluid temperatures.Figure 8
Two operation strategies for the BHE array: (a) SHL mode and (b) SIT mode.
(a)(b)According to the definition of the two working scenarios, the total heat load for the pipe array with 9 BHEs is shown in the following equations:(4)QSIT=∑i=19Acv−fufTi−out−Tin,QSHL=∑i=19Acv−fufTi−out−Ti−in.In the investigation, the total heat load for the BHE array in the two scenarios was kept constant, indicatingQSIT = QSHL. In this way, the two working scenarios are compared by the variation of ground temperature and carrying fluid outlet temperature.


### 4.1. Ground Temperature Variation

Figure9 shows the temperature field of the two working scenarios after 2400 days (7th year, day 210) of heat extraction. The results show that the surrounding temperature of the BHEs is abnormally low due to the heat exchange between the pipe and the surrounding medium. The closer the center of the BHE array, the lower the soil temperature. The middle BHEs have lower temperature in the same HER scenario, and the temperature could reach −7 °C, which is 2 °C lower than the same inlet temperature scenario.Figure 9
Two different working scenarios on day 2400 (7th year, day 210): (a) same HER and (b) same inlet fluid temperature.
(a)(b)Figure10 shows the ground temperature contour in the vicinity of the BHEs during the service period of the BHEs. The result shows that the surrounding temperature decreases by year, and a larger area would be influenced. Specifically, the surrounding temperature on day 210 of the 1st year is 8 °C higher than that in the 7th year.Figure 10
Ground temperature contour around the BHEs for the inlet temperature scenario after (a)t = 210 d (1st year, day 210); (b) t = 1305 d (4th year, day 210); (c) t = 2400 d (7th year, day 210).
(a)(b)(c)Figure11 shows the temperature field of the surrounding soil after the recovery season for the 2nd, 5th, and 8th years. The result shows that the surrounding temperature around the BHEs is not able to fully recover. The surrounding temperature even decreases by year after the temperature recovery season (from spring to autumn). Specifically, the temperature can recover to 10 °C after the recovery season in the second year (2nd year, day 90), while the surrounding temperature can only be recovered to −1.5 °C in the 8th year (8th year, day 90).Figure 11
Ground temperature during the recovering season: (a)t = 455 d (2nd year, day 90); (b) t = 1550 d (5th year, day 90); (c) t = 2645 d (8th year, day 90).
(a)(b)(c)


### 4.2. Outlet Temperature of the Carrying Fluid

The outlet temperature for the two scenarios during the service season of the BHE array is shown in Figure12. The results show that the outlet temperature decreases during the service season of a year, and the difference between the beginning and the end of the service season can be as large as 7 °C. Since the service of the BHE array causes the cold accumulation surrounding the BHEs, the outlet temperature decreases by year. For the two scenarios, the outlet temperature difference for the two scenarios enlarges by year.Figure 12
Average outlet temperature during the service period (8 years) of the BHEs for the 2 scenarios.The annual average fluid outlet temperature for the two scenarios during the 8-year operation of the BHE array is given in Figure13. The result shows that the annual average carrying fluid outlet temperatures for the two scenarios decrease by year. Specifically, the average carrying fluid outlet temperatures in the 1st and the 8th year for the same inlet temperature scenario are 5.08 and −5.13 °C, respectively, while the average fluid outlet temperatures for the same HER scenario are 5.07 and −5.33 °C in the first and the last year. The result shows also that the SIT mode has higher annual outlet temperature than the SHL mode, indicating a better performance.Figure 13
Annual average carrying fluid outlet temperature during the 8-year operation of the BHEs.


## 4.1. Ground Temperature Variation

Figure9 shows the temperature field of the two working scenarios after 2400 days (7th year, day 210) of heat extraction. The results show that the surrounding temperature of the BHEs is abnormally low due to the heat exchange between the pipe and the surrounding medium. The closer the center of the BHE array, the lower the soil temperature. The middle BHEs have lower temperature in the same HER scenario, and the temperature could reach −7 °C, which is 2 °C lower than the same inlet temperature scenario.Figure 9
Two different working scenarios on day 2400 (7th year, day 210): (a) same HER and (b) same inlet fluid temperature.
(a)(b)Figure10 shows the ground temperature contour in the vicinity of the BHEs during the service period of the BHEs. The result shows that the surrounding temperature decreases by year, and a larger area would be influenced. Specifically, the surrounding temperature on day 210 of the 1st year is 8 °C higher than that in the 7th year.Figure 10
Ground temperature contour around the BHEs for the inlet temperature scenario after (a)t = 210 d (1st year, day 210); (b) t = 1305 d (4th year, day 210); (c) t = 2400 d (7th year, day 210).
(a)(b)(c)Figure11 shows the temperature field of the surrounding soil after the recovery season for the 2nd, 5th, and 8th years. The result shows that the surrounding temperature around the BHEs is not able to fully recover. The surrounding temperature even decreases by year after the temperature recovery season (from spring to autumn). Specifically, the temperature can recover to 10 °C after the recovery season in the second year (2nd year, day 90), while the surrounding temperature can only be recovered to −1.5 °C in the 8th year (8th year, day 90).Figure 11
Ground temperature during the recovering season: (a)t = 455 d (2nd year, day 90); (b) t = 1550 d (5th year, day 90); (c) t = 2645 d (8th year, day 90).
(a)(b)(c)


## 4.2. Outlet Temperature of the Carrying Fluid

The outlet temperature for the two scenarios during the service season of the BHE array is shown in Figure12. The results show that the outlet temperature decreases during the service season of a year, and the difference between the beginning and the end of the service season can be as large as 7 °C. Since the service of the BHE array causes the cold accumulation surrounding the BHEs, the outlet temperature decreases by year. For the two scenarios, the outlet temperature difference for the two scenarios enlarges by year.Figure 12
Average outlet temperature during the service period (8 years) of the BHEs for the 2 scenarios.The annual average fluid outlet temperature for the two scenarios during the 8-year operation of the BHE array is given in Figure13. The result shows that the annual average carrying fluid outlet temperatures for the two scenarios decrease by year. Specifically, the average carrying fluid outlet temperatures in the 1st and the 8th year for the same inlet temperature scenario are 5.08 and −5.13 °C, respectively, while the average fluid outlet temperatures for the same HER scenario are 5.07 and −5.33 °C in the first and the last year. The result shows also that the SIT mode has higher annual outlet temperature than the SHL mode, indicating a better performance.Figure 13
Annual average carrying fluid outlet temperature during the 8-year operation of the BHEs.


## 5. Conclusion

The investigation aims to study the long-term performance of BHE array considering two working scenarios (same inlet temperature and same heat load).A TRT was conducted firstly in-situ to get the local soil thermal properties, and the BHE array installed in the same location was studied numerically. The results show that both ground temperature and BHE performance decrease with year in the heat extraction during the service time of the BHE array, and ground temperature could not be fully compensated after the temperature recovery season. Ground temperature in the vicinity of the BHE in the SHL scenario during the service period of the BHEs is lower than that in the SIT scenario. The difference is even larger with the longer service period of the BHEs. Specifically, the ground temperature around the center BHEs in the SHL scenario is around 2 °C lower than the SIT scenario in the 7th year. The outlet temperature in the SIT scenario is higher than that in the SHL scenario, which is also getting more obvious with year. The difference can be 0.2 °C in the 8th year.The results show that SIT strategy is better than the SHL strategy. It should be noted that there are also other similar operation strategies like zoning operation strategy [31]. The combination of the SIT and these strategies has the potential to enhance the performance of BHE array.



---

*Source: 3232062-2023-02-28.xml*