Abstract
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 COMSOL Multiphysics software, a thermo-hydromechanical coupling model 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. The 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. The 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. The melting of the frozen layer is bidirectional, so the duration of slope melting is shorter than that of the freezing process. The 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.
1. Introduction
Seasonally frozen soil is a special soil-water system wherein ice and water coexist. The 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].
The 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. Throughout 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 thermo-hydromechanical (THM) coupling have been performed. In particular, with the development of computer technology, many studies [12–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. The 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.
2. 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.
2.1. 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 is the volume content of unfrozen water in the soil; is the volume content of ice in the soil; 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; is the density of ice (kg/m3); is the density of water (kg/m3); is the hydraulic conductivity of unsaturated soil (m/s), and is the diffusivity of seasonally frozen soil (m2/s).
2.2. Temperature Field Governing Equation
The heat conduction equation for the latent heat of the phase transition as the internal heat source is expressed as follows [17]:where is the temperature of the soil (°C), C is the volume heat capacity of the soil (J/m3/°C), is the soil thermal conductivity (W/m/°C), and is the phase-change latent heat (J/kg).
2.3. Stress Field Governing Equation
The equilibrium equation of isotropic linear thermoelastic materials can be expressed as follows:
The Cauchy equation is expressed as
The stress continuity equation can be expressed as follows [18]:where is the shear modulus (Pa), is the Lamé constant (Pa, ), µ is Poisson’s ratio, is the volume strain, and is the comprehensive coefficient of body thermal expansion and frost heave (°C−1). When the soil temperature is lower than the freezing temperature (), is the frost heaving coefficient of the soil (°C−1). When the soil temperature is higher than the freezing temperature (), is the thermal expansion coefficient of the soil (°C−1). is the Kronecker delta notation, is the displacement component, and is the volume force component for the plane problem i = (x, z), and 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
2.4. Phase Transition Dynamic Equilibrium Relation
The content of unfrozen water in the pores of the soil can be expressed aswhere and represent the volume of water and ice, respectively.
The relation between the content of unfrozen water and the negative temperature is always in a dynamic equilibrium [19], which can be expressed aswhere is the empirical constant related to the properties of the soil and can be selected according to the type of soil: sand (0.61), silt (0.47), and clay (0.56).
The water in soil is composed of pore ice and pore water, the volume content of unfrozen water , the volume content of ice , and the total-volume water content , in which is the initial water content of the soil slope, taking .
2.5. 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 coefficient I is introduced to describe the impedance of ice to the formation of water flow. The hydraulic conductivity and water diffusivity of frozen soils are reduced to 1/I of the thawed soil with the same water content of unfrozen soil.
The impedance coefficient is mainly determined by experience, Taylor and Luthin [5] and Jame and Norm [20] suggested that the impedance coefficient 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 coefficient is very arbitrary. Chen et al. [21] suggested that the impedance coefficient can be expressed aswhere is the unsaturated hydraulic conductivity of thawed soil under the condition of saturated unfrozen water content (), is the saturated hydraulic conductivity of thawed soil, and is the saturated water content of the soil.
The saturated unfrozen water content in the frozen soil is
Equation (10) can effectively avoid the arbitrariness of the value of the impedance coefficient.
2.6. Heaviside Step Function
The two-dimensional Heaviside step function is introduced to characterise the ice-water transition process during freezing. The expression of the Heaviside step function is as follows:where is the slope temperature (°C); is the temperature of the phase transition point (°C), taking = −0.3; and dT is the transition gap (°C), taking dT = 0.3.
The Heaviside step function is shown in Figure 1.
3. Computational Model
3.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. The geometry of the slope is shown in Figure 2. The diffusivity and hydraulic conductivity of the soil are generally measured via experiments. Herein, the diffusivity and hydraulic conductivity of unsaturated soils in the unfrozen zone are obtained by fitting the experimental data as follows:
The diffusivity and hydraulic conductivity of unsaturated soil in the frozen area are divided by the impedance coefficient 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 affecting 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 influence 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. The effect 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 defined as constants are adopted. The rest of the main calculation parameters are shown in Table 1; the subscripts f and u represent frozen soil and unfrozen soil, respectively.
3.2. Multiphysical Field Coupling Technique and Boundary Conditions
COMSOL Multiphysics is a finite-element numerical simulation software that can couple multiphysical fields and solve nonlinear differential equations. Herein, on the basis of the module of the coefficient partial differential equation in the COMSOL software, the THM control equation of frozen soil is transformed into a unified general differential equation group provided by COMSOL. Then, the temperature, humidity, and displacement fields are obtained by applying certain boundary and initial conditions.
The 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 is the air temperature (°C), is the annual average air temperature (°C), is the half of the difference in the annual air temperature (°C), 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:
The 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. The 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. The displacement boundary of AF and DE is a horizontal constraint, and the EF edge is a vertical constraint, whereas the other boundary is free.
3.3. Verification of Solution
The calculated values are compared with the experimental data to verify the calculation model. The test data are obtained from a temperature-monitoring sensor placed in the subgrade slope of the the Monghua Railway test section. The temperature and humidity probe adopts the method of excavation and layered embedding, as shown in Figure 3. The 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.
(a)
(b)
4. Simulation Results and Analysis
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 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). The 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. The 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. The slope is concave at the slope foot and is least affected by ambient temperature; thus, the maximum frozen depth is smallest.
(a)
(b)
(c)
4.1.2. Distribution of the Temperature Field in the Slope at Different Times
The 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). The 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. The frozen layer is simultaneously melting at the maximum frozen depth affected 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 Figures 7(e) and 7(f).
(a)
(b)
(c)
(d)
(e)
(f)
The 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. The 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. The 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. The slope surface temperature is basically consistent with the ambient temperature. In contrast, the deep soil temperature is mainly affected by the geothermal temperature, and the influence depth of the ambient temperature is about 6 m from the slope surface.
The temperature of point J on the slope surface varies periodically and is completely consistent with the ambient temperature (as shown in Figure 8). The 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 interior points of the slope lags behind the ambient temperature. The lowest temperature of point J is first reached at 380 d, and the lowest temperature of points K, L, and M below the surface appears for the first 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. Thus, the number of freeze-thaw 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 influence of the ambient temperature on the slope is gradually weakened. The 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.
4.2. Slope Moisture Field Analysis
The isoline of the slope moisture field is similar to that of the temperature field, as shown in Figure 9. The isoline is denser, and the humidity changes considerably in the shallow area of the slope. The 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. The 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.
(a)
(b)
(c)
(d)
(e)
(f)
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. Thus, 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.
4.3. Slope Displacement Analysis
Figure 11 shows the variation and distribution of slope surface displacement at different 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. The 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. The displacement caused by freezing and thawing is obviously larger than that caused by thermal expansion and contraction. Therefore, in the seasonal freezing zone, the effect of stress caused by frost heaving and thaw settlement on the slope is stronger than that caused by expansion with heat and contraction with cold.
(a)
(b)
5. Conclusions
(1)The 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. The freezing front moves inward continuously, and the freezing occurs only in the shallow range of the slope soil. The 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. The 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 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)The variation in slope displacement is closely related to the change in the ambient temperature. The 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. The effect of frost heaving and thaw settlement is obviously greater than that of thermal expansion and contraction.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This research is supported by the Youth Innovation Promotion Association CAS (2015270), Outstanding Youth Fund of Hubei Province (2017CFA056), Jilin Transportation Science and Technology Project (2016-1-8), Natural Science Foundation of China (Nos. 41472286, 41472290, and 41672312), and Science and Technology Service Network Initiative (KFJ-STS-ZDTP-037). These financial supports are gratefully acknowledged.