Long-Term Thermal Performance of Group of Energy Piles in Unsaturated Soils under Cyclic Thermal Loading

: Geothermal energy piles (GEPs) are an environmentally friendly heat exchange technology that dualizes the role of the structural foundation pile for load support and in meeting the building heating/cooling need. Energy loops made from high-density polyethylene, which allow heat carrier ﬂuid circulation, are ﬁtted into the pile foundation elements to extract or inject and store heat energy in the soil surrounding the pile. This paper reports the results of a numerical study investigating the long-term behaviour of a group of energy piles embedded in unsaturated soils (sand and clay) under continuous cyclic heating and cooling load. Additionally, two scenarios were investigated where: (1) the whole GEPs were heated and cooled collectively; (2) alternate piles were heated and cooled. It was found that the trend of temperature magnitude at all the observed locations decreases with time as a result of the continuous heating and cooling cycles. Furthermore, subjecting alternate GEPs to the heating and cooling cycles result in lower temperature development in comparison to thermally activating all the GEPs in the group. This is attributed to the applied thermal load, which is 0.5 times that considered in the ﬁrst case. However, this might not be the case where equal thermal load is applied on the GEPs in the two cases investigated. to understand the extent of the soil domain that could be inﬂuenced by the cyclic heating and cooling of the group of GEPs under different degrees of saturation condition.


Introduction
The use of concrete pile foundation elements coupled with a heat pump unit has proven to be an environmentally effective and sustainable approach towards achieving space heating and cooling demand of the overlying structure. The coupling process uses energy loops made from high-density poly-ethylene (HDPE) plastic pipes incorporated into the pile to connect the heat pump unit with the foundation element. Within the pipes, a heat carrier fluid (HCF) is circulated to exchange heat energy with the shallow earth surface. The combined system is often referred to as a ground source heat pump (GSHP) or geothermal energy piles (GEP) system [1,2].
The system operates by extracting low-grade heat energy from the ground and transfers it to the building to achieve space heating in winter. However, in summer, the process is reversed, i.e., heat is removed from the building and injected into the ground for storage whilst achieving space cooling.
This cyclic heating and cooling process alters the mechanical behaviour of a typical structural foundation pile by imposing additional thermal stresses and strains on the pile and its surrounding soil. Nevertheless, several studies have shown that these effects behave in a thermoelastic manner; meaning, they are reversible in nature [3,4]. Similarly, the type of soil surrounding the pile and the restraint at the pile-head/toe have a significant effect on the axial load induced in the pile [5][6][7].
Equally, the cyclic heating and cooling process has an influence on the thermal performance of the system. In situations where the heating and cooling demand are equal, the system transfers back the heat energy harvested in winter while achieving space cooling in summer. However, in an unbalanced system, other approaches should be used to ensure the heat energy lost in winter is recovered back [4,[7][8][9] for use in the subsequent heating season.
Several studies were carried out investigating the usage of the soil domain for the purpose of heat storage. Hesaraki et al. [10] reported that utilising a collective regional development (comprising of piles from several housing units) for the purpose of seasonal heat storage results in higher efficiency in comparison to that of a single residential unit. Ochsner et al. [11] reported that factors such as reference storage depth and heat capacity of the storage domain should be accurately determined to maximise the thermal storage capacity of the soil domain. Mccartney and Baser [12] have shown that the presence of vapour in the soil pores has a positive impact on the heat storage efficiency of the system. However, in the situation where GEPs are installed in unsaturated soils, extreme heat injection was shown to drive moisture away from the location of the pile due to temperature gradient. This causes the reduction of the soil moisture content in the region where higher temperature exists, i.e., pile surface [13]. As a result, the thermal conductivity and the heat exchange rate of the soil next to the pile reduces. Eventually, thereby causing temperature build-up in that region owing to the heat injection process. However, this can be positively beneficial because the injected/stored heat can be harnessed for heating use in winter times.
Dupray et al. [14] reported a study investigating the long-term performance of a group of GEPs under cyclic thermal loading. They found that the increase in the heat injection rate and the injection and extraction ratio has an influence on the overall efficiency of the global energy storage system.
Rees and Van Lysebetten [15] developed a numerical model to investigate the longterm thermal response of a group of GEPs using the Dynamic Thermal Network (DTN) approach. Their choice for DTN was because of its suitability to represent complex geometries and heterogenous thermal properties. The validation exercise was carried out using the experimental field results of [16,17]. They reported that it is crucial to account for the differences in thermal load conditions over a long term, particularly depending on the position of the GEP in the array and the effect of the long-term heat injection into the soil. These have implications on the maximum HCF temperatures for a given energy demand, thereby influencing the long-term energy efficiency and sustainability of the system.
To further understand the sustainability of a group of GEPs over a long-term period, Sutman et al. [18] conducted a life cycle analysis (LCA) of a group of energy piles under three different climatic conditions. The results obtained using a GEP heating scheme were compared with that of a conventional heating and cooling system. The LCA looked at the effect of each system on climate change, resource consumption, human health, and ecosystem quality. Their findings show that the use of the GEP scheme indicated a reduction in environmental impacts in the majority of the cases considered.
Ferrantelli et al. [19] investigated the performance of a group of energy piles for a commercial hall-type building in a cold climate. They reported that a nonlinear relationship exists between the dimensioning of the pile group (length and spacing of the GEPs) and the heat pump extraction power. However, among the limitations of their study, the findings could only be applicable to a commercial hall-type building and with GEPs installed in clay.
The arrangement of GEPs in various configurations-ranging from a single GEP unit, 2 × 2, 3 × 3, 4 × 4, and 5 × 5 rectangular grids-on the long-term performance of the system was reported by Olgun et al. [20]. In all the cases, the temperature magnitude at the GEP centre, GEP wall, 1 2 -diameter away from the GEP wall, and 1-diameter away from the GEP wall significantly increases in the 5 to 8 years of the simulations. Afterwards, the temperature increase becomes insignificant depending on the climatic conditions employed during the simulations. However, Abdelaziz [21] reported that the temperature magnitude in the GEP (pile centre and periphery) and in the soil domain ( 1 2 -and 1-diameter) becomes negligible after the first 2 years, and thus, the thermal efficiency of the system is equally affected.
Hepburn [22] showed that a period of 3 years was sufficient for a horizontally laid ground loop GSHP scheme to reach a steady state. This finding is useful considering that horizontal ground loop systems are laid at a shallower depth (1.5-2.0 m) [23], thereby making them more susceptible to atmospheric interference. However, GEPs are less likely to be greatly influenced by the ambient temperature because they are installed at greater depth and are often covered by the building structure. Thus, based on the findings from literature, a period of 2-8 years can be considered sufficient for the thermal changes to reach a steady state when investigating the performance of a group of GEPs.
Furthermore, one of the main issues associated with numerical investigations is related to the strenuous effort required in simulating complex and intricate models. Thus, to reduce the complexity of the modelling process, Abdelaziz [21] modelled a quarter of a pile group. The whole grid comprises 9 energy piles (3 × 3 in a grid). This greatly decreases the computation time required for simulating the whole pile grid. In addition, Abdelaziz [24] showed that the 2D modelling approach of GEPs was found to be adequate provided the top ground surface is insulated and that the approach can effectively account for the heat transfer through walls of the energy loops, in the pile material (concrete), and within the soil [21]. Similarly, via the use of 2D and 3D approaches, Ferrantelli [25] showed that the HCF and energy loops could be neglected.
All the above-cited studies were carried out by superimposing either heating load (heat injection), cooling load (heat extraction), or cyclic heating and cooling loads on all the GEPs in the group. However, none of the studies investigated thermally activating some selected piles in the group. Similarly, it is highly likely for the GEPs to be installed in soil that is either dry, unsaturated, or saturated conditions depending on the depth of pile installation and the water table depth.
Hence, this paper investigates, via the use of finite element modelling (FEM), the longterm thermo-hydraulic (TH) behaviour of a group of GEPs under cyclic heat injection and extraction process. Two scenarios are studied, namely: (1) the whole piles in the group are thermally activated; (2) where alternate piles are thermally activated in the group in order to reduce the applied thermal load by half. This is especially important where excessive temperature development is anticipated during the operation of the system. Similarly, the behaviour of the group of GEPs in sand and clay soils is investigated and how they vary with an increase in soil moisture content. Findings from this investigation will be useful in making informed decisions during the preliminary and actual design stages of the GEP system. However, the findings here are limited to the 2D modelling approach and provides an opportunity for future studies in this area.

Finite Element Modelling
The finite element modelling and analyses were carried out in COMPASS (COde for Modelling PArtially Saturated Soils). It has the capability of numerically solving Thermo-Hydraulic-Mechanical and Chemical (THM-C) processes in a partially saturated porous media. The extensive details of its theoretical formulations are given in various publications [26][27][28][29] and dissertations [30][31][32][33].
To achieve the aim of this work, the TH capability of the COMPASS code was used. The governing equations for solving the TH problem in unsaturated soil media are reported here in terms of primary variables. The equations solve for the coupled heat and moisture transfer processes via the use of pore water pressure (u l ), pore air pressure (u a ) and temperature (T), respectively. The numerical solution for the TH problem is solved via the use of the finite element approach to achieve spatial discretization, whilst temporal discretization is achieved by the use of an implicit finite difference algorithm. The detailed equations are reported in Ewen and Thomas [34], Thomas and Sansom [29] and Thomas and Li [28].

Soil Moisture Transfer Mechanism
The transfer of moisture in unsaturated soil domain occurs in two forms: liquid water and vapour water. According to the law of conservation of mass, it can be expressed mathematically as: where n is the soil porosity, S a and S l are the degree of saturation of pore air and pore water, respectively. The terms v l , v v, and v a are the velocities of liquid, vapour, and air, respectively. The terms ρ v and ρ l represent the densities of water vapour and liquid water, respectively. ∇ is gradient operator, and t is time.
The flow of liquid water through an unsaturated media is defined mathematically using Darcy's law, expressed as: where k l is the intrinsic permeability, µ l is the absolute viscosity of pore liquid, K l is the unsaturated hydraulic conductivity, γ l is the unit weight of liquid, and z is the elevation. The vapour transfer occurs due to diffusive and pressure flows. The diffusive flow may be solved using the expression proposed by Philip and De Vries [35] and Ewen and Thomas [34] and extended by Cleall et al. [36,37], given as: where D atms is the molecular diffusivity of vapour through air, and v v is a mass flow factor. (∇T) a /∇T is the microscopic pore temperature gradient factor, h is the relative humidity, s is the suction, ρ o is the saturated vapour density, and u a is the pore air pressure.

Dry Air Transfer Mechanism
The dry air present in an unsaturated porous soil media can exist in two forms, namely bulk and dissolved air. The latter is driven by the gradient of air pressure and can be determined using Darcy's law. Meanwhile, dissolved air is transported advectively with the pore liquid. The proportion of dry air contained within the pore liquid can be determined using Henry's law.
In addition, 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, mathematically defined as: where θ a and θ l are the volumetric air and liquid content, respectively, H s is Henry's volumetric coefficient of solubility, ρ da is the density of dry air, and ∂V is the incremental volume.

Soil Heat Transfer Mechanism
Heat transfer in soils occurs via conduction, convection, and radiation. Equally, heat is added due to phase change from liquid to vapour as latent heat of vaporization. However, heat transfer via radiation in soils is insignificant due to the smaller soil pore sizes and how close the soil skeletons are packed next to each other. Farouki [38] reported that the radiation effect accounts for only 5% of the coupled heat transfer process in gravels with a particle size of 20 mm. Thus, the contribution of heat transfer via radiation is neglected here. The law of conservation of energy for heat flow dictates that the temporal derivative of the heat content, Ω, is equal to the spatial derivative of the heat flux, Q. This is mathematically defined as: The heat flux per unit area, Q, is defined as: where λ s is the coefficient of thermal conductivity of unsaturated soil; C ps , C pl , C pv , and C pda are the specific heat capacities of solid particles, liquid, vapour, and dry air, respectively; L is the latent heat of vaporisation; ρ s is the density of solid particles; T i is the reference temperature in Kelvin.

Geometry Description
The modelling exercise was carried out considering the geometry of a group of GEPs shown in Figure 1a comprised of 49 piles. The piles are characterised with a diameter of 600 mm and a length of 30 m, spaced at 7 m apart and arranged in a 7 × 7 grid. The geometry is further simplified by considering one of the seven grids, as shown in Figure 1b. Although this might not be a true representative of the whole group of 49 piles, it has, however, reduced the complexity of the modelling exercise. Equally, the results obtained from the 2D model will pave the way and be useful for future research of the full 3D domain.
During the model development, half of the GEPs shown in Figure 1b were considered and used for the modelling. Thus, a quasi-2D axisymmetric model was developed and set up in COMPASS (shown in Figure 1c). In the setup, only the concrete and the soil domain were modelled, i.e., the HDPE pipes and HCF were neglected. Hence, the heat flux at the perimeter of the pile can be easily calculated, which depends on factors, including thermal conductivity of HCF, HDPE, concrete, HCF flowrate, HCF type, soil's thermal and hydraulic properties, GEP diameter and length, etc.
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.
Additionally, in the soil region, a domain size of 45 m wide and 50 m in height was chosen for the group of GEPs to be installed. The domain was made to be large enough to ensure there were no boundary effects. In addition, all thermal evolution and temperature changes, which might arise due to cyclic heating and cooling, occurs within the soil domain.
To discretize the GEPs and soil domain, several mesh sizes were used to carry out a sensitivity analysis in order to arrive at an optimum mesh size. Numerous cases, A-G, with a different total number of elements ranging from 2350 to 24,320 investigated are indicated in Table 1. Similar to the Table, the elements' sizes ranging from 0.05 m to 1 m were used to discretize the GEPs. Conversely, in the case of the soil domain, the elements' sizes vary along the radial distance and along the depth and range between 0.2 m to 1 m, as shown in the Table. The sensitivity analysis was conducted by applying 25 W/m 2 to the surfaces of the GEPs for a duration of 10 days of continuous heat injection. The results of temperature evolution with time were obtained at mid-depth of the GEP-soil interface (i.e., 0.3 m from the line of symmetry and 15 m from top surface) and shown in Figure 2a. In addition, Figure 2b shows the result of temperature magnitude observed at the end of the heat injection period for the seven cases investigated. Case G was chosen as the optimum mesh size due to very satisfactory temperature results and fewer number of elements, hence, minimising the computation time in comparison to cases C, D, E, and F. Henceforth, case G mesh information was employed and used to discretized the compass model. 600 mm and a length of 30 m, spaced at 7 m apart and arranged in a 7 × 7 grid. The geom etry is further simplified by considering one of the seven grids, as shown in Figure 1 Although this might not be a true representative of the whole group of 49 piles, it ha however, reduced the complexity of the modelling exercise. Equally, the results obtain from the 2D model will pave the way and be useful for future research of the full 3 domain. During the model development, half of the GEPs shown in Figure 1b were considered and used for the modelling. Thus, a quasi-2D axisymmetric model was developed and set up in COMPASS (shown in Figure 1c). In the setup, only the concrete and the soil domain were modelled, i.e., the HDPE pipes and HCF were neglected. Hence, the heat flux at the perimeter of the pile can be easily calculated, which depends on factors, including thermal conductivity of HCF, HDPE, concrete, HCF flowrate, HCF type, soil's thermal and hydraulic properties, GEP diameter and length, etc.
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.
Additionally, in the soil region, a domain size of 45 m wide and 50 m in height was chosen for the group of GEPs to be installed. The domain was made to be large enough to ensure there were no boundary effects. In addition, all thermal evolution and temperature changes, which might arise due to cyclic heating and cooling, occurs within the soil domain.
To discretize the GEPs and soil domain, several mesh sizes were used to carry out a sensitivity analysis in order to arrive at an optimum mesh size. Numerous cases, A-G, with a different total number of elements ranging from 2350 to 24,320 investigated are  was reported by Loveridge et al. [39]. 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 [40]. The values for the initial temperature and suction conditions are given in Table 2.

Initial Conditions
The initial temperature adopted and used for the FEM simulations corresponds to 286.55 K (13.4 • C). A value measured during a thermal response test in East London and was reported by Loveridge et al. [39].
Furthermore, the initial degree of saturation (S l ) was varied from 0%, 30%, 60%, and 100%, respectively. The initial suction (s) corresponding to the respective S l percentages were determined using the soil-water characteristic curve fitting proposed by van Genuchten [40]. The values for the initial temperature and suction conditions are given in Table 2.

Boundary Conditions
A heat flux boundary condition of 25 (W/m 2 ) was applied at the perimeter of the GEPs. This corresponds to the axial distance of 0.3 m from the GEP central axis. The heat flux, q, was determined using the data for heat injection rate reported by Gawecka et al. [41]. The q value was back-calculated for the GEP of 30 m length used in this study. In addition, the sides, top and bottom of the domain were considered as adiabatic and impermeable boundaries to ensure that no heat and moisture losses occur outside of the domain.

Material Parameters
The different materials used for carrying out the numerical analyses comprise concrete and soils (sand and clay). The properties of the concrete and Leighton Buzzard sand were measured, and the type of tests conducted are given in Sani [7]. However, for the Speswhite kaolin clay soil, its properties were adopted from the work of Singh [33]. Parameters, including saturated hydraulic conductivity (k sat ), unsaturated hydraulic conductivity (K l ), degree of saturation (S l ), thermal conductivity (λ), soil specific heat capacity (C ps ), the specific heat capacity of vapour (C pv ), the specific heat capacity of water (C pw ), dry density (ρ d ), porosity (n), latent heat of vaporisation (L), and Henry's volumetric coefficient of solubility (H s ), are given in Table 3.
van Genuchten [40] fitting parameters The thermal conductivity coefficient (λ s ) is an important parameter for estimating the thermal properties of a material. It relates the rate at which heat disperses through a material per unit time. In soils, the λ s depends on factors such as soil type, degree of saturation, porosity, soil particle shape, soil granularity distribution, and mineralogy. The λ s can be mathematically expressed as a function of saturation:

Thermal Parameters
To capture the variation in λ s with S l for the sand, the expression obtained from Ewen and Thomas [34] shown in Equation (8) was used.
where θ is the volumetric moisture content. In addition, the expression for the λ s of clay was adopted from the work of Melhuish [42], which was obtained based on the linear interpolation of the laboratory experimental data of Börgesson and Hernelind [43]. The expression is mathematically given as: where λ dry and λ sat are the dry and saturated thermal conductivity of the clay soil. The variation in thermal conductivity for the sand and clay against the degree of saturation is shown in Figure 3. where λdry and λsat are the dry and saturated thermal conductivity of the clay s The variation in thermal conductivity for the sand and clay against the de uration is shown in Figure 3.

Soil-Water Characteristic Curve and Unsaturated Hydraulic Conductivi
The soil-water characteristic curve (SWCC), also referred to as the wate curve, provides a relationship between volumetric/gravimetric soil moisture c suction or water potential. The full range of the SWCC curve is often plotted on mic scale Fredlund and Xing [44]. In the current work, the van Genuchten c technique, given in Equation (10), was used to establish the SWCC of the soils: where θl is the volumetric liquid water content, θlr is the residual volumetric l content, θls is the saturated volumetric liquid water content (i.e., taken equal t

Soil-Water Characteristic Curve and Unsaturated Hydraulic Conductivity
The soil-water characteristic curve (SWCC), also referred to as the water retention curve, provides a relationship between volumetric/gravimetric soil moisture content and suction or water potential. The full range of the SWCC curve is often plotted on a logarithmic scale Fredlund and Xing [44]. In the current work, the van Genuchten curve fitting technique, given in Equation (10), was used to establish the SWCC of the soils: (10) 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 φ, ε, and ϕ are constant fitting parameters. The SWCC parameters for the respective soils are given in Table 3.
The relationship between the soil degree of saturation and the matric suction is shown in Figure 4. The relationship of the SWCC for sand was measured, and the details of the tests are given in Sani [7], whereas that for clay was obtained from the work of Singh [33].
In addition to the SWCC, another important parameter needed for the numerical simulation is the unsaturated and saturated hydraulic conductivity. The saturated hydraulic conductivity was measured, and the details of the tests are also given in Sani [7]. To determine the unsaturated hydraulic conductivity (K l ), the expression proposed by Melhuish [42] was used. It relates the K l as a function of S l , void ratio and temperature: where k sat is the saturated hydraulic conductivity (m/s), and δ is a parameter ranging between 3 and 10 according to Börgesson and Hernelind [43].
[42] was used. It relates the Kl as a function of Sl, void ratio and temperature: where ksat is the saturated hydraulic conductivity (m/s), and δ is a parameter r tween 3 and 10 according to Bӧrgesson and Hernelind [43].

Numerical Simulation
Thermo-hydraulic numerical simulations were carried out by applying cyc and cooling load to the surfaces of the GEPs to investigate the long-term beha group of energy piles in unsaturated soils. The heating and cooling loads were by applying 25 W/m 2 and −25 W/m 2 to the GEPs, respectively. The heating cyc plied for 6 months, followed by 6 months of the cooling cycle. In total, the num ulations of the thermal cycles were carried out for a duration of 10 years for eac and degree of saturation condition. Two cases were investigated during the simulations: (1) where all the GEPs in the group were heated and cooled collec (2) where alternate GEPs were heated.

Case 1: Heating and Cooling of All the GEPs in the Group
In these simulations, transient cyclic heating and cooling loads were app the GEPs in the group, shown in Figure 1. The GEPs were subjected to cyclic h cooling loads, shown in Figure 5, and the changes due to temperature and mois surrounding soil of the GEPs are investigated. The numerical simulations were for a duration of 10 years, i.e., 10 heating and cooling cycles.  (2007)).

Numerical Simulation
Thermo-hydraulic numerical simulations were carried out by applying cyclic heating and cooling load to the surfaces of the GEPs to investigate the long-term behaviour of a group of energy piles in unsaturated soils. The heating and cooling loads were simulated by applying 25 W/m 2 and −25 W/m 2 to the GEPs, respectively. The heating cycle was applied for 6 months, followed by 6 months of the cooling cycle. In total, the numerical simulations of the thermal cycles were carried out for a duration of 10 years for each soil type and degree of saturation condition. Two cases were investigated during the numerical simulations: (1) where all the GEPs in the group were heated and cooled collectively, and (2) where alternate GEPs were heated.

Case 1: Heating and Cooling of All the GEPs in the Group
In these simulations, transient cyclic heating and cooling loads were applied on all the GEPs in the group, shown in Figure 1. The GEPs were subjected to cyclic heating and cooling loads, shown in Figure 5, and the changes due to temperature and moisture in the surrounding soil of the GEPs are investigated. The numerical simulations were conducted for a duration of 10 years, i.e., 10 heating and cooling cycles.

Case 2: Heating and Cooling of Alternate GEPs in the Group
Generally, the spacing between pile foundation elements is chosen based on the structural requirement that provides an economical and safe design for the overlying and surrounding structures. However, in the case of geothermal energy design, the thermal load has no direct control over the spacing. Equally, where excessive temperature development is anticipated during the operation of the system, the thermal load can be reduced by applying half the total thermal load unto the GEPs. Thus, in these simulations, the cyclic thermal load shown in Figure 5 was applied to GEP-2 and GEP-4 (shown in Figure 1) for a duration of 10 years. The changes to temperature and soil moisture as a result of the cyclic thermal loading were investigated for the sand and clay soils with 0-100% initial degree of saturation.

Numerical Simulation
This section discusses the results of the numerical investigations of the long-term thermal performance of a group of energy piles in unsaturated soils. Firstly, a validation exercise had to be carried out prior to the start of the numerical simulations. The purpose of the validation exercise was to check and validate the numerical constants and curve fitting parameters that are required for carrying out the numerical simulations. These parameters are the van Genuchten [40] fitting parameters for defining the hydraulic properties of the Leighton Buzzard sand. The validation exercise was successfully carried out against an experimental study conducted on dry and unsaturated sand.
The test was conducted by circulating water at a temperature of about 35 and 30 • C through the inlet and outlet PVC pipes, respectively, installed in a concrete pile. The test was conducted for a duration of 6 days of a continuous heat injection process. The full details of the experiment setup and test sequence are reported in Sani [7]. The test schematic is shown in Figure 6. In addition, Section 4.1 presents and discusses the results of the numerical analyses in which all the GEPs are subjected to cyclic thermal loading. Moreover, the results of the numerical simulations where the cyclic thermal loads were applied on alternate GEPs are presented and discussed in Section 4.2.  The results of changes in temperature and degree of soil saturation (S l ) with time were recorded throughout the duration of the tests and compared with the results of numerical modelling, as shown in Figures 7-9, respectively. It can be seen that the results of temperature and S l for the dry and unsaturated sands were in very good agreement with the results obtained using COMPASS code. The maximum deviation between the experimental and numerical results of temperature obtained for the case with dry and unsaturated sand were 6% and 2%, respectively. For the case of the degree of saturation shown in Figure 9, a maximum difference of 7% was observed when the experimental and numerical results were compared. In addition, Section 4.1 presents and discusses the results of the numerical a in which all the GEPs are subjected to cyclic thermal loading. Moreover, the resul numerical simulations where the cyclic thermal loads were applied on alternate G presented and discussed in Section 4.2.   The results of temperature evolution versus time for the cyclic heating and cooling of all the GEPs in the group are shown in Figure 10. The investigation was carried out by applying a heat flux value of 25 W/m 2 and −25 W/m 2 to simulate the heating and cooling loads, respectively. The results were obtained at points A, B, and C, as depicted in Figure 1.   From the results shown in Figure 10, the maximum temperature observed a       Additionally, from Figure 10, it can be clearly observed that the greatest temperature changes at all the locations were observed in the sand than in the clay soil during the cyclic heating and cooling simulations. This is owing to the higher thermal conductivity of the sand, which significantly increases with saturation in comparison to that of clay ( Figure  3). The lower thermal conductivity values in clays permit the development of greater temperature magnitude around the GEPs and prevent heat from radiating away to greater From the results shown in Figure 10, the maximum temperature observed across all the locations, i.e., points A, B, and C, occurred in the case with dry soils (i.e., 0% S l ). The average temperature at the observed locations increases uniformly to a maximum value of 55.5, 60, and 56.6 • C at locations A, B, and C, respectively, for the dry sand. The reason why higher temperature magnitude was observed at location B could be attributed to the contribution of heat from the GEPs located in grids 1 and 4 towards the centre of the GEPs. Additionally, as the degree of saturation of the sand increases to 100%, the temperature at A, B, and C decreases to a minimum value of 39.55, 43.65, and 42.55 • C, respectively. This signifies an average drop in temperature with a difference of about 15.95, 16.35, and 14.05 • C at locations A, B, and C, respectively, as the soil S l value increases from 0% to 100%. The drop in temperature due to the soil saturation increases as a result of the contribution of moisture towards heat dissipation. This phenomenon can be explained using Figure 3, where the soil thermal conductivity significantly increases and reaches a maximum as the soil degree of saturation approaches 100%. In addition, it should be noted that the maximum temperature magnitude was witnessed at the end of the heat injection period during the 1st cycle (i.e., the end of 1 year).
The magnitude of the maximum temperature observed decreases nonlinearly with time, at a rate of about 2 • C per year, during the 10-year period to a value of 34.5, 38.2, and 39.4 • C at locations A, B, and C, respectively, in the sand with 0% S l . Similar trends were observed in the soil with 30%, 60%, and 100% S l .
However, during the cooling phases of the tests, the minimum temperature magnitude was observed during the 10th year of the heating and cooling cycles, i.e., the end of the thermal loading. An average lowest temperature value of about −2.5, −10, and −9.8 • C were witnessed at locations A, B, and C, respectively, in the sand. This is directly opposite to the phenomenon that was observed during the heating cycle, where the maximum magnitude of temperature change was witnessed in the 1st thermal cycle. Furthermore, the average temperature magnitude during the cooling cycles decreases nonlinearly with time at a rate of about 1.67 • C per year when the 1st and the 10th cycle of the cooling phases were compared.
Conversely, in the clay soil, the maximum magnitude of temperature change observed at the end of the first heating cycle were 43.3, 51.1, and 51.45 • C at locations A, B, and C, respectively, for the 0% S l . The temperature magnitude at these locations decreases to a value of 37.95, 45.3, and 44.1 • C as the soil degree of saturation increases from 0% to 100%. This shows an average temperature drop of about 5.35, 5.8, and 7.35 • C as the soil saturation increases from 0 to 100% at locations A, B, and C, respectively. Similar to the trend observed in the sand, the greatest temperature change in the clay soil during the heat injection phase was witnessed in the 1st thermal cycle. The temperature magnitude during the 1st heating cycle decreases at an average rate of about 1.4 • C per year to a value of 33.72, 33.4, and 34.65 • C at locations A, B, and C, respectively, during the 10th year.
Moreover, during the cooling phases of the simulations, the major drop in temperature was observed in the 10th cycle of the 10 years simulations for the 0% S l condition. At the end of the 10th year heat extraction phase, an average minimum temperature value of about 13.1, 10.95, and 8.3 • C were observed at locations A, B, and C, respectively. In comparison to the 1st cycle of the cooling phases, an average temperature change of about 0.8 • C per year was observed at all the locations.
Additionally, from Figure 10, it can be clearly observed that the greatest temperature changes at all the locations were observed in the sand than in the clay soil during the cyclic heating and cooling simulations. This is owing to the higher thermal conductivity of the sand, which significantly increases with saturation in comparison to that of clay ( Figure 3). The lower thermal conductivity values in clays permit the development of greater temperature magnitude around the GEPs and prevent heat from radiating away to greater distances in the soil domain. Conversely, in the case of sand, its thermal conductivity was about twice that of the clay as the two soils approach full saturation. This allows the sand to easily transfer heat to greater distances within the soil whilst preventing the development of greater temperature magnitudes at the regions close to the GEPs.

Radial Temperature Distribution
The results were obtained for the radial temperature distribution at the end of the 10th heating and 10th cooling cycles. The results shown in Figure 11 were obtained at the pile mid-depth, from the pile labelled number 4 (or grid 4) to the farthest boundary to the right, located at 21.9 m, from the group of GEPs. The main purpose of obtaining these results is to understand the extent of the soil domain that could be influenced by the cyclic heating and cooling of the group of GEPs under different degrees of saturation condition.

Radial Temperature Distribution
The results were obtained for the radial temperature distribution at the end of the 10th heating and 10th cooling cycles. The results shown in Figure 11 were obtained at the pile mid-depth, from the pile labelled number 4 (or grid 4) to the farthest boundary to the right, located at 21.9 m, from the group of GEPs. The main purpose of obtaining these results is to understand the extent of the soil domain that could be influenced by the cyclic heating and cooling of the group of GEPs under different degrees of saturation condition. The results of the maximum radial temperature distribution with normalised distance at the end of the 10th cycle are shown in Figure 11a for the sand and clay soils with 0 to 100% Sl conditions. The maximum temperature witnessed at the end of the 10th heat injection cycle for the sand and clay soils ranges between 32 and 42 °C, respectively. The temperature magnitude decreases nonlinearly with distance, away from the group of GEPs, towards the far-field boundary to the right, located at a distance of 21.9 m (or 36.5 when normalised with the GEP diameter of 0.6 m). Equally, the magnitude of the temperature developed with distance drops with the increase in soil Sl values from 0% to 100%.
In addition, Figure 11b shows the minimum temperature distribution obtained at the end of the 10th cycle for the two different soil types. At the end of the heating and cooling cycle, the greatest temperature change between −2 °C to −10 °C was witnessed at the pile surface, in the sand and clay soils having 0%-100% Sl values. The minimum temperature associated with the cyclic heating/cooling effect reduces with increasing distance away from the pile. The largest region of temperature influence in the soil domain is about 10 and 15 times the pile diameter for the clay and sand, respectively.
Between the two soils, greater temperature dissipation effects were witnessed in the sand owing to its higher thermal conductivity. This is an intrinsic property of the sand due to its mineralogical composition being quartz [45], in comparison to that of clay soil considered in this study [33].
Thus, the higher thermal conductivity value in the sand would lead to a larger soil area/volume that could be influenced by the cyclic heating and cooling effects of the group of GEPs. Similarly, another important factor that significantly contributes towards heat The results of the maximum radial temperature distribution with normalised distance at the end of the 10th cycle are shown in Figure 11a for the sand and clay soils with 0 to 100% S l conditions. The maximum temperature witnessed at the end of the 10th heat injection cycle for the sand and clay soils ranges between 32 and 42 • C, respectively. The temperature magnitude decreases nonlinearly with distance, away from the group of GEPs, towards the far-field boundary to the right, located at a distance of 21.9 m (or 36.5 when normalised with the GEP diameter of 0.6 m). Equally, the magnitude of the temperature developed with distance drops with the increase in soil S l values from 0% to 100%.
In addition, Figure 11b shows the minimum temperature distribution obtained at the end of the 10th cycle for the two different soil types. At the end of the heating and cooling cycle, the greatest temperature change between −2 • C to −10 • C was witnessed at the pile surface, in the sand and clay soils having 0%-100% S l values. The minimum temperature associated with the cyclic heating/cooling effect reduces with increasing distance away from the pile. The largest region of temperature influence in the soil domain is about 10 and 15 times the pile diameter for the clay and sand, respectively.
Between the two soils, greater temperature dissipation effects were witnessed in the sand owing to its higher thermal conductivity. This is an intrinsic property of the sand due to its mineralogical composition being quartz [45], in comparison to that of clay soil considered in this study [33].
Thus, the higher thermal conductivity value in the sand would lead to a larger soil area/volume that could be influenced by the cyclic heating and cooling effects of the group of GEPs. Similarly, another important factor that significantly contributes towards heat dissipation is the increase in soil granularity and S l . These two important factors, in addition to the mineral composition, are the main reasons that easily allow heat dissipation to a larger area in the sand.

Temperature Distribution with Depth
The results of temperature with depth obtained at half the distance between the GEP-2 and GEP-3 (presented in Figure 1) are shown in Figure 12 for the sand and clay soils, obtained at the end of the 10th heating and 10th cooling cycle. ture magnitude decreases by about 10.65 and 5.4 °C for the sand and clay as the soil satu ration increases from 0% to 100%. The decrease in temperature owing to an increase in soil moisture content is a result of the effect of heat transfer via convection. This is even more prominent as the soil granularity increases, thus allowing moisture in the form o liquid water and vapour to easily flow through the soil voids, as shown in the case o sand. Similarly, at the end of the 10th cooling cycle or heat extraction process, a change in temperature of about −10.3 and −0.25 °C along the pile length was observed in the sand and clay soils, respectively, with a degree of saturation value of 0%.
In addition, during the cooling phases or heat extraction process, the minimum temperature magnitude in the soils with 30% Sl value was observed to be about −8.5 and +2 °C for the sand and clay. Additionally, in soils with 60% Sl, the changes to the temperature magnitude measured were about −7 and +4 °C for the sand and clay, respectively. This temperature range, below 0 °C, is likely to cause ground freezing. Thus, it is advisable tha this should be avoided by either spacing out the GEPs or thermally activating alternate GEPs, which is later discussed in Section 4. At the end of the 10th heating cycle or heat injection process, a temperature magnitude of about 39.3 and 31 • C was observed along the pile length at the equidistant of the GEP-2 and GEP-3 for the sand and clay, respectively, having 0% S l value. The temperature magnitude decreases by about 10.65 and 5.4 • C for the sand and clay as the soil saturation increases from 0% to 100%. The decrease in temperature owing to an increase in soil moisture content is a result of the effect of heat transfer via convection. This is even more prominent as the soil granularity increases, thus allowing moisture in the form of liquid water and vapour to easily flow through the soil voids, as shown in the case of sand.
Similarly, at the end of the 10th cooling cycle or heat extraction process, a change in temperature of about −10.3 and −0.25 • C along the pile length was observed in the sand and clay soils, respectively, with a degree of saturation value of 0%.
In addition, during the cooling phases or heat extraction process, the minimum temperature magnitude in the soils with 30% S l value was observed to be about −8.5 and +2 • C for the sand and clay. Additionally, in soils with 60% S l , the changes to the temperature magnitude measured were about −7 and +4 • C for the sand and clay, respectively. This temperature range, below 0 • C, is likely to cause ground freezing. Thus, it is advisable that this should be avoided by either spacing out the GEPs or thermally activating alternate GEPs, which is later discussed in Section 4.2 of this study. Similarly, this temperature range can be avoided by reducing the amount of heat that is been extracted from the ground. In such a case, other alternative source(s) of energy can be used to compliment and deliver the remaining energy required for space heating operation. Figure 13 presents the results of the degree of saturation versus depth, obtained at half the distance between GEP-2 and GEP-3, for the sand and clay soils with 30 and 60% S l conditions. This location was chosen because the greatest changes in temperature were witnessed at this position. The results were taken at the end of the 10th heating and cooling cycle. range can be avoided by reducing the amount of heat that is been extracted from ground. In such a case, other alternative source(s) of energy can be used to complim and deliver the remaining energy required for space heating operation. Figure 13 presents the results of the degree of saturation versus depth, obtain half the distance between GEP-2 and GEP-3, for the sand and clay soils with 30 and Sl conditions. This location was chosen because the greatest changes in temperature witnessed at this position. The results were taken at the end of the 10th heating and ing cycle. During the heat injection process or heating phase, the changes in soil degree o uration along the GEP-depth remain relatively constant, i.e., within 1.5% for the sand clay soils. This could be attributed to the range of temperature that was witnessed du the cyclic heating operation. In the soil with 30% Sl values, the maximum tempera witnessed was 36 and 30 °C for the sand and clay soils. Similarly, with soils having Sl, the maximum temperature magnitude observed was about 32.5 and 28.3 °C fo sand and clay, respectively ( Figure 12). This temperature range is not expected to c any significant drying of the soil at the observed location.

Temperature Evolution
The results of temperature evolution with time for the cyclic heating and coolin alternate GEPs in the group are shown in Figure 14 for the sand and clay soils with 100% Sl. In the course of these analyses, a heat flux of +25 W/m 2 and −25 W/m 2 wer plied at the surface of the GEPs in grids 2 and 4 to numerically simulate the effect of c heat injection and extraction associated with the energy piles. The main rationale be During the heat injection process or heating phase, the changes in soil degree of saturation along the GEP-depth remain relatively constant, i.e., within 1.5% for the sand and clay soils. This could be attributed to the range of temperature that was witnessed during the cyclic heating operation. In the soil with 30% S l values, the maximum temperature witnessed was 36 and 30 • C for the sand and clay soils. Similarly, with soils having 60% S l , the maximum temperature magnitude observed was about 32.5 and 28.3 • C for the sand and clay, respectively ( Figure 12). This temperature range is not expected to cause any significant drying of the soil at the observed location.

Temperature Evolution
The results of temperature evolution with time for the cyclic heating and cooling of alternate GEPs in the group are shown in Figure 14 for the sand and clay soils with 0 to 100% S l . In the course of these analyses, a heat flux of +25 W/m 2 and −25 W/m 2 were applied at the surface of the GEPs in grids 2 and 4 to numerically simulate the effect of cyclic heat injection and extraction associated with the energy piles. The main rationale behind carrying out this analysis is to understand the heat flow behaviour as a result of the increase in the spacing between GEPs in a group. However, the main limitation is that the applied total thermal load on the GEPs is not equal in the two cases investigated. This is to ensure that excessive temperature development is avoided. Nonetheless, it provides an opportunity for future work in this area. Thus, it was observed that the phenomenon of cyclic heat injection and extraction results in a nonlinear temperature decrease within the sandy soil domain, with the rate of temperature decrease declining with increasing thermal cycle. However, in the clay soil, the phenomenon was observed to follow a linear pattern. Perhaps, this could be associated with the inherent properties of the soils, such as mineralogy, grain size distribution, and other thermal and hydraulic properties, which govern the flow of heat from the heating 0% to 100% for the sand. The drop in temperature is associated with the contribution of soil moisture towards heat dissipation away from the thermally activated GEPs.
The temperature changes in the 1st heating cycle at locations A, B, and C drops nonlinearly at a rate of about 1.2 • C per year to about 33.95, 22.55, and 26 • C in the 10th cycle of the heat injection and extraction process. Equally, as the soil saturation increases to 100%, the magnitude of temperature changes at the observed locations further decreases at a rate of about 0.9 • C per year to a value of 25.1, 18.6, and 20.8 • C, respectively.
At the end of the 1st cooling cycle of the simulations, the minimum temperature observed at locations A, B, and C was found to be 5.1, 14.05, and 8.3 • C, and the effect of the temperature changes on initial temperature decreases to about 9.35, 14.65, and 11.8 • C as the degree of saturation of the sand increases from 0% to 100%. However, at the end of the 10th cycle, temperatures of about −8.15, 4.15, and 0.45 • C were observed in the sand with 0% S l . As the soil S l increases to 100%, the changes in temperature magnitude caused by the cooling load reduce to about 1.45, 8.8, and 5.1 • C at the observed locations, respectively. Figure 14 also presents the results of temperature evolution at locations A, B, and C for the clay soil with S l values ranging from 0% to 100%. The maximum temperature magnitude occurred during the 1st thermal cycle. Similar to the sand, the highest temperatures were observed in soil with 0% S l values. During the 1st cycle, a temperature of about 40.4, 29.8, and 34.55 • C were witnessed at locations A, B, and C, respectively. The temperature magnitude decreases to 32.8, 24.45, and 27.25 • C as the clay soil saturation increases to 100%. In comparison to the 10th cycle, the temperature at the observed locations decreases nonlinearly to about 27.85, 24.55, and 25.7 • C, respectively. The values further decrease down to about 22.85, 19.65, and 20.8 • C as the soil reaches full saturation. However, during the 1st cooling cycle, the minimum temperature observed in the clay soil were 20.85, 16.6, and 15.95 • C, respectively, with the magnitude of temperature change observed decreasing to about 14.45, 15.9, and 14.1 • C as the clay reaches 100% saturation. In the 10th cooling operation or heat extraction process, the temperature values decrease linearly at a rate of about 0.7 • C per year due to the cumulative heat extraction to a value of about 8.4, 11.2, and 11 • C, respectively. In addition, the temperature observed at the 10th cooling cycle in the clay with 0% S l further decreases to 9, 11.5, and 8.05 • C as the S l value increases to 100%.
Thus, it was observed that the phenomenon of cyclic heat injection and extraction results in a nonlinear temperature decrease within the sandy soil domain, with the rate of temperature decrease declining with increasing thermal cycle. However, in the clay soil, the phenomenon was observed to follow a linear pattern. Perhaps, this could be associated with the inherent properties of the soils, such as mineralogy, grain size distribution, and other thermal and hydraulic properties, which govern the flow of heat from the heating source. In the sand, the higher thermal and hydraulic properties help in easily dissipating the thermal load from the GEP. Conversely, in the clay, these properties are quite moderate and thus allowing a very slow heat dissipation process.
In addition, it was observed that the soil degree of saturation results in a positive influence on the thermal performance of the GEPs. In heat injection mode (i.e., space cooling mode), the soil moisture helps to easily dissipate the heat being injected by the heat pump system to the surrounding soil, thereby, preventing excessive temperature build-up in the soil surrounding the piles. This effect rises as the soil saturation increases and reaches a maximum at an S l value of 100%. In contrast, during heat extraction from the ground (i.e., space heating), the soil temperature decreases depending on the rate of the cooling load superimposed on to the GEP elements. The presence of soil moisture helps in ensuring that the excessive heat extraction does not cause ground freezing. Thus, as observed in the myriads of numerical simulations carried out, the minimum sub-zero temperature observed was in the soil with 0% S l . However, as the initial soil S l value increases, the lowest temperature magnitude remains above 0 • C throughout the 10 years of numerical simulations. The effect further decreases as the soil S l increases towards 100%. Furthermore, lower temperature values were observed at all the locations in the simulations where alternate GEPs were thermally activated in comparison to the heating or cooling of all the GEPs in the group. A temperature difference of about 10 • C or more was observed, especially in the drier soils and that with lower S l values. However, the temperature difference decreases as the soil saturation increases for the sand and clay soils, respectively. In addition, this approach is useful, especially where piles are closely situated next to each other. In such a case, alternate GEPs can be thermally activated to maximise the thermal performance of the GEP system. Figure 15a presents the results of maximum radial temperature distribution versus normalised distance obtained at the end of the 10th heating cycle for the sand and clay soils. The maximum temperature distribution at the observed location (distance between GEPs to the far-field boundary to the right of the GEPs, located at a distance of 21.9 m or 36.5 when normalised with the GEP diameter of 0.6 m) ranges between +32 and 43 • C for the sand and clay soils. The magnitude of temperature distribution decreases nonlinearly with distance away from the group of GEPs. The maximum region of temperature influence within the clay and sand soils ranges at a maximum normalised distance of about 10 to 15 times the pile diameter, respectively. A similar range of temperature magnitudes was observed when the results of radial temperature distribution for thermally activating all the GEPs shown in Figure 11 are compared with the results where alternate GEPs were thermally activated in the group as shown in Figure 15. This could be as a result of the fact that in all the cases, the GEP in the 4th column or grid was thermally activated. Thus, indicating that the radial temperature distribution was a result of the heat that is majorly emanating from the GEP in grid 4. Figure 16 shows the results of temperature magnitude versus depth obtained at half the distance between GEP-2 and GEP-3 for the sand and clay soils. Figure 16a shows the results obtained at the end of the 10th heating cycle for the two soils with Sl values ranging from 0% to 100%. In the soil with a 0% Sl value, the maximum temperature magnitude observed at this location was about 24.15 and 20.95 °C for the sand and the clay soils, respectively. The temperature magnitude decreases to about 19.3 °C for the two soils as the soil saturation reaches maximum.

Temperature Distribution with Depth
However, at the end of the 10th cycle during the cooling period, the minimum temperature magnitude witnessed along the pile depth was about 4 and 8.7 °C for the sand and clay soils, respectively, in a fully dry state, as shown in Figure 16b. However, as the In addition, Figure 15b shows the results of the minimum temperature distribution with normalised distance at the end of the 10th cooling cycle. The observed changes in radial temperature distribution range between −2 and −12 • C for the sand and clay soils, respectively. The temperature change magnitude decreases with an increase in soil saturation. Similar to the case of the heat injection cycle, the region of thermal influence around the group of GEPs lies at a distance of about 10-15 times the diameter of the GEP for the clay and sand, respectively.
A similar range of temperature magnitudes was observed when the results of radial temperature distribution for thermally activating all the GEPs shown in Figure 11 are compared with the results where alternate GEPs were thermally activated in the group as shown in Figure 15. This could be as a result of the fact that in all the cases, the GEP in the 4th column or grid was thermally activated. Thus, indicating that the radial temperature distribution was a result of the heat that is majorly emanating from the GEP in grid 4. Figure 16 shows the results of temperature magnitude versus depth obtained at half the distance between GEP-2 and GEP-3 for the sand and clay soils. Figure 16a shows the results obtained at the end of the 10th heating cycle for the two soils with S l values ranging from 0% to 100%. In the soil with a 0% S l value, the maximum temperature magnitude observed at this location was about 24.15 and 20.95 • C for the sand and the clay soils, respectively. The temperature magnitude decreases to about 19.3 • C for the two soils as the soil saturation reaches maximum.

Degree of Saturation with Depth
The results of the degree of saturation with depth, obtained at half the distance between GEP-2 and GEP-3, for the sand and clay soils having 30 and 60% Sl conditions are shown in Figure 17. The results were obtained at the end of the 10th heating and cooling cycle of the numerical analyses carried out by thermally activating alternate GEPs in grids 2 and 4.
At the end of the 10th heating cycle, the changes in the degree of saturation of the sand and clay soils along the GEP-depth were found to be within 1% for both the 30 and 60% Sl values. Thus, the changes are considered negligible at this particular location. Perhaps, this is associated with the maximum temperature magnitude witnessed at this location, which was in the region of about 22 and 20 °C for the 30% and 60% Sl, for both the two soils at the end of the 10th heat injection process.
However, at the end of the 10th cycle, during the cooling period for the 30% Sl value, a minimum temperature of about 4.7 and 9.6 °C were observed for the sand and clay soils. In addition, the soils with 60% degree of saturation had a temperature of about 5.6 and 10.3 °C for the sand and clay, respectively. Thus, the temperature range is unable to cause any major changes to the initial degree of saturation. However, at the end of the 10th cycle during the cooling period, the minimum temperature magnitude witnessed along the pile depth was about 4 and 8.7 • C for the sand and clay soils, respectively, in a fully dry state, as shown in Figure 16b. However, as the soil saturation increases towards 100%, the minimum temperature magnitude observed in the sand and clay soils were found to be about 7.75 and 11.1 • C, respectively. Furthermore, when the temperature distribution versus depth for the two cases: (1) where all the four GEPs were heated and cooled together and (2) thermally activating alternate GEPs, investigated in this study are compared; it was observed that thermally activating alternate GEPs results in lower temperature distribution at the location where the results were obtained, i.e., half the distance between GEP-2 and GEP-3. This could be attributed to the fact that the GEP-1 and GEP-3 were not thermally activated. Thus, the associated contribution of thermal changes from the thermally inactive GEPs is absent. In contrast, in the case where all the GEPs were thermally activated, the GEP-2 and GEP-3 equally contributed to the thermal changes that took place at the centre of the two piles.

Degree of Saturation with Depth
The results of the degree of saturation with depth, obtained at half the distance between GEP-2 and GEP-3, for the sand and clay soils having 30 and 60% S l conditions are shown in Figure 17. The results were obtained at the end of the 10th heating and cooling cycle of the numerical analyses carried out by thermally activating alternate GEPs in grids 2 and 4. Furthermore, the results of the degree of saturation involving the case where all the GEPs were thermally activated and that of heating and cooling of alternate GEPs are compared and discussed here. In the case where all the GEPs were heated, the range of temperature results in the maximum temperature changes within the soil domain for the clay and sand. Similarly, the heat extraction mode resulted in lower temperatures in the subzero range. A minimum temperature magnitude of about −10 °C was observed, which could result in ground freezing. This temperature is expected to keep decreasing with increasing thermal cycles. Thus, this should be avoided by controlling the heat that is extracted from the system to ensure that the overall thermal performance and geotechnical behaviour of the GEP are not jeopardised.

Conclusions
This paper presents the results of series of numerical investigations investigating the long-term thermal and hydraulic behaviour of a group of energy piles installed in unsaturated sand and clay. The simulations were carried out by applying cyclic transient heating and cooling loads on the GEPs for a total duration of 10 years. The thermal load pattern comprises 6 months of heating followed by 6 months of cooling each year. In addition, the numerical simulations investigated two cases, including (1) a case where all the GEPs in At the end of the 10th heating cycle, the changes in the degree of saturation of the sand and clay soils along the GEP-depth were found to be within 1% for both the 30 and 60% S l values. Thus, the changes are considered negligible at this particular location. Perhaps, this is associated with the maximum temperature magnitude witnessed at this location, which was in the region of about 22 and 20 • C for the 30% and 60% S l , for both the two soils at the end of the 10th heat injection process.
However, at the end of the 10th cycle, during the cooling period for the 30% S l value, a minimum temperature of about 4.7 and 9.6 • C were observed for the sand and clay soils. In addition, the soils with 60% degree of saturation had a temperature of about 5.6 and 10.3 • C for the sand and clay, respectively. Thus, the temperature range is unable to cause any major changes to the initial degree of saturation.
Furthermore, the results of the degree of saturation involving the case where all the GEPs were thermally activated and that of heating and cooling of alternate GEPs are compared and discussed here. In the case where all the GEPs were heated, the range of temperature results in the maximum temperature changes within the soil domain for the clay and sand. Similarly, the heat extraction mode resulted in lower temperatures in the sub-zero range. A minimum temperature magnitude of about −10 • C was observed, which Energies 2021, 14, 4122 24 of 28 could result in ground freezing. This temperature is expected to keep decreasing with increasing thermal cycles. Thus, this should be avoided by controlling the heat that is extracted from the system to ensure that the overall thermal performance and geotechnical behaviour of the GEP are not jeopardised.

Conclusions
This paper presents the results of series of numerical investigations investigating the long-term thermal and hydraulic behaviour of a group of energy piles installed in unsaturated sand and clay. The simulations were carried out by applying cyclic transient heating and cooling loads on the GEPs for a total duration of 10 years. The thermal load pattern comprises 6 months of heating followed by 6 months of cooling each year. In addition, the numerical simulations investigated two cases, including (1) a case where all the GEPs in the group are thermally activated and (2) where alternate GEPs are activated thermally. It was found that:

•
The maximum temperature magnitude in all the simulations were observed at point B, i.e., the mid-point for all the group of GEPs. This could be as a result of the contribution of heat radiating from all the four GEPs towards the centre (point B). In addition, the temperature magnitude at this location and all the other locations decreases as the initial degree of saturation of the soil increases from 0% to 100% in the sand and the clay soils.

•
Similarly, the temperature magnitudes at locations A, B, and C were observed to be higher in the sand than in clays. This phenomenon can be described based on the thermal conductivity of the two soils: sand and clay. The lower thermal conductivity values in clays permit the development of greater temperature magnitude around the GEPs and prevent heat from radiating away to greater distances in the soil domain. In contrast, in the case of sand, its thermal conductivity was about twice that of the clay as the two soils approach full saturation. This allows the sand to easily transfer heat to greater distances within the soil whilst preventing the development of greater temperature magnitudes at the regions close to the GEPs. Similarly, soil permeability is another factor that is greater in sands and allows heat to be easily transferred via convection due to its high porosity.

•
The continuous heating and cooling cycles of the GEPs, during the 10 years of numerical simulations, resulted in a continuous decrease in temperature trend at all the observed locations (i.e., A, B, and C, respectively). 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 intermittent heating and cooling mode of operation permits the soil domain to naturally recover and prevent excessive heat build-up or deficit in the soil [46]. This is important in safeguarding the long-term thermal performance of the system.

•
In addition, thermally activating some selected and alternate GEPs in the group results in lower temperature magnitude at all the observed locations in the soil domain, and the temperature change decreases with an increase in soil saturation. This is a result of the increase in the spacing between the GEPs that are thermally activated, thus allowing greater soil volume for heat rejection and solicitation. This is especially important where the structural pile foundation elements are situated close to each other and/or where they are installed in problematic clays. This could lead to the heaving or contraction of the clays and the GEPs due to excessive heat injection or extraction. Thus, thermally activating alternate GEPs offers an alternative approach towards installing an energy-efficient GEP system, and ensure the longevity of the system. • Furthermore, temperatures that are well below freezing were observed at some point during the simulations. It is strongly advised that the application of any temperature magnitude that could result in extreme changes to the soil and pile(s) should be avoided. Because the superimposition of extreme heating and cooling loads results