Numerical analysis of energy piles in a hypoplastic soft clay under cyclic thermal loading

This work is a numerical investigation of the effect of thermally induced volumetric collapse of normally consolidated clays on the performance of energy piles. A series of coupled thermo‐hydro‐mechanical Finite Element simulations were carried out using the commercial software ABAQUS. These examined a single free‐head energy pile embedded in a normally consolidated clay layer subjected to a constant mechanical load and to a number of heating/cooling cycles, to reproduce operating conditions. The soil behaviour was described with two advanced hypoplastic constitutive models for clays, one of which incorporates the thermally induced volumetric collapse using an ad‐hoc algorithm developed by the authors. Both models predict a cyclic accumulation of settlement and excess pore water pressure, especially when the thermal collapse effect is considered. While the excess pore pressure distribution stabilises within a few cycles, the rate of settlement of the pile head does not show any tendency to decrease from one cycle to another. These results are in agreement with data from small scale tests on an isolated energy pile in normally consolidated clay, indicating that the numerical model developed in this study can be used to investigate the complex soil/pile/raft interaction processes occurring in real piled foundations incorporating energy piles.

(see e.g. 1 ). This is possible because the soil beneath a relatively shallow boundary layer close to the ground surface can be considered as a reservoir at an almost constant temperature, regardless of the seasonal variation of the atmospheric temperature. During the hot seasons, the heat pump extracts the heat from the spaces to be cooled to transfer it to the soil at a significantly lower temperature than the outside air. The opposite will occur during the cold seasons. From an energy point of view, the basic operational principles of energy piles and other energy geo-structures are well known, and their adoption has been growing almost exponentially over the last two decades, due to their environmentally friendly features. Several review works on this topic are available in the literature. [2][3][4][5][6][7][8][9][10] From the mechanical point of view, energy piles are subjected to cyclic thermally induced internal forces that arise due to the interaction with the surrounding soil, and that add to the mechanical load transferred by the structure. This can significantly change the performance of the soil/foundation system and may have adverse effects on the serviceability of the supported structure.
With reference to a single pile, designers should focus on two main aspects: i) the evaluation of the thermally induced axial force distribution along the pile axis; ii) the assessment of the additional permanent settlement of the pile head caused by thermal cyclic loading.
The first aspect has been widely investigated. [11][12][13][14][15][16][17][18][19] Recently, Iodice et al. 20 proposed a set of analytical solutions to evaluate the design axial force induced by thermal loading to be considered in the geotechnical and structural Ultimate Limit State assessments. A method for defining the stiffness profile to be used in the calculation as well as a practical rule for the evaluation of the degree of restraint to be considered at pile head have been subsequently proposed by the same authors. 21 The behaviour of the soil-energy pile-structure system under normal operating conditions-that is, subjected to thermal cycles-has been mainly investigated experimentally by laboratory or field tests on instrumented energy piles. [22][23][24][25][26][27][28][29][30][31][32][33][34][35] An experimental study conducted on small scale models tested in a geotechnical centrifuge at 40g 36 has shown that an energy pile installed in lightly overconsolidated (OC) clay accumulates larger irreversible settlements if compared to a pile in heavily overconsolidated clay. In a series of small-scale pile loading tests at 1 g, Wu et al. 37 observed that a single energy pile subjected to heating and cooling cycles undergoes an accumulation of irreversible axial displacements at the pile head, progressing at an almost constant rate, while the soil experiences an increase and decrease of pore water pressures during heating and cooling stages, respectively, with progressively smaller increments of excess pore water pressure as the number of thermal cycles increases. These two studies represent the main experimental contributions concerning the cyclic behaviour of energy piles in soft clay in normal operating conditions, while there is a dearth of corresponding numerical studies 38,39 as the main contributions are mainly related to overconsolidated clay and sandy soils. 27,[40][41][42][43] The present work is a numerical study of the main phenomena controlling the cyclic behaviour of energy piles installed in soft clays, with particular reference to the pile performance in terms of accumulated settlements. Non-linear, fully coupled thermo-hydro-mechanical Finite Element simulations were performed on a single ideal energy pile installed in a normally consolidated (NC) clay layer, whose mechanical behaviour was described using an advanced inelastic constitutive model capable of capturing the main aspects of the coupled thermo-mechanical behaviour of the soil under cyclic mechanical and thermal loading conditions.

SUMMARY OF AVAILABLE EXPERIMENTAL EVIDENCE ON THE THERMO-MECHANICAL BEHAVIOUR OF CLAYS
A major feature of the thermo-mechanical behaviour of clays is the strong dependence of the volumetric response to thermal loading on the previous loading history, quantified by the overconsolidation ratio, in both monotonic and cyclic loading conditions. [44][45][46][47][48][49][50][51][52][53][54][55] While heavily overconsolidated clays dilate with increasing temperature (similarly to materials such as steel or concrete), soft NC or slightly OC clays experience volumetric contraction when heated in drained conditions, see Figure 1(A). Campanella & Mitchell 44 attributed this phenomenon to a decrease in the bond strength at the interparticle contacts, which in turn induces a rearrangement of the soil fabric with a net irreversible decrease of void ratio. Whatever the physical explanation, the experimental observations indicate that the thermally induced contraction stabilizes after a few cycles, possibly because additional contacts between the particles are established, see Figure 1(B).
The number of cycles required for the stabilization of the irreversible volumetric contraction increases with increasing plasticity index and voids ratio. 52,55 Another important experimentally observed effect of thermo-mechanical coupling is that both the preconsolidation pressure of a clayey soil and the position of its Normal Compression Line (NCL) in the void ratio-mean effective stress plane are affected by temperature changes. 44,49,50,51,55,56 In particular, the preconsolidation pressure decreases with increasing F I G U R E 1 Volumetric strains of NC and OC clay subjected to temperature variations in drained conditions. temperature and the NCL translates downward. On the contrary, temperature variations have a relatively small effect on the critical state strength of the soil. 47,48,50,[56][57] In undrained conditions, Campanella and Mitchell 44 have shown that the pore water pressure-temperature relationship has a tendency to be hysteretic from the early temperature cycles. A similar behaviour has been observed by Plum & Esrig, 58 who performed cyclic heating and cooling tests in a triaxial apparatus in undrained conditions on a low plasticity clay.
Finally, an increase in both the coefficient of consolidation and the hydraulic conductivity with temperature has been reported by Towhata et al., 59 Burghignoli et al., 48 Abuel-Naga et al., 60 Di Donna and Laloui. 55 All the features of the observed thermo-mechanical behaviour of clay soils can be reproduced qualitatively and, in most cases, quantitatively by the advanced constitutive models adopted in the numerical simulations of the energy pile behaviour presented in the following sections.

Details of the pile model
The problem considered in the numerical simulations consists in the mechanical and thermal loading of a single energy pile, 25 m long and with a diameter of 0.5 m, installed in a 50 m deep homogeneous soft NC clay layer (see Figure 2). As the focus of the numerical investigation is the response of the pile to thermal cycles, the pile installation effects were not considered, assuming the pile as wished in place.
In the first stage of the numerical simulations, the pile was kept at constant ambient temperature and was subjected to a monotonically increasing axial load, applied in drained conditions. In the simulations, two different values of the axial load, Q, were considered, namely 0.3R t and 0.6R t , where R t is the bearing capacity of the pile, evaluated using the β method for the shaft resistance (β = 0.25, see, e.g., 61 ) and by using classical bearing capacity factors for pile base resistance (N q = 5). 62 In the second stage of the simulations, the axial load was kept constant, while the pile temperature was varied cyclically, for a total of five full thermal cycles. Each cycle consisted of a 6-month phase in which the pile temperature was rapidly increased (heating, h) and then maintained constant for 5 months (continuous heating, hf), followed by a thermal rest period of 1 month (hr) in which only the temperature at the pile head was rapidly reduced back to the ambient temperature and then kept constant. The subsequent 6-month phase is characterized by a rapid decrease of the pile temperature (cooling, c), subsequently kept constant for 5 months (continuous cooling, cf), followed by a second thermal rest period of 1 month (cr), in which the temperature at the pile head was rapidly increased back to the ambient temperature and then kept constant. The details of the thermal cycle are provided in the inset of Figure 2.

F I G U R E 2
Problem under investigation (T 0 = initial temperature; pp = pore water pressure; γ w = unit weight of water; x and y = radial and vertical coordinates; u x and u y = displacements along x and y).
The time and space evolution of the displacement, pore water pressure and temperature in the soil-pile system during the mechanical and thermal loading stages were evaluated by axisymmetric, non-linear, fully-coupled thermo-hydromechanical FE simulations performed with the FE code ABAQUS Standard 6.14. 63 In the FE model, the pile was discretized using 1000 thermo-mechanical quadrilateral elements (CAX4T), while the soil was modelled using 27956 three-field elements with nodal displacement, pore water pressure and temperature degrees of freedom (CAX4PT). Further details concerning the FE discretisation of the domain, the pile properties, the thermal and mechanical load paths along with the selection of parameters for the advanced constitutive models can be found in Iodice et al. 21

3.2
Constitutive models adopted for the soil Over the last decades, a relatively large number of constitutive models have been proposed to describe the non-isothermal behaviour of soils within the framework of the classical theory of thermo-elastoplasticity, see, for example, the works by Cui et al., 64 Laloui & Cekerevac, 65 Abuel-Naga et al., 53 Laloui and Francois. 66 These isotropic hardening plasticity models are capable of reproducing the non-linear, coupled response of soils under combined mechanical and thermal monotonic loading conditions, but fail to capture the accumulation of irreversible deformations and pore water pressure build-up under cyclic loading conditions. The adoption of tensor-like internal variables capable of providing enough memory of the past loading history, as in kinematic hardening plasticity, or the extension of the constitutive framework to Bounding Surface plasticity 67 or Generalized Plasticity 68 may remove this limitation. An alternative approach, pursued in this work, consists in formulating the constitutive equations in rate-form within the framework of the theory of hypoplasticity with internal state variables, 69-71 the main features of which are: the absence of any elastic domain in stress space; the absence of any kinematic decomposition of the strain rate tensor into an elastic, reversible part and a plastic, irreversible part; the incrementally non-linear character of the constitutive equations, in which the tangent stiffness tensor is a continuous function of the strain rate direction.
Within this framework, Mašín 72,73 developed a constitutive equation for fine grained soils incorporating the main features of Critical State Soil Mechanics and capable of reproducing the cyclic behaviour of clays with the introduction of the intergranular strain, 74 that is, a second-order tensorial internal variable, which keeps track of the recent strain history of the material. The model was subsequently improved to incorporate explicitly defined asymptotic states 75 and anisotropic small-strain stiffness 76 . This last version of the model-from now on referred to as the "Hypo" model-was selected as the first option to model the NC clay layer in the numerical simulations of the energy pile. The Hypo model can reproduce, both from a qualitative and a quantitative point of view, the main features of the monotonic response of clay soils-such as pressure dependence, history dependence, and the existence of critical state conditions for extreme shear deformations-as well as their cyclic/dynamic response to strain levels of different amplitudes (see, e.g. 77 ). The constitutive equations of the Hypo model are summarized in Appendix A.
In the Hypo model, the thermo-mechanical coupling is limited to the thermal expansion/contraction of the solid skeleton associated to a temperature increase/decrease. Recently, the model has been further enhanced by Ma et al. 78 to improve the description of thermo-mechanical coupling effects and allow the accumulation of irreversible contractive volumetric strains in NC conditions, using the shakedown concept, originally proposed for cyclic mechanical loading by Mašín and Khalili. 79 This model-from now on referred to as the "Hypo-T" model-was selected as the second option to model the thermo-mechanical response of the clay in the present study. The relevant constitutive equations for the Hypo-T model are summarized in Appendix B.
The Hypo-T model response to thermal loading is based on the following key assumptions: i) a Thermal Stabilization Line (TSL)-equations (B4) and (B5)-exists in the bilogarithmic compression plane whose slope is smaller than the slope * of the Normal Compression Line (NCL, see Figure 3); ii) the position of both the NCL and TSL in the ln(1 + ) ∶ ln( ) plane changes with temperature according to equations (B3)-(B5). In particular, they translate downward for increasing temperature, . The accumulation of irreversible volumetric contraction during thermal cycles continues with decreasing rate until the soil state reaches the TSL; iii) during cooling stages, there is no irreversible change in void ratio and therefore the volumetric strains are completely recoverable.
A consequence of the last assumption is that during cooling the soil skeleton contracts at constant void ratio, since the volume change of voids and of the solid skeleton are proportional. 79,80 Figure 3 reports the qualitative behaviour predicted by the Hypo-T model for a NC soil specimen subjected to heating and cooling cycles at constant mean effective stress.
In the initial isothermal condition ( = in ), the NCL and TSL intersect at the state 0. When the soil is heated to its final temperature ( = f in ), the TSL shifts downwards more than the NCL and the soil state remains on the NCL (point 1). Therefore, a contractive volumetric deformation is accumulated which corresponds to the reduction in void ratio Δ ln (1 + ) 1 . If the temperature is then decreased back to the initial value (point 2), according to assumption (iii), no additional strain is accumulated and both the NCL and the TSL return to their original position. In the subsequent heating cycle, the TSL and the NCL shift downward again, and the behaviour remains elastic until the TSL crosses the soil state at point 2. After that, additional contractive volumetric strains accumulate, quantified by Δ ln (1 + ) 2 < Δ ln (1 + ) 1 . The soil state is dragged downward and its distance from the TSL reduces (point 3). In the following second cooling stage, no further strains develop (point 4). During the subsequent heating, volumetric contraction continues to build-up but at a decreasing rate. After a few cycles, the soil state reaches the TSL and the response stabilizes with no additional permanent volumetric strain accumulation. It follows that, in a cyclic thermal loading with constant amplitude, the final temperature controls the position of the TSL and thus the final steady-state void ratio.

Model calibration
The isotropic version of the Hypo model detailed in Appendix A is fully characterized by the 13 material constants listed in Table 1 The second group of constants are related to the incorporation of the intergranular strain concept to provide the material with a certain amount of memory of the previous deformation history and of its effects on the soil stiffness at small strains. In particular: constant defines the size of the 'pseudo-elastic' range in the strain space; and control the shear modulus of the soil at very small strains, 0 , and its dependence on mean effective stress, respectively; , and are exponents controlling the evolution of the intergranular strain tensor with accumulated strain-equation (A4)-and the effect of the ratio = ∕ on the linear and non-linear parts of the stress rate-equations (A2) and (A3); finally, defines the ratio of the small strain stiffness for orthogonal loading conditions (strain rate orthogonal to ) and for full strain reversal (strain rate and pointing in opposite directions). Constants and can be calibrated from direct measurements of the small strain shear modulus 0 , for example by means of resonant column or bender element tests, and its dependence on mean effective stress (see, e.g. 81 ). Constants , and can be determined by trial and error, fitting available modulus decay curves obtained in static conditions from TX-CU tests performed with accurate local strain measurements or from resonant column tests. Constant controls the size of the 'linear' very small strain region, while and control the shape of the ∕ 0 curve. In principle, calibration of constant would require access to smallstrain stiffness measurements performed upon a full stress path reversal and for a 90 • rotation of the stress path direction starting from the same state, as in the experiments by Atkinson et al. 82 In most practical situations, 78 this kind of data is not available; however, usually varies in a relatively narrow range, typically between 0.5 and 0.7 71 . Constant influences the amount of strain/pore water pressure accumulation in drained/undrained cyclic mechanical loading and prevents the model predictions from displaying excessive ratcheting. Again, the values of vary in a relatively narrow range, typically from 8 to 10. 71,83 TA B L E 1 Material constants adopted for the Hypo and Hypo-T models in the numerical simulations.

Constant
Set The Hypo-T model described in Appendix B is fully characterized by the same constants of the first two groups plus a third group that characterizes the thermo-mechanical response of the clay in non-isothermal conditions, see Table 1, Set #2. In particular, and control the translation and rotation of the NCL with temperature; is the slope of the TSL, which should be larger than the slope of the unloading-reloading lines, * 78 ; and control the accumulation of irreversible volumetric contraction during heating/cooling cycles. Ma et al. (2017) noted that the effect of temperature on the slope of the NCL is usually negligible, that is constant T can be set to zero. Based on experimental results by Campanella and Mitchell, 44 Vega and McCartney, 54 Di Donna and Laloui, 55 typical values of lie in the range between 0.4 and 0.5. 78 Finally, a sensitivity analysis performed by Ma et al. 78 has shown that the accumulation of irreversible strain stabilizes within five thermal cycles when is set to 0.1. This is in accordance with the experimental results by Vega and McCartney. 54 For the numerical simulations of the energy pile performance, the two sets of constants #1 (Hypo model) and #2 (Hypo-T model) were adopted. Note that Set #1 is derived from a calibration performed by Mašín 71 for London clay. Constant , not available in this source, was assigned according to the indications provided by Wegener & Herle, 83 while and were chosen based on the results of laboratory tests on reconstituted NC London clay by Viggiani and Atkinson. 81 As no non-isothermal tests were available for London clay, the constants of the third group of the Hypo-T model were assigned within the ranges suggested by Ma et al. 78 The complete characterisation of the numerical model discussed in Section 3.1 requires the definition of additional soil properties not related to the constitutive models, namely: the saturated unit weight , the soil hydraulic conductivity , the coefficients of thermal expansion and for the solid grains and the pore water, respectively, the soil thermal conductivity th , the soil specific heat capacity , and the ambient temperature T 0 .
The values adopted in the FE simulations for these constants are provided in Table 2, where the data for Set #1 (Hypo model) and Set #2 (Hypo-T model) are the same. In principle, the thermal expansion coefficient should be considered a known function of temperature. This would introduce an additional source of strong non-linearity in the problem formulation. In view of the relatively small temperature variations expected during the energy pile operations and to reduce the computational effort, in all the analyses the thermal expansion coefficient of water was kept constant and equal to the value at T 0 . Finally, the physical, mechanical and thermal properties for the pile material are provided in Table 3, where γ CLS is the unit weight, E CLS the Young's modulus, ν CLS the Poisson's ratio, α CLS the linear coefficient of thermal expansion, λ th,CLS the thermal conductivity, and c p,CLS the specific heat capacity. TA B L E 2 Physical, hydraulic and thermal soil properties adopted in the numerical simulations.

TA B L E 3
Physical, mechanical and thermal properties adopted for the linear thermo-elastic pile.

Numerical implementation of the Hypo and Hypo-T models
For the implementation of the two models into the chosen FE platform, a stress-point algorithm consistent with the evolution equations provided by the constitutive model must be defined and incorporated into the code to: i) update the stress and additional state variables from their known values at the beginning of the generic time step [ , +1 ] to their final values at the end of the step ( = +1 ) for given increments of strain, Δ +1 ≔ +1 − , and temperature, Δ +1 ≔ +1 − ; ii) evaluate the consistent Jacobian matrices: required in the iterative solution of the discrete equilibrium equations cast in monolithic form via Newton's method. ABAQUS Standard allows to incorporate a stress-point algorithm using a standard user interface for user-defined constitutive models, that is a UMAT routine. A UMAT routine is already available for the isothermal version of the Hypo model, on the SoilModels.com website. 84 This is based on an explicit adaptive Runge-Kutta-Fehlberg algorithm of order 2/3 (RKF-23 algorithm), with adaptive substepping and error control, similar to the ones employed by Sloan 85 for classical plasticity and Tamagnini et al. 86 for sand hypoplasticity. In this algorithm, the main idea is to estimate the sub-step size that would provide the update sought with the desired level of accuracy by comparing two solutions obtained with two Runge-Kutta algorithms of different order (in this case, the second and third orders).
For the thermo-hypoplastic model Hypo-T, the standard RKF-23 algorithm had to be modified to account for the higher level of thermo-mechanical coupling existing in the constitutive equations. In particular, in this case, the evolution equation for consists of two distinct contributions, the first one depending on the mechanical strain rate and the second on the temperature rate, as immediately apparent from Equations (B1), rewritten here as Equations (2) for convenience:̇=̇+̇( 2a) and (̇) is the Heaviside step function, equal to 1 iḟ> 0 and equal to 0 otherwise. Let us now consider the problem of integrating the evolution equations by partitioning the time axis in several given time steps, of size Δ +1 = +1 − . The rate-independent nature of the constitutive equations allows a rescaling of the time axis that is useful for the setting of the adaptive sub-stepping integration scheme. Let: be a non-dimensional time factor taking values = 0 at = and = 1 at = +1 . It is then easy to show that Equations (2) can be recast in the following alternative format: with initial conditions: In Equations (5), the quantities: that is the (prescribed) mechanical strain and temperature increments over the time step, represent the two driving mechanisms of the evolution process. Note that in the evaluation of the two tensors and of Equation (5a), the identitieŝ ⋅ (̇−̇) =̂⋅ Δ +1 ∕Δ +1 and  (̇) = (Δ +1 ) must be used to replace time rates with the prescribed increments over the step.
By defining the vectors: the evolution problem can be finally recast in the standard first order ODE format: with initial conditions: Box 1. RKF-32 algorithm.
In the explicit, adaptive RKF-23 method, the normalized time step is split into smaller sub-steps [ , +1 ] with size: and two distinct explicit Runge-Kutta methods of order 2 and 3, respectively, are used to update the solution within the generic sub-step, from its initial known value at time to its final value at time +1 . In particular, the second and third order updates are given by:̃+ in which:̃1 = 0,̃2 = 1,̂1 = 1∕6,̂2 = 2∕3,̂3 = 1∕6 and: A suitable norm of the difference in the two solutions provided by Equations (12) and (13), +1 , can be compared to a prescribed tolerance TOL. If +1 < TOL, the sub-step is accepted and the ratio TOL∕ +1 can be used to extrapolate a larger sub-step size for the next sub-step; otherwise, the two updates are rejected and a smaller sub-step size is evaluated using a different extrapolation formula based again on the ratio TOL∕ +1 . The algorithmic procedure is illustrated in Box 1. The integration procedure stops when the elapsed time within the sub-step equals the step size or when the maximum required number of sub-steps is attained. At the end of each completed step, the consistent Jacobian matrices of Equation (1) are approximated with the corresponding continuum tangent operators: in which the tensors and are provided by Equations (A6) and (A3). For relatively large step sizes, this approximation may produce the loss of the quadratic convergence of the global Newton iterations. In our case, preliminary numerical simulations, confirmed by the results of the simulation program discussed in the next Section 4, have shown that the convergence of the global iterative solution was still sufficiently fast, and therefore the approximations of Equation (15) were acceptable.

Testing of the implementation of the Hypo-T model
The implementation of the model was validated against the predictions made by Ma et al. 78 for a drained isotropic consolidation test with thermal cycles on a NC reconstituted Illite clay, performed by Campanella and Mitchell. 44 In the first part of the test, the specimen had been subjected to an isotropic consolidation stage up to a mean effective stress of 200 kPa, at constant ambient temperature of about 18 • C. After the mechanical compression stage, the specimen temperature had been varied over a wide range, from about 4 • C to about 60 • C with a total of three heating and cooling cycles, in drained conditions. The numerical simulation of this test was carried out using the constants of Set #3 of Tables 1 and 2. These constants correspond to the ones adopted by Ma et al. 78 in their simulation of the laboratory test results.
The results of the numerical simulations are compared with the experimental data and the predictions by Ma et al. 78 in Figure 4, in terms of volumetric strain as a function of temperature .
The results show that the proposed implementation matches reasonably well the predictions made by Ma et al. 78 The slight differences in the two predictions may be due to fine details of the integration algorithms adopted to match the solution in time for the prescribed loading path.

RESULTS OF FE SIMULATIONS OF ENERGY PILE PERFORMANCE
A selection of the results obtained from the FE simulations program detailed in Section 3 is reported here focusing first on the overall performance of the energy pile in terms of pile head displacements, and then looking in detail to the effects F I G U R E 5 Dimensionless global settlement after five thermal cycles.
that the evolution of temperature with time has on the volumetric deformations, pore water pressures and effective stress state at several points along the soil-pile interface. Figure 5 shows the evolution with time of the vertical displacement at the pile head over a period of 60 months (five full cycles of the loading program of Figure 2). Figure 5(A) refers to simulations performed with a stationary load at the pile head = 0.3 ; Figure 5(B) to simulations performed with a stationary load at the pile head = 0.6 . In isothermal conditions, the response of the pile predicted by the two models is clearly the same. During the cyclic thermal loading stage, the overall downward displacements predicted by the Hypo-T are larger than those predicted with the Hypo model. This is a consequence of the permanent positive volumetric deformations (thermal collapse) accumulated in the NC clay over the thermal cycles. For both load levels considered, the accumulation of irreversible displacements increases almost linearly with the number of cycles, without showing any tendency to stabilize-that is without any significant reduction in the accumulated settlement rate.

Predicted behaviour of the soil-pile system in terms of pile head displacements
As expected, in the simulations with the Hypo model the largest settlement increment occurs during the first cycle. In particular, for the smaller static load at the pile head ( = 0.3 ) the ratio between the settlement after five thermal cycles and the one at the end of the first thermal cycle are between 4.3 (Hypo model) and 7 (Hypo-T model). A similar trend is obtained for the larger static load at the pile head ( = 0.6 ). In this case, due to the higher stress level in the soil and the stronger non-linearity in its mechanical response, the magnitude of the thermally induced displacements is larger, and this is particularly evident for the Hypo-T model.

Predicted response of the soil at the material point scale
It is interesting to note that at the very beginning of the thermal loading stage the response predicted with the two models is almost identical. To understand this result, it is useful to clarify some fundamental aspects of the behaviour of a saturated soil subjected to a temperature increase under both drained and undrained conditions. In the absence of volumetric collapse, solid grains experience a volumetric deformation equal to Δ . If the thermal expansion coefficient of the water were the same as that of the solid grains (i.e. = ), the water would merely expand in the same manner as the soil grains and therefore there would be no thermally induced changes in pore water pressures. On the contrary, since > , the water expands more than the void volume due to the same temperature increase. In drained conditions, this would produce a net water flow leaving the heated soil volume, quantifiable as ( − )Δ per unit soil volume, where is the soil porosity. However, if the soil hydraulic conductivity is so low that, at least at the beginning of the heating process the relative motion between the solid skeleton and the pore water can be considered F I G U R E 6 Volumetric strain versus temperature (Q = 30% R t ).
negligible (undrained conditions), to guarantee the volumetric deformation compatibility between the two phases, the pore water would experience an increase in pore pressure proportional to − . In the presence of thermally induced volumetric collapse, the soil skeleton contracts for a positive increase in temperature. Let this additional (positive) volumetric deformation be quantified by the quantity Δ , where the coefficient measures the apparent volumetric contraction of the solid skeleton. Volumetric collapse is due to a reduction in porosity resulting from the rearrangement of the grains, still experiencing an expansion equal to α s ΔT. It follows that, in drained conditions more water will flow out of the heated soil volume as compared to the case in which no thermal volumetric collapse occur, while in undrained conditions a larger excess pore water pressure will develop. As a consequence, the effect of volumetric collapse on the pore water flow or pore water pressure is the same as that of a larger in both drained and undrained conditions. It can be shown that the apparent water thermal expansion coefficient, which would take implicitly into account the effect of volumetric collapse is given by ,app = + ∕ . Figure 6 reports the computed volumetric strain as a function of temperature at four selected locations along the pilesoil interface (namely at = 6.25 m, 12.50 m, 18.75 m and 25.0 m-this last being at the pile tip), for the pile loaded at = 0.3 . The plots on the left provide the results obtained with the standard Hypo model, while those on the right provide the results obtained with the Hypo-T model. As the temperature is increased, only the soil in the vicinity of the pile experiences a temperature variation. As the heating stage is of very short duration as compared to the characteristic time of the consolidation process, the soil deformations in this stage occur in undrained condition. As a consequence, volumetric expansion is observed during the heating stage in the predictions obtained with the Hypo model, which superimposes to F I G U R E 7 Excess pore pressure (Q = 30% R t ).
the pre-existing contraction of the mechanical loading phase. The same response is observed in the predictions obtained with the Hypo-T model, the only difference being in the larger excess pore water pressure developed during the heating stages. This is clearly visible in Figure 7, which provides, for the same couple of simulations, the evolution with time of the excess pore water pressures computed at the same locations along the pile-soil interface, with the Hypo (plots on the left column of Figure 7) and Hypo-T models (plots on the right column of Figure 7).
During the subsequent continuous heating stage, the increase of temperature extends to larger portions of the soil around the pile, while the excess pore water pressures developed at the soil-pile interface-which now is at constant temperature-start dissipating because of the simultaneously occurring hydrodynamic consolidation process. In this stage, larger volumetric compactions are predicted by the Hypo-T model at the soil-pile interface ( Figure 6) due to the higher excess pore water pressures accumulated during the heating stage.
A different response is observed in the two simulations during the thermal rest stage after the first cooling (i.e. the last step before the end of the first cycle). Indeed, the simulation with the Hypo model predicts soil expansion due to the temperature increase, contrasted by a contraction of smaller magnitude caused by the concurrent dissipation of the excess pore pressure. As a result, the overall volumetric deformation is characterized by a slight expansion at any depth. In the simulation with the Hypo-T model, the dissipation of the larger excess pore water pressure together with the volumetric collapse induced by the thermal cycle result in an overall contraction of the soil.
The other stages of the first cycle are characterized by contraction. In particular, during the continuous cooling stage, in which all the elements are at constant temperature, negative excess pore water pressures are developed which increase in absolute terms. Figure 8 shows the distribution of pore water pressure in a region of soil close to the pile at a depth of 6.25 m, at different times during the continuous cooling stage of the first cycle (cooling times equal to 1 day, 5 days, 18 days, 2 months and 5 months, respectively). It can be observed that, in the region close to the pile, the pore water pressure tends to decrease with increasing distance to the pile-soil interface, up to a distance of about 2 pile diameters from the interface. As a consequence, water is forced to flow outwards of the soil-pile interface.
As far as the response of the pile-soil system during the thermal cyclic loading is concerned, both models predict a sequence of open loops in the volumetric strain versus temperature plots (Figure 6), with a continuous accumulation of contractive volumetric strains at any depth, in parallel with a progressive accumulation of excess pore water pressure at a decreasing rate. The pore water pressure accumulation stabilizes within about 3 cycles in the simulation with the Hypo-T model (Figure 7).
In addition, the difference in the thermal expansion coefficients of the pile material and of the soil causes the development of shear stresses at the pile-soil interface which, in turn, could determine significant modifications in the distribution of the axial forces within the pile shaft, as well as in the stress state of the surrounding soil. Figure 9 shows the stress paths experienced by the same soil elements at the soil-pile interface considered in Figures 6  and 7, as predicted with the Hypo model (plots on the left column of Figure 9) and the Hypo-T model (plots on the right column of Figure 9) in the simulations performed with = 0.3 . A common feature of the stress paths at all points located on the pile shaft (z = 6.25 , 12.50 and 18.75 m) is the progressive reduction of both mean effective stress and deviator stress as the number of thermal cycles increase, for both models considered. In the simulation with the Hypo model, the largest reduction in is observed in the first cycle, while the variations in mean effective stress tend to reduce in size as the number of cycles increases. In the simulation with the Hypo-T model, the activation of the thermal collapse term in the model formulation leads to a larger reduction of mean stress in the first two cycles. In agreement with the observations made on the time evolution of excess pore water pressures (Figure 7), the maximum reduction in takes place during the first rest stage after cooling.
At the soil element located at the pile tip, a net increase of both mean effective stress and deviatoric stress with respect to the geostatic stress condition is predicted in the simulation with the Hypo model at the end of the fifth cycle. A slight decrease in is observed at the same time station in the simulation with the Hypo-T model.
The observations made at the element scale highlights that the stress paths which the soil experiences at all the material points considered are largely different to those imposed to soil specimens in triaxial experimental campaigns such as the one reported by Campanella and Mitchell 44 in which the soil sample was isotropically consolidated in a triaxial apparatus and then subjected to cyclic temperature variations under controlled drainage conditions. These perfectly drained or undrained conditions imposed in laboratory tests do not apply, if not for a very limited time span, to the soil surrounding the energy pile, which is subjected to a complex, non-stationary coupled deformation and flow (of mass and energy) process during the entire operational life of the foundation.
These observations may explain why the results presented in previous Section 4.1, in terms of the evolution of pile head displacements with time, do not show any tendency for the pile head settlements to stabilize, while in laboratory tests F I G U R E 9 Stress paths (Q = 30% R t ).
the stabilization of volumetric deformations or excess pore water pressure occurs within a few thermal cycles. In this respect, it is worth noting that the predicted pile performance at the global scale is consistent with available experimental observations on energy pile prototypes. This aspect is explored in detail in the following Section 5.

Long-term pile head displacement
To explore the long-term pile behaviour, the numerical simulations with Q = 30% R t were extended to consider 50 operational cycles. Figure 10 shows that the settlement accumulation progresses without stabilization after 50 years of regular operation of the energy pile. The settlements predicted using the Hypo model continue to increase at the same constant rate of the first 5 cycles (about 0.3 %/year). On the contrary, the rate of settlement obtained using the Hypo-T model is characterized by an initially higher value (about 0.6 %/year in the first 5 cycles), which then progressively reduces during the first 15 cycles to a final value of about 0.4 %/year. The thermal collapse effect is responsible for the differences observed in the initial part of the simulations, as well as for the different final settlement rate computed in the two cases.
In order to assess which model better captures the real behaviour of energy piles, data from long duration full scale or laboratory tests would be needed. Unfortunately, due to the difficulties in performing experiments involving many thermal cycles, most of the data available are related to only few years of operation. For the short duration test performed by Wu et al., 37 the numerical predictions are in good qualitative agreement with the experimental data (see Figure 11).

COMPARISON OF THE GLOBAL PERFORMANCE OF THE ENERGY PILE WITH AVAILABLE LITERATURE DATA
Wu et al. 37 performed a series of pilot tests on small-scale floating energy piles (of diameter = 23 mm and length = 450 mm) installed in a NC saturated clay layer under natural gravity. The Authors considered three different testing layouts: i) single energy pile with an adjacent non-energy pile without a cap (EP-F); ii) energy pile with an adjacent non-energy pile connected by a cap (EP-R); and iii) isolated energy pile (EP-S).
In the tests, the energy pile was first axially loaded up to about 40% of the estimated pile axial bearing capacity, and then subjected to five thermal cycles between 14 • C and -13 • C at constant axial load. Although the test results cannot be considered representative of the behaviour of energy piles at full scale, due to the impossibility to achieve the condition of complete similarity between the model pile and a prototype under normal gravity, they are nonetheless very useful from a qualitative point of view. Figure 11 reports the measured pile head displacements, normalized with respect to the pile diameter , as well as the measured excess pore pressures at four different points within the soil, as indicated in the inset of Figure 11(B).
The data reported in Figure 11(A) show that each thermal cycle gives rise to an accumulation of permanent axial displacements which continues even after the end of the fifth thermal cycle. This response is qualitatively very close to the results of the numerical predictions reported in Figure 5. The final accumulated displacement at the end of the thermal cyclic loading is about 1.5% of the pile diameter. The experimental results in Figure 11 show also that the presence of a thermally inactive pile (EP-R test) causes a strong reduction (about 40%) of the accumulated permanent settlement. This behaviour is also confirmed by Rotta Loria and Laloui 87 who reported an increase of about 158% in the pile head settlement in case all the piles in the foundation system are thermally activated.
The accumulation of cyclic settlements is also confirmed by the results of centrifuge tests carried out by Ng et al. 36 for a single pile installed in a lightly overconsolidated clay (OCR = 1.7). On the other hand, negligible settlements after the first cycle were obtained by Di Donna and Laloui 38 from numerical analyses employing an advanced constitutive model capable of modelling the cyclic response of the soil. In the analyses by Di Donna and Laloui, 38 the pile was first subjected to mechanical load and then to 10 years of thermal cycles to explore its long-term performance. The maximum and minimum temperature variations applied to the pile were very similar to that employed in the present study, about +14 • C and-13 • C, respectively. The main difference between the analyses reported here and by Di Donna & Laloui 38 is in the adopted values of soil hydraulic conductivity, k w = 10 −7 m/s, and in the applied heating and cooling rates, 0.47 • C/day and 0.43 • C/day, respectively. During the thermal cycles, the pile showed an average settlement slightly larger than that experienced at the end of the mechanical loading phase, with no accumulation. A possible explanation of this apparent inconsistency could be that, because of the relatively higher soil hydraulic conductivity, the deformation of the solid skeleton occurs in drained conditions, due to the very slow heating rate adopted in the simulations.
From the evolution of the pore pressures for EP-F test by Wu et al. 37 reported in Figure 11(B), it can be noted that no significant accumulation of excess pore pressure is observed as the number of cycles increases. The results of the numerical simulations performed with an axial load = 0.3 ( Figure 12) show a somewhat different trend for the points located in the upper portion of the pile shaft (down to = 12.50 m). This could be due to: i) the different drainage conditions adopted in the small scale model tests and in the numerical simulations; ii) the different temperature-time histories considered in the model tests and in the numerical simulations; iii) the presence of the thermally inactive pile in the vicinity of the energy pile in the test EP-F, which is not included in the FE models adopted in this work; and, iv) the adoption of a constant coefficient of water expansion in the numerical simulations, in place of a more realistic temperature-dependent function, which would allow to better reproduce the hysteretic response in the excess pore pressure development from the beginning of the thermal cyclic loading stage.
An important role in the development of excess pore water pressures in the soil surrounding the energy pile is played by the combination of soil hydraulic conductivity which controls the duration of the consolidation process-and the imposed heating rate. In the simulations performed in this work, a heating rate of 14 • C/day was considered, in a soil with a hydraulic conductivity of 10 −10 m/s. To study the effect of heating and cooling rates on the response of the pilesoil system, additional analyses were carried out, with the Hypo-T model, using the following heating/cooling rates: (1) 0.47 • C/day, that is the pile temperature was increased in 30 days followed by a continuous heating phase in which the temperature was maintained constant for 4 months, and by a thermal rest period of 1 month, after which the path was reversed; (2) 1 • C/day, that is the pile temperature was increased in 14 days followed by a continuous heating phase in which the temperature was maintained constant for 4.5 months, and by a thermal rest period of 1 month, after which the path was reversed; (3) 28 • C/day, that is the pile temperature was increased in 12 h followed by a continuous heating phase in which the temperature was maintained constant for about 5 months, and by a thermal rest period of 1 month, after which the path was reversed. The results reported in Figure 13 suggest that the excess pore pressure predicted using such a wide range for the heating/cooling rates are essentially the same. Also, the displacement profiles reported in Figure 14 indicates a slight increase of the settlement when higher heating/cooling rates are employed, yet their trend is quite similar.
It is worth mentioning that, if in the anlyses a much higher hydraulic conductivity (i.e. = 10 −7 m/s) is used, the thermal loading at low heating/cooling rates results in almost negligible predicted excess pore water pressures, and the soil deformation can be considered to occur in drained conditions. F I G U R E 1 2 Excess pore pressure versus temperature (Q = 30% R t ).

CONCLUSION
To achieve a better understanding of the performance of energy piles installed in soft fine-grained soils under normal working conditions-that is subjected to thermal cyclic loading-a numerical study was carried out on an ideal, but realistic, floating energy pile installed in a deep NC clay layer. Two advanced hypoplasticity models for fine-grained soils were considered, one of which (the Hypo-T model) incorporating explicitly the thermo-mechanical coupling observed experimentally upon thermal cyclic loading. An explicit algorithm with adaptive sub-stepping and error control was developed for the implementation of this model in ABAQUS Standard.
The main outcomes of this work can be summarised as follows: -the overall performance of the energy pile predicted by the two hypoplastic models is qualitatively similar, as both models are capable of reproducing the progressive accumulation of pile head displacements as the number of thermal cycles increases both in short and long-term conditions; -the temperature changes in the soil surrounding the pile generate excess pore water pressures that tend to accumulate with increasing number of thermal cycles; at the same time, excess pore water pressure dissipate due to the concurrent consolidation process; F I G U R E 1 3 Excess pore pressure for different heating/cooling rates using the Hypo-T model.

F I G U R E 1 4
Dimensionless global settlement after five thermal cycles for different heating/cooling rates using the Hypo-T model.
-in the simulations with the Hypo-T model, due to the low soil hydraulic conductivity considered, thermal volumetric collapse of the NC clay is prevented during the heating stages of each cycle; as a consequence, larger excess pore water pressures are predicted than in the simulations performed with the standard Hypo model; -when used to simulate laboratory test results on clay specimens subjected to thermal cycles, the Hypo-T model predicts volumetric strains that tend to stabilise after a few cycles in drained condition, while in undrained conditions the model predicts an accumulation of excess pore water pressure during the first few cycles followed by stabilization. This behaviour is not reproduced at the boundary value problem scale, as the predicted permanent vertical displacements at the pile head tend to increase with the number of cycles with no sign of stabilisation. This feature of the global response of the soil-pile system might be attributed to the following two factors. First, during the thermal cycling loading, some excess pore pressure dissipation occurs due to the concurrent consolidation process which, in turn, generates soil deformations and pile head settlements; second, in the energy pile problem, soil elements are subjected to stress changes generated by the complex soil-pile interaction processes which are significantly different from those imposed in the laboratory tests reported in the literature; -in view of the previous considerations, it would be desirable to perform model calibration based on experimental results reproducing more closely the expected stress paths for the specific problem at hand (in this case, the response of the energy pile to thermal cyclic loading); -the tendency of the single energy pile to accumulate permanent vertical settlements at the pile head, which could be as large as about 9% of the pile diameter over five thermal cycles only (in the simulation performed with the Hypo-T model at an axial load level of 60% of the pile bearing capacity), is probably the most important result of the present study and has important implications for the assessment of the foundation safety with respect to Service ability Limit State; -the lack of stabilisation of accumulated permanent displacements predicted during cyclic thermal loading is not related to the specific choice of the constitutive model for the soil; both models considered provide qualitatively similar responses, although the capability of the Hypo-T model of reproducing the effects of thermal volumetric collapse is responsible for a significant increase of the accumulated displacement rate; -as far as pile head displacements are concerned, the results obtained from the numerical simulations appear qualitatively comparable to the experimental observations made by Wu et al. 37 on small scale models of energy piles subjected to thermal cyclic loading. Note that this comparison can only be qualitative due to the impossibility of 1 g small scale experiments on model pile foundations to properly reproduce the stress level in the soil.
Although the numerical simulations performed in this work refer to a single free-head energy pile, they can be considered relevant for the behaviour of a pile group in which all the piles are thermally activated and sufficiently spaced to avoid thermal interference. The presence of inactive piles connected to the same cap may significantly reduce the magnitude of the irreversible pile head displacements. In this case, the development of an additional axial force at the head of the energy piles would be balanced by the development of axial forces of opposite sign at the head of the inactive piles. As a consequence of the above findings and previous literature studies, the latter configuration is the only design option when the settlements due to the cyclic temperature variations threaten the serviceability of the supported structure. In this case, a careful design of the foundation is necessary to satisfy both ULS and SLS safety requirements.

A C K N O W L E D G M E N T S
Open Access Funding provided by Universita degli Studi della Campania Luigi Vanvitelli within the CRUI-CARE Agreement.

D ATA AVA I L A B I L I T Y S TAT E M E N T
The data that support the findings of this study are available from the corresponding author upon reasonable request. where the two constitutive tensors (tangent stiffness tensor) and (hardening tensor) are given by:

O R C I D
In equations (A2)-(A4):̇∶ = (̇)1 is the strain rate due to thermal expansion, with 1 the second-order identity tensor and the linear thermal expansion coefficient of the solid particles; the quantitieŝ≔ ∕ and ≔ ∕ , with a material constant, denote the direction of the intergranular strain tensor and its scaled norm, respectively, and is the fourth-order identity tensor, with components = (1 1 + 1 1 ). In equation (A4), is a material constant controlling the rate of evolution of the intergranular strain with accumulated strains. Note that the second equation (A1) implies that thermal expansion/contraction of solid particles alone does not induce any volume change in the solid skeleton, as demonstrated by Khalili et al. 80 .
Let ≔ − 1 be the translated stress tensor, accounting for the true cohesion of the material quantified by the material constant , and let: In equation (B2), and are material constants controlling the variation of the intercept and slope of the NCL with . The additional term in equation (B1) for the stress rate is the product of a scalar function , known as collapse potential factor, and a second-order tensor , depending on the temperature ratė. To define these two terms, the existence of a Thermal Stabilization Line (TSL) is postulated with equation: the scalar 0 is a reference temperature and and ′ are the void ratios at points O and O' on the NCL and on the TSL lines, respectively, at the reference pressure , see Figure B1. When the temperature is lower than the reference temperature, ′ is fixed equal to . The constants is the slope of the TSL line, while controls the accumulated volumetric contraction after the first thermal cycle.
The expression for the collapse potential factor is given by: where the symbol ⟨ ⟩≔( + | |)∕2 denotes the positive part of , while and * are the values of the void ratios determined on the NCL and on the TSL at the current effective pressure (see Figure B1) and is a material constant controlling the rate of irreversible volumetric contraction.
Finally, the second-order tensor function is defined as: