Numerical investigation of the performance of group of geothermal energy piles in unsaturated sand

The use of foundation structures (piles) coupled to a heat pump system, commonly referred to as geothermal energy pile (GEP) system, provides a renewable energy solution of achieving space heating and cooling in buildings; whilst also being utilised for the structural stability of the overlying structures. The system operates by exchanging the low-grade heat energy within the shallow earth surface with the building, via the circulation of heat carrier fluid enclosed in a high-density polyethylene plastic pipes. In summer, heat energy is extracted from the building and transferred into the ground to achieve space cooling. While in winter, the ground heat energy is harnessed and transferred to the building to achieve sustainable space heating. This paper investigates the thermal performance of group of GEP system under the effects of different initial soil pore water content. Through the five-year simulation’s period, it was found that the increase in soil pore water content decreases the possibility of thermal interaction between the GEPs in the group. Also, it was observed that the trend in maximum temperature witnessed within the soil domain decreases nonlinearly during the five years period.


Introduction
Geothermal energy pile (GEP) system comprises the process of coupling ground heat exchangers with heat pump units to harness the low-grade heat energy within the shallow earth surface for the purpose of space heating and/or cooling in buildings. The system installation involves incorporating ground energy loops, made from high density polyethylene (HDPE) plastic, within the foundation element(s) often referred to as ground heat exchangers or pile heat exchangers. Inside the HDPE pipes, heat carrier fluid (HCF) is circulated to exchange heat energy with the ground. In winter, the system extracts heat from the ground and transfers it to the building for space heating, whereas the process is reversed in summer to achieve space cooling. The working principle of the system relies on the difference between the ambient air temperature and the ground temperature to transfer the energy between the two mediums. Thus, to ensure the longevity and sustainability of the system, researchers advised that the system should not be subjected to extreme operating conditions such as extreme monotonic or monodirectional thermal loading [1,2]. Conversely, operating the system for bidirectional use minimises or negates the danger of subjecting the system to extreme loading conditions. For example, where the heating and cooling thermal demand of the building are equal, the heat energy harvested in winter for heating is transferred back to the ground for storage whilst achieving space-cooling operation in summer. However, in a situation where the heating or cooling demand is greater, other means of heating and cooling methods should be used to supplement for the resultant overall energy demand [2]. This ensures that excessive ground temperature deficit, due to the space heating operation; or ground temperature build-up, owing to the corresponding space cooling operation is avoided to safeguard the efficiency of the heat exchanger unit, and consequently the GEP system. Factors such as soil initial temperature, soil thermal conductivity and diffusivity, soil initial moisture content, groundwater flow, soil type, soil chemical and minerals composition etc. do have a significant influence on preventing the excessive temperature changes within the soil domain surrounding the heat exchanger unit [1][2][3][4][5].
This study presents the results of a series of numerical analyses carried out to investigate the thermal performance of group of heat exchanger units installed in Leighton Buzzard sand. In addition, the initial soil pore water content was varied from fully dry, unsaturated and saturated conditions during the different numerical simulations.

Finite element modelling
The finite element modelling and simulations in this study were carried out in COMPASS (Code for Modelling PArtially Saturated Soils) code. The code was developed in Cardiff, based on mechanistic approach. Its theoretical and numerical formulations are given in various publications including Thomas and Sansom [6] and Thomas et al. [7]. The numerical code is capable of carrying out Thermo-Hydro-Chemo-Mechanical (THCM) analyses in a partially saturated porous media. The current study used the TH capability of the code, which describes the governing equations of the code in terms of primary variables. The TH problem is solved by computing the equations of heat and moisture flow through the use of pore water pressure (ul), pore air pressure (ua) and temperature (T), respectively. The numerical solution is solved using finite element method to achieve spatial discretisation, while implicit finite difference algorithm is used to achieve temporal discretization.

Moisture transfer
Moisture transfer in a porous media occurs in liquid and vapour form. Based on the law of conservation of mass, it can be mathematically expressed as: where n is porosity, Sa and Sl are the degree of saturation of pore air and water. vl, vv and va represent the velocity of liquid, vapour and air. ρv and ρl are the densities of water vapour and liquid water respectively. ∇ is the gradient operator and t is time.
According to Darcy's law, the flow of liquid water through an unsaturated soil media is given by: where kl is the intrinsic permeability, μl is the pore liquid absolute viscosity, Kl is the unsaturated hydraulic conductivity, γl is the liquid unit weight and z is the elevation. The transfer of vapour occurs due to diffusive and pressure flows. The diffusive flow may be solved using the expression proposed by Philip and De Vries [8]; Ewen and Thomas [9] and extended by Cleall et al. [10,11], given as: where Datms is the molecular diffusivity of vapour through air, vv is mass flow factor. (∇Τ)a/∇Τ is the microscopic pore temperature gradient factor, h is the relative humidity, S is the suction and ρo is the saturated vapour density.

Dry air transfer
In an unsaturated porous media, the dry air can be expressed as bulk and dissolved air. The latter is driven by the gradient of air pressure, and can be determined using Darcy's law. While, dissolved air is transported via advection along with the pore liquid. The proportion of dry air contained within the pore liquid can be defined using Henry's law.
The law of conservation of mass dictates that the temporal derivative of the dry air content is equal to the spatial derivative of the dry air flux, given as: where θa and θl are the volumetric air and liquid content respectively, Hs is Henry's volumetric coefficient of solubility, ρda is the density of dry air and ∂V is the incremental volume.

Heat transfer
Heat transfer in soil occurs via conduction, convection, and radiation. Heat is added due to phase change from liquid to vapour as the latent heat of vaporization. In this study, the radiation effect was neglected, because its contribution towards heat transfer is insignificant in soils with smaller grain size [12,13].
The law of conservation of energy for heat dictates that the temporal derivative of the heat content, Ω, is equal to the spatial derivative of the heat flux, Q, as: The heat flux per unit area, Q, is defined as: where, λT is the unsaturated soil thermal conductivity coefficient, Cps, Cpl, Cpv and Cpda are the specific heat capacities of solid, liquid, vapour and dry air respectively, L is the latent heat of vaporisation, ρs is the density of solid particles and T is temperature.

Model description
The modelling approach was based on the scenario where piles exist in a grid of 7 by 7 group of piles (shown in Fig  1a), spaced at 7 m in both longitudinal and lateral directions. In total, 49 piles exists in the group with each having a diameter of 600 mm, and a length of 30 m. As regards to the transitioning from the 3D problem into 2D axisymmetric, the findings of [14,15] were adopted and used. They reported that heat predominantly disperses in the radial plane of a GEP. Thus, to minimise the complexity of the modelling process and study the groups' thermal behaviour, one of the 7 grids was chosen for the modelling process (Fig. 1a). Likewise, since equal amount of heat is expected to radiate from each pile, half of the grid can be modelled into a quasi 2D axisymmetric, which serves as a representative of the whole grid.
During the modelling process in COMPASS (Fig. 1b), only the concrete and the soil domain were modelled, i.e. the HDPE pipes and HCF were neglected. Most importantly, neglecting the effect of the HCF and HDPE in the geometry reduces the total number of required elements in the model, and consequently, minimising the computational time required for the numerical analyses.
The GEPs are installed in a soil domain of 40 m wide and 50 m in height. The domain was made large enough to ensure there were no boundary effects. Also, that all thermal evolution and temperature changes, that might arise due to the cyclic heating and cooling, occurs within the soil domain. The GEPs and the soil domain were discretized in COMPASS using 8-node quadrilateral elements. The elements sizes increases from 0.3 m at the location of the GEPs to about 1 m at the farthest boundaries from the group of piles. The model is characterised by having 2350 total number of elements and 2448 total number of nodes.

Initial & boundary conditions
The initial temperature used for the FEM simulations corresponds to 286.55 K (13.4°C), a value measured during a thermal response test in East London as reported by Loveridge et al. [16]. Furthermore, the initial degree of saturation (Sl) was varied from 0, 30, 60 and 100% respectively. The initial suction (s) corresponding to the respective Sl percentages were determined using the soil water characteristic curve fitting proposed by van Genuchten [17]. The values for the initial temperature and suction conditions are given in Table 1.
In addition, heat flux boundary condition of 25 (W/m2) was applied at the perimeter of the GEPs. The heat flux value was back-calculated for a GEP of 30 m length, based on the data for heat injection rate reported by Gawecka et al. [18]. Similarly, the sides, the top and bottom of the domain were considered as adiabatic and impermeable boundaries to ensure that no heat and moisture losses occurs outside of the domain.

Material parameters
The materials used for the modelling comprises Leighton Buzzard sand and concrete. Their thermal properties, and the hydraulic properties of the sand were measured and presented in Table 2. The detailed description of the different laboratory tests carried out on the materials can be found in Sani [2].

Thermal conductivity
The coefficient of thermal conductivity of unsaturated soils (λT) is described as a function of the degree of saturation (Sl), mathematically expressed as: The variation in the λT of the sand was computed using the relationship obtained from Ewen and Thomas [9] shown in Equation (8) where θ is the volumetric moisture content. The λT values ranges between 0.35 to 2.85 W(mK) -1 for the Sl values between 0 to 100%. The full curve for the λT relationship obtained using Equation (8) compared to the lab results of Sani [2] is shown in Fig. 2.

Soil water retention curve (SWRC) and unsaturated hydraulic conductivity (Kl)
The soil water characteristic curve provides a relationship between the volumetric/gravimetric water content of soils and suction. In the current work, van Genuchten curve fitting technique was used to establish the SWRC of the sand, given in Equation (9). and, where θl is the volumetric liquid water content, θlr is the residual volumetric liquid water content, θls is the saturated volumetric liquid water content (i.e. taken equal to porosity), S is the matric suction, ϕ, ε and φ are fitting parameters. In addition, Melhuish [19] reported that the unsaturated hydraulic conductivity, Kl, is a function of the void ratio (e), Sl and the temperature (T), expressed as: where ksat is the saturated hydraulic conductivity (m/s), δ is a parameter ranging between 3 and 10 [20].

Numerical simulations
Firstly, a numerical validation was carried out prior to using the numerical model for the analyses. The detailed description of the results for the validation is reported in Sani [2]. In addition, transient thermo-hydraulic numerical simulations were carried out by applying cyclic heating and cooling loads of +25 and -25 W/m 2 , respectively to the surfaces of the GEPs (Fig. 3). In total, the simulations were conducted for 5 years i.e. five thermal cycles. Within each thermal cycle, a heating load was applied for 6 months followed by superimposing cooling loads on to the GEPs for 6 months. In addition, the effect of the variation of the soil initial pore water content was also investigated by varying the initial degree of saturation between 0, 30, 60 and 100%Sl. The temperature at the observed points increases uniformly until it reached a maximum value of 55.5 and 60°C at points A and B respectively for the sand with 0% Sl. As the degree of saturation (Sl) increases from 0% to 100%, the maximum temperature observed at points A and B decreases to temperature values of 39.6 and 43.7°C respectively. Signifying a temperature drop of about 16 and 16.4°C when the results witnessed at the end of the heating cycles for the 0% Sl and 100% Sl were compared. In addition, the maximum temperature witnessed at these points occurred at the end of the 1 st heating cycle and decreases nonlinearly in the subsequent heating cycles. The temperature trend declines at a rate of about 3.5°C/year, to a value of 36.9 and 40.9°C at A and B, at the end of the 5 years cyclic heating and cooling period of the numerical simulations.
However, during the cooling phase of the simulations, maximum temperature magnitudes of about 11.3 and 6.2°C were observed at points A and B respectively after imposing a heating power of -25 W/m 2 to the surfaces of the GEPs embedded in dry soil. Similarly, as the soil initial saturation increases to 100%, the temperature changes imposed by the application of the cooling load decreases at these to 14.3 and 10.9°C, respectively. Similarly, as the pore water increases from 0 to 100%, the temperature magnitude at the observed points decreases over the 5 year cycle. In both the heating and cooling simulations, greater temperature magnitudes were observed in dry soils, with the effect decreasing with increase soil pore water content. This is largely owing to the contribution of the pore water towards heat dissipation. As the water content rises from 0% to 20% Sl, the thermal conductivity significantly increases from 0.35 W(mK) -1 to about 2.4 W(mK) -1 . Thereby, preventing excessive temperature build-up or deficit in the soil during heat injection or extraction process.
In addition, Fig. 5 presents the results of radial temperature distribution with normalised distance at the end of the 5 th heating and cooling cycles. The results were obtained at the pile mid-depth from the pile labelled no. 4 (GEP-4) to farthest boundary to the right; located at 21.9 m from the group of GEPs. Maximum and minimum temperature magnitudes of 32 and -2°C were observed at the surface of GEP-4 at the end of the heating and cooling simulations. The temperature decreases linearly with distance away from the pile group. The maximum temperature changes were observed in the dry soil, with the effect decreasing with increase in soil pore water. Furthermore, the region that was influenced by the effect of the cyclic thermal load application was found to be about 15 times the pile diameter.

Conclusion
This paper presented the results of numerical simulations investigating the behaviour of group of GEPs, and how that is influenced by the increase in soil pore water content. It was found that: • the increase in soil pore water content significantly increases the heat transfer in soils, thereby preventing extreme heat build-up during heat injection or heat deficit during heat extraction process. In the context of GEP systems, this ensures that the ground temperature is not excessively disturbed, beyond recoverable limits, thus ensuring the longevity of the system. • the continuous heating and cooling cycles of the GEPs, during the 5 years of numerical simulations, resulted in continuous decrease in temperature trend at the observed locations (i.e. A and B). This could be as a result of the continuous heat injection and extraction mode without allowing the system any opportunity to recover throughout the period of the numerical simulations. However, adopting daily transient heating and cooling mode of operation permits the soil domain to naturally recover and prevent excessive heat build-up or deficit in the soil • lastly, this paper provides an understanding of the behaviour of group of GEPs when subjected to thermal cycles. It provides information on the extent of the soil domain that is influenced by the heating/cooling operation of GEPs. Thereby, can serve as a guide to ensure piles that are thermally activated are situated at a certain distance to prevent GEP thermal interaction. However, it should be noted that the findings here are not applicable where groundwater flow exists.