Volume change behavior of unsaturated soils under non-isothermal conditions

With the development of modern geotechnical engineering practices such as the construction of high level radioactive waste repositories, exploitation and utilization of geothermal resources, energy-saving buildings and underground storage of CO2, research into the influence of temperature on the basic mechanical properties of unsaturated soils has become an important issue internationally. By using the work expression and considering the influence of temperature on the basic properties of unsaturated soils, the average soil skeleton stress, modified suction and temperature were selected as state variables of generalized forces in thermodynamics and the soil skeleton strain, saturation and entropy were chosen as state variables of generalized flows conjugate to the variables of generalized forces. Based on the nonlinear multi-field coupled model and by using existing experimental results, an elastic-plastic constitutive model of unsaturated soils under non-isothermal conditions was developed to analyze the influence of temperature on the deformation properties of unsaturated soils. The model was used to predict and analyze the influence of suction and temperature on the deformation properties of unsaturated soils under isotropic conditions, and was successfully verified using experimental results.

Soil layers on the earth's surface are affected by rainfall and evaporation and most are in an unsaturated state. Because unsaturated soil is composed of solid, liquid and gas phases, both its basic properties and engineering characteristics are quite different from those of a saturated soil. There are many problems of geotechnical analysis or design with unsaturated soils to be solved. Therefore, since the 1990s, study of the basic mechanical properties of unsaturated soils has become important in international geotechnical engineering [1,2], while predicting deformation and strength is one of the most important research subjects in unsaturated soil mechanics.
The development of modern geotechnical engineering brings new challenges to traditional soil mechanics. For example, in the construction of high level radioactive waste repositories, unsaturated bentonite, as a buffer/backfill material, is an important man-made barrier. When subjected to heat from the nuclear waste, its temperature will increase and its water content gradually decrease. Eventually, as the groundwater level rises, the overlying soil wets and then becomes saturated. During this process, under the joint influences of temperature change, capillary effect and exterior load, the deformation properties of the unsaturated soil will change. In many engineering construction practices such as the exploitation and utilization of geothermal resources, CO 2 underground storage, highway sub-grade stability analysis, high-voltage cable burial and petroleum exploitation, research on the constitutive model of unsaturated soils under non-isothermal conditions is of important theoretical significance and application value.
Unsaturated soil is always found in certain natural environments and its deformation behavior will be influenced by moisture, temperature and other environmental factors. Unsaturated soil mechanics is the result of considering moisture change of a saturated soil and is therefore based on saturated soil mechanics. The influences of temperature on the deformation and strength properties of unsaturated soils are also very significant [3]. As the temperature changes, interactions among the three phases within unsaturated soils change accordingly, resulting in differences in their preconsolidation pressure and further affecting the evolution law of the yield surface during the deformation process, making the deformation properties of unsaturated soils dependent on temperature to a certain extent. The influence of temperature on soil mechanical properties has long drawn attention [4], but because of experimental limitations, progress has been very slow.
Hueckel and Borsetto [5,6] proposed an early elasticplastic constitutive model of soil, considering heat effects. Laloui et al. [7] and Cui et al. [8] have also studied the effect of temperature on soil deformation. Current research on the effect of temperature on soil properties has focused largely on saturated soils, but seldom on unsaturated soils, which is the more general condition. Furthermore, the existing constitutive models of unsaturated soils also rarely consider the impact of temperature, so that the interaction of temperature and suction on soil deformation properties are difficult to describe. Therefore, it is particularly urgent and important to develop a constitutive model of unsaturated soils taking account of temperature influences. Using advanced experimental techniques, Tang and Cui [9], François et al. [10], and Uchaipichat et al. [11] carried out a large number of experiments on the effect of temperature on the fundamental properties of unsaturated soils, which have laid the foundations for the development of constitutive models.
Multi-phase porous media theory is used to describe the complex interactions between different phases in unsaturated soils and their responses to external stress or chemical actions. It appears able to unify the theoretical basis of unsaturated soil mechanics [12][13][14][15][16][17]. Recently, multi-phase porous media theory was used by our research team as the basis for some useful insights into unsaturated soil mechanics including the principle of generalized effective stress and constitutive relations for unsaturated soils [18][19][20], continuum porous medium soil mechanics and its application in constitutive relationships of unsaturated soils [21], and the nonlinear multi-field coupled model for multi-constituent three-phase soils [22]. This paper is based on the above research and is centered on assessment of the theory framework proposed by Cai et al. [22]. It presents, from a constitutive point of view, an elastic-plastic constitutive model of unsaturated soils under non-isothermal conditions. Two coupled constitutive aspects are used to fully describe non-isothermal behavior, including both the mechanical constitutive part of the soil skeleton and an elasto-plastic model for water retention characteristics. In this paper, we focus on the influence of temperature on the deformation properties of unsaturated soils. The model is used to predict and analyze the influence of suction and temperature on the deformation properties of unsaturated soils in isotropic conditions. Validation of the model is based on existing experimental results.

Selection of stress state variables
According to thermodynamic and multi-phase porous media theory, the selection of stress and strain variables should meet the requirements of the deformation work in unsaturated soils. Based on the porous media theory, Zhao et al. [18] derived an expression for the total deformation work: where W is the total deformation work,  is the total stress, p l , and p g are the pressure of fluid and gas respectively, S r is saturation,  is a unit tensor, n is porosity, s = p g   p l is matric suction, K g is the volume compression modulus, and v s is the velocity of the solid phase. Eq. (1) shows that there are three components contributing to the total deformation work for unsaturated soils under isothermal conditions. The first component denotes the work produced by the deformation of the soil skeleton, the second is produced by the change of water content, and the final contribution is from pressure variation of the gas phase. Unsaturated soils are composed of solid, liquid and gas, and there are interactions among each of these phases. Usually, the interaction between liquid and gas is established through the soil water characteristic curve, describing the dependency of the degree of saturation on suction. In this paper, we follow convention in unsaturated soil mechanics and neglect any changes in gas pressure, while the effect of temperature, both in deformation of the soil skeleton and in changes in the liquid phase, is considered. Within this context, the generalized stress and strain variables expressed by matrix takes the form: where  is the effective stress, which can be defined as are the axial stress, radial stress and gas pressure, respectively;  a and  r are the axial and radial strains.

Elastic constitutive equations
Elastic-plastic coupling is not considered in this paper; that is to say, the elastic modulus is not affected by the plastic strain. Furthermore, we assume that the elastic modulus and the elastic deviatoric strain are not affected by temperature, so the elastic stress-strain relationships for the solid phase are written as where e w K is the elastic modulus for liquid, and its specific form is determined by the selected model for the soil-water characteristic curve. Here we also neglect temperature effects on elastic deformation of the liquid phase.

Elastic-plastic constitutive equations
A nonlinear multi-field coupled model for multi-constituent three-phase soils has been derived by using the hybrid mixture theory in [22]. Based on this model, this paper considers the deformation of soil skeleton and the change of liquid content in unsaturated soils under non-isothermal conditions, and neglects other effects on the soils. At the same time, we assume that temperature is homogeneous in space, so heat conduction is not considered. Then, the nonlinear multi-field coupled model in [22] is simplified as where X  is the effective generalized thermodynamic force, and Y is the generalized thermodynamic flux coupled with X  . Generally, stress is taken to be a controlled variable when we describe the stress-strain constitutive relationship for soils. So, taking generalized forces as independent variables is more convenient, and generalized fluxes are the function of generalized forces. Hence the generalized fluxes can be derived from the dissipative function * ( )  D X , with the generalized forces  X as independent variables, as sions for the coupling effect among different fields of soils can be derived in the same way as in [22], just needing to exchange the generalized thermodynamic forces and generalized thermodynamic fluxes. The expression for generalized thermodynamic forces with generalized thermodynamic fluxes as independent variables can then be obtained as in the following equations (6) and (7). When building specific constitutive models, corresponding equations for the respective solid and liquid phases should be established. For the solid phase, the generalized thermodynamic flux that should be considered is the plastic strain d p (for triaxial tests, it can be expressed as the sum of the plastic volumetric strain d p v and plastic shear strain d p s ), and the generalized thermodynamic forces include effective stress d , modified suction d s and the change of temperature dT. Therefore, the coupled constitutive relationship for the solid phase is given by where A, B, C are coupling coefficients. For the liquid phase, the generalized thermodynamic flux that should be considered is the plastic change of the degree of saturation d p r S . From the microscopic view of deformation mechanism, the elasto-plastic phenomenon of the liquid phase is the hydraulic process of water inflow and outflow to individual voids. For a single pore channel, when the water content in soil changes, the radius of curvature of the interfaces between gas and water changes first. This kind of movement of the interfaces is a reversible (elastic) process and so the corresponding saturation change is also reversible and can be expressed by d e r S . If the water content changes beyond a critical value, where the radius of the interface is just sufficient to bridge the largest entry route to the void, gases will break through into the void during drying and during water movement into the void during wetting. The breakthrough of gas or water into a void during drying or wetting is an irreversible process that produces plastic changes in the degree of saturation d p r S , because the void will not immediately refill with water or gas if the suction changes.
Based on the above deformation mechanism, Wheeler et al. [2], Sheng et al. [23], Wei et al. [24], Sun et al. [25] established constitutive models for unsaturated soils. In this paper, the generalized thermodynamic forces, which causing the change of d p r S , include effective stress d , modified suction d s and the change of temperature dT. Therefore, the coupled constitutive relationship for the liquid phase is given by where A′, B′, C′ are coupling coefficients. When considering the temperature effect on the plastic strain of unsaturated soils and the interactions between solid and fluid phases, the following method is adopted: for the solid phase, the effect on its plastic strain from both the liquid phase and temperature is expressed by the function for yield stress for the solid phase c p  being dependent on both the degree of saturation p r S and temperature T, while for the liquid phase, the effect of both solid phase and temperature to the plastic change of degree of saturation is expressed by the function for yield stress for the liquid phase y s  being dependent on both the plastic volume strain p v  and temperature. Under non-isothermal conditions, the incremental constitutive equations coupling the capillary effect and elasto-plastic deformation for unsaturated soils are derived in the next section, based on the yield condition, hardening law, flow rule and consistency condition.

Yield conditions
(1) Yield condition in the p q   plane. The yield surface in the p q   plane reflects the yield condition for the solid phase in unsaturated soils. In soil mechanics, the use of thermodynamics to develop constitutive relation for the solid phase has been studied by several researchers [16,17,21,26]. Compared with traditional methods, the thermodynamic method is more rigorous in theory and can effectively avoid the situation that equations do not satisfy the fundamental laws of thermodynamics.
Zhao et al. [21] derived the increment equation for Helmholtz free energy. Based on the increment equations for dissipative and stored energy proposed by Collins [26], the yield surface equation in the dissipative stress space can be derived: p  is the preconsolidation pressure, M is the slope of the critical state line, and  is a material parameter that ranges from 0 to 1; parameter , which has a value between 0 and 1 with = 1 for clays, is related to the amount of the stress ratio and reflects the stored energy in the proportion of the total plastic work;  and  are the average effective stress and deviatoric stress in the dissipative stress space. Then, the yield condition in the dissipative stress space is transformed into true stress space by the shift stress, so the yield surface equation in the true stress space is as follows: (2) Yield condition in the p s    plane I-LC yield curve.
The LC (loading collapse) yield surface in the p q s     space describes the slip between particles or between aggregates composed of particles. Under isotropic compression conditions, it can be simplified to an LC yield curve in the p s    plane. The LC yield curve can describe not only the plastic shrinkage when wetting under a high level confining pressure condition, but also the increment of preconsolidation pressure resulting from an increment of suction. As the preconsolidation pressure represents the limit of the elastic region, the LC yield curve reflects the strengthening effect of suction on the soil to some extent. According to Loret and Khalili [27], we adopt the following equation for the LC yield curve: (3) Yield condition in the p s    plane II-SI/SD yield curve. The yield equation for the liquid phase is described by the soil-water characteristic curve, and the linear model proposed by Wheeler et al. [2] is adopted here, in which the slope of the main drying/wetting curve and the canning curve is denoted by  w and  w , respectively. On the canning curve, the change of degree of saturation is elastic, while plastic change occurs on the main drying/wetting curve. The effect from the elastic deformation created by the increase of p  on the modified suction is not considered in this paper. Based on this assumption, the SI (suction increase) and SD (suction decrease) yield curve in the p s    plane is horizontal, and the yield equation for the liquid phase is given by , where y s  is the yield stress for the liquid phase when y=I, and I s  is the yield stress created by a suction increase; when y=D, D s  is the yield stress created by a suction decrease. The plastic change for the degree of saturation caused by suction changes can be expressed as is the plastic modulus of the liquid phase.
(4) Yield condition in the Tp plane I-LY yield curve. A new plastic mechanism is proposed to consider the temperature effect on the elasto-plastic deformation for unsaturated soils. For simplicity, and in relation with the loading collapse (LC) yield curve, this new yield curve is called the loading yield (LY). Below LY, the behavior of a soil subjected to increments of temperature and mean effective stress is purely elastic. The LY yield curve can be determined by the preconsolidation pressure changes with soil temperature at a constant plastic condition. Based on the experimental results of Hueckle & Baldi (1990) [6], the fitted exponential expression for the LY yield curve is adopted here: where ( ) c p T  and 0 c p  are the preconsolidation pressure at the tested temperature and at reference temperature, respectively, and  Tp is a model parameter.
(5) Yield condition in the Tp plane II-TY yield curve. Experimental results reported from different types of soils shows that the thermally induced volumetric strain at a constant elevated temperature depends primarily on the stress history (OCR) [8,28]. As the soil changes from a normally consolidated state to an over-consolidated state, the thermally induced contractive volumetric strain keeps decreasing and beyond a certain OCR value it exhibits dilative behavior. In fact, there is another kind of yield mechanism caused by temperature for over-consolidated soil, and it is believed that any temperature increase above the maximum temperature ever supported by the soil can generate a plastic strain. This yield curve, called the thermal yield (TY), has been proposed by Cui et al. [8] for saturated soils where T CT is yield temperature, T 0 and T c are the initial and reference temperatures, respectively, and  T defines the effect of p/p 0 (the ratio of net mean stress at the current state to the net mean stress at the initial state) on the yield temperature.
Based on the research of Tang et al. [9], the temperature yield curve for unsaturated soils can be expressed as: where s 0 and s are the initial suction and the present suction, respectively, and parameter Ts  defines the suction effect on the yield temperature.

Hardening law
The hardening law is proposed to describe evolution of the yield surface.
(1) Hardening for the soil skeleton.In this paper, the unsaturated soil is assumed to be under non-isothermal condi-tions, and the yield conditions for soil at different degrees of saturation and at different temperatures should be considered. In this case, the yield stress of the soil skeleton c p  depends not only on the volumetric plastic strain, but also on the degree of saturation and temperature. Therefore, the hardening law for the soil skeleton consists of three parts: volume strain hardening term d cs p for the soil skeleton, d cw p for the liquid phase and d cT p caused by the change of temperature. The hardening law for the soil skeleton is assumed to be the linear superposition of these three parts: For the hardening term for the soil skeleton ( d cs p ), the plastic volume strain hardening in the Cam-Clay model is adopted here. For the hardening term for the liquid phase ( d cw p ), suction is selected as the hardening parameter in many models, but the capillary hysteresis from the drying/wetting cycle cannot be considered in this way. In contrast, the plastic change of degree of saturation, which corresponds to filling or emptying of voids with water, is much more important than the suction in producing a change in the stabilizing effect of meniscus water lenses. So it is more reasonable to choose the plastic change of degree of saturation as the hardening parameter when considering the effect from liquid deformation to the yield stress of the soil skeleton. For the hardening term caused by the change of temperature, the experimental result of eq. (13) is adopted. Therefore, the hardening equation for the soil skeleton can be written as where , (2) Hardening for the liquid phase.When developing the hardening law for the liquid phase in unsaturated soil under a non-isothermal condition, the effects on the SWCC from soil skeleton deformation and temperature should be considered. So, the hardening equation for the liquid phase is where 0  y s is the initial yield stress, and k ws and  Ts are coupling coefficients. Wheeler et al. [2] established a hardening equation for the liquid phase that is similar to eq. (18), but it did not reflect the impact of temperature.
(3) Hardening caused by movement of the TY yield surface.When the thermo-mechanical path crosses the TY yield surface, a hardening phenomenon moves the TY surface upwards. Cui et al. [8] supposed that this movement is governed by the hardening parameter  T , and the following expression for the plastic volume strain caused by a change of temperature was derived based on experimental results: where a is a material parameter, and 0 CT T is the initial yield temperature. Following Cui et al. [8,9], we assume that increases in net mean stress and suction have no effect on the movement of the TY yield surface. So the hardening law caused by the movement of TY yield surface is given by where the parameter /1

Flow rule and consistency condition
A flow rule is used to determine the magnitude of plastic strain increment. When the unsaturated soil is under non-isothermal conditions, there is interaction between the solid and liquid phases, and temperature affects the two phases. For the solid phase, the plastic strain is given by the sum of any direct movement caused by yielding on the LC or LY curve itself and any coupled movement produced by yielding on the SI/SD or TY curves. Similarly, for the liquid phase, the plastic change of degree of saturation is given by the sum of any direct movement caused by yielding on the SI/SD curve itself and any coupled movement produced by yielding on the LC, LY or TY curves. So, the general expression for the plastic strain increment is S are the plastic change of the degree of saturation caused by yielding on the LC and LY curve, the SI/SD curve and the TY curve, respectively. For the non-expansive clay in this paper, the plastic strain created by the SI/SD yield is not caused by the direct movement of the SI/SD curve but by the coupled movement of the LC curve, so . Similarly, the plastic change of the degree of saturation created by the LC yield is not caused by direct movement of the LC curve but by the coupled movement of the SI/SD curve, so . The plastic change of the degree of saturation and the shear strain for solid phase produced by the temperature yield are neglected here, so Based on the above analysis, eq. (21) can be simplified to where g s , g w are the plastic potentials for the solid and liquid phase, respectively,  s and  w are the corresponding plastic factors, and d p vT is the plastic strain caused by temperature change expressed by eq. (19).
For the solid phase, the non-associated flow rule is adopted, and the plastic potential g s is given by eq. (8), while for the liquid phase, the associated flow rule is adopted, so According to the yield eqs. (9) and (11) for the solid and liquid phases, the consistency condition can be written as Substituting eqs. (17), (18) and (22) into eq.(24) yields variables in eq. (29) can be calculated based on the stress incremental variables. There are 17 model parameters in the derived model, including the elasto-plastic parameters of saturated soils (, , M, v, u, ), parameters of the LC yield curve (k m ,   ), coupling parameters (k sw , k ws ,  Tp ,  Ts ), elasto-plastic parameters of the liquid phase ( w ,  w ), the thermal expansion coefficient (′ s ) and parameters of temperature yield (  T , a). The triaxial test can be used to determine the elasto-plastic parameters of saturated soils (, , M, v, u, )); the isotropic consolidation test for unsaturated soil is used to determine the parameters of the LC yield curve (k m ,   ) and the coupling parameters (k sw , k ws ); the drying-wetting cycle test is used to determine the elasto-plastic parameters of the liquid phase ( w ,  w ); under different temperatures, the consolidation test, the drying-wetting cycle test and thermal expansion test are used to determine the thermal expansion coefficient (′ s ) and the parameters of temperature yield ( T , a).

Validation of the constitutive model
The proposed model has been validated using the results of existing experiments. In this section, numerical simulations based on the program compiled by Fortran language are compared with the isotropic compression test data for the MX80 clay from Tang and Cui [9].  Figures 1 and 2 reproduce the numerical simulations of isotropic compression tests at two temperatures (T=25°C and T=60°C) and at different suctions. The simulations show good agreement with experimental data. When T=25°C, the yield stresses are 1.3, 7.2 and 18 MPa corresponding to suctions of 9, 39 and 110 MPa, respectively, and when T=60°C, the yield stresses are 5.6 and 11 MPa corresponding to suctions of 39 and 110 MPa, respectively. At a given temperature, an increase in suction induces an increase in yield stress, which means that the suction increases the yield strength of the soil. When the temperature increases, the suction effect on the soil reinforcement remains. Figures 3-5 reproduce the numerical simulations of isotropic compression tests at three suctions (9, 39 and 110 MPa) and at different temperatures. The simulations show good agreement with experimental data. When s=9 MPa, the yield stresses are 1.3 and 1.1 MPa corresponding to temperatures of 25°C and 80°C, respectively, when s=39 MPa, the yield stresses are 7.2, 5.6 and 3.8 MPa corresponding to temperatures of 25°C, 40°C and 60°C, respectively, and when s=110 MPa, the yield stresses are 18 and 11 MPa corresponding to temperatures of 25°C and 60°C, respectively. Thus, the yield stress decreases with an increase of temperature at a fixed suction, which shows that temperature decreases the yield stress of the soil and makes the soil easier to yield. The change of temperature has no effect on the slope of the compression line of the soils.

Conclusions
Based on experimental evidence, a constitutive model coupling the capillary effect and elasto-plastic deformation for unsaturated soils under non-isothermal conditions was developed. The effect of temperature on volume change behavior was focused on. Existing experimental results were used to determine the model parameters and the proposed model was validated by predictions of experimental results on MX80 clay. The main conclusions are as follows: (1) According to thermodynamic and multi-phase porous media theory, the selection of stress and strain variables should meet the requirements of deformation work in unsaturated soils. By using the work expression and considering the influence of temperature on the basic properties of unsaturated soils, the average soil skeleton stress, the modified suction and temperature were selected as state variables of generalized forces in thermodynamics and the soil skeleton strain, saturation and entropy were chosen as state variables of generalized flows conjugate to the variables of generalized forces.
(2) Under non-isothermal conditions, there are five different yield mechanisms arising from different stress paths in the p q s T     stress space for unsaturated soils. The yield in the p q   plane describes the yield conditions for the solid phase, and the yield equations based on thermodynamics were given by eqs. (8) and (9). The LC yield in the p s    plane reflects the strengthening effect of suction on the soil and the corresponding yield condition is expressed by eq. (10). The SI/SD yield in the p s    plane reflects the yield conditions for the liquid phase, which can be described by the Soil Water Characteristic Curve, and eq. (11) is the corresponding yield equation. The LY yield in the Tp plane reflects the weakening action arising from a temperature increase in the soil and is expressed by yield equation (13). The TY yield in the Tp plane reflects another yield mechanism, which is different from the LY yield, and the corresponding yield equation considering the effects of the current stress state and suction is eq. (15).
(3) In the proposed model, the effects of temperature and the other phases are considered in the hardening law for both the solid and liquid phase, expressed by eqs. (17) and (18).
(4) For friction materials, the increment dissipation function depends on the true stress. So, the non-associated flow rule is adopted for the solid phase in the true stress space, while the associated flow rule is adopted for the liquid phase and the yield caused by temperature change.
(5) The constitutive model coupling the capillary effect and elasto-plastic deformation for unsaturated soils under a non-isothermal condition is presented as eq. (29). The 17 model parameters can be fitted using experimental data. Comparisons between numerical simulations and experimental results show that the proposed model can reproduce the volume change behavior of unsaturated soils subject to any coupled suction, temperature and external stress variations under isotropic conditions.