A Thermal Conductivity Model for Deformable and Unsaturated Soils to Assess the Thermal Behaviour of Energy Piles

This paper presents an empirical model that predicts the thermal conductivity of soils by accounting for the effect of both the degree of saturation and void ratio on the heat exchange capacity of shallow geothermal reservoirs. The model is generated by the product of two terms: the former accounts for the effect of void ratio on the dry thermal conductivity, whereas the latter describes the influence of the degree of saturation on the moisture-dependent thermal conductivity. The model is a function of three parameters, which are easy to calibrate based on their physical meaning. Model predictions are validated against five different sets of experimental data from literature by means of two alternative approaches: blind prediction of thermal conductivity measurements not employed during calibration and numerical simulations of thermal tests performed on energy piles. Results show that the proposed model is capable of accurately predicting both the thermal conductivity of deformable unsaturated soils as well as reproducing the thermal behaviour of energy piles.


Introduction
Shallow geothermal systems are composed of heat exchangers that exploit ground temperatures which remain approximately constant throughout the year, at depths below the ground surface of about 10-15 m.Heat exchangers are embedded within piled foundations to create thermo-active energy piles that use the constant ground temperature as a heat source or sink to warm up or cool down the indoor environment of buildings [1].Energy piles can reduce energyrelated costs by at least 30% compared with conventional heating and cooling systems [2,3].
An accurate analysis of the thermal performance of energy piles requires a reliable prediction of the thermal conductivity of the surrounding soil as well as its dependency on both the degree of saturation and porosity of the geothermal reservoir.Several moisture-dependent thermal conductivity models have been proposed in the literature, ranging from physical models that account for the thermal conductivities of the volume fractions of soil grains, water and air [4][5][6][7] to empirical relationships that relate the thermal conductivity with physical properties of the soil, such as the water content (either gravimetric or volumetric) and porosity [8][9][10].
This paper proposes a simple relationship to predict the moisture-dependent thermal conductivity of deformable and unsaturated soil by improving the expression proposed by Bruno and Alamoudi [10].In particular, the proposed model postulates that the thermal conductivity of a deformable and unsaturated soil depends on two interacting terms accounting for (a) the effect of the degree of saturation on the moisture-dependent thermal conductivity and (b) the effect of void ratio on the dry thermal conductivity of the soil.This represents a significant improvement over the model of Bruno and Alamoudi [10], which can only predict a constant thermal conductivity under dry conditions.The proposed model has been validated against five different sets of experimental data either via blind prediction of experimental measurements of thermal conductivity not considered during calibration or via numerical simulations of the thermal behaviour of energy piles to reproduce the thermal tests performed by Akrouch et al. [11].These simulations follow a similar approach to that adopted by several other studies not only for energy piles [12][13][14] but also for other thermal active geotechnical constructions [15][16][17][18][19][20][21][22].The comparison between the experimental and predicted behaviour of the energy piles highlights the excellent capability of the proposed thermal conductivity model to predict the heat exchanges between the energy piles and the surrounding soil.

Thermal Conductivity Model
Bruno and Alamoudi [10] proposed a simple empirical relationship between the ratio dry (where is the moisturedependent thermal conductivity and dry is the dry thermal conductivity) and both the degree of saturation S r and the dry density d of the shallow geothermal reservoir.This relationship is expressed as follows: where w is the gravimetric water content, s is the density of the solid particles, m f and m d are two model parameters named moisture factor and dry density factor, respectively.The relationship proposed by Bruno and Alamoudi [10] successfully predicted the thermal conductivity of a broad range of soil granularities from fine clays to coarse sands, compacted at various levels of dry density, ranging from 1465 kg/m 3 to 2310 kg/m 3 , and spanning the whole range of the degree of saturation (i.e. from dry to fully saturated conditions).However, the model proposed by Bruno and Alamoudi [10] can only predict a constant dry thermal conductivity regardless of the level of dry density.To overcome this limitation, the present paper introduces a novel relationship between the moisture-dependent thermal conductivity and both the degree of saturation S r and the void ratio e of the soil.This new relationship is given by the following expression: where dry,e=1 is the thermal conductivity of the dry soil at a constant and unitary void ratio (i.e. at a porosity of 0.5) and m e is the void ratio factor.Note that the thermal behaviour of deformable and unsaturated soils not only depends on the degree of saturation and porosity but also on a number of other factors (e.g.material mineralogy, grain size (1) e m e distribution, and soil fabric).Moreover, the thermal conductivity also depends on the actual distribution of solid, water and air phases.This is particularly important when the saturation state evolves from funicular (i.e.water is continuous whilst air is discontinuous) to pendular (i.e.water is discontinuous whilst air becomes continuous) conditions.In fact, the continuity of the water fluid in the funicular state creates preferential paths for heat exchanges that progressively vanish as the soil desaturates and the less-conductive air fluid becomes continuous.Finally, soils are porous media characterised by complex natural structures and significant variation of thermal properties in different points (i.e.heterogeneity) and along different directions (i.e.anisotropy).A more accurate description of the heat exchanges in deformable and unsaturated soils should account for all these aspects.However, the main scope of the present work is to capture the macroscopic thermal behaviour of soils whereas a more sophisticated modelling is deemed to be outside the scope of the present paper and it constitutes matter for future research.
Inspection of Eq. ( 2) also provides a physical interpretation of the moisture factor m f , which represents the variation rate of the moisture-dependent thermal conductivity with a varying degree of saturation.The moisture factor m f can also be interpreted as the difference between the ratio e=1 dry,e=1 (with e=1 being the moisture-dependent thermal conductivity at a unitary void ratio) calculated under fully saturated conditions, and the same ratio determined under dry conditions, thus resulting in the following relationship: The dry thermal conductivity dry can be formulated as a function of the void ratio of the geothermal reservoir as follows: which can be re-casted by taking the logarithms as follows: Hence, Eqs.(4,5) introduce the dependency of the dry thermal conductivity on the volumetric state of the soil, thus overcoming the limitation of the relationship proposed by Bruno and Alamoudi [10].Inspection of Eqs.(5) shows that the void ratio factor m e represents the slope of the straight line in the double logarithmic plane log dry dry,e=1 − loge .Note also that Eqs.(4,5) would predict an infinite thermal conductivity when the void ratio becomes equal to zero (i.e.all =−m e loge porosity is erased).It is therefore recommended to limit the validity of Eqs.(4,5) only to void ratios which are higher than a minimum value, which can be determined by means of standard testing procedures for measuring the relative density in coarse soils or by means of cyclic minimum void ratio tests in finer soils, as recommended by Kim et al. [23].The dry thermal conductivity of the soil at minimum values can be therefore expressed by means of the following equation: with e min being the minimum void ratio.Figure 1 shows a schematic representation of Eq. ( 5) together with the geometrical meaning of the parameters dry,e=1 and m e as well as the recommended limitation for the minimum value of void ratio.Note that the experimental data sets considered in the present work are characterised by values of void ratio sensibly higher than the minimum ones, thus avoiding the need of considering both the above-mentioned limitation and the calibration of the thermal conductivity of the solid grains.
By combining Eqs. ( 2) and ( 4), the moisture-dependent thermal conductivity of the geothermal reservoir can therefore be expressed as follows: Equation ( 7) formally recovers the expression proposed by Liuzzi et al. [9], with the difference that the effect of ( 6) dry,e=e min = s 1 + e min (7) the volumetric behaviour on the thermal conductivity of the soil is now taken into account via the dry thermal conductivity dry , as expressed by Eqs.(4,5).
The proposed thermal conductivity model, as expressed by Eq. ( 7), has been calibrated and then validated against five different sets of experimental data by Tang et al. [24], Hall and Allinson [25], Akrouch et al. [11], Mansour et al. [26] and Song et al. [27].Two distinct strategies have been adopted for validating the proposed thermal conductivity model.The blind prediction method consists in predicting the experimental measurements of moisture-dependent thermal conductivity not considered during calibration.Blind prediction was employed to validate the model against the data sets of Tang et al. [24], Hall and Allinson [25], Mansour et al. [26] and Song et al. [27].Instead, numerical simulations, performed by means of the commercial software CODE_BRIGHT, were developed to validate the thermal conductivity model against the data set from Akrouch et al. [11], who performed thermal tests on small portions of energy piles, as described in the following section.

Experimental Data Sets
This section summarises the main physical properties of the different soils considered in this work, including the grain size distribution, plasticity properties, specific gravity of solid particles and both the water content and dry density after compaction.Tang et al. [24] measured the moisture-dependent thermal conductivity of an expansive MX80 bentonite by means of the sensor KD2 from Decagon Devices Inc.The MX80 bentonite exhibited a specific gravity of solid particles G s of 2.76, a liquid limit of 520% and a plastic index of 478%, thus attaining a very high plasticity compared with the other data sets considered in the present study.Tang et al. [24] compacted the MX80 bentonitic samples at a dry density that ranged from 1465 kg/m 3 to 1801 kg/m 3 and with a degree of saturation that varied from 0.29 to 0.86.

Data Set by Hall and Allinson [25]
Three earth mixes composed of 14-6.3 mm rounded pea gravel, 5 mm-down medium grade grit sand and a silty clay were tested by Hall and Allinson [25].In agreement with the nomenclature proposed by the authors, also in this paper, the three mixes are named hereafter as 433, 613 and 703 based on the mass proportion in units of ten of pea gravel (1st digit), grit sand (2nd digit) and silty clay (3rd digit).All soil mixes were stabilised with the addition of 6% by total dry mass of ordinary CEM IIa Portland cement, mixed with the optimum water content of 8% and then compacted according to the standard Proctor procedure as prescribed by the norm BS 1377-4 [28].The thermal conductivity was measured by means of a computer-controlled P.A. Hilton B480 heat flow metre apparatus with descending vertical heat flow, as prescribed by the standards ISO 8301 [29].These measurements were taken on earth samples made by the three earth mixes at a dry density ranging from 1939 kg/m 3 to 2230 kg/ m 3 and at a degree of saturation varying from 0 up to 0.6.

Data Set by Mansour et al. [26]
Mansour et al. [26] performed hot wire tests to measure the moisture-dependent thermal conductivity of soil samples composed by 59.9% of clay and silt, 39.1% of sand and 1.0% of gravel with a liquid limit of 22.4% and a plastic index of 6.6%.The soil was mixed with a water content of 13% and then compacted at different dry densities ranging from 1577 kg/m 3 to 2147 kg/m 3 , which resulted in a degree of saturation varying from 0.076 to 0.256.

Data Set by Song et al. [27]
Song et al. [27] measured the thermal conductivity of a Xyanyang clay by means of the Hot Disc Transient Plane Source technique.The Xyanyang clay presented a liquid limit of 37.0%, a plasticity index of 13.9% and a specific gravity of the solid particles G s equal to 2.74.This soil was compacted at different dry densities ranging from 1550 kg/ m 3 to 2070 kg/m 3 and at various levels of water content varying from 3 to 24% (i.e.corresponding to a variation of degree of saturation from 0.107 to 0.980).

Data Set by Akrouch et al. [11]
Akrouch et al. [11] performed laboratory tests to investigate the effect of the soil saturation and porosity on the thermal behaviour of energy piles.In particular, these laboratory tests were conducted inside sand boxes of dimensions 1200 mm × 1200 mm × 25 mm in which a circular energy pile of a diameter of 300 mm was built.The energy pile hosted two thermal pipes made of PVC materials and with an inner diameter of 38 mm and an outer diameter of 42 mm.Water was circulated inside the thermal pipes at a constant temperature of 37 °C whilst the sand box was stored at a constant room temperature of 21 °C.The temperature gradient between the thermal pipes and the sand box was maintained constant for a minimum time of 48 h during which temperature was monitored by six probes positioned in six different points along the two perpendicular directions in the horizontal cross section.Temperature probes B and E were placed at the pile-soil interface whereas the couples of probes C, F and D, G were respectively placed at a distance of 7.5 cm and 15 cm from the pile-soil interface, always along the same two perpendicular directions.The soil surrounding the energy pile was compacted at a dry density ranging from 1416 kg/m 3 to 1488 kg/m. 3and equalised to various levels of degree of saturation ranging from 0.015 to 1.Note that these ranges of both dry density and degree of saturation are solely referred to the Geometry I tested by Akrouch et al. [11] and this is because only this geometry was numerically modelled in the present work.Further details on the testing set-up can be found in Akrouch et al. [11].
Inspection of the five data sets indicates that different experimental protocols have been followed to measure the thermal conductivity of deformable and unsaturated soils.This is a particularly challenging task because of the complex interactions amongst the three soil phases (solids, water, and air) as well as the heat exchanges occurring between the soil and the atmosphere.However, for the purpose of the present work, the experimental measurements of thermal conductivity taken from the five data sets are assumed reliable.In summary, the five different data sets, considered in the present work, cover a broad range of soil granularities varying from bentonitic clay to coarse sands.Also, soil samples were compacted at various levels of dry density ranging from 1418 kg/m 3 to 2230 kg/m 3 and spanning the whole spectrum of degree of saturation (i.e. from dry to fully saturated conditions).Note that the proposed thermal conductivity model could be used for modelling heat exchanges in shallow geothermal or other geotechnical and geo-environmental applications as long as the variations of granularity, dry density and degree of saturation are within the above-mentioned ranges.Future research will be directed at further extending the validation of the proposed analytical formulation to additional experimental data in view of broadening its field of application.

Model Parameters Calibration
The proposed thermal conductivity model is therefore dependent on a total of three model parameters, i.e. dry,e=1 , m f and m e , which were calibrated by means of a least-square regression of some of the available experimental measurements of thermal conductivity for each of the five data sets.
The experimental measurements of thermal conductivity used for calibration were selected so that the corresponding levels of degree of saturation fall within either the upper or the lower 25% of the whole saturation range investigated by a given data set.This choice is dictated by the need of making the validation stage more demanding as the model is required to extrapolate the values of thermal conductivity at levels of degree of saturation not considered during calibration.Moreover, the selection of the experimental thermal conductivity for model calibration was based on the degree of saturation and not on the void ratio because (a) soil samples exhibited proportionally larger variation of the saturation state in comparison with the changes in volumetric state and (b) the variation of thermal conductivity exhibits a stronger dependency on the changing degree of saturation rather than the void ratio, as observed in all the data sets considered in the present study.
Figure 2 shows both the experimental and calibrated thermal conductivity for all data sets considered in the present work.Inspection of Fig. 2 indicates the excellent capacity of the calibrated model to reproduce the experimental thermal conductivity of the soils from all the data sets, regardless of the soil granularity or both the saturation and volumetric state of the soil.Table 1 lists the numerical values of all the model parameters, which will be used for validation  purposes (i.e. both blind prediction and numerical simulations) in the following section.

Model Validation-Blind Prediction
Model predictions have been validated against experimental values of thermal conductivity not considered during the calibration stage and results are reported in Fig. 3.Note that the model extrapolates the values of thermal conductivity at levels of degree of saturation different than those used for calibration and this choice is dictated by the need of imposing a more demanding validation procedure.Despite the imposed calibration condition, the quality of the model predictions is excellent considering the high values of correlation factor r and the corresponding low mean relative error MRE for all the considered data sets, as reported in Fig. 3. Once again, the proposed model provides an accurate prediction of the thermal conductivity regardless of the soil granularity and both the saturation and volumetric states.Figure 4 compares the experimental values of thermal conductivity against those predicted by the model by Bruno and Alamoudi [10].Inspection of Fig. suggests that model Fig. 3 Predicted and experimental moisture-dependent thermal conductivity for all data sets Fig. 4 Experimental moisturedependent thermal conductivity for all data sets compared with predictions from Bruno and Alamoudi [10] predictions are comparable to those of the present formulation as long as the void ratio remains approximately constant, i.e. data set by Hall and Allinson [25].Model predictions from Bruno and Alamoudi [10] tend to diverge from the experimental measurements for data sets characterised by more significant variations of the void ratio (e.g.Tang et al. [24]; Mansour et al. [26]; Song et al. [27]).This highlights the improved consideration of the effect of soil porosity on the variation of the thermal conductivity proposed in the present model.

Model Validation-Numerical Simulations
This section first presents the main features of the numerical simulations developed by means of the finite element software CODE_BRIGHT v.8.4 [30][31][32].In the second part of the section, outcomes from the numerical simulations are then compared with the experimental thermal behaviour of the energy pile as measured by Akrouch et al. [11].
First, the geometry of the circular energy pile was built by reproducing the experimental set-up adopted by Akrouch et al. [11], in a 2D cross section.An unstructured mesh was assigned, with 2506 nodes constituting 2433 quadrilateraltype elements, with mesh refinement near the pipe-pile and pile-soil interfaces.Numerical simulations were conducted by assuming that thermal properties of both soil and pile are isotropic and not dependent on temperature.At the beginning of all the simulations, the whole domain is at the same initial temperature of 21 °C before the two thermal pipes are heated up to a temperature of 37 °C.Then, the temperature of the thermal pipes is maintained constant for a minimum period of 48 h similar to the experimental protocol adopted by Akrouch et al. [11].The thermal pipes are therefore modelled as cylindrical heat source.Note that conduction is the only heat transport mechanism that is modelled in the FEM simulations whilst the contribution from heat convection is considered negligible due to lack of heat transport originating from advective fluxes in the porous medium.
The theoretical framework of CODE_BRIGHT consists of multiphase (solid, liquid, and gas) and multispecies (solid, water, and dry air) approach.The state variables are the liquid pressure and the soil temperature.The formulation therefore consists of defining the mass balance equations per species following a compositional approach, which overcome the need of modelling the heat transfer due to phase changes. is because it does not appear explicitly in the balance equations.The solid mass balance equation is expressed in terms of the porosity, whilst the liquid mass balance equation is written in terms of the unknown liquid pressure as a state variable.The sum of the advective j E (where stands for the phases: solid, liquid, and gas) and non-advective i phase fluxes constitute the mass movement fluxes.Assuming thermal equilibrium between phases at a given node [33], the internal energy balance equation can be expressed in terms of the changing temperature (due to fluctuating heat fluxes) and the internal energy of each of the phases E (J m −2 s −1 ).Darcy's law represents the advective single-phase flow in isotropic porous media, where flow is assumed proportional to the pressure drop whilst being inversely proportional to the viscosity.
where g is the gravity vector, n is the porosity at a given timestep, n 0 is the initial porosity, (Pa.s) is the dynamic viscosity of the phase .The intrinsic permeability tensor k 0 (m 2 ) is computed using Kozeny's model for a continuum medium from the hydraulic conductivity whilst the phase relative permeability k r is separately defined as a function of the degree of saturation (van Genuchten [34]).Nonadvective fluxes generally embrace molecular diffusion and mechanical dispersion and are calculated using Fick's law in terms of gas saturation, tortuosity and dispersion coefficients.Heat conduction is expressed by Fourier's law in terms of the moisture-dependent thermal conductivity (Wm −1 K −1 ), as follows: where the thermal conductivity is calculated by means of Eq. (7).To replicate the experimental conditions, the numerical boundary conditions did not allow evaporation, and the soil-atmosphere interaction was limited to the thermal component.The general numerical boundary condition form, applied to external soil boundaries and thermal pipes, is obtained by adding energy fluxes at the interface nodes, which are reduced to imposed temperature: where T is the soil temperature, T 0 is the imposed tempera- ture, and e is the heat transfer coefficient acting as a thermal diffusion leakage coefficient.According to the numerical assumptions made by Cuadrado et al. [35], a high numerical value of e = 10 3 W/m 2 K was assigned at the thermal pipe-concrete interface indicating a high heat transfer gradient at the thermal pipes circulating the heat carrying fluid.A much lower value of e = 10 −3 W/m 2 K was assigned at the soil-atmosphere interface, indicating very low heat transfer occurring between the soil and the external atmosphere, to numerically replicate the heat insulation conditions imposed by Akrouch et al. [11] in the experimental set-up.
Table 2 summarises the constitutive equations employed to relate the state variables with the dependent variables (e.g.link between the liquid pressure and the degree of saturation).Instead, Table 3 lists the numerical values of all the (8) During the simulations, the variation of temperature over time is monitored by the temperature probes B, C, D, E, F and G, positioned as described above and reported by Akrouch et al. [11].Results from these simulations are presented in Fig. 5, which shows the excellent ability of the numerical model, coupled with the proposed thermal conductivity function, to reproduce the heat propagation and the consequent variation of temperature over time, as experimentally measured by Akrouch et al. [11] with the temperature probes B, C and D (Fig. 5a) and the probes E, F and G (Fig. 5b) in the soil equalised at a degree of saturation equal to 0.015.
To further validate the proposed numerical model, simulations were run at six different levels of degree of saturation ranging from almost dry (i.e. S r = 0.015) to fully saturated (i.e. S r = 1.0) conditions.For each simulation, the tempera- ture gradient between the thermal pipes at 37 °C and the boundaries of the sand box at 21 °C maintained constant for a minimum time of 48 h.The temperature calculated at the positions of the six temperature probes after the equalisation time of 48 h is compared against the experimental measurements taken by Akrouch et al. [11] with the probes B, C and D (Fig. 6a) and the probes E, F and G (Fig. 6b).Inspection of Fig. 6 suggests the good ability of the proposed numerical model in reproducing the experimental behaviour of the thermal piles tested by Akrouch et al. [11].
Finally, the present thermal conductivity model of deformable unsaturated soils can easily be coupled with a mechanical law [37, 38] and a soil-water retention law [39] to enable thermo-hydro-mechanical modelling of soils.Interestingly, this will allow for the consideration of timedependent behaviour (e.g.soil consolidation) and thermohydro-mechanical hysteresis as long as the present analytical model is coupled with a suitable thermos-hydro-mechanical law.This is however considered outside the scope of the present work and it constitutes matter for future research.

Conclusions
The present paper presented an empirical moisture-dependent thermal conductivity model for unsaturated and deformable soils.The model is generated by the product of two terms accounting for a) the effect of the degree of saturation on the moisture-dependent thermal conductivity at a unitary and constant void ratio and b) the effect of void ratio on the dry thermal conductivity.The proposed model has been calibrated and then validated against five different sets of experimental data obtained on soils with different granularities and equalised at various levels of both saturation and volumetric state.The validation of the model has been performed by two distinct methods: blind prediction and numerical simulations.Both methods provided a good agreement between the experimental measurements and the predicted thermal behaviour.The main outcomes of the present work can be summarised as follows: -The proposed model well predicts the thermal conductivity of various types of soils ranging from fine clays to coarse sands compacted at different dry densities ranging from 1465 kg/m 3 to 2310 kg/m 3 and equalised at different degrees of saturation (i.e.spanning from dry to fully saturated conditions).-The adopted validation procedure indicates that the proposed model is capable of extrapolating the thermal conductivity of soils for levels of degree of saturation not considered during the calibration range, as observed from both validation methods (i.e.blind prediction and numerical simulations).Future work will be directed to the assessment of the thermal of energy piles under different conditions of degree of saturation and porosity of the shallow geothermal piles as well as at different thermal conductivities of the energy piles.Finally, numerical simulations will be developed to assess the energetic performance of energy piles with an optimised geometry in view of maximising the heat exchanges with the surrounding soil.
Funding Open access funding provided by Università degli Studi di Genova within the CRUI-CARE Agreement.This research did not receive any specific grant from funding agencies in the public, commercial or not-for-profit sectors.

Fig. 1
Fig. 1 Relationship between the ratio λ dry λ dry,e=1 and the void ratio e

Fig. 2
Fig. 2 Calibrated and experimental moisture-dependent thermal conductivity for all data sets ) c =− ∇T (10) j e = e T 0 − T model parameters required to run the numerical simulations.Note that the value of the intrinsic permeability of 10 -12 m 2 has been assumed as a typical value for a sandy clay soil, as suggested by Gens et al. [36].

Fig. 5
Fig. 5 Model prediction against experimental data by Akrouch et al. [11]: probes B, C and D a and probes E, F and G b at a degree of saturation of 0.015

Fig. 6
Fig. 6 Model prediction against experimental data by Akrouch et al. [11]: probes B, C and D a and probes E, F and G b after the 48 h equalisation stage

Table 1
Numerical values of the model parameters for the different experimental data sets

Table 3
Numerical values of model parameters [11]ibration taken from Vanapalli et al.[40]as indicated by Akrouch et al.[11]for sandy clay P (J/kg K)