A Coupled Thermo-Hydromechanical Model of Soil Slope in Seasonally Frozen Regions under Freeze-Thaw Action

Soil slope diseases in seasonally frozen regions are mostly related to water migration and frost heave deformation of the soil. Based on the partial differential equation defined using the COMSOLMultiphysics software, a thermo-hydromechanical couplingmodel considering water migration, ice-water phase change, ice impedance, and frost heave is constructed, and the variations in the temperature field, migration of liquid water, accumulation of solid ice, and deformation of frost heave in frozen soil slopes are analysed. .e results show that the ambient temperature has a significant effect on the temperature and moisture field of the slope in the shallow area. In addition, the degree of influence gradually weakens from the outside to the inside of the slope, and the number of freeze-thaw cycles in deep soil is less than that in shallow soil. During the freezing period, water in the unfrozen area rapidly migrates to the frozen area, and the total moisture content abruptly changes at the vicinity of the freezing front. .e maximum frozen depth is the largest at the slope top and the smallest at the slope foot. During the melting period, water is enriched at the melting front with the frozen layer melting; the slope is prone to shallow instability at this stage..emelting of the frozen layer is bidirectional, so the duration of slope melting is shorter than that of the freezing process. .e slope displacement is closely related to the change in temperature—a relation that is in agreement with the phenomenon of thermal expansion and contraction in unfrozen areas and reflects the phenomenon of frost heave and thaw settlement in frozen areas.


Introduction
Seasonally frozen soil is a special soil-water system wherein ice and water coexist.e seasonally frozen soil region in China accounts for 53.5% of the total area of China [1].Given the development of China's western region, the implementation of strategy for revitalizing the northeast, and the development strategy known as the Belt and Road Initiative, a large amount of infrastructure-related construction will be required in the seasonally frozen soil regions, e.g., the Sino-Russian oil pipeline project in Northeast Asia, the Harbin-Dalian high-speed railway project in China, and China's National Highway 301.More than two-thirds traffic trunk lines in China are located in seasonally frozen soil regions, and freeze-thaw slope diseases along the lines are very serious.Most of these diseases are related to water migration and frost heave deformation of the soil [2,3].e coupling of water, heat, and force is a key problem in the study of frozen soil in seasonally frozen regions and is also at the frontier of international research in this field.roughout the years, many studies have proposed various frozen soil models.Harlan [4] first established the water-heat coupling model based on the principles of unfrozen water dynamics and energy conservation, and the equations in this model had a clear physical meaning.Taylor and Luthin [5] calculated the water field and temperature field using the implicit finite-difference method based on the Harlan model and compared these calculations with the experimental data.O'Nell and Miller [6] studied the formation and growth of ice lenses in frozen soil and proposed and developed a rigid ice model.An et al. [7] analysed the coupling of water and heat during soil freezing.Shen and Ladanyi [8] addressed the coupling problem of water, heat, and force in three fields by proposing a simplified coupling model based on the Harlan hydrodynamic model.Lai et al. [9] first derived the mathematical mechanical model and solution for the coupled problem of temperature, seepage, and stress fields based on phase changes predicted by heat transfer theory, seepage theory, and frozen soil mechanics.Li et al. [10] established a heat, moisture, and deformation coupling model based on the theory of equilibrium, continuity, and energy principles of a multiphase porous medium.Gens et al. [11] established a new mechanical model that encompasses frozen and unfrozen behaviours within a unified framework based on effective stress.In recent years, numerous studies on thermohydromechanical (THM) coupling have been performed.In particular, with the development of computer technology, many studies [12][13][14][15][16] have addressed the multi-field coupling problem using the finite-element method.
In this study, a THM coupling model is constructed considering water migration, ice-water phase change, ice impedance, and frost heave.
e coupling simulation is realized using the finite-element analysis software COMSOL Multiphysics, and the variations in the slope temperature field, migration of liquid water, accumulation of solid ice, and deformation of frost heave are analysed under freezing and thawing environments.

THM Coupling Theory
Based on the main physical processes in seasonally frozen soil, several hypotheses have been proposed: the freeze-thaw soil medium of the slope is incompressible, homogeneous, and isotropic; there is only movement of liquid water in frozen soil, while the ice remains stationary; the effect of water vapour migration on unfrozen water and heat flow migration is ignored; the effect of water migration caused by temperature gradients and convection heat transfer is ignored; the unfrozen water content in frozen soil is in dynamic equilibrium with the negative temperature of the soil.Based on these assumptions, the THM-coupled model in seasonally unsaturated frozen soils is established.

Water Field Governing Equation.
For plane problems, the law of water migration in unsaturated frozen soils can be expressed by Richards' equation with a phase transition [17]: where θ l is the volume content of unfrozen water in the soil; θ i is the volume content of ice in the soil; t is the time coordinate (s); x, z are the space coordinates (m); x is the horizontal direction, and z is the vertical direction, taking z as positive upward; ρ i is the density of ice (kg/m 3 ); ρ l is the density of water (kg/m 3 ); K(θ l ) is the hydraulic conductivity of unsaturated soil (m/s), and D(θ l ) is the diffusivity of seasonally frozen soil (m 2 /s).

Temperature Field Governing Equation.
e heat conduction equation for the latent heat of the phase transition as the internal heat source is expressed as follows [17]: where T is the temperature of the soil ( °C), C is the volume heat capacity of the soil (J/m 3 / °C), λ is the soil thermal conductivity (W/m/ °C), and L is the phase-change latent heat (J/kg).

Stress Field Governing Equation.
e equilibrium equation of isotropic linear thermoelastic materials can be expressed as follows: e Cauchy equation is expressed as e stress continuity equation can be expressed as follows [18]: where G is the shear modulus (Pa), λ m is the Lamé constant (Pa, λ m � 2Gμ/(1 − 2μ)), µ is Poisson's ratio, ε V is the volume strain, and α T is the comprehensive coefficient of body thermal expansion and frost heave ( °C−1 ).When the soil temperature is lower than the freezing temperature (α T � αθ i ), α is the frost heaving coefficient of the soil ( °C−1 ).When the soil temperature is higher than the freezing temperature (α T � α s ), α s is the thermal expansion coefficient of the soil ( °C−1 ).δ ij is the Kronecker delta notation, u i is the displacement component, and f i is the volume force component for the plane problem i � (x, z), and K � (2G(1 + μ))/(3(1 − 2μ)) is the bulk elastic modulus of the soil (Pa).By using the above static equilibrium equation and Cauchy equation, the THM-coupled stress-control equation can be expressed as

Phase Transition Dynamic Equilibrium
Relation.e content of unfrozen water in the pores of the soil can be expressed as where V l and V i represent the volume of water and ice, respectively.e relation between the content of unfrozen water and the negative temperature is always in a dynamic equilibrium [19], which can be expressed as where b is the empirical constant related to the properties of the soil and b can be selected according to the type of soil: sand (0.61), silt (0.47), and clay (0.56).

Advances in Civil Engineering
e water in soil is composed of pore ice and pore water, the volume content of unfrozen water θ l θ 0 χ(T), the volume content of ice θ i (θ 0 − θ l ) • (ρ l /ρ i ), and the totalvolume water content θ θ l + θ i • (ρ i /ρ l ), in which θ 0 is the initial water content of the soil slope, taking θ 0 25%.

Impedance of Ice to the Formation of Water Flow in Frozen
Soils.For frozen zones, because of the ice in the soil slope, the impedance coe cient I is introduced to describe the impedance of ice to the formation of water ow.e hydraulic conductivity and water di usivity of frozen soils are reduced to 1/I of the thawed soil with the same water content of unfrozen soil.
e impedance coe cient is mainly determined by experience, Taylor and Luthin [5] and Jame and Norm [20] suggested that the impedance coe cient depends on the volume content of ice in the frozen soil and can be expressed as follows: Using this treatment method, the value of the impedance coe cient is very arbitrary.Chen et al. [21] suggested that the impedance coe cient can be expressed as where K us is the unsaturated hydraulic conductivity of thawed soil under the condition of saturated unfrozen water content (θ us ), K s is the saturated hydraulic conductivity of thawed soil, and θ s is the saturated water content of the soil.e saturated unfrozen water content in the frozen soil is Equation ( 10) can e ectively avoid the arbitrariness of the value of the impedance coe cient.

Heaviside Step Function.
e two-dimensional Heaviside step function is introduced to characterise the ice-water transition process during freezing.
e expression of the Heaviside step function is as follows: where T is the slope temperature ( °C); T trans is the temperature of the phase transition point ( °C), taking T trans = −0.3;and dT is the transition gap ( °C), taking dT = 0.3.e Heaviside step function is shown in Figure 1.

Physical Parameter Value.
Taking a typical soil slope of the Monghua Railway in a seasonally frozen region as an example, the height of the slope is 15 m and the angle of the slope is 30 °; the slope lithology is mainly silty clay.e geometry of the slope is shown in Figure 2. e di usivity and hydraulic conductivity of the soil are generally measured via experiments.Herein, the di usivity and hydraulic conductivity of unsaturated soils in the unfrozen zone are obtained by tting the experimental data as follows: e di usivity and hydraulic conductivity of unsaturated soil in the frozen area are divided by the impedance coe cient I.When the freezing state of the soil changes, the thermal conductivity and volume heat capacity of the soil also change.Herein, the thermal conductivity and volume heat capacity of the frozen soil and the unfrozen soil are determined.Under certain conditions of soil texture and porosity, the water content is an important factor a ecting the thermal conductivity and volume heat capacity of the slope soil.When the water content of the unfrozen soil exceeds the liquid limit [22], the in uence of water content on the thermal conductivity and volume heat capacity of the soil is weakened.When the soil is frozen, water molecules are not as active as in the melting state.e e ect of water content on frozen soil is not as obvious as that on unfrozen soil; therefore, the thermal conductivity and volume heat capacity of the high water content slope de ned as constants are adopted.e rest of the main calculation parameters are shown in Table 1; the subscripts f and u represent frozen soil and unfrozen soil, respectively.

Multiphysical Field Coupling Technique and Boundary
Conditions.COMSOL Multiphysics is a nite-element numerical simulation software that can couple multiphysical elds and solve nonlinear di erential equations.Herein, on the basis of the module of the coe cient partial di erential equation in the COMSOL software, the THM control equation of frozen soil is transformed into a uni ed general Advances in Civil Engineering differential equation group provided by COMSOL.en, the temperature, humidity, and displacement fields are obtained by applying certain boundary and initial conditions.e surface layer of a slope is mainly affected by the external ambient temperature; a periodic change in the ambient temperature causes a change in the slope surface temperature, which can be estimated using a sinusoidal curve [23]: where T is the air temperature ( °C), T 0 is the annual average air temperature ( °C), A is the half of the difference in the annual air temperature ( °C), t is the time (d), and ϕ is the initial phase.For example, the variation in the average annual air temperature in Harbin, China, from 1971 to 2000, can be expressed as follows: e sinusoidal temperature function is applied to the upper surface of the slope, i.e., the temperature load is applied simultaneously on the three sides of AB, BC, and CD in Figure 2. e annual geothermal temperature in this region is about 8 °C, so the initial temperature of the slope is set to 8 °C.Because rainfall and evaporation are not considered, there is no infiltration or exudation at AB, BC, CD, DE, EF, and AF boundaries in Figure 2. e displacement boundary of AF and DE is a horizontal constraint, and the EF edge is a vertical constraint, whereas the other boundary is free.

Verification of Solution.
e calculated values are compared with the experimental data to verify the calculation model.
e test data are obtained from a temperaturemonitoring sensor placed in the subgrade slope of the the Monghua Railway test section.e temperature and humidity probe adopts the method of excavation and layered embedding, as shown in Figure 3. e calculated and measured soil temperature and moisture content are compared in Figure 4. Figure 4 shows that the temperature calculated by the model is consistent with the magnitude and trend of the measured values, whereas the moisture content is only different in the surface layer mainly because the calculation model does not consider rainfall and evaporation on the surface.Figure 5 shows the temperature distribution at the top, waist, and foot of the slope for a vertical range of 15 m from the slope surface for about 2 years (a temperature curve was drawn every day).e maximum frozen depth is 2.80 m at the slope top, 2.14 m at the slope waist, and 1.80 m at the slope foot.e slope is convex at the slope top and is affected by the bidirectional freezing of the slope surface and slope top; thus, the maximum frozen depth of the slope top is the largest.e slope is concave at the slope foot and is least affected by ambient temperature; thus, the maximum frozen depth is smallest.

Distribution of the Temperature Field in the Slope at Different Times.
e isoline of the slope temperature field is shown in Figure 6; the slope temperature fluctuates periodically with the change of seasons.In the shallow area of the slope, the isoline is denser and the temperature gradient is larger, which is strongly influenced by the ambient temperature.When the ambient temperature is higher than the slope surface temperature, heat is transferred from the external higher temperature to the internal lower temperature, as shown in Figure 6(a).When the ambient temperature is lower than the slope surface temperature, the temperature of the slope surface initially begins to decrease and heat is transferred from the inside to the outside of the slope, as shown in Figure 6(b).When the ambient temperature is below freezing point, the slope surface begins to freeze.Over time, the freezing front keeps moving inward, but the freezing occurs only in the shallow range of the slope soil, as shown in Figures 6(c) and 6(e).e total duration of the freezing process of the slope is about 170 d.When the ambient temperature is higher than the freezing point and enters the stage of spring melting, the temperature of the slope surface begins to rise first; then, the frozen zone of the surface begins to melt, and this formed melting front keeps moving inward.e frozen layer is simultaneously melting at the maximum frozen    Advances in Civil Engineering e temperature distribution of the slope waist along vertical distance from slope surface is shown in Figure 7 at t 193, 290, 375, 430, 459, 480, and 560 d. e slope is in an unfrozen state when t 193 and 290 d.It heats up from the initial state to t 193 d; at this time, the ambient and slope temperatures are at their highest.When the ambient temperature changes from high to low and then drops to its freezing point at t 290 d, the slope surface will be frozen, and the slope will enter the freezing stage with a further decrease in the ambient temperature.
e ambient temperature reached its lowest point at t 375 d, as did the slope temperature.After that time, the ambient temperature begins to rise until t 459 d; the freezing period ends, and the slope begins to enter the stage of spring thaw.When the ambient temperature reaches the next-highest temperature at t 560 d, the slope experiences a complete freeze-thaw process.e slope surface temperature is basically consistent with the ambient temperature.In contrast, the deep soil temperature is mainly a ected by the geothermal temperature, and the in uence depth of the ambient temperature is about 6 m from the slope surface.
e temperature of point J on the slope surface varies periodically and is completely consistent with the ambient temperature (as shown in Figure 8).e response temperature also varies periodically with the ambient temperature at points K, L, and M, whose vertical distances from the slope surface are 0.5, 1.0, and 3.0 m, respectively.However, the time of the highest and lowest temperatures at the    6 Advances in Civil Engineering interior points of the slope lags behind the ambient temperature.e lowest temperature of point J is rst reached at 380 d, and the lowest temperature of points K, L, and M below the surface appears for the rst time at 390, 410, and 540 d, respectively.It means that the shallow soil has been frozen and thawed for some time, and the deep soil is just beginning to freeze and thaw.us, the number of freezethaw cycles in deep soil is less than that in shallow soil, and the degree of freeze-thaw damage will also decrease.With an increase in the depth of the slope, the in uence of the ambient temperature on the slope is gradually weakened.e peak temperature range of points K and L changes from 23.0 °C to −18.3 °C of ambient temperature from 19.1 °C to −11.4 °C and 16.1 °C to −6.3 °C.Point M has a peak temperature of 10.9 °C -1.9 °C, which is already in the nonfreeze-thaw area and is close to the initial temperature of the slope for most of the time.

Slope Moisture Field Analysis.
e isoline of the slope moisture eld is similar to that of the temperature eld, as shown in Figure 9. e isoline is denser, and the humidity changes considerably in the shallow area of the slope.e distribution of the total water content of the slope waist along vertical distance from slope surface at t 193, 290, 375, 430, 459, 480, and 560 d is also studied, as shown in Figure 10.In the initial unfrozen state, the total water content of each region of the slope is equal to the initial water content.When the ambient temperature is below the freezing point, the slope begins to freeze from the surface, the unfrozen water content in the frozen area decreases, and the suction increases.Under the action of suction, the water in the unfrozen area of the slope moves rapidly toward the frozen zone.Part of the migrating water freezes near the freezing front; the ice content increases, and the frozen front moves inward.Part of the migrating water is blocked by ice.Water migration leads to an increase in the total moisture content in the frozen area and a decrease in the moisture content in the unfrozen area.e total water content abruptly changes at the vicinity of the freezing front.When the ambient temperature is negative for a long time, the shallow layer of the slope is frozen and stable.Because of the existence of ice, water migration in the stable frozen area is almost stagnant, and the total water content in the stable frozen zone remains almost unchanged.
However, with the melting of the surface layer of the slope, water that accumulated in the surface layer during the freezing period is enriched at the melting front.us, there is a peak value of the water content of the slope at the melting front.Because the frozen layer still exists in the lower part of the thawing layer for a relatively long time, it is very likely that the slope will lose its stability at the freeze-thaw interface during this stage.

Slope Displacement Analysis.
Figure 11 shows the variation and distribution of slope surface displacement at di erent times.In the unfrozen area, the slope displacement is in agreement with the phenomena of expansion with heat and contraction with cold (the curves of t 193 and 290 d in Figure 11).In contrast, in the frozen area, the displacement of the slope is the result of frost heaving and thaw settlement (the curves of t 375, 459, and 480 d in Figure 11).When the temperature drops below the freezing point, the unfrozen water phase becomes ice and the volume expands.e vertical and horizontal displacements of the slope is large at the slope top and small at the slope foot, and the slope is inclined outward, so frost heaving and tensioning is more likely to occur at the slope top in winter.When the temperature increases from negative to positive, the ice in the frozen soil melts and the slope experiences thaw settlement.Simultaneously, because of the existence of the freeze-thaw interface, the slope is prone to shallow sliding failure at the interface.e displacement caused by freezing and thawing is obviously larger than that caused by thermal expansion and contraction.erefore, in the seasonal freezing zone, the e ect of stress caused by frost heaving and thaw settlement on the slope is stronger

Advances in Civil Engineering
than that caused by expansion with heat and contraction with cold.

Conclusions
(1) e slope temperature fluctuates periodically with the change in seasonal ambient temperature, and the influence of ambient temperature weakens gradually from the outside to the inside of the slope.During the freezing period, the slope gradually freezes from outside to inside.e freezing front moves inward continuously, and the freezing occurs only in the shallow range of the slope soil.e maximum frozen depth of the slope is the largest at the slope top and the smallest at the slope foot.During the melting period, the melting of the frozen layer is bidirectional.Compared with the freezing process of the slope, the melting process is shorter and faster.e number of freeze-thaw cycles in deep soil is less than that in shallow soil.
(2) In the shallow area of the slope, the humidity is considerably affected by the ambient temperature.
During the period of freezing, the water content of the unfrozen zone moves rapidly toward the frozen zone and the total moisture content abruptly changes at the vicinity of the freezing front.In the spring thaw Advances in Civil Engineering period, with the shallow layer melting, the water is enriched at the melting front and the slope is prone to shallow instability at the freeze-thaw interface.(3) e variation in slope displacement is closely related to the change in the ambient temperature.e displacement corresponds to the phenomena of expansion with heat and contraction with cold in the unfrozen area and presents the phenomena of frost heaving and thaw settlement in the frozen area.e effect of frost heaving and thaw settlement is obviously greater than that of thermal expansion and contraction.

4. 1 .
Analysis of Temperature Field on Slope 4.1.1.Maximum Frozen Depth.Based on the boundary conditions, the annual minimum temperature of the external environment is −18.3 °C.We analysed the maximum frozen depth of the slope at three key positions: top, waist, and foot.

Figure 3 :
Figure 3: Temperature and humidity sensor embedded in the test section.
depth a ected by the geothermal temperature, so the melting of the frozen layer is bidirectional.When the melting process lasts 20 d, there is almost no frozen area in the slope.It can be seen that the melting time is shorter and faster compared with the freezing time, as shown in Figures7(e) and 7(f ).

Figure 4 :Figure 5 :
Figure 4: Comparison of calculated values with test values.(a) Temperature distribution with depth; (b) water content distribution with depth.

tFigure 7 :Figure 8 :
Figure 7: Temperature distribution of the slope waist along vertical distance from slope surface at di erent times.