A thermo-hydro-mechanical model for energy piles under cyclic thermal loading

This paper introduces a thermo-hydro-mechanical finite element model for energy piles subjected to cyclic thermal loading. We address four particular features pertaining to the physics of energy piles: three-di-mensionality, embedded heat exchangers, soil constitutive modeling and pile–soil interface. The model is de- signed to capture the strong coupling between all important physical and thermomechanical processes occurring in a concrete pile embedding U-tubes heat exchangers and surrounded by a saturated soil mass. It encompasses solid and fluid compressibility, fluid and heat flow, thermoplastic deformation of soil, buoyancy, phase change, volume change, pore expansion, melting point depression, cryogenic suction and permeability reduction due to ice formation. The model is distinct from existing energy pile models in at least two features: (1) it can simulate the detailed convection-conduction heat flow in the heat exchanger and the associated unsymmetrical thermal interactions with concrete and soil mass; and (2) it can simulate cyclic freezing and thawing in the system and the associated changes in physical and mechanical properties of the soil mass that likely lead to thermoplasticity and deterioration of pile shaft resistance. The performance of the model is demonstrated through a numerical experiment addressing all its features.


Introduction
An energy pile (also known as thermal pile) is a dual-purpose structural element, functioning as a structural foundation and a ground source heat exchanger. This technology is appealing because it makes use of shallow depths to extract renewable energy suitable for heating and cooling of buildings. It reduces the use of fossil fuels and mitigates CO 2 emissions, and, compared to the conventional borehole heat exchanger, reduces installation costs and saves space. Nonetheless, construction companies are reluctant to apply this technology in daily practice because engineers are yet striving for more insight into the consequences of adding heat exchangers inside the pile on its structural functionality and the integrity of the pile-soil interface. Furthermore, so far, there is no compelling understanding of the long term behavior of energy piles when subjected to extreme operational scenarios leading to cyclic freezing and thawing. Current design methods do not account for the detailed thermal interaction between the heat exchangers and concrete, neither for the cyclic thermal loading and its effect on the thermo-hydro-mechanical behavior of soil and pile-soil interface.
Despite the increasing use of this technology, there are relatively few numerical models designed to study the thermo-hydro-mechanical (THM) behavior of energy piles and their surrounding soil mass. Several notable THM studies have been introduced, including those by Yavari et al. [35], Di Donna and Laloui [12], Gawecka et al. [16] and Anongphouth et al. [3]. They model energy piles with different levels of physical complexity, material constitutive relationships and pile-soil interaction.
A common feature of these studies is that they simulate the heat exchanger as a line heat source with a constant heat flux or a prescribed temperature. Normally, the heat exchanger is a U-tube, ~25 mm in diameter, made of high-grade polyethylene, fixed on the reinforcement cage inside the concrete pile. The heating (cooling) system works by circulating a fluid through the U-tube that collects (rejects) heat arising from a series of thermal interactions between the fluid, pipe wall, concrete and surrounding soil mass. These distinct geometrical and physical features are ignored in the line heat source approach, eliciting three main shortcomings: (1) it ignores the conductive-convective heat flow in the U-tube that varies following the daily and seasonal thermal load demands; (2) it ignores the three-dimensionality of the problem which results from the U-tubes configuration and their associated unsymmetrical heat flow and thermal stresses; and (3) it does not allow assessing the energy efficiency of the energy pile, which constitutes the main goal of using this technology.
The effect of temperature on the soil material parameters and its structural behavior constitutes the main disparity of energy piles from the conventional concrete piles, and hence given a special treatment in these studies. Yavari et al. [35], Gawecka et al. [16], and Anongphouth et al. [3] adopt constitutive relationships similar to those used in conventional soil mechanics; i.e. without specifically taking the temperature effects on the soil material parameters into account. However, the thermo-hydro-mechanical behavior of the soil mass is considered via the balance equations. 24,14,25], on the other hand, do consider the temperature effect on the soil material parameters. They introduced a phenomenological constitutive modeling approach based on continuum thermoplasticity. We briefly discuss this approach in Section 2.
The pile-soil interaction is another important aspect in energy pile design, currently treated with different levels of complexity. Gawecka et al. [16] and Yavari et al. [35] do not model the interface between the pile and soil. The pile finite elements share the same nodes with the soil elements, inhibiting thus any relative displacement between them. Anongphouth et al. [3], on the other hand, model the interface between the pile and soil, but adopt a constitutive model for the interface element similar to the surrounding soil mass, i.e. the soil in the interface has the same strength as the soil mass, but can exhibit relative displacement. Suryatriyastuti et al. [32] employ frictional interface elements based on the Mohr-Coulomb failure criterion. Sutman et al. [33] use thin-layer elements with elastic properties to model the pile-soil interface. Di Donna and Laloui [12] use a thin-layer element based on their ACMEG-T model (Advanced Constitutive Model for Environmental Geomechanics), but the same constitutive parameters as those used for the surrounding soil. Rotta Loria et al. [29], Rotta Loria et al. [30], and Di Donna et al. [13] have employed more advanced thin-layer elements based on the thermo-elastoplastic Mohr-Coulomb failure criterion with reduced mechanical properties compared to the surrounding soil.
In this paper, we address these three aspects of energy piles: heat exchangers, soil constitutive behavior and pile-soil interaction, together with the three-dimensionality of the problem. We introduce a 3D finite element model describing the THM behavior of an energy pile system subjected to cyclic thermal loads. Details of the modeling approach are given hereafter.

Modeling approach
In a previous work, the authors introduced a detailed thermo-hydromechanical computational model for freezing and thawing in porous media [5]. It is formulated based on the averaging theory and discretized using axial-symmetric finite elements. The model is capable of capturing the strong coupling between all important phenomena and processes occurring in the porous domain, including solid and fluid compressibility, fluid and heat flow, buoyancy, phase change, volume change, pore expansion, melting point depression, cryogenic suction and permeability reduction due to ice formation. Here, this model is customized to be used for energy piles, and tailored to consider four particular features: three-dimensionality, embedded heat exchangers, soil constitutive modeling, and pile-soil interface.
The energy pile is by definition three-dimensional (3D), pertaining to the configuration of the U-tube heat exchanger inside it (Fig. 1a). The temperature in pipe-in at cross sections along the pile is different than that in pipe-out (Fig. 1b). As the pile diameter is relatively large compared to the U-tube, the temperature gradient between pipe-in and pipe-out gives rise to radially unsymmetrical thermal stresses and strains in the pile and the surrounding soil mass. Furthermore, piles are usually installed in groups with arbitrary configurations, and thus, detailed modeling of pile groups must be three-dimensional. In this work, 3D brick and wedge finite elements are utilized for this purpose.
Heat flow in the heat exchanger is conductive-convective arising from the thermal interaction between the circulating fluid in the Utubes and the surrounding concrete and soil mass. The U-tube heat exchanger is relatively small in diameter ( 25 mm) which entails that the heat flow inside it is nearly one-dimensional (1D). This allows modeling heat flow in the U-tube by a line element with its thermal interaction with the pile expressed explicitly in the governing heat equation [1]. Here, the embedded finite element method is utilized to discretize the heat exchanger. Elements of this kind, can go through the host 3D elements, avoiding the need for fine meshes and keeping the mesh structured. The soil is by definition a multiphase material constituting, basically, a solid matrix phase and a water phase. Within the temperature range of this application, the solid matrix, except for thermal expansion and contraction and its Young's modulus, is not significantly affected by the temperature variation. Rather, it is the water phase that is being affected: the temperature alters its mass density, heat capacity, thermal conductivity and viscosity. Additionally, the water might change phase, leading to changing soil strength, permeability and volume. Based on this, we adopt a multiphase mixture approach to formulate the constitutive relationships for the soil material. This approach is different than the phenomenological approach of 24,14,25], who consider the soil as a homogeneous porous matrix. The behavior of the soil in their approach is described by a thermoelastoplastic model with a yield function expressed in the stress-temperature space, as shown schematically in Fig. 2a. This kind of models is basically appropriate for the soils under investigation, and any other type of soils must be examined experimentally to establish its parameters dependency on temperature. In our multiphase mixture approach, however, the soil is described by its constituents: a solid matrix and a water phase (Fig. 2b). Following this, the conventional constitutive relationships for soil can be utilized to model the solid matrix, and the equations of state (EOS) can be utilized to describe the water phase. We model the solid matrix using the modified Cam-Clay yield function, and the water using equations of state describing its pressuretemperature-volume relationships and the dependency of its physical parameters on temperature. The thermal dependency of soil is thus included fundamentally via the water EOS, together with the thermal expansion coefficient and temperature dependent Young's modulus of the solid matrix.
Understanding the behavior of the pile-soil interface is essential to prevent the deterioration of the pile shaft resistance and bearing capacity under the effect of cyclic thermal loads. We adopt the Desai thinlayer approach [11], but extending it to a multiphase mixture exhibiting sliding, debonding and volume change pertained to the induced thermo-hydro-mechanical forces. The interface material behavior is expressed using the Drucker-Prager yield function and the water EOS.

Governing equations
The physical domain consists of a concrete energy pile embedded in a soil mass and subjected to mechanical and cyclic thermal loads. The pile embeds U-tube heat exchangers exhibiting conduction-convection heat transfer due to the thermal interaction between the circulating fluid in the U-tube and the surrounding concrete. The concrete pile exhibits thermo-mechanical behavior due to the applied mechanical loads, the thermal interaction with the heat exchangers and the thermomechanical interaction with the surrounding soil mass. The soil mass exhibits thermo-hydro-mechanical behavior arising from the thermomechanical interaction with the concrete pile and the fluid flow due to consolidation, buoyancy and cryogenic suction in the soil body. The field equations governing these processes in the soil mass, concrete, pile-soil interface and heat exchanger are given hereafter.

Soil mass
The soil mass is considered saturated, isotropic and non-isothermal with local thermal equilibrium. It constitutes a multiphase mixture composing a solid matrix and water, which can change phase between liquid water and ice. The soil exhibits solid and fluid compressibility, heat and fluid flow, thermo-elastoplasticity, buoyancy, volume change, pore expansion, permeability change, cryogenic suction and melting point depression. The three phases (solid matrix, liquid water and ice) interact physically with each other and exchange mass, momentum and energy. The governing balance equations describing the conserved quantities in a multiphase porous medium domain are given in detail in . Here, we state the balance equations together with the constitutive relationships relevant to the energy pile system. They are expressed in terms of the primary state variables: solid matrix displacement u, water mixture pressure p m , water mixture specific enthalpy h m , solid matrix temperature T s , and cryogenic suction s c .

Conservation of momentum
The averaged macroscopic linear momentum balance equation of a multiphase mixture constituting a solid matrix, a liquid water phase and an ice phase, and subjected to thermo-hydro-mechanical forces can be expressed in an incremental form as where = + m p ' s is the effective stress, with being the total stress; is Biot's coefficient; = m [1, 1, 1, 0, 0, 0] T ; g is the gravitational vector;p s is the pressure exerted by the water phase on the solid matrix, defined as where p lw and p ice are the liquid water and ice pressures, and S lw and S ice are their degrees of saturation; and eff is the effective mass density, defined as in which is the porosity, s , lw and ice are the mass density of solid matrix, liquid water and ice, respectively.
Eq. (5) is the equilibrium equation of the volume-averaged, compressible multiphase mixture of solid matrix, liquid water and ice.

Constitutive relationships
The effective stress increment in Eq. (5) is described as in which s is the volumetric thermal expansion coefficient of the solid matrix, L is the differential operator, and D ep is the tangential elastoplastic stiffness matrix of the solid matrix, described, as where F is a yield function, Q is a potential plastic function, p is the plastic strain tensor, and D e is the tangential elastic stiffness matrix, expressed, for a three-dimensional solid matrix, as in which is Poisson's ratio, and E T ( ) is a temperature-dependent elastic modulus, defined here as where E 0 is the Young's modulus at a reference temperature, T 0 , and b is a material parameter.
For the relatively small thermo-hydro-mechanical strain level generated by the cyclic thermal loading, it is reasonable to assume that the soil behavior is governed by an isotropic hardening/softening rule with an associated flow rule, giving = Q F in Eq. (7). The modified Cam-Clay yield function is utilized for this purpose. In the p q stress space, it is expressed as in which p is the mean normal stress, q is the deviatoric stress, M cs is the slope of the critical limit state line, p c is the instantaneous consolidation pressure, defined as [26] where p c0 is the pre-consolidation pressure, n is a material constant related to the hardening and softening of the material, v p 0 is the initial volumetric plastic strain, and = m v p T p is the volumetric plastic strain. The yield surface is smooth with an initial size governed by the magnitude of the pre-consolidation pressure (p c0 ). The yield surface expands in size due to strain hardening arising from compaction, and shrinks in size due to strain softening resulting from volume increase. More details can be found in Chen and Baladi [9] and Lewis and Schrefler [26].
The consistency condition of the yield surface, Eq. (10), is where which is the plastic strain increment defining an associated flow rule, with The other constitutive relationships in Eq. (5)

Momentum field equation
Substituting the constitutive equations into the momentum balance equation, Eq. (5), gives the momentum field equation of the solid matrix.

Conservation of mass
The averaged macroscopic mass balance equation for the multiphase soil mixture is formulated by summing the mass conservation equations of the individual phases.

Solid matrix
The mass balance equation for a homogeneous solid matrix phase can be described as Following Lewis and Schrefler [26], the solid matrix mass density can be described as This equation reveals that the porosity is not constant and can be altered by the thermo-hydro-mechanical forces in the soil.

Liquid water phase
The mass balance equation for the liquid water phase can be expressed as in which v lw is the extrinsic averaged velocity of liquid water, and m lw ice is the mass exchange rate arising from the phase change between the liquid water and ice. Inserting Eq. (18) into Eq. (19) and expressing the liquid water density and saturation in terms of the state variables: lw m m lw m m , yields Ice phase Similar to the liquid water phase in Eq. (20), the mass balance equation of the ice phase can readily be derived to give in which v ice is the extrinsic averaged velocity of ice, considered to be zero in this application (

Water mixture (liquid water and ice)
Considering the following identities: and by summing Eqs. (20) and (21), assuming that the liquid water flow rate is governed by Darcy's law, and defining the water mixture pressure as = + p p s m lw c , gives [5]: where µ lw is the dynamic viscosity of liquid water, k is the absolute permeability, and k rlw is the relative permeability. Eq. (23) is the continuity equation of the volume-averaged, compressible multiphase mixture of solid matrix, liquid water and ice with phase change.

Water equation of state (EOS)
The thermodynamic state variables and properties of the liquid water, ice and the water mixture, m , T m , ice , lw , S ice , S lw , µ lw are obtained from the equation of state of water, adopted from IAPWS [21] and other relevant literature, given in Appendix A. Kurylyk and Watanabe [23] presented an interesting review describing different forms of the Clapeyron equations and empirical relationships for soil freezing curves (SFC). Here, we adopt an exponential function of the form:

Melting point depression
in which S is the residual unfrozen water content at a relatively cold condition, T f is the bulk freezing temperature, and a is a material constant.

Cryogenic suction
The cryogenic suction, s c , exhibits a substantial change for each degree Celsius below zero. As a consequence, the cryogenic suction is considered here a primary state variable, to have it directly computed from solving the finite element equations, rather than being calculated in the post processing (see  for more details). The computed quantity has to satisfy the Clausius-Clapeyron relation [26]: where L f is the latent heat of fusion of water. To satisfy this condition, the following constraint is imposed:

Relative permeability
Even though the domain is fully saturated, the water exhibits a phase change during freezing and thawing, giving rise to a quasi-partially saturated condition within the water phase. As for the partially saturated conditions, the relative permeability of liquid water is described based on the Brooks and Corey relationship [8]: where is a material constant.

Mass field equation
Substituting the constitutive relationships from Eqs. (24)- (27) and Appendices A and C into the macroscopic mass balance equation, Eq. (23), gives the mass field equation.

Conservation of energy
The averaged macroscopic energy balance equation for a multiphase mixture in local equilibrium is formulated by summing the heat equations of the solid matrix and the water mixture.

Solid matrix
in which T sm is the solid matrix temperature and s is its thermal conductivity.

Water mixture
Summing Eq. (28) to Eq. (29) and considering a thermal local equilibrium, is the effective thermal conductivity of the porous domain, with s denoting the thermal conductivity of the solid matrix, and ice ice , with lw and ice denoting the thermal conductivity of liquid water and ice, respectively. (17) and (18) where K T is the bulk modulus of the solid skeleton.

Energy field equation:
Substituting the involved constitutive equations of water from Appendix A into Eq. (32) gives the energy field equation.

Concrete
In this study, the concrete of the pile is considered non-porous, linear elastic material. Though, it is worth mentioning that concrete is a porous material and might exhibit damage due to freezing [22,17], but the focus in this paper is on the effect of freezing on the soil mass and the pile-soil interaction.
The relevant conservation laws are the momentum and energy equations.

Momentum field equation
The momentum field equation for concrete can be expressed as where T c is the concrete temperature, c is the volumetric thermal expansion coefficient, c is the mass density, and D c is the linear elastic stiffness matrix, similar to that for soil, Eq. (8), but its elastic modulus E c is constant (not function of temperature).

Energy field equation
The energy field equation for concrete can be expressed as where c c denotes the specific heat capacity of concrete, and c is its thermal conductivity.

Heat exchanger
The heat exchanger is considered a 1D, non-deforming conductiveconvective heat source. It is governed only by the energy equation.

Energy field equation
The energy field equation of the heat exchanger is expressed as where T r , r , c r , r , and v r denote the temperature, mass density, specific heat capacity, thermal conductivity, and velocity of the circulating fluid, respectively, and b cr is the thermal interaction coefficient between the circulating fluid and concrete, described in Appendix B. It is worth noting that the inclusion of the thermal interaction term in the differential equation allows embedding the heat exchanger element in the 3D concrete element.

Pile-soil interface
Accurate estimate of the load-deformation behavior of the pile-soil interface is essential for determining the integrity of the interaction between the pile and the soil. Most current numerical models for pile-soil interaction use interface elements based on the joint element, first introduced by Goodman et al. [36]. This element is of zero thickness, and formulated based on the relative displacements of solid elements surrounding the interface, Fig. 3. It is expressed in terms of stressdisplacement relationship, as in which d is the incremental stress, C is the stiffness matrix and u d is the increment of the relative displacements between the two sides of the interface.
This kind of elements can describe slip/no-slip conditions at the boundary between solid elements. They are basically suitable for piles with mainly axial loading and exhibiting nearly no volume deformation at the interface with the soil mass. However, for energy piles, the induced thermo-mechanical forces give rise to thermal expansion and contraction in and around the pile, leading to changes in volume. Furthermore, cyclic thermal loading can lead to slip, debonding and rebonding at the pile-soil interface. The volume change becomes more significant with thermal loads leading to cycles of freezing and thawing. Hence, for energy piles, it is more appropriate to utilize a solid element which exhibits volume change. This kind of elements was first proposed by Zienkiewicz et al. [37] who utilized an isoparametric finite element to model the interface of jointed rock systems. Based on this, Desai et al.
[11] introduced a thin-layer element for modeling interfaces with different deformation modes, including slip, no-slip, debonding and rebonding. This element describes the soil at the interface as a homogeneous solid matrix. Here, we extend this element to a multiphase porous mixture.
The thin-layer element is basically formulated similar to the solid elements, and its constitutive relationship is expressed in the usual stress-strain space, as Similar to the soil momentum balance equations, Eqs. (5) and (6), the momentum balance equation for the thin-layer interface can be expressed as in which C ep is the elastoplastic stiffness matrix of the interface material, i is the volumetric thermal expansion coefficient, T i is the interface temperature, and eff i is the effective mass density, defined as where i is the porosity of the interface, and si is the mass density of the interface solid matrix. The porosity change of the soil at the interface is governed by Eq. (18).
The stiffness matrix, C ep , is described in Eq. (7) with an associated yield and potential functions, F i , governed by the Drucker-Prager elastoplastic criterion: where is the friction angle and c is the cohesion, described as [26] = c c e n 0 1 ( ) in which c 0 is the initial cohesion and n i is an empirical constant for

Initial and boundary conditions
At = t 0, the primary state variables are described by their initial values, as where q is the Neumann boundary; t is a prescribed traction; q lw is a prescribed mass flow rate of liquid water; and Q soil air is the convective heat transfer between the ground surface and air, with T air the air temperature, and b sa the associated thermal interaction coefficient, given in Appendix B. Pile: Pile-heat exchanger: where b cr is the thermal interaction coefficient between the heat exchanger and the pile concrete, given in Appendix B. This equation indicates that this heat flux is equal in magnitude and opposite in direction to Q rc (right-hand side of Eq. (35)). Pile-soil: in which Q cs is the heat transfer between the soil and concrete pile, with b sc the associated thermal interaction coefficient, given in Appendix B.

Finite element discretization
The governing equations are linearized using the modified Newton-Raphson method and discretized using the Galerkin finite element method. The linearization and discretization procedures are given in detail in . Here, we present the finite element equation, tailored to the energy pile system, and we discuss in detail the discretization of the embedded 1D element for the heat exchanger.
The global finite element equation for the energy pile can, symbolically, be described as where the subscripts/superscripts s c i , , denote the soil matrix, concrete and interface, respectively, and the subscripts m r , represent the pore water mixture and circulating water, respectively. The vector on the right-hand side represents the THM force vector, including heat fluxes M.M. Arzanfudi, et al. Computers and Geotechnics 125 (2020) 103560 due to the thermal interactions at the interface boundaries between the heat exchanger, pile and soil mass. The matrix entries of Eq. (49) represent 3D and 1D element matrices, expressed hereafter.

3D element
The matrix and vector entries for 3D soil, concrete and interface elements, including C K ,

Embedded 1D element
The heat exchanger is discretized using a linear embedded 1D element. This element can go through the host 3D elements without affecting the underlying structure of the mesh, Fig. 4. The essence of this kind of elements is that the host and embedded elements are discretized in the usual manner, except for their interaction terms; f cr and f rc (Eq. (49)). This enables a simpler discretization of the problem.
The are the finite element weighting functions for the host element and embedded element, respectively, and Q cr is the coupling term between the two elements, Eq. (47), expressed as where N eh is embedded-host coupled shape function. For a two-node 1D element embedded in an 8-node 3D element, Fig. 4, using the Gauss quadrature, Eq. (51) can be expressed as where n is the number of Gaussian integration points in the 1D element ( = n 3 in Fig. 4), w k is the Gaussian quadrature weight, and 1, 2 ej are the weighting and shape functions of the host and embedded elements at sampling point , respectively.

Model verification against London energy pile experiment
There are relatively few well-documented full-scale energy pile experiments in the literature. Here, the results of the well-known London energy pile experiment presented by , and Amatya et al. [2] are utilized to verify the proposed model.

Physical domain
The London energy pile experiment is conducted using a single energy pile, 23 m in length, in which a double U-tube heat exchanger is embedded. The lower 18 m of the pile is 0.55 m in diameter, and the   Arzanfudi, et al. Computers and Geotechnics 125 (2020) 103560 upper 5 m is cased in two diameters: 0.61 m for 4.1 m and 0.9 m for 0.9 m, Fig. 5. The configuration of the U-tubes inside the pile is shown in Fig. 1a. The energy pile is embedded in a soil mass consisting of three layers: 1.5 m of made ground, 2.5 m of sandy gravel, and the rest is London clay. The physical and thermal parameters of the system are given in Table 1.

Initial and boundary conditions
Initially, the temperature in the domain, as reported in , is 19.5 °C.
Mechanical boundary condition: The pile is subjected to two mechanical vertical loading/unloading forces: 1200 kN and 1800 kN.
Thermal boundary conditions: The system is subjected to 31 days of cooling, followed by 12 days of heating. Fig. 6 shows the input temperature which ranges between −6°C and 40 °C, applied at the inlet of the U-tube. According to the literature of this experiment, a power failure occurred after 36 days and lasted for 4 days.

Computational domain
The physical processes in the system are three-dimensional, but the geometry is planar symmetric with respect to the plane that passes through the middle of the pile, Fig. 7. Accordingly, only half of the geometry is considered. The computational domain is 100 m in diameter and 70 m in depth, divided into three soil layers resembling made ground, sandy gravel and London clay (Table 1). An energy concrete pile, embedded in which double U-tubes, is located at the axis of symmetry.
The mesh consists of 1548, 3D linear hexahedron and wedge-shaped elements for the pile and the soil mass; 84, 3D hexahedron elements for the pile-soil interface; and 26, embedded 1D linear finite elements for the heat exchanger. The thickness of interface is 0.02 m. A view of the finite element mesh at the proximity of the pile is given in Fig. 7.
The prescribed initial and boundary conditions are: where z is the depth.

Boundary conditions:
Thermal boundary conditions: Porous:

Hydraulic boundary conditions:
Top/side boundaries: drained Bottom boundary: undrained Mechanical boundary conditions: = F pilehead varying, see Fig. 6b Top boundary: no constraints (free to move) Bottom boundary: fully constrained Side boundaries: horizontally constrained

Experiment-computation comparison
The computed and measured temperatures are compared at a point in the pile 9 m below the surface and 0.2 m from the axis of symmetry; and two points in the soil mass 9 m in depth, and 0.5 m and 2 m radially. The angular coordinates of the sensors are not indicated in the literatures, and hence, the computational results are presented at five angular positions , =°°°°°0 , 45 , 90 , 135 , and180 , Fig. 8a. Fig. 8b shows the measured and computed temperatures in the pile, and Fig. 8c and Fig. 8d show those in the soil mass. In general, the figures display a good match between the two results. Though, due to power failure at time, = t 36 40days, there is some mismatch between the computed and measured data in the pile. The power failure has been modeled as switching off of the heat pump (see Fig. 6), but the exact consequences of this failure on the measurement devices are not clear. Therefore, the comparison is not reliable in this period. Fig. 9 shows the computed thermo-mechanical pile head displacement as compared to the measured data. The figure shows a good match between the two results, especially before the power failure. Fig. 10 shows the computed mechanical and thermo-mechanical axial strain along the pile, together with the measured data. Fig. 10a shows that there is a good agreement in the axial strain due to mechanical loading. Fig. 10b displays the thermal strain at the end of cooling ( = t 35 days). This figure shows a reasonable matching between the computed and measured results, though, there is some deviation between the two results, which might be attributed to the mobilization of the pile-soil interface due to the thermal contraction of the concrete and soil. The behavior of the pile-soil interface of the London energy pile is not reported in the literature. Fig. 10c shows the strain at the end of heating ( = t 47 days). Even though the computed results are not matching the measured data, they exhibit a similar trend. The mismatch can be attributed to the power failure and the uncertainty at the pile-soil interface. Other literatures reveal a similar mismatch, see [35,16,3].
It is worth mentioning that the observed reversibility in the axial strain in Fig. 10 occurs due to a combination between the reversibility of the elastic strain and the change of stress direction between the heating mode and the cooling mode. However, its magnitude is influenced by the configuration of the heat pipes inside the pile and their temperature gradient along the pile; the pile-soil interaction; and the type of the surrounding soil mass and its mechanical, thermal and hydraulic properties. Therefore, detailed description of the physics is of paramount importance to describe the behaviour of energy pile systems. Mimouni and Laloui [27] and Rotta Loria and Laloui [31], among others, highlighted the issue of reversibility of the axial and radial strains in energy piles.

Numerical energy pile experiment
The verification example given in the previous section indicates that the proposed multiphase thermo-hydro-mechanical modeling approach can provide reasonably accurate computational results describing the behavior of typical energy pile set-ups. However, this is not exceptional as many other computational approaches, including Yavari et al. [35], Di Donna and Laloui [12], Gawecka et al. [16], and Anongphouth et al. [3], are capable of simulating full-scale experiments to a practically accepted level of accuracy. This is feasible since the boundary conditions adopted in most full-scale experiments are restrained to a limited range. The heating system in these experiments is usually operated with relatively low heat extraction rates to make sure that the soil does not freeze. This makes the simulation of the thermo-hydro-mechanical processes feasible via standard constitutive models and numerical discretization procedures. However, in practice, such a constraint can lead to limiting the energy efficiency of the system in at least two important engineering scenarios: 1. Limiting the operation of the system to relatively low circulating flow rates would lead to low heat harvest, and 2. Limiting the system operation to above freezing level would lead to shorter operation time, and hence less efficiency.
The proposed model, on the other hand, is designed to simulate conditions that go beyond the current operation constraints to allow for higher energy extractions and longer periods of operation. To examine these capabilities, we extended the London energy pile experiment numerically to attain cycles of freezing and thawing conditions in the system.

Initial and boundary conditions
The initial and boundary conditions are the same as those for the verification example in Section 5, except the following:

Initial condition:
The temperature in the system is assumed 10 °C. This relatively low temperature was chosen to accelerate the freezing condition in the system.

Mechanical boundary condition:
A mechanical load of 1200 kN is applied at the pile head prior to the thermal loading.

Thermal boundary conditions:
The upper boundary of the soil mass is subjected to a constant air temperature of 10 °C, and the bottom boundary is prescribed at 10 °C. The heat exchanger is subjected to a cyclic thermal load, ranging between 0 and 7 kW. Eight cycles of switching on and off for 2 h are applied at the inlet of the heat exchanger for 16 h, followed by 8 h switching off, Fig. 11a. The cycles are repeated for 30 days, after which the heat pump is switched off for 5 months, as shown in Fig. 11b. According to current specifications this thermal load is not realistic, but adopted to study the feasibility of using the proposed model for high energy extraction rates and cyclic freezing and thawing conditions.
The cyclic heat pump power is converted to a prescribed inlet temperature, T in , utilizing the heat pump-heat exchanger set-up schematically given in Fig. 12. Depending on the HVAC design, the heat pump extracts a required amount of heat from the fluid coming out of the heat exchanger outlet, and returns it to its inlet. The amount of power extracted from the heat pump can be calculated by in which m is the mass flow rate of the circulating fluid, c r is the specific heat capacity of the circulating fluid, T p in is the fluid temperature entering the heat pump and T p out is the fluid temperature leaving the heat pump. Ignoring the heat loss in the pipes, which carry the circulating fluid, T p in is equal to T out (temperature leaving the heat exchanger); and T p out is equal to T in (temperature entering the heat exchanger). Having this heat exchange set-up between the heat pump and the heat exchanger, T in can be calculated from Eq. (55).

Physical and computational domains
The geometry and material parameters of the physical domain, together with the types of the finite elements and the finite element mesh size are identical to those of the verification example in Section 5.

Results and discussion
The behavior of the system is presented in terms of deformation at the head and toe of the pile, temperature at the inlet and outlet of the heat exchanger, and temperature in the soil mass, 9 m below the surface, 0.3 m in the radial direction, with angular coordinates: =°°°°°0 , 45 , 90 , 135 , and 180 . To highlight the cyclic loading, the diagrams are presented in a semi-log format. Fig. 13a shows the temperature variation with time in the soil. It shows that freezing starts after 14 days of operation and the temperature reaches −33 °C after 30 days. A closer look at the figure reveals that in the period between day 28 and day 30 the temperature exhibits a sudden drop. This occurs due to the accumulation of ice in the soil. Fig. 13b shows the ice volume fraction during this period. It reveals that in day 14, with the onset of freezing, the ice volume fraction in the pores was 1%. According to Eqs. (A.5) and (A.8) in Appendix A, the thermal conductivity for this water mixture is 0.56 W/m⋅K. However, in day 30, the ice content became 9% and the thermal conductivity of the pore water rose to 0.70 W/m⋅K. This increase in the thermal conductivity allows more heat to transfer between the soil and the pile leading to this drop in temperature. Fig. 13c shows the temperature profiles at the inlet and outlet of the heat exchanger. It reveals that with this relatively large heat pump power, the temperature in the circulating fluid dropped to below 0 °C from the first day, and in 30 days the temperature dropped to nearly −20 °C. The figure also reveals that the sudden drop in temperature in the soil in the period between day 28 and 30 has been reflected on the circulating fluid. Fig. 13d shows the pile head and toe displacements. The displacement at the head generated by the mechanical load is 2.2 mm and at the toe is 0.9 mm. Upon the start of the cyclic thermal loading, the displacement in the pile exhibits ups and downs due to the thermal expansion and contraction of concrete and its interaction with the  surrounding soil mass. This can be clearly seen in the sketch given in Fig. 14. This figure shows the pile at its initial condition, upon the application of the mechanical load and during the thermal load cycles. The ups and downs of the pile are scaled to the amount of displacement at the head and toe of the pile, and the arrows indicate the movement direction. The following can be observed: = t 0day: upon mechanical loading, both head and toe moved in the direction of the applied force. = t 18days: four days after the onset of soil freezing on day 14, the pile exhibits contraction. = t 28days: with the increase of ice content in the soil (see Fig. 13b), the pile exhibits lifting up due to the soil heave at the toe. = t 30days: the pile temperature drops due to its thermal interaction with the heat exchanger and soil, causing the pile to contract again. = t 40days: during thawing, the pile exhibits expansion. = t 55days: with the increase of plastic deformation in the soil (see Fig. 15f), the soil exhibits settlement around the pile, causing the pile to lift up. = t 180days: at the end of thawing, the pile exhibits expansion again. Fig. 15 presents three-dimensional plots for the temperature distribution and fluid flow, the effective plastic strain, and the deformation in the proximity of the pile on day 30 (maximum freezing) and day 180 (end of thawing). Fig. 15a shows the 3D temperature distribution in the pile and the soil around it, together with the fluid flow velocity vectors. The figure shows that the concrete and approximately 0.1 m of the surrounding soil have been frozen on day 30. It also shows the water flow towards the frozen zone due to cryogenic suction. (Note that the flow pattern is not symmetric because the temperature distribution in Fig. 13. Temperature, ice fraction and pile displacement: (a) temperature in soil, (b) ice fraction in soil, (c) temperature at inlet and outlet of heat exchanger, and (d) pile head and toe displacements. and around the pile is not symmetric due to the unsymmetric configuration of the heat exchanger inside it (Fig. 1a). Fig. 15b reveals that upon switching off the heat pump for 5 months, the temperature in the system has recovered as a result of the heat flow between the soil mass and its boundaries. Fig. 15c to f show the deformed mesh and their associated effective plastic strains. Fig. 15e shows that the plastic strain has already began in the first 30 days and increased considerably after thawing, Fig. 15f. This can be noticed from the deformed mesh on day 180 (Fig. 15d) as compared to day 30 (Fig. 15c). The relatively large deformation in the soil can be attributed to the fact that, upon freezing, the soil stiffness increases due to ice formation, accompanied by expansion of pores. Upon thawing, the soil stiffness reduces significantly in regions where excessive pore expansion has occurred, causing the soil above it to collapse under gravity.

Conclusions
Current energy pile design procedures put a limit to the amount of extracted energy to make sure that operating the energy pile for heating would not lead to soil freezing. This constraint is good for safety of buildings, but short on the economy. It lessens the energy efficiency of the system because it requires low flow rates and shorter operation time. One of the reasons for this constraint is attributed to lack of understanding of the coupled thermo-hydro-mechanical forces arising from freezing and thawing of porous media and their effects on the soil properties and the pile-soil interface. Freezing and thawing in a soil mass gives rise to solid and fluid compressibility, fluid and heat flow, thermo-elastoplasticity, buoyancy, phase change, volume change, pore expansion, melting point depression, cryogenic suction and permeability reduction due to ice formation. None of the existing numerical tools used for energy piles design and analysis is capable of capturing the coupling of these processes. This paper addresses this issue and demonstrates that the proposed model is capable of capturing these processes and their coupling effects on the integrity of the energy pile systems.
The paper addresses four particular features pertaining to the physics of energy piles: three-dimensionality, embedded heat exchangers, soil constitutive modeling and pile-soil interface. It highlights several aspects related to the physics of energy piles and their operation that need to be considered in energy pile design.

A concrete pile exhibits deformation and movement under cyclic
thermal loading with magnitudes related to its thermo-mechanical properties and the thermo-hydro-mechanical properties of the pile-soil interface and the soil mass. Fig. 14 clearly illustrates how the pile can expand and contract locally or move as a rigid body during cyclic thermal loads, leading ultimately to deteriorating the integrity of the pile-soil interaction. This signifies the fact that the thermo-hydro-mechanical properties of the system can alter during daily and seasonal operation of the energy pile, requiring special attention during design. 2. Freezing in soil leads to pore expansion which, upon thawing, might lead to weakening the strength of the pile-soil interface. Fig. 15d shows how the thawing has led to weakening part of the interface that caused the collapse of the soil above it. This emphasizes the need of taking the likely occurrence of freezing and thawing into consideration during the design. Even if current design procedures prohibit reaching to the freezing point, freezing might occur due to malfunctioning in the mechanical parts, misuse of the system or unexpected change in weather. 3. Heat conduction-convection in the heat exchanger gives rise to axial and unsymmetrical radial thermal strain gradients in the pile and the surrounding soil mass. The numerical experiment has shown that such gradients can lead to localized damage that can affect the pile-soil shaft integrity. This effect cannot be captured by the line heat source model, requiring thus revaluating the validity of this model in designing energy piles.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. where r , c r , µ r , r and v r , denote the density, specific heat capacity, dynamic viscosity, thermal conductivity and velocity of the circulating fluid, respectively.

Pile-Soil
The thermal interaction coefficient the pile-soil interface is: in which r r 2 1 represents the thickness of interface.

Soil-Air
The soil-air thermal interaction can be described as