Changes in shaft resistance and pore water pressures during heating of an energy foundation

. This study focuses on the evolution of shaft resistance during operation of a geothermal energy foundation installed in a saturated glacial till layer. Energy foundations are a sustainable alternative to traditional space heating and cooling approaches for buildings. Despite efficient operational performance, there are still valid concerns regarding the effects of heating on the structural performance of foundations. To investigate the effect of heating at the soil-pile interface, four drilled shafts are utilized as energy foundations on the Urbana-Champaign campus of the University of Illinois and instrumented. Although the energy foundations are not yet operational, a theoretical investigation is possible to understand the effects of heating on the evolution of thermally induced pore water pressures and the shaft resistance of an energy foundation. A thermo-poroelastic numerical model is validated against an analytical solution, then is used to analyze the thermo-mechanical response of the soil-structure system under different conditions. The results indicate that the evolution of pore water pressure is affected by the rate of heating and the hydraulic conductivity of the surrounding soil, as expected. Higher pore water pressures are generated in the case of low hydraulic conductivity and higher rates of heating. Prior to the dissipation of excess pore pressures, the changes in shaft resistance are variable and influenced by the thermally-induced deformation of the foundation and the surrounding soil.


Introduction
Energy foundations are used to exchange heat between structural elements of foundations and surrounding soil for heating and cooling of buildings, and they are proven to provide efficient and sustainable heating and cooling [1]. They function by circulation of a fluid through pipes installed with different configurations in the foundations. One advantage of geothermal energy foundations compared to other shallow geothermal energy technologies, such as thermal energy storage systems and ground source heat pumps where horizontal or vertical boreholes are used, is that they don't require extra drilling or excavation [2,3].
One factor impeding increased numbers of installation of energy foundations is the impact of heating and cooling cycles on the structural performance of the foundation. Potential issues include thermally induced axial and radial stresses and strains within the foundation and adjacent soil, as well as the generation of thermally induced excess pore water pressures in soils having low hydraulic conductivity. These stresses and pore water pressures could affect the mobilized or ultimate foundation shaft resistance and lead to instability in the overlying structure. For an accurate design of energy foundations, a better understanding of the behavior of the foundation-soil system under combined mechanical and thermal loads is required.
To investigate the impact of heating on the shaft resistance, four energy foundations are installed and instrumented to provide heating and cooling to the new Ven Te Chow Hydrosystems Laboratory building on the Urbana-Champaign campus of the University of Illinois. Energy foundations are not yet operational but a theoretical investigation is possible to study thermohydro-mechanical behavior of energy foundations and surrounding soil during heating. In this study, a coupled thermo-hydro-mechanical model was developed based on poroelastic theory and solved with finite element method. The model is first validated against results from two laboratory experiments in the literature that allow to perform parametric study to investigate the factors affecting excessive pore water pressure development and consequent changes in shaft resistances.

Background
When a soil layer is heated, volumetric strains are generated due to the differential expansion of fluid and solid phases [4][5][6]. The magnitude of the strains depends on the rate of heating and drainage conditions based on type of loading. If the soil has a sufficiently low hydraulic conductivity, excess pore water pressures are immediately generated upon heating [6]. To investigate the thermally induced strains in undrained conditions, Biot's [7] generalized theory of poroelasticity was extended by Booker and Savvidou [8] to include thermal effects in an elastic, saturated soil layer. Britto et al. [9] and Lewis et al. [10] applied similar formulations to simulate pore water pressure generation and consolidation around buried heat sources with the finite element method. In recent years, many investigators have favored more advanced and complex semi-empirical constitutive models based on Modified Cam-Clay theory (e.g., [5,11]). These models can estimate plastic deformations but require laboratory characterization of several nonisothermal material parameters, including parameters that describe changes in preconsolidation stress, elastic moduli, and thermal expansion as functions of temperature. Thermoplastic flow rules include various soil-specific coefficients to accurately describe both isotropic and deviatoric yielding. On the other hand, poroelastic theory provides a simple framework for simulating the thermal response of overconsolidated soils with few input parameters.
Excess pore water pressures reduce the shaft resistance of a pile foundation under undrained loading conditions. The shaft resistance can be expressed as follows [12]: where fs is the shaft resistance (kPa), ca is the adhesion between the pile and the soil (kPa), σ ' h is the average effective horizontal stress on the pile (kPa), and δ is the angle of friction between the pile and the soil. The effective horizontal stress is calculated as follows: where K is an earth pressure coefficient, σ ' v is the effective vertical stress (kPa), σv is the total vertical stress (kPa), and u is the pore water pressure (kPa). Shaft resistance decreases as excess pore water pressures are generated. Thermally induced deformations in a geothermal energy foundation may affect shaft resistance during heating. Laloui et al. [13] and Bourne-Webb et al. [14] instrumented field-scale energy foundations and observed elastic axial strains and changes in shaft resistance induced by heating loads. Bourne-Webb et al. [15] introduced a qualitative framework to describe thermallyinduced deformations in energy foundations and their effect on axial stress and shaft resistance. Thermal loads cause the pile to expand or contract around a null point and mobilize axial stresses and shaft resistance in proportion to the magnitude of the thermal load, the stiffness of the end restraints, and the stiffness of the surrounding soil. Since the direction of thermally induced deformations above the null point is upwards, opposite to the direction of the weight of the foundation or mechanical load, shaft resistance decreases. Below the null point, where thermally induced deformation is in the same direction as mechanical loads, the shaft resistance increases.
Bodas-Freitas et al. [16] performed numerical simulations to investigate the effect of heating on the shaft resistance in an energy foundation. They developed a simplified formulation assuming both the pile and the soil were elastic and found that the shaft resistance decreased above the null point and increased below the null point. When the coefficient of thermal expansion of the soil was increased, the change in shaft resistance decreased. An opposite effect was observed when the Young's modulus of the soil was increased. One drawback of Bodas-Freitas et al. [16] is that they did not consider pore water pressures in their model. Holzbecher et al. [17] later extended the thermo-elastic model of Bodas Freitas et al. [16] to couple pore water pressure, stress, and strain, but did not consider the separate thermal expansion coefficients of the pore fluid and the soil matrix. Thus, their model cannot simulate thermally-induced excess pore water pressures.
Fuentes et al. [18] used a one-dimensional finite difference approach to study the effects of soil hydraulic conductivity and compressibility on pore water pressure generation due to heating. They found that an increase in pore water pressure is generated in soils with low permeability and high compressibility. However, Fuentes et al. [18] didn't consider thermally induced deformations.

Model description and formulation
In this study, poroelastic theory is extended to include thermally induced excessive pore water pressure and thermo-elastic strains in the soil skeleton to investigate changes in shaft resistance during heating of an energy foundation installed in a glacial till layer. The poroelastic theory is valid when the expected thermally induced strains are recoverable. The mass balance is coupled with conservation of energy. The 2D axisymmetric numerical model is implemented with COMSOL Multiphysics v5.4b and solved using the finite element method. In poroelastic theory stress, strain, and pore pressure are related as follows: where σ' is the effective stress tensor (kPa), ε is the strain tensor, C is the elasticity matrix from Young's Modulus E (kPa) and Poisson's ratio ν, αB is the Biot-Willis coefficient, pf is the pore water pressure (kPa), and I is the identity matrix. The thermo-elastic strains in the soil skeleton can be expressed as follows: where εel is the elastic strain from mechanical loads, εth is the thermal strain, βs is the soil thermal expansion coefficient (1/°C), T is the temperature (°C), and Tref is a reference temperature (°C). Since the displacement response is much faster than the fluid response, quasistatic equilibrium is assumed, and it is expressed as follows: where ρ (kg/m 3 ) is the density of the matrix, g (m/s 2 ) is gravitational acceleration, and F (kN) includes all external forces. After the conservation of mass is applied, the governing equation to include thermally induced volumetric strains is expressed as follows: where ρf is the density of the pore fluid (kg/m 3 ), S is a storage coefficient, u is the fluid velocity vector (m/s), εvol is the volumetric strain, and Qβ describes the thermal expansion of the soil constituents. Darcy's law is used to describe the fluid velocity: where K is the hydraulic conductivity (m/s). The storage coefficient from Equation 6 is given as follows: where εp is the porosity, Kf is the bulk modulus of the fluid (kPa), and Kd is the bulk modulus of the matrix under drained conditions (kPa). The thermal expansion source term from Equation 6 is given as follows: where βf is the fluid thermal expansion coefficient (1/°C). Finally, conservation of energy is applied: where Cp is the specific heat capacity of the soil (J/kgK), λ is the thermal conductivity of the soil (W/mK), and QTh is the heat source term (W/m 3 ).

Validation
To validate the numerical model, two studies focused on the investigation of energy foundations from the literature were selected and the results were replicated using the numerical model developed in this study. The first study focuses on prediction of thermally induced pore water pressures reported in Booker & Savvidou [8] while the second focuses on changes in shaft resistance reported in Bodas Freitas et al. [16].

Validation of pore water pressure
Booker & Savvidou [8] developed an approximate analytical solution for the evolution of pore water pressure around a buried canister with constant volumetric heat flux. The analytical solution was derived assuming the pore fluid and soil solid are incompressible and both gravity and convective heat flux are ignored. The geometry and mesh used to predict thermally induced excessive pore water pressures in [8] are shown in Figures  1(a) and 1(b). An initial temperature of 20°C and initial hydrostatic pressure distribution with groundwater at the surface were applied over the entire domain. A constant volumetric heat flux of 1,000 W/m 3 was applied to the buried canister. For heat transfer, a Neumann-type boundary condition was applied to all external boundaries. The top surface of the domain was free to displace in any direction. The lateral exterior boundary was constrained in the horizontal direction, and the bottom exterior boundary was constrained in the vertical direction. For subsurface flow, a Neumann-type boundary was applied to the exterior bottom and lateral boundaries. A fixed Dirichlet-type atmospheric pressure was applied to the top boundary. The input parameters for the simulation are shown in Table 1. The pressures from both numerical analyses and the original study were normalized by the theoretical maximum pressure that would occur under the same thermal load if the soil were completely impermeable and are plotted in Figure 2. This maximum pressure was calculated directly from the equations presented by Booker & Savvidou in their analytical formulation [8].
The pressures from the model are higher than the analytical solution by approximately 0.5% of the normalization pressure. The difference between the two solutions increases at the end of the analysis. The difference can be attributed to the temporal changes in intrinsic material properties with heating. Overall, the results from the analytical solution and numerical model are in good agreement. An initial temperature of 20°C and initial stress state at equilibrium were applied to the entire domain. A temperature of 20°C was maintained on the exterior boundaries of the soil domain. For heat transfer, the top of the foundation was constrained by a Neumann-type boundary. The displacement and subsurface flow boundary conditions at the external boundaries were equivalent to those used to replicate Booker & Savvidou [8]. A constant temperature of 50°C was applied to the outer boundaries of the foundation. This value is representative of regions experiencing very hot summers. A mechanical load of 6 MPa was applied to the top of the pile. The inputs for this simulation are shown in Table 2. The results from the numerical analyses and the original study are plotted in Figure 4. The two solutions are within 5 kPa of each other.

Parametric study
The numerical model was used to perform a parametric study by using different values of hydraulic conductivities and the rate of heating. Input parameters were equivalent to those in Table 2, except the porosity εp was used as 0.5 and the compressibility of the pore fluid Kf was set equal to 2x10 9 Pa. The hydraulic conductivity values used in the analyses are representative of those observed in glacial tills in Baser et al. [19]. The Biot-Willis coefficient αB was used as 0.99 based on reported values for the compressibility of dense sand and overconsolidated clay by Mitchell and Soga [20]. A temperature-dependent coefficient of thermal expansion for liquid water βf from COMSOL's material library was also used.

Geometry
The geometry used in the model is representative of one of the energy foundations at the University of Illinois, as shown during construction in Figure 5.

Initial and boundary conditions
An initial temperature of 12 °C, an initial hydrostatic pressure distribution with groundwater at the surface, and an initial stress distribution at equilibrium were applied to the entire domain. A fixed temperature was applied to the boundaries of the heat source inside the foundation to represent the fluid temperature. For heat transfer, the top of the foundation was constrained by a Neumann-type boundary. Constant pressure applied at the top of the soil domain. Boundary conditions for the displacement, mechanical load, and flow were identical to those used in the verification against [16].

Analysis and results
Pressures at the mid depth of the soil-foundation interface are shown in Figure 7 with different hydraulic conductivities and a constant fluid temperature. An increase in maximum pore water pressure from 105 kPa to 170 kPa was observed as the hydraulic conductivity was decreased from 1x10 -9 m/s to 1x10 -12 m/s.  Vertical profiles of pore water pressure at the foundation radius are shown in Figure 9. Excess pressures exceeded the total vertical stress in top 3 meters of the foundation. Profiles of the change in shaft resistance at the foundation radius are shown in Figure 10. After excess pore water pressures dissipate, the steady state condition is in agreement with the framework of Bourne-Webb et al. [15]. Before excess pressures dissipate, the magnitude and direction of changes in shaft resistance are variable. Long term changes are approximately 10-20 kPa, but short term changes are up to 5-10 kPa higher in the positive and negative directions after five days.

Discussion and conclusions
Based on the results of this study, the risk of failure in structures with energy foundations is higher in soils having low hydraulic conductivity and at relatively higher amounts of heat injection. Pore pressures exceeding the total stress on portions of the foundation were observed at hydraulic conductivities on the order of 5x10 -11 m/s. The magnitude of the changes in shaft friction as excess pressures dissipate are higher than the steady-state condition. The combined effects of soil deformation, foundation deformation, and pore water pressure dissipation contributed to the transient response.
The magnitude of excess pore water pressure increases with decrease in hydraulic conductivity and increases in rates of heating. Maintaining a constant fluid temperature of 35°C and changing the hydraulic conductivity from 1x10 -9 m/s to 1x10 -12 m/s resulted in maximum pressures of 105 kPa to 170 kPa, respectively. With a constant hydraulic conductivity of 5x10 -11 m/s, varying the fluid temperature between 20°C and 50°C resulted in maximum pressures of 106 kPa to 167 kPa, respectively.
After the dissipation of excess pressures, the shaft resistance decreased by 11 kPa above the null point and increased by 15 kPa below the null point. As pore water pressures increased over the first five days, shaft resistance decreased by 14 kPa above the null point and increased by 20 kPa below the null point. As the excess pressures dissipated, the change in shaft resistance changed direction until ultimately reaching the steadystate condition. Neglecting the effects of thermallyinduced excess pore water pressures and the combined deformation of the soil-foundation system could lead to unconservative designs.