Thermo-Hydromechanical Coupling Responses Driven by a Central Heat Source in Unsaturated Soils

Based on the thermodynamic and thermoelastic theory, the coupled governing equations of deformation, heat transfer, and moisture migration in unsaturated soils were given. &e coupled calculation process was realized by the finite element method. &e hydrothermal coupling characteristics of two cases of unsaturated kaolin clay were studied by using a self-developed test device, and the test results were compared with the numerical results.&e results showed that when the initial saturation of the soil is high, the volumetric water content of the measured point increases, and then it is in a stable state during the heating process. When the initial saturation is low, the volumetric water content near the heat source increases firstly and then decreases during the heating process. &e rise and fall of temperature make the volumetric water content of the soil irreversible. &e volumetric water content of each measured point is lower than the initial state. &e closer it is to the internal measured heat source, the more obvious this phenomenon is; the lower the initial saturation is, the more obvious it is.


Introduction
e study on the unsteady responses of internal temperature transfer and moisture migration induced by the coupling effect of suction potential, gravitational potential, and temperature potential in unsaturated soils has very important significance. e application respects are involved in geothermal resource development [1], thermal energy storage [2], nuclear waste disposal [3], and heating pipeline design [4]. e study mainly relates to the coupling of different physical fields, for instance, hydraulic and mechanical thermal responses.
In order to analyse the responses of various physical fields, scholars have studied them from different angles. Chen et al. [5] constructed a complicated multifield coupling model for unsaturated porous media considering chemical and hydration swelling effects based on mixture coupling theory, nonequilibrium dynamics, and Biot's elastic theory. Chen et al. [6] proposed a method for analysing thermo-hydromechanical processes in porous media based on the continuum mechanics of the conservation of momentum, energy, and mass and the averaging method in the three-phase mixture theory (solid-liquid-gas) system. is model considered the effects of stress-strain, water-air-vapour phase flow, thermal transfer, and porosity evolution on the coupled response. Rutqvist et al. [7] gave basic governing equations of the thermo-hydromechanical coupling process of unsaturated fractured rock and soil and compared four finite element programs. Chen et al. [8] conducted experimental evaluations on several models describing the relationship between permeability and capillary pressure-saturation, which showed that LNM, VGM, GDM, and BRB models can better describe the flow of two-phase fluids. Wu et al. [9] established a coupled thermo-hydromechanical stress-strain relationship of unsaturated soils, which can describe the influence of heat on the hydraulic characteristics of soils and can quantitatively simulate thermal softening phenomena. Finally, the test results were compared with the numerical simulation results. Akrouch et al. [10] discussed the coupled hydrothermal characteristics of energy piles to unsaturated soils by experiments, theory, and numerical simulation and compared the experimental results with analytical and finite element solutions. Ghorbani et al. [11] established a theory that can describe the hydraulic interaction behavior for typical unsaturated soils under static and dynamic loads on the basis of previous studies. By solving the dynamic problems of unsaturated soils with elastic frameworks and nonlinear frameworks, the static problem of the saturated soil verified the correctness of the established model.
Bai [12] discussed the semianalytical solutions of porous and saturated spheres with cavities under several boundary conditions and compared the effects of thermal response on field variables such as radial displacement, pore pressure, and tangential or radial stress through numerical analysis, which is a pioneering work. Caulk et al. [13] proposed a pore-scale numerical method for simulating the thermohydromechanical coupling behavior of granular media on the basis of the existing fluid mechanics model. For unsaturated soils, Takayama et al. [14] established a constitutive model based on the Cambridge model to describe the changes of the mechanical properties in view of bentonite in the nuclear waste repository due to the hydration process. Zeng et al. [15] established a vertical one-dimensional, mass flow heat, and two-phase model considering diffusion, advection, and dispersion and used the model to discuss the advection effects of high-permeability soil as well as lowpermeability soil. For saturated porous media, Bai [16] established a normalized coupled governing equation and discussed the theoretical response under cyclic thermal load for one-dimensional cyclic thermal consolidation. Among these, the evolution process of displacement, pore pressure, and temperature from an unsteady state (i.e., instantaneous state) to a quasi-steady state was analysed. An et al. [17] established a two-dimensional hydrothermal coupling model and used meteorological data to estimate the surface heat flow and water flow boundary. On this basis, the variation rules of temperature as well as the volumetric water content of embankment were also studied. Zhou and Ng [18] established a new coupled thermomechanical theory of saturated soil based on the boundary surface plasticity theory. e shape and size of the boundary surface of this theoretical model can change with temperature, which can predict the elastic-plastic response of saturated soil under small deformation. For unsaturated porous media, the multifield coupling process is greatly affected by saturation degree. e effect of saturation on the coupled thermohydromechanical process in porous media is of great significance to understand the properties of the multifield coupling response, even including the influence of chemical factors [19][20][21]. However, all the constitutive models described above are derived based on the general elastic-plastic theory.
In this paper, the relationship between the degree of saturation and pore gas pressure as well as pore water pressure in the soil-water characteristic curve [22] is utilized to derive the seepage equation of unsaturated porous media through the continuity equation of pore water and pore gas. Based on the combination of thermodynamics and thermoelastic theory, the governing equations with coupled thermo-hydromechanical characteristics as well as state equations of the water-gas two-phase flow in porous media are further given. e coupling calculation process is realized by the finite element method. en, a self-developed test device was utilized to analyse the hydrothermal coupling characteristics of two different unsaturated kaolin clay under the same temperature load, and the obtained experimental results verified the rationality of the proposed theoretical model. is work provides further insight into the hydrothermal effects of kaolin soil from the multifield coupling point of view.

Seepage Differential Equation.
Under the condition of ignoring the evaporation of water in porous media, the mass conservation equation of pore water, pore gas, and solid matrix can be expressed as [7,8] z nS w ρ w zt where ∇·is the divergence calculation; n is the porosity of the porous medium; ρ w , ρ s , and ρ g are the density of pore water, solid matrix, and pore gas (ML ), respectively; S w and S g are pore water and pore gas saturation; and v w , v s , and v g are the velocity vectors of pore water, solid substrate, and pore gas (LT −1 ), respectively. Considering the effect of gravity, the velocity of the pore water as well as pore gas phase to the solid substrate can be expressed as [23][24][25][26] According to [27][28][29][30][31], where K is the porous medium permeability (L 2 ); μ w and μ g are the pore water and pore gas viscosity (ML −1 T −1 ), respectively; K rw and K rg are the relative permeability of the liquid and gas phase, respectively; P w and P g are the pore water and pore gas pressure (ML −1 T −2 ); v rw and v rg are the relative velocity of the liquid and gas phase, respectively; z is the height relative to the reference plane; g is the gravitational acceleration; u g and u w are the average relative flow 2 Advances in Materials Science and Engineering rates of pore gas and pore water (LT −1 ) in the cross section of the porous medium, respectively. According to [32][33][34][35][36], it can be known that pore water saturation should be a function of temperature and matrix suction, which has the following relationship: where T is the current absolute temperature and P c is the matrix suction (ML −1 T −2 ). en, the derivative of the saturation of pore water S w with respect to time is where C P and C T are the water-holding characteristic coefficients of capillary action (M −1 LT 2 ) and temperature action (θ −1 ), respectively, which are determined by the soilwater characteristic curve. Substituting equations (4)-(8) into equations (1) and (2), the seepage differential equations of pore water as well as pore gas can be obtained as follows:

Momentum Conservation Equation.
For unsaturated porous media, the momentum conservation equation without considering the inertial force can be expressed as [7,8] Considering the skeleton deformation, thermal strain, wet strain, pore pressure deformation, and equivalent stress principles of porous media, the stress-strain relationship is where where β T and β M are the linear thermal expansion coefficient and linear wet expansion coefficient of the porous solid matrix, respectively; D is the tangential stiffness matrix, which is the weighted average pore pressure (ML −1 T −2 ); S w0 and T 0 are the initial saturation and initial temperature, respectively; and K m is the bulk modulus of the matrix of the porous medium (ML −1 T −2 ).

Energy Conservation Equation.
Considering the change rate of thermal strain and wet strain energy, for unsaturated porous media, the energy conservation equation can be expressed as [7,8] z zt where e w , e g , and e s are the internal energy per unit mass of pore water, pore gas, and solid matrix (L 2 T −2 ), respectively, e � c v ·(T − T 0 ); h w , h g , and h s are the heat content of pore water, pore gas, and solid matrix per unit mass, respectively (L 2 T −2 ); Q is the total heat source intensity (L −1 MT −3 ); k t is the total thermal conductivity of porous media (MT −3 θ −1 ); c v and c p are the specific heat of constant volume and pressure, respectively (L 2 T −2 θ −1 ); and c is the specific heat rate.
On the left side of equation (13), the first term is the change rate of the internal energy of the system, the second term is the change rate of reversible elastic energy caused by the solid phase change [37][38][39][40], the third term is the heat conduction, the fourth term is the energy change produced by pore water pressure [41][42][43][44], and the fifth term is the difference between the incoming and outgoing energy per unit time.

Advances in Materials Science and Engineering
According to previous studies [45][46][47][48][49], considering the effects of saturation and porosity degree on the effectivity of the thermal conductivity, the total thermal conductivity can be expressed as where k s , k g , and k w are the thermal conductivity of the porous media matrix, pore gas, and pore water (MT −3 θ −1 ) and a 1 and a 2 both are experimental fitting parameters, respectively.

Water Retention Characteristics.
In view of the VG model [50][51][52][53], the relationship between matrix suction and relative permeability and relative effective saturation can be expressed as where S ew is the relative effective saturation; m and n are the fitting parameters, where m is related to the asymmetry of the soil-water characteristic curve and n is related to the slope of the curve; and H c is the matrix suction head. Pore water saturation degree can be written as where S wr and S ws are the residual saturation and maximum saturation, respectively. According to equation (8), it is assumed that the parameters α, m, and L have approximately linear relationships with temperature, which can be commonly written as When H c > 0, according to equations (8), (14), (15), and (16), one has

State Equations.
In order to simulate the real situation effectively, each state quantity adopts a nonlinear calculation formula related to temperature. According to [54][55][56], the density ρ of pore gas, pore water, and solid matrix, the viscosity coefficient of pore gas and pore water μ (Pa·s), the heat capacity at constant pressure c p (J/ (kg·K)), and thermal conductivity k (W/(m·K)) can be expressed as, respectively, 4 Advances in Materials Science and Engineering where where R g is actually the gas constant of air (L 2 T −2 θ −1 ); σ ′ is the effective stress acting on the matrix in the porous medium (L −1 MT −2 ), in which σ 0 ′ is the initial value of the effective stress (L −1 MT −2 ); ρ s0 is the initial density of the porous medium matrix (ML −3 ); and tr(·) is the trace operation in mathematics. According to equations (3), (11), and (19), the porosity can be expressed as Advances in Materials Science and Engineering 5

Test Apparatus Description.
A test device with a cylindrical chamber is used to study the water-heat coupling effect in unsaturated soils, as shown in Figure 1. e upper and lower ends of the test device are fixedly connected with the stainless steel plate through movable bolts, and an O-shaped rubber ring is set for sealing treatment. e test device is made of stainless steel (304 stainless steel), with an inner radius of R � 50 cm, a cylinder of 1 cm wall thickness, and a height of H � 50 cm. In order to meet the axisymmetric temperature, hydraulic, deformation, and migration boundary conditions, the stainless steel bottom plate and top cover are inlaid with a heat-insulating and impervious PC board. A layer of asbestos material is placed between the top cover and the bottom plate and the PC board. e centre of the test device is a cylindrical heat source T in with a radius of r � 2 cm, whose temperature is applied by using a constant temperature water bath circulation device, which can provide a long-term constant temperature possessing an accuracy of ±0.1°C. Besides, the outer temperature boundary of the test device also circulates through the constant temperature water bath. e device is applied and set to T out � 25°C.
In order to analyse the heat-water coupling effect in unsaturated soils, temperature sensors (K-type thermocouples) and volumetric water content sensors (EC-5) were embedded at different radial distances on the middle plane of the soil during the filling process of the soil sample, as shown in Figure 1. Four measuring points of temperature and moisture were arranged from the edge of the test tube to the centre, with an interval of 10 cm. By the temperature and moisture collector, the moisture and temperature changes at different locations were recorded in real time.

Test Schemes.
In the experiment, washed kaolin clay from Lingshou County, Hebei Province, China, was used as the porous medium, and the relevant property parameters are shown in Table 1. A mortar mixer is used to mix kaolin clay with a particle size of 37.5 μm into two different saturation degrees. rough experiments, the changes in temperature and seepage field of kaolin clay with two different saturation degrees are studied during the heating and cooling process. e two cases include (1) the kaolin sample with the initial saturation of 0.47 (mass moisture content 0.130), the internal temperature boundary of 80°C, and the outer temperature boundary of 25°C. e coupled properties of the seepage and temperature field during the heating process are studied. Both internal and external temperature boundaries are at 25°C. Meanwhile, the coupling response of seepage and temperature field is considered during the cooling process; (2) the initial saturation of the kaolin sample is 0.29 (mass moisture content 0.08), the internally measured temperature boundary is 80°C, and the outer temperature boundary is 25°C. e schemes are emerged in Table 2.

Physical Model and Boundary Conditions.
is paper uses "COMSOL Multiphysics 5.6" finite element software to solve the governing equations in Section 2 and to simulate the experimental process. e calculation model adopts an axisymmetric model. e inner heat source and the outer cylindrical wall were set as impervious and incapable of displacement boundary conditions; the top cover and bottom plate were set as thermally insulated with impervious boundary conditions. e temperature boundary of the inner heat source was set to 80°C and 25°C during the rising and falling process. e temperature boundary of the outer cylindrical wall was set to 25°C, and the initial temperature of the soil was 25°C. us, the optimization results of relevant parameters required in the calculation are shown in Table 3. e model mesh was divided into a free triangle mesh, the maximum element size was 2.65 cm, and the minimum element size was 0.15 mm (Figure 2). e complete grid contains a total of 884 domain elements and 74 boundary elements.
In the calculation process, both the internal heat source temperature boundary T in and the outer temperature boundary T out were set as the first type of boundary conditions. As a result, the internal heat source temperature boundary is within the first 5.15 d, and the temperature was maintained at about 80°C. During the test, the temperature was basically maintained at about 25°C. It should be noted that, during the 5.15 d cooling process, the water bath needs to be replaced with cold water to achieve the cooling effect. e constant temperature water bath had no cooling function, so the internal temperature of the heat source was measured. After the temperature drop starts, it first rises and then gradually decreases to 25°C, as shown in Figure 3(a), while the outer temperature boundary has been maintained at about 25°C (Figure 3(b)), which was arranged on the outer cylindrical wall of the test device (Figure 3(b)). e temperature measured by sensors changes over time at four different locations.

Discussion.
According to the test results, it can be found that the volumetric water content of each monitoring position in the soil has changed after the temperature rising and falling process. Due to different soilwater characteristics of unsaturated soils in the process of moisture absorption and dedrying, the VG model parameters are different during the heating and cooling processes. After trial calculations, the parameter α 1 � 1.89 m −1 during the heating process; it is taken as α 2 � 1.93 m −1 during the cooling process. In the calculation process, the initial pore water pressure as well as pore gas pressure was 0.94 × 105 Pa and 1.01 × 105 Pa, respectively. Figure 3 shows the comparison between the calculated results and the experimental results of the volumetric water content at different monitoring positions in the soil for case 1. By comparison, it can be observed that the calculation results are consistent with the experimental results. During the heating process, the calculation results of positions B 1 , C 1 , and D 1 are in good agreement with the test results, while the calculation results of position A 1 are slightly lower than the test results, but the trend of change is still consistent. It can be 6 Advances in Materials Science and Engineering observed that when the parameter α does not change during the heating and cooling process (α 1 � α 2 � 1.89 m −1 ) by comparing Figures 4(a) and 4(b), the calculation results present that the content of volumetric water at each location after the heating and cooling process has little change from the initial state and is far from the test result (Figure 4(a)). When the parameter α is set differently during the heating and cooling process (e.g., α 1 � 1.89 m −1 and α 2 � 1.93 m −1 ), the calculation results show that the volumetric water content at each location after the heating and cooling process is lower than the initial state and closer to the test results (Figure 4(b)). Figure 5 shows the comparison between the test results and the calculated results of the temperature at different monitoring positions in the soil for cases 1 and 2. In the calculation process, the parameter α of case 1 is set to α 1 � 1.89 m −1 and α 2 � 1.93 m −1 , and the parameter α of case 2 is set to α 1 � 1.89 m −1 and α 2 � 1.94 m −1 . By comparison, it can be seen that the numerical calculation results of the temperature changes at each measuring point during the heating process of the two cases are in good agreement with the test results, but there is a certain gap during the cooling process. is is due to the errors in the parameters and boundary conditions of the laboratory experiments and numerical calculations. Collector Figure 1: Schematic diagram of the test device.   Advances in Materials Science and Engineering Figure 6 shows the comparison between the test results and the calculated results of the volumetric water content at different monitoring positions in the soil for case 2. rough comparison, it can be observed that the calculation results are consistent with the experimental results. In the calculation process, the parameter α is set to α 1 � 1.89 m −1 and α 2 � 1.94 m −1 , and the initial pore water pressure and pore gas pressure are set to 1.56 × 105 Pa and 1.46 × 105 Pa, respectively. e permeability K in saturated soil is set to 4 × 10 −19 m 2 . It can be observed from Figure 6 that the numerical calculation results in the heating process are in good agreement with the test results, while the calculation results of the measuring points B 1 , C 1 , and D 1 in the cooling process have a certain deviation compared with the test results. Figure 7 shows the distribution of saturation degree in the soil for case 2 when t is equal to 1 d (heating process) and when t is equal to 5.5 d (cooling process). It can be seen from Figure 7 that the bottom saturation degree of the soil, due to gravity effect, is higher than the top saturation degree; that is, the pore water pressure at the bottom is greater. e saturation degree near the inner heat source is almost lower than that of other locations, indicating that pore water migrates from the inside to the outside due to the thermal driving  effect. During the heating process, the saturation degree of the soil showed a cone-shaped distribution with a narrow top and a wide bottom, while during the cooling process, the saturation degree of the soil showed a funnel-shaped distribution with a wide top and a narrow bottom. Besides, and the direction of the pore water pressure gradient during the cooling process is from the bottom close to the boundary to the top close to the heat source.

Conclusion
For unsaturated porous media, the governing equations of the three-field coupling effect of thermo-hydromechanical were established based on the theory of thermodynamics and thermoelasticity. e corresponding state equations were given, and each state variable was adopted according to the temperature-dependent nonlinear expressions.
A self-developed test device was utilized to discuss the hydrothermal coupling effects of two cases of unsaturated kaolin under the same temperature load. e test results were compared with the numerical calculation results. e calculation results at different positions are in good agreement with the test results, which shows that the established theoretical model is able to simulate the water-heat coupling process of unsaturated soils. Different initial saturation degrees have little effect on the change of temperature over time. Compared with the process of moisture migration, the temperature reaches a steady state faster. e characteristics of moisture migration during the temperature rise and fall process are more complicated. When the initial saturation degree of the soil is high, the volumetric water content of the measuring point during the heating process increases, and then it is in a stable state; when the initial saturation is low, the volumetric water content near the heat source increases firstly and then decreases. e volumetric water content of the soil undergoes an irreversible process during the temperature rising and falling process. e test results indicate that the volumetric water content of the measuring point after the temperature rising and falling process is lower than the initial state.  Data Availability e data are available from the corresponding author upon reasonable request.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this manuscript.

Authors' Contributions
Wei Wang conceptualized the study and provided the methodology. Jingjing Liu and Jing Chen investigated the study.