Simulation of the thermo-hydraulic response of energy piles in unsaturated soils

. This paper focuses on the simulation of the coupled heat transfer and water flow in unsaturated soil layers surrounding a solitary energy pile undergoing heating and cooling cycles typical of a field-scale energy pile. The results indicate that heating leads to drying of the soil surrounding the energy pile, which has been shown in previous studies to result in an increase in axial capacity. During cooling, the degree of saturation was observed to recover to the value present before the start of heating initially, however, it will not recover in the following years. Which will lead to a cumulative effect after several cycles of heating and cooling. Heating and cooling cycles lead to an overall reduction in the thermal conductivity of the subsurface, reducing the heat transfer from the energy pile but also leading to greater storage of heat in the subsurface surrounding the pile.


Introduction
Energy piles are deep foundations used for supporting buildings and exchanging heat between the ground and an overlying structure. They can help improve the energy efficiency of buildings and can be used as part of heat storage systems. Although energy piles are often installed in unsaturated soil deposits, questions remain on the impacts of heating and cooling on the response of surrounding unsaturated soil profiles and soil-structure interaction. Başer et al. (2018) found that a heatingcooling cycle applied to a geothermal heat exchanger in an unsaturated soil layer led to a permanent decrease in degree of saturation around the heat exchanger, as well as a corresponding decrease in thermal conductivity and volumetric heat capacity. It is well established that the thermal conductivity of soil has a major effect on the heat transfer and efficiency of geothermal heat exchangers, along with the heat exchanger dimensions, heat exchanger spacing, inlet heat exchange fluid temperature, and heat exchange fluid flow rate (Bourne-Webb 2013; Catolico et al. 2016). Unsaturated conditions also affect soil-structure interaction in energy piles due to changes in soil stress state ( Behbehani and McCartney (2020) showed that the increases in axial capacity of energy piles in unsaturated silt tested by Goode and McCartney (2015) were related to changes in degree of saturation. The objective of this study is to simulate the transient response of an unsaturated soil layer surrounding a solitary energy pile during heating and cooling cycles. Specifically, simulations of coupled heat transfer and water flow are used to understand the evolution in the degree of saturation around the pile.

Background
Many models have been developed to study nonisothermal multiphase flow in unsaturated soil. Philip and de Vries (1957) were the first to provide a comprehensive model to predict the nonisothermal fluid movement in porous materials under combined suction and temperature gradients. Most models for coupled heat transfer and water flow are based on their approach that considers enhanced vapor diffusion and equilibrium phase change. Their formulation was improved by de Vries (1958), who combined water vapor and liquid flow in the mass balance equation. Cass et al. (1984) included the temperature effect on the thermal conductivity of the soil and corrected the vapor diffusive flux by adding a phenomenological enhancement factor for vapor diffusion. Thomas and He (1995) accounted for the effects of elastic volume changes that may occur due to thermal expansion and wetting/drying effects. Other studies focused on the coupled effects of temperature and degree of saturation on the thermal and hydraulic properties of unsaturated soils and their effects on the coupled heat and fluid flow in unsaturated soils (e.g., Campbell et al. 1985

Numerical Model
In this study the model of Smits et al. (2011) is extended to investigate the behaviour of unsaturated soils under the boundary conditions associated with heat exchanger fluid flow through a closed-loop pipe loop in an energy pile. The model governing equations, material properties, and model geometry are given in the following sections.

Governing equations
The momentum and continuity equations for fluid flow in a pipe (Barnard et al. 1966) are given by: and where is the density of a glycol-water mixture; is the averaged velocity across the cross-section of a heat exchanger pipe; is time; is the fluid pressure in the pipe; is the Darcy friction factor; = 4 / is the mean hydraulic diameter; A is the pipe cross section area; and is the wetted perimeter. The energy equation for an incompressible fluid flowing in a pipe is (Lurie 2008): where is the heat capacity of water; is the temperature; is the thermal conductivity of the glycolwater; and is a heat source. The nonisothermal water mass balance equation in porous medium (Bear 1972 where is the porosity; is the degree of water saturation; is the density of water as function of temperature (Hillel 1980); is the capillary pressure (Grant and Salehzad 1996); is the temperaturedependent relative permeability function for water; is the intrinsic permeability; is the water dynamic viscosity as function of temperature; is the pore water pressure, is the acceleration due to gravity; and is the phase change rate between water and vapor due to evaporation or condensation and it is calculated as (Bixler 1985 where is a fitting parameter; is the universal gas constant; is the water molecular weight; is the equilibrium vapor density; is the vapor density; and is the residual degree of saturation as a function of temperature given by where is the degree of gas saturation; is the temperature-and pressure-dependent density of gas; is the relative permeability function for gas in soil; and is the gas dynamic viscosity is the pore gas pressure. The mass balance equation of water vapor is given by (Smits et al. 2011): where is the mass fraction of water vapor in the gas phase; and = is the effective vapor diffusion coefficient; is the diffusion coefficient of water vapor in air, and = ⁄ ⁄ is the tortuosity (Millington and Quirk 1961). The enhancement factor for vapor diffusion is calculated using the form of the model of Cass et al.
where a fitting parameter; and is the fines content. The heat transfer energy balance equations for porous media containing water and gas (Whitaker 1977;Moradi et al. 2016) is given by: where is the total density of soil; is the heat capacity of soil; is the heat capacity of gas; is the thermal conductivity of soil; , are the velocity of water and gas, respectively; is the latent heat of vaporization (Monteith and Unworth 1990); and is a heat source.

Material Properties
The energy pile investigated in this study is a concrete pile with embedded high-density polyethylene (HDPE) heat exchanger tubing filled with glycol-water. The concrete pile has a density of 2400 kg/m 3 , a heat capacity of 960 J/(kg·K), and a thermal conductivity of 1.4 W/(m·K). The HDPE heat exchanger tubing has a thermal conductivity of 0.48 W/(m K). The glycol-water mixture has a density of 1008 kg/m 3 , a heat capacity of 3267 J/(kg·K), and a thermal conductivity of 0.58 W/(m·K).
The soil profile investigated in this study consists of Bonny silt, which has a specific gravity G of 2.65. For a porosity of 0.46 the soil has a saturated hydraulic conductivity of 1.24×10 −7 m/s. The SWRC model of van Genuchten (1980) was used to represent the soil with parameters n = 1.6 and α = 0.2 kPa . The enhancement factor was calculated with = 8. Temperature dependency of the SWRC is considered through the temperature effect on the matric suction defined by Grant and Salehzadeh (1996) and the effects of temperature on the residual saturation defined by She and Sleep (1998) and with constant c = 0.08. The soil thermal conductivity was calculated using a model of Lu and Dong (2015): where and are the thermal conductivity of dry and saturated soil, respectively; is a soil saturation parameter; and is a pore fluid conductivity parameter. The volumetric heat capacity of unsaturated soil was calculated using the following model: where and are the volumetric heat capacity values of dry and saturated soil, respectively.

Model Geometry
The 3D model geometry consists of a single driven concrete pile with three closed-loop U-tubes that is installed in a uniform unsaturated soil layer as shown in

Initial and Boundary Conditions
The soil layer was simulated with initial uniform temperatures of 10, 14 and 18 °C for comparison and is thermally insulated with zero water and air flux at the boundaries. These assumptions were made as the energy pile investigated is beneath an insulated building where there are negligible atmospheric interaction effects. The energy pile boundary conditions were defined based on observations from McCartney and Murphy (2017) to study the effects of heating and cooling cycles on the unsaturated soil layer. A continuous relationship was fitted to the measured inlet heat exchanger fluid temperatures for an energy pile at their site, which were then applied to the simulated energy pile for a duration of three years as shown in Figure 2. The calculated exit temperatures during the heat transfer process are also shown in this figure.

Results
Different simulations were performed for the purpose of understanding the long-term effect of energy piles on an unsaturated soil deposit. The simulations included a sensitivity analysis of the soil for different initial temperatures. Then the model was analyzed for three different initial SWRCs representative of Bonny silt with n = 1.6 and the α values shown in Figure 3 and the three initial soil temperatures of 10, 14 and 18 °C. The different scenarios were analyzed using COMSOL Multiphysics software. A coarse mesh was used at the far boundaries, and a fine mesh was used near the pile boundary where the flow processes are likely to have the greatest effect on soil-structure interaction. Time series of the simulation results in terms of the changes in soil temperatures and degrees of saturation at different radial distances from the pile at a depth of 5 m from the surface are shown in Figures 4 and 5. Although the soil temperatures show a sinusoidal trend with time, the degree of saturation shows a slight downward inclined trend reflecting gradual permanent drying similar to that observed by Başer et al. (2018), who applied much higher temperatures than investigated in this study and observed greater drying. The results indicate that the changes in the degree of saturation relate to changes in soil temperature, the shape of the SWRC via the α parameter, and the duration over which the thermal load is applied. Although the degree of saturation gradually decreases with time, it had a maximum change of about 4.3% for the case with the lowest initial degree of saturation (α = 0.8 kPa ). The initial soil temperature on the other hand had an inverse correlation with the change of the degree of saturation. For α = 0.2 kPa the change in the degree of saturation was the highest (3.7%) for the case where T = 10 °C. The maximum changes (during the summer season) for all the cases occurred after three years, reflecting an accumulated drying of the soil. Due to the similar durations of heating and cooling, the pore water will condense, and the degree of saturation will nearly recover its initial state after each heating-cooling. However, in the long run, the evaporation rate will become higher than the condensation rate and as the soil temperature at the interface will increase, the soil will gradually dry. The degree of saturation could reach lower values for longer thermal loading durations, higher inlet flow rates and temperatures, and different pile spacings. The initial degree of saturation had a key effect on the changes in temperature of the soil during energy pile operation. The soil experienced maximum changes in temperature of ΔT = 12.7, 13.3 and 13.7 °C for α = 0.2, 0.4, and 0.8 kPa , respectively, as shown in Figure 4. This occurs because the higher α value leads to a lower initial degree of saturation and based on the model of Lu and Dong (2015), decreases in degree of saturation will lead to a corresponding decrease in thermal conductivity. This means that the temperature of the soil will increase by a greater amount for initially dryer soil.
At a radial distance of 1 m from the pile, the change in the soil temperature and the degree of saturation will reach its peak one month after the value at the soilstructure interface reaches its peak. The timing of the temperature changes seems to be closely linked with the degree of saturation changes, despite the different thermal and hydraulic properties governing heat transfer and water flow. The factors affecting the soil temperature are the heat capacity and thermal conductivity which depend on degree of saturation (Equations 9 and 10). The main factors that affect degree of saturation changes include the hydraulic conductivity and the enhanced vapor diffusion parameter a (Equation 7).
The annual extreme soil temperature profiles at distances of 0.05 m from the energy pile for T = 10 and 14 °C are shown in Figure 6. The temperature profiles indicate that the soil temperature will decrease in winter if the initial soil temperature is much higher than the minimum inlet temperature. However, this decrease in the temperature will be less in the following years even when applying the same low inlet temperature. This can be observed at the depths just below the pile where the temperature remained high during cooling. The distributions in degree of saturation at distances of 0.05 m from the energy pile for Tinitial = 14 °C with α = 0.2 and 0.8 kPa -1 are shown in in Figure 7. The degree of saturation remains constant during the first winter then it will slightly decrease the following year. The decrease in the degree of saturation may be higher if surface evaporation was permitted as a boundary condition.

Conclusion
A model formulation for coupled thermal and fluids transport in unsaturated soil was used to analyze the longterm effects of thermal heating and cooling cycles on unsaturated soil for a single energy pile. The results indicate that the cyclic heating and cooling will gradually dry the soil surrounding the energy pile. The magnitude of drying is not as significant as that observed in previous studies due to the lower temperature ranges experienced in energy piles. However, the cumulative drying of soil will lead to more gas flow in soil and so more water evaporation and will increase the soil temperature at larger radial distances from the pile. Because of the decrease of the degree of saturation, the soil layer thermal conductivity will decrease. In the long term, the soil surrounding the pile will have greater heat storage which will lead to a decrease in heat transfer from the energy pile to the surrounding soil and better pile performance. An important practical implication of the results is that measurement of temperature or water content further from the pile interface may not be representative of those at the soil-structure interface. The interface values will likely have the greatest impact on soil-structure interaction. Further study is needed to understand the effects of longer durations of seasonal heating and cooling cycles, different boundary conditions including the ambient temperature effects, the soil surface water evaporation, the water table level, and applying uneven spans of thermal loads.