Influencing Factors of Soil Mechanical Properties in Seasonally Frozen Regions and Thermo-Hydro-Mechanical Coupled Simulation of Embankment Stability

1. Introduction

River embankments serve as critical safeguards for flood protection, and their instability is governed by complex interactions among hydraulic processes, soil properties, and thermal conditions. In high-latitude and high-altitude regions, freeze-thaw cycles play a critical role in controlling riverbank stability. The stability assessment of river embankments becomes particularly challenging under the coupled effects of freeze-thaw processes and seepage. Such coupling effects can significantly alter soil physical properties and soil structure, thereby increasing the risk of embankment failure.

Freeze-thaw processes induce changes in soil structure and pore characteristics, thereby altering its physical and mechanical properties [1,2]. Aoyama and Ogata reported that freeze-thaw cycling caused mechanical degradation in various soils, characterized by reduced cohesion, relatively unchanged internal friction angle, and decreased shear strength [3,4]. Laboratory freeze-thaw tests further demonstrated that cohesion and internal friction angle experienced the most significant changes after the first freeze-thaw cycle and gradually approached a stable state after approximately 8–9 cycles [5,6]. The degradation characteristics of soil mechanical properties during freeze-thaw processes depend strongly on freezing temperature, cycle number, soil type, and water content. For gravelly soils, shear strength decreases markedly with increasing water content, whereas cohesive soils exhibit an opposite trend [7]. Yu et al. developed a low-temperature freeze-thaw triaxial testing apparatus and systematically investigated silty soils from the Qinghai–Tibet Plateau. Their study evaluated the effects of confining pressure, freeze-thaw cycles, and freezing temperature on soil physical and mechanical properties, and proposed a strength evolution model during freeze-thaw processes [8,9]. During freezing, water migrates toward the surface, increasing the water content of near-surface soil layers toward saturation. During thawing, melting of surface frozen soil generates pore water, while the underlying frozen layer restricts drainage and increases pore water pressure. This process may induce gradual deformation at the frozen–unfrozen soil interface and eventually trigger solifluction (freeze-thaw-induced soil flow) and slope instability [10,11,12]. Failure modes of frozen slopes can generally be classified into landslides, collapses, and solifluction [13]. The limit equilibrium method (LEM), finite element method (FEM), and strength reduction method (SRM) have been widely applied in frozen slope stability analyses [14,15]. Leshchinsky et al. demonstrated the feasibility of finite element analysis for three-dimensional slope stability simulations by comparing rigorous three-dimensional limit equilibrium analysis with finite element methods [16]. Xu et al. incorporated temperature-dependent shear strength parameters into a strength reduction finite element framework to investigate the influence of temperature on slope stability [17]. With advances in multiphysics coupling theories, temperature, moisture, and stress fields have been increasingly incorporated into freeze-thaw slope stability assessments. Previous studies have demonstrated that the combined effects of freeze-thaw cycles and seepage can modify effective stress distributions, accelerate shear strength degradation, and reduce slope stability reserves, highlighting the necessity of considering coupled thermo-hydro-mechanical processes in freeze-thaw-induced instability assessments [18].

River ice processes can significantly influence water-level fluctuations in the Inner Mongolia section of the Yellow River, making riverbank stability assessments under freeze-thaw conditions more complex. In this study, a representative river embankment section was selected for investigation. Field sampling and laboratory tests were conducted to characterize the physical and mechanical properties of soils before and after freeze-thaw cycles. By integrating experimental characterization with a thermo-hydro-mechanical (THM)-coupled finite element analysis, the stress–strain responses and temporal evolution of bank slope displacement were investigated under freeze-thaw conditions. The results provide insights into freeze-thaw-induced riverbank failure mechanisms and support stability assessment of cold-region river embankments.

2. Study Area and Experimental Design

2.1. The Study Area Overview

The Shisifenzi Bend, located approximately 4 km upstream of the Toudaoguai Hydrological Station along the Inner Mongolia section of the Yellow River, was selected as the study area. The geographical location of the study area is shown in Figure 1. The channel slope is approximately 0.15‰, with an average river width of about 300 m and an average water depth of approximately 2.85 m [19,20]. The river flows from northeast to southwest, and the bend exhibits an approximately 180° curvature with a curvature coefficient of 3.15. The channel width at the upstream entrance of the bend is approximately 620 m, whereas the downstream exit section narrows substantially to only 215 m, representing a reduction of nearly two-thirds in river width. The Shisifenzi Bend is a key location for ice-jam-induced river closure along the Inner Mongolia section of the Yellow River. Since 2000, five initial ice-jam events have occurred at this site. Following ice-jam formation, the river water level increases substantially, reaching a maximum elevation of 988.4 m. After river breakup, the minimum water level decreases to 985.5 m, resulting in a water-level difference of 2.9 m during the ice period. The bend exhibits distinct hydraulic characteristics, including deposition along the convex bank and erosion along the concave bank. Therefore, it is considered a critical river reach along the Inner Mongolia section of the Yellow River. Due to ice-period water-level fluctuations and freeze-thaw cycles, instability and collapse frequently occur along the concave bank during the ice-melting period.

2.2. Experimental Design

In accordance with the Standard for Soil Test Methods [21], soil samples were collected from the concave bank of the Shisifenzi Bend. Soil samples were collected vertically at 60 cm intervals, with sampling depths of 0, 60, 120, 180, and 240 cm below the ground surface. Particle size distribution analysis, dry density measurements, and water content measurements were conducted to characterize the basic physical properties of the riverbank soils.

To investigate the effects of different experimental conditions on soil mechanical properties, triaxial shear tests were conducted under four water content levels (12.0%, 15.0%, 18.0%, and 21.0%), three confining pressure levels (50, 100, and 150 kPa), and different numbers of freeze-thaw cycles. Based on the recorded minimum winter temperature of −22.6 °C and maximum temperature of 15 °C in recent years [22], the freezing and thawing temperatures were set at −20 °C and 15 °C, respectively. The detailed experimental program is summarized in Table 1.

Prior to the freeze-thaw triaxial shear tests, soils with four target water contents were sealed in plastic bags and cured for 24 h to ensure uniform moisture distribution. The required soil masses for triaxial and compression specimens with different initial dry densities and water contents were calculated. The specimens were prepared using a three-piece mold and compacted in five layers. Each layer was scarified before the subsequent layer was compacted to prevent density nonuniformity between layers. Remolded cylindrical triaxial specimens with a height of 80 mm and a diameter of 39.1 mm were obtained. The prepared specimens were subsequently placed in a freeze-thaw chamber for freeze-thaw cycling. Freeze-thaw cycles were performed using a walk-in high-low temperature test chamber (DHT(H)R). Triaxial shear tests were conducted using a strain-controlled triaxial testing apparatus (TSZ10-1.0). The unconsolidated undrained (UU) triaxial shear test was adopted with a strain rate of 0.08 mm/min. The test was terminated when the axial strain reached 15%, and the maximum deviator stress was adopted as the failure criterion. The failure principal stress lines under different confining pressures were plotted with σ₁ as the vertical axis and σ₃ as the horizontal axis. The inclination angle of the principal stress line was defined as α, and the intercept on the vertical axis was defined as d. The internal friction angle and cohesion of the soil were determined according to Eqs. (1) and (2), and the shear strength was subsequently calculated using Eqs. (3) and (4).

$$\phi ={\sin }^{-1}\left(\tan \alpha \right)$$

(1)

$$c={\displaystyle \frac{d}{\cos \phi }}$$

(2)

$$\tau =c+\sigma \tan \phi$$

(3)

$$\sigma ={\displaystyle \frac{{\sigma }_1+{\sigma }_3}{2}}+{\displaystyle \frac{{\sigma }_1-{\sigma }_3}{2}}\cos 2\left({45}^{\circ }+{\displaystyle \frac{\phi }{2}}\right)$$

(4)

where φ is the internal friction angle of the soil (°); c is the cohesion of the soil (kPa); τ is the shear strength at failure of the soil (kPa); and σ is the normal stress acting on the failure plane of the soil (kPa).

Figure 2. Schematic diagram of principal stress lines at soil failure.

Figure 2. Schematic diagram of principal stress lines at soil failure.

3. Tests Results and Analysis

3.1. Soil Particle Size Distribution

A BT-9300LD laser particle size analyzer (LPSA) was used to characterize the particle size distribution of soil samples. The particle size distribution curves of the surface, middle, and bottom layers are presented in Figure 3. The riverbank soils exhibit evident vertical variations in particle composition. From the surface layer (0–60 cm) to the middle layer (120–180 cm), all soil samples were classified as silty sand (SM). The content of fine particles (<75 μm) exhibited a fluctuating trend with increasing sampling depth, characterized by a slight increase, followed by a rapid increase and subsequent decrease, with corresponding proportions of 35.60%, 36.60%, 47.64%, and 43.87%. The particle size distribution characteristics of the riverbank soils were generally consistent within the study area, although the proportions of different soil fractions varied among sampling depths. Overall, the soils were characterized by high fine particle contents and relatively low coarse particle fractions. Fine components, mainly silt and clay particles, dominated the soil composition. The absence of abrupt changes in the particle size distribution curves suggests relatively uniform particle size characteristics within the study area.

3.2. Mechanical Properties of Soils

3.2.1. Stress–Strain Behavior of Soils Under Different Water Contents

Figure 4 presents the measured shear stress–strain relationships of soils under different confining pressures with varying initial water contents. The results indicate that initial water content has a clear negative relationship with soil shear strength. At all confining pressure levels, the deviator stress of the soil samples exhibited an overall decreasing trend as the initial water content increased from 12% to 21%. Under low water content conditions, the soils exhibited the highest initial tangent modulus and peak strength, accompanied by relatively small axial strains. Subsequently, a certain degree of strain softening was observed after peak strength. As the water content increased to 21%, the peak strength decreased substantially. The stress–strain curves became flatter and smoother, exhibiting more pronounced strain-hardening behavior and ductile failure characteristics. The ultimate strength was also considerably lower than that of specimens with lower water contents. In addition, although increasing confining pressure substantially increased the deviator stress of specimens under all water content conditions by enhancing lateral confinement, it did not alter the fundamental trend of decreasing soil strength with increasing water content.

3.2.2. Evolution of Soil Cohesion

Figure 5 illustrates the evolution of cohesion in riverbank soils with different initial dry densities under four initial water content conditions with increasing numbers of freeze-thaw cycles. Regardless of the combinations of water content and dry density, soil cohesion exhibited a monotonic degradation trend with increasing freeze-thaw cycles, and the degradation process showed distinct stage-dependent characteristics. During the early stage of freeze-thaw cycling (0–3 cycles), the curves exhibited relatively steep declines, indicating that frost heaving pressure caused rapid structural deterioration of the soil. As the number of freeze-thaw cycles increased from 3 to 10, the curves gradually flattened, and the rate of cohesion loss decreased, indicating that the accumulated damage gradually approached stabilization.

Initial dry density was the primary factor controlling the absolute level of soil cohesion during freeze-thaw cycling. Under all four water content conditions, the five curves were consistently distributed according to initial dry density, with higher-density soils maintaining substantially higher residual cohesion than lower-density soils after 10 freeze-thaw cycles. Meanwhile, water content affected not only the initial cohesion level but also the freeze-thaw degradation characteristics of soils. At a water content of 15%, soils with different dry densities exhibited the highest initial cohesion. However, when water content increased to 21%, the initial cohesion of soils at all dry density levels decreased substantially. After 10 freeze-thaw cycles, the difference in cohesion among soils with different dry densities became smaller, indicating that high water content weakened the strengthening effect of dry density.

The degradation of cohesion during freeze-thaw cycling is primarily associated with repeated frost-induced stresses generated by pore water phase changes, which disrupt particle bonding and interlocking structures. Soils with higher initial dry density possess smaller void ratios and limited freezable water content, while their denser skeletons provide stronger resistance to frost-induced deformation. In contrast, under high water content conditions, increased pore water content amplifies frost-induced stresses and weakens particle bonding, leading to irreversible structural loosening and cohesion loss during freeze-thaw cycles.

3.2.3. Evolution of Soil Internal Friction Angle

Figure 6 presents the variations in the internal friction angle of riverbank soils with freeze-thaw cycles and initial dry density under different initial water contents. Under different water contents, initial dry density consistently dominated the absolute level of internal friction angle. In all panels, the five curves exhibited a clear vertical arrangement according to initial dry density, with denser soils maintaining higher internal friction angles throughout the 10 freeze-thaw cycles. Regarding the effect of freeze-thaw cycling, the internal friction angle generally exhibited a slight decreasing trend with increasing freeze-thaw cycles. However, Figure 6(b) and Figure 6(c) show that the internal friction angle of soils with low initial dry density remained stable or even slightly increased locally during the early freeze-thaw stage (0–3 cycles). This phenomenon may be attributed to the temporary enhancement of particle interlocking in loose soils caused by ice crystal compression during the initial freeze-thaw stage. However, with continued cycling, the degradation rate increased, indicating the irreversible accumulation of freeze-thaw damage. Meanwhile, water content had a considerable influence on the internal friction angle. Within the water content range of 15%–18%, soils with different dry densities exhibited relatively high initial internal friction angles, which was consistent with the existence of an optimum water content (approximately 15%) for the tested soils. When the water content increased to 21%, the overall internal friction angle decreased substantially, and high-density soils (e.g., 1.49 g/cm³) also exhibited greater strength losses under freeze-thaw conditions.

3.3. Sensitivity Analysis of Factors Influencing Soil Shear Strength

To quantitatively determine the relative importance of initial dry density, water content, and freeze-thaw cycles on soil shear strength, analysis of variance (ANOVA) was employed to evaluate the statistical significance of these key influencing factors. The effects of initial dry density, water content, and freeze-thaw cycle number on soil shear strength parameters were systematically examined, providing statistical evidence for understanding the mechanical evolution of riverbank soils under multi-field coupling conditions.

Considering the potential interactions among different factors, a statistical model containing second-order interaction terms was established:

$$Y=\mu +A_w+B_{\rho }+C_N+A_w{B}_{\rho }+A_w{C}_N+B_p{C}_N+\epsilon$$

(5)

where Y represents the soil shear strength parameters, including cohesion and internal friction angle under different normal stress conditions. A_w, $B_{\rho }$ , and C_N represent the effects of water content, dry density, and freeze-thaw cycle number, respectively. A_w$B_{\rho }$, A_wC_N, and B_ρC_N denote the interaction effects between water content and dry density, water content and freeze-thaw cycles, and dry density and freeze-thaw cycles, respectively. ε represents the random error.

Four water contents, five dry densities, and six freeze-thaw cycle levels were considered in this experiment, resulting in a total of 120 test combinations. The effects of different factors on soil shear strength were evaluated using the sum of squares (SS), mean square (MS), F-statistic, and significance level (p-value). In addition, the contribution ratio was used to quantify the contribution of each factor to the total variation:

$${\eta }_i={\displaystyle \frac{S{S}_i}{S{S}_T}}\times 100\%$$

(6)

where η_i represents the contribution ratio of the ith influencing factor, SS_i is the corresponding sum of squares of the ith factor, and SS_T represents the total sum of squares.

Table 2 presents the ANOVA results for cohesion. The results indicate that dry density was the dominant factor affecting cohesion, with a contribution ratio of 60.199%. Increasing dry density enlarged the contact area among soil particles and enhanced particle bonding, thereby improving soil structural stability. Water content was the second most influential factor, with a contribution ratio of 19.618%. Variations in water content altered the pore conditions between soil particles. During freeze expansion and thaw contraction processes, changes in water content weakened particle interlocking, thereby affecting soil cohesion. Freeze-thaw cycle number also exhibited a considerable influence on cohesion, with a contribution ratio of 15.940%. This effect was mainly attributed to the repeated destruction of internal bonding structures caused by freeze-thaw processes, resulting in cohesion degradation. Furthermore, the three interaction terms accounted for 4.244% of the total contribution, indicating that interactions among water content, dry density, and freeze-thaw cycles also influenced the cohesive structure of the soil.

Table 3 presents the ANOVA results for the internal friction angle $\phi$. Similar to cohesion, the internal friction angle was predominantly controlled by dry density, which accounted for 88.556% of the total contribution, far exceeding the contributions of water content and freeze-thaw cycles. This is because the internal friction angle primarily reflects friction and interlocking among soil particles, whereas changes in dry density directly modify particle arrangement and interparticle contact conditions. The contribution of freeze-thaw cycles to the internal friction angle was relatively small, at only 3.215%, consistent with the experimental observations. This suggests that freeze-thaw action primarily affects the structural integrity of the soil, while its influence on interparticle friction is comparatively limited.

Overall, the three-factor ANOVA results indicate that dry density, water content, and freeze-thaw cycle number ranked in descending order of influence on both cohesion and internal friction angle. Dry density was the dominant factor controlling the initial level of soil shear strength, whereas freeze-thaw cycling played an important role in cohesion degradation by progressively disrupting the internal bonding structure of the soil. The shear behavior of freeze-thaw-affected soils therefore reflects the combined effects of three processes: dry density governs the initial strength level, water content regulates the extent of freeze-thaw-induced deterioration, and repeated freeze-thaw cycles promote progressive structural degradation.

4. Development and Application of a Thermo-Hydro-Mechanical Coupled Model for Riverbank Stability

4.1. Governing Equations

4.1.1. Governing Equation for the Temperature Field

Based on the principle of energy conservation, the net heat inflow into an infinitesimal soil element per unit time is equal to the rate of increase in its internal energy plus the internal heat-source term. Accordingly, the heat-conduction equation for frozen soil, accounting for the latent heat associated with ice–water phase change, can be expressed as:

$$\lambda {\nabla }^{2}T\cdot dxdydzdt+L{\rho }_i{\displaystyle \frac{\partial {\theta }_i}{\partial t}}dt\cdot dxdydz=\rho C\left(\theta \right){\displaystyle \frac{\partial T}{\partial t}}dt\cdot dxdydz$$

(7)

where $\lambda$ is the effective thermal conductivity of the soil (W·m^-1K^-1); $\rho$ is the dry density of the soil (kg·m^-3); $T$ is the temperature (°C); $t$ is time (s); ${\rho }_i$ is the density of ice (kg·m^-3); $L$ is the latent heat of phase change, taken as $3.34\times {10}^5$J·kg^-1; $\theta$ is the porosity, with$\theta ={\theta }_w+{\theta }_i$,where ${\theta }_w$ and ${\theta }_i$ are the volumetric contents of unfrozen water and ice, respectively (m³·m^-3); and $C\left(\theta \right)$ is the specific heat capacity of the soil (J·kg^-1K^-1).

Equation (7) can be simplified to obtain the governing equation for the temperature field of frozen soil [23]:

$$C_v{\displaystyle \frac{\partial T}{\partial t}}=\lambda \left(\theta \right){\nabla }^{2}T+L{\rho }_i{\displaystyle \frac{\partial {\theta }_i}{\partial t}}$$

(8)

Where, C_v=ρC(θ) is the volumetric heat capacity of the soil (J·m^-3K^-1).

4.1.2. Governing equation for the Moisture Field

According to the law of mass conservation and Darcy’s law for unsaturated soils, the change in the mass of liquid water within an infinitesimal element of frozen soil per unit time is equal to the sum of the net inflow of liquid water and the mass change induced by the ice–water phase transition. The flow of liquid water is driven by the gradient of the total hydraulic head. Accordingly, the governing equation for moisture migration during freeze-thaw processes can be expressed as [24]:

$${\displaystyle \frac{\partial {\theta }_w}{\partial t}}+{\displaystyle \frac{{\rho }_i}{{\rho }_w}}{\displaystyle \frac{\partial \theta i}{{\partial }_t}}=\nabla [D(\theta )\nabla {\theta }_w+K({\theta }_w)]$$

(9)

where ${\rho }_w$ is the density of water (kg·m^-3); $K\left({\theta }_w\right)$ is the unsaturated hydraulic conductivity; and $D\left(\theta \right)$ is the moisture diffusivity. The unsaturated hydraulic conductivity $K\left({\theta }_w\right)$ and moisture diffusivity $D\left(\theta \right)$ of frozen soil can be calculated as follows:

$$D\left({\theta }_w\right)={\displaystyle \frac{K\left({\theta }_w\right)}{c\left({\theta }_w\right)}}\cdot I$$

(10)

$$K\left({\theta }_w\right)=k_s{S}^l{\left[1-{\left(1-S^{1/m}\right)}^m\right]}^2$$

(11)

$$c\left({\theta }_w\right)={\displaystyle \frac{1}{a}}{\left(S^{-1/m}-1\right)}^{1-m}$$

(12)

where $k_s$ is the saturated hydraulic conductivity (m·d^-1); $l$is the pore-connectivity parameter, commonly taken as 0.5; and $m$is an empirical parameter related to the pore-size distribution. Here, $c\left({\theta }_w\right)$ is the specific moisture capacity (m^-1), and $I$ is the impedance factor representing the inhibitory effect of ice content on moisture migration. The impedance factor is calculated as [25]:

$$I={10}^{-10{\theta }_i}$$

(13)

4.1.3. Governing Equation for the Stress Field

Under freeze-thaw conditions, the total strain of the soil consists of the elastic strain induced by external loading and the volumetric strain induced by moisture phase change and migration [26], expressed as:

$$\epsilon ={\epsilon }_{ij}^e+{\epsilon }_{ij}^v$$

(14)

where ${\epsilon }_{ij}^e$ represents the elastic strain induced by external loading, and ${\epsilon }_{ij}^v$ represents the volumetric strain induced by moisture phase change and migration. The strain associated with moisture phase change and migration can be calculated as:

$${\epsilon }_{ij}^v=0.09\left({\theta }_o+{\Delta }_{\theta }-{\theta }_u\right)+\Delta \theta +\left({\theta }_o-n\right)$$

(15)

A linear elastic constitutive relation accounting for the initial strain was adopted and combined with the Mohr–Coulomb strength criterion. The resulting constitutive relations are expressed as:

$$\left\{\sigma \right\}=\left[C\right]\left(\left\{\epsilon \right\}-\left\{{\epsilon }_0\right\}\right)$$

(16)

$$\tau =c+\sigma \tan \phi$$

(17)

where $\left\{\sigma \right\}$ is the stress vector; $\left[C\right]$ is the material stiffness matrix; $\left\{\epsilon \right\}$ is the total strain vector; and $\left\{{\epsilon }_0\right\}$ is the initial strain vector. In Eq. (17), $c$ is the soil cohesion, $\phi$ is the internal friction angle, $\sigma$ is the normal stress, and $\tau$is the shear stress. By combining the above governing equations with the constitutive relations and discretizing the system using the finite element method (FEM), the displacement field, stress field, and distribution of plastic zones within the bank slope can be obtained.

4.1.4. Coupling Theory for the Moisture and Temperature Fields

During soil freezing, moisture migration and temperature distribution are strongly coupled through bidirectional interactions. On the one hand, the temperature gradient not only indirectly affects matric potential and moisture transport parameters by altering the viscosity and surface tension of water, but can also directly drive moisture migration. On the other hand, the latent heat associated with moisture phase change and variations in water content significantly affect the temperature distribution. To enable the coupled solution of the temperature and moisture fields, the ratio of the volumetric ice content to the volumetric unfrozen water content, $B_i$, is introduced as the key coupling variable, thereby linking the dependent variables of the two fields through the phase-equilibrium relationship. The resulting hydrothermal coupling relationship is expressed as [27,28]:

$$B_i={\displaystyle \frac{{\theta }_i}{{\theta }_u}}\left\{\begin{array}{c}\begin{array}{cc}1.1{\left({\displaystyle \frac{T}{T_f}}\right)}^b-1，& T<T_f\end{array}\\ \begin{array}{ccc}0,& & T\ge T_f\end{array}\end{array}\right.$$

(18)

where $T_f$ is the soil freezing temperature (°C), and $b$is a constant that varies with soil type and salt content, with values of 0.61 for sandy soil, 0.47 for silty soil, and 0.56 for clayey soil.

Based on the temperature-field and moisture-field formulations described above, the coupled hydrothermal governing equations for frozen soil can be expressed as:

$$\theta ={\theta }_w+{\displaystyle \frac{{\rho }_i}{{\rho }_w}}{\theta }_i$$

(19)

$$B_i={\displaystyle \frac{{\theta }_i}{{\theta }_w}}=f\left(T\right)$$

(20)

The above system of equations couples the temperature and moisture fields through the dynamic evolution of unfrozen water and ice contents, thereby describing the coupled evolution of moisture migration, latent heat associated with phase change, and temperature distribution during freeze-thaw processes. By explicitly incorporating the water–ice phase-transition relationship and the latent-heat term, bidirectional feedback between the two fields is represented.

4.2. Geometric Model Development and Validation

4.2.1. Model Development

COMSOL Multiphysics was employed for the numerical simulations. The temperature and moisture fields were implemented using the Coefficient Form PDE interface, whereas the displacement field was modeled using the built-in Solid Mechanics interface. In this way, a fully coupled thermo-hydro-mechanical (THM) model was established.

The general form of the coefficient-form partial differential equation in COMSOL Multiphysics is expressed as [29]:

$$e_a{\displaystyle \frac{{\partial }^{2}u}{\partial t^2}}+d_a{\displaystyle \frac{\partial u}{\partial t}}+\nabla \left(-c\nabla u-\alpha u+\gamma \right)+\beta \cdot \nabla u+au=f$$

(21)

where $e_a$ is the mass coefficient (s²); $d_a$ is the damping coefficient [J·m^-3·K^-1]; $c$ is the diffusion coefficient [W·m^-1·K^-1]; $\alpha$is the conservative-flux convection coefficient (m^-1); $\gamma$ is the conservative-flux source term (W·m^-2); $\beta$is the convection coefficient [W·m^-2·K^-1]; $a$is the absorption coefficient [W·m^-3·K^-1]; and $f$ is the source term (W·m^-3).

The governing equation for the temperature field was then combined with the hydrothermal coupling relationship, yielding:

$$\rho C\left(\theta \right){\displaystyle \frac{\partial T}{\partial t}}=\lambda \left(\theta \right){\nabla }^{2}T+L{\rho }_i\left({\theta }_s-{\theta }_r\right)\left[S{\displaystyle \frac{\partial B_i}{\partial t}}+B_i{\displaystyle \frac{\partial S}{\partial t}}\right]+L{\rho }_i{\theta }_r{\displaystyle \frac{\partial B_i}{\partial t}}$$

(22)

Rearranging Eq. (22) into the coefficient form required by the PDE interface gives:

$$\left\{\begin{array}{c}{d}_a=\rho C-L{\rho }_i\left[\left({\theta }_s-{\theta }_r\right)S+{\theta }_r\right]{\displaystyle \frac{\partial B\left(T\right)}{\partial T}}\\ c=\lambda \left(\theta \right)\\ f=L{\rho }_i\left({\theta }_s-{\theta }_r\right)B\left(T\right){\displaystyle \frac{\partial S}{\partial t}}\end{array}\right.$$

(23)

Because the remaining terms in the general coefficient-form equation have no corresponding terms in the temperature-field governing equation, the coefficients $a$, $\gamma$, $\beta$, $\alpha$, and $e_a$ were all set to zero.

Similarly, the coefficient inputs for the moisture field were obtained as:

$$\left\{\begin{array}{c}{d}_a=1+B_i{\displaystyle \frac{{\rho }_i}{{\rho }_w}}\\ c=D\left(S\right)\\ \gamma =-K\left(S\right)\\ a={\displaystyle \frac{{\rho }_i}{{\rho }_w}}{\displaystyle \frac{\partial B_i}{{\partial }_t}}\\ f={\displaystyle \frac{{\rho }_i}{{\rho }_w}}{\displaystyle \frac{{\theta }_r}{{\theta }_s-{\theta }_r}}{\displaystyle \frac{\partial B_i}{{\partial }_t}}\end{array}\right.$$

(24)

For the moisture-field equation, $e_a$, $\alpha$, and $\beta$ in the general coefficient-form equation were set to zero.

For the simulation of frost-heave deformation, the Thermal Expansion feature in the Solid Mechanics interface was used to represent changes in ice volume as an equivalent thermal strain. Based on previous studies, the threshold ice content for the onset of frost heave in loam was taken as 2.7% [30].

The relationship between the frost-heave coefficient and ice content was defined as:

$$\eta \left(x,y\right)=\left\{\begin{array}{ll}0.0156{\omega }_{\theta i}-0.0042,& {\omega }_{\theta i}>0.027\\ 0,& {\omega }_{\theta i}<0.027\end{array}\left(loam\right)\right.$$

(25)

where $w_{\theta i}$is the mass fraction corresponding to the ice content ${\theta }_i$, calculated as:

$${\omega }_{\theta i}={\displaystyle \frac{{\rho }_i}{{\rho }_{\omega }}}\cdot {\displaystyle \frac{{\theta }_i}{{\rho }_s}}$$

(26)

4.2.2. Model Validation

To ensure the reliability and numerical accuracy of the developed thermo-hydro-mechanical (THM) coupled model in reproducing temperature evolution, moisture migration and accumulation, and frost-heave deformation induced by water–ice phase transition in seasonally frozen soils, a laboratory one-dimensional freezing soil-column test conducted by Lu was selected for model validation. The hydrothermal coupling component and the Solid Mechanics component of the model were validated independently.

In the closed-system one-dimensional freezing soil-column experiment, monitoring points were arranged at depths of 10, 20, 30, and 40 cm below the top of the soil column, and temporal variations in temperature and water content were recorded at each depth [31]. Figure 7 compares the simulated and measured temperature and water-content values at different depths and freezing times. Overall, the simulated results reproduced the measured trends well at all four monitoring depths throughout the 0–80 h freezing period.

During the initial stage of freezing, the simulated rate of temperature decrease agreed well with the experimental observations, and the temperature curves were nearly coincident. The maximum absolute temperature error was less than 2 °C, with the relative error remaining within 8% over the corresponding range. For water content, the maximum absolute error was less than 3%, and the relative error remained within 10%. As the freezing front migrated downward, pronounced moisture migration and accumulation occurred near the freezing interface. At this stage, the mean absolute temperature error at each monitoring depth remained below 3 °C, while the relative error over the full range was within 12%.

As the phase-transition process approached completion, the rate of decrease in unfrozen water content gradually slowed and tended toward a stable state. The mean absolute error in unfrozen water content at each depth was less than 4%. A slightly larger deviation occurred at the shallow 10 cm monitoring point during the later stage, mainly because of the stronger influence of boundary heat exchange. Overall, the relative error of the hydrothermal model over the entire freezing period was within 12% for temperature and within 15% for unfrozen water content, while the relative errors during the active phase-transition stage were both maintained within 10%.

Figure 8 compares the measured and simulated displacement at the top of the soil column over time. Overall, the simulated and measured curves exhibited consistent temporal trends, both showing a distinct three-stage evolution characterized by rapid initial growth, a subsequent reduction in growth rate, and eventual stabilization. During the initial stage (0–60 h), both curves increased steeply, corresponding to the period of the most rapid frost-heave development over the entire freezing process. The simulated displacement was generally slightly lower than the measured value during this stage, with a maximum absolute error of approximately 0.014 cm and a corresponding relative error of approximately 28%. The two curves intersected at approximately 60 h. During the subsequent stage from 60 to 120 h, the displacement growth rates decreased simultaneously and gradually approached a stable state. The simulated values became slightly higher than the measured values, and the deviation increased slowly with time. At 120 h, the final absolute error was approximately 0.004 cm, corresponding to a relative error of approximately 6.3%. Throughout the freezing process, the simulated and measured displacements exhibited consistent evolutionary behavior and generally good agreement, indicating that the model can reasonably reproduce frost-heave deformation.

4.3. Simulation Results and Analysis for the Typical Cross-Section

4.3.1. Geometric Model and Mesh Discretization

Figure 9(a) shows the selected typical cross-section of the riverbank at the Shisifenzi Bend along the Inner Mongolia reach of the Yellow River. The riverbank is approximately 6 m high. The middle portion of the bank slope consists of an inclined section with a slope angle of approximately 45°, while the bank toe comprises a 2 m-high vertical bank face and a 0.2 m-high stepped gentle-slope transition section.

The computational domain was discretized using an unstructured triangular mesh, as shown in Figure 9(b). The maximum and minimum element sizes were set to 0.2 m and $7.5\times {10}^{-4}$m, respectively, resulting in a total of 2,941 mesh elements.

For temporal discretization, the time step was set to 1 day over the complete 130-day freeze-thaw period. The temperature field, moisture field, and Solid Mechanics field were incorporated into a unified transient solution framework. Through temporal discretization, the continuous transient process was divided into 130 discrete time steps, at which the coupled multi-physics equations were solved simultaneously. This scheme enabled the coupled freeze-thaw simulation to proceed on a daily basis from the prescribed initial conditions throughout the entire freeze-thaw period.

4.3.2. Boundary Conditions

A time-dependent surface temperature boundary was prescribed for the temperature field to reproduce the environmental temperature variation over a complete 130-day freeze-thaw period in the study area. The air temperature along the riverbank generally exhibits a sinusoidal seasonal variation. The daily mean temperature falls below 0 °C in November and decreases progressively to its minimum in late January, with a mean January temperature of approximately −15 °C. It then gradually rises and exceeds 0 °C around early April, corresponding to a complete freeze-thaw period of approximately 130 days.

For the moisture field, a closed-system boundary condition was adopted, with no external water supply or drainage considered. The analysis therefore focused on moisture migration and redistribution within the bank slope driven by temperature gradients and gravity. Zero-flux conditions were imposed on the top, bottom, and lateral boundaries to ensure conservation of water within the system. This configuration was used to isolate the dominant sequence of moisture redistribution during freezing and thawing, including temperature-gradient-driven migration toward the freezing front, moisture accumulation near the freezing front, and subsequent redistribution toward deeper soil layers during thawing.

For the mechanical field, the bottom boundary was fixed to constrain both vertical and horizontal displacements, while roller constraints were applied to the two lateral boundaries to restrict only the normal displacement. The slope surface was treated as a free boundary [32]. This combination of mechanical boundary conditions allows the deformation response of the bank slope under the combined effects of self-weight and freeze-thaw action to be represented while minimizing artificial disturbance of the internal stress–strain evolution caused by boundary constraints.

4.3.3. Simulation Results and Analysis

Figure 10 presents the spatial evolution of the stress and strain fields within the bank slope during the freeze-thaw process. As shown in Figure 10(a) and Figure 10(b), during the earlier freezing stage, relatively high stresses were mainly concentrated along the slope surface and in the upper part of the bank toe, whereas the stress contours within the interior of the slope remained comparatively smooth and the internal stress levels were generally low. By Day 90, corresponding to the deep-freezing stage, the stress field had undergone pronounced redistribution under the combined influence of sustained low-temperature boundary conditions and extensive freezing of the near-surface soil. Lateral compression and volumetric expansion induced by frost heave generated substantial additional stresses in the shallow slope and bank-toe regions. The originally low-stress region at the bank toe developed into a pronounced stress-concentration zone, where the peak stress exceeded 303.50 kPa.

This strong stress concentration associated with water–ice phase transition and frost-heave deformation indicates a marked increase in the mechanical loading of the bank toe during the deep-freezing stage. The bank toe therefore represents a potentially vulnerable location in the subsequent warming and thawing period, during which stress redistribution and release may promote shear deformation and reduce slope stability.

Consistent with the evolution of the stress field, the strain field exhibited clear spatial migration and accumulation between Days 45 and 90, as shown in Figure 10(c) and Figure 10(d). At Day 45, during the early freezing stage, the exposed slope surface and crest responded first to the low-temperature boundary, producing frost-heave-induced tensile strains of approximately 0.011. Meanwhile, a slight compressive strain of approximately −0.001 developed within the slope interior, indicating that deformation was still mainly concentrated in the shallow near-surface soil.

By Day 90, corresponding to the peak frost-heave stage, the freezing front had propagated substantially downward into the soil mass. The zones of frost-heave-induced tensile strain along the slope surface and around the bank toe progressively expanded and became interconnected, while strain contours exceeding 0.006 extended further into the slope interior. The high-strain zone near the bank toe also expanded further. Although the maximum surface strain decreased slightly from 0.011 at Day 45 to 0.009 at Day 90, this reduction does not necessarily indicate an attenuation of the overall frost-heave deformation. Rather, together with the spatial expansion of the high-strain region, it suggests that deformation became progressively redistributed from a localized near-surface response toward a broader zone within the slope as freezing advanced.

Figure 11 presents the time histories of cumulative displacement at different depths within the bank slope, together with the spatial distribution of displacement over the slope profile on Day 90. The displacement evolution curves indicate that the monitoring point at 1 m remained relatively stable throughout the simulation period, with only minor fluctuations, whereas the points at 2 m, 3 m, and 4 m exhibited pronounced dynamic responses. During the freezing period from Day 40 to Day 90, the displacements at 2 m, 3 m, and 4 m underwent successive sharp decreases or uplift-like fluctuations, reflecting structural adjustment of the soil induced by frost-heave-related volumetric change.

After Day 90, when the slope entered the thawing stage, especially during the interval from Day 90 to Day 120, the displacement curves at 2 m, 3 m, and 4 m all showed a distinct synchronous turning point followed by rapid increase, indicating that the thawing stage was associated with the most pronounced displacement variation during the entire freeze-thaw cycle. This intense displacement fluctuation suggests that, as air temperature rises and the near-surface frozen soil begins to thaw, the shallow bank soil undergoes a strong mechanical response under the combined influence of ice–water phase transition and soil-structure readjustment.

Because displacement variation was most pronounced during the thawing stage, the displacement field on Day 90 was selected for further spatial analysis. As shown in Figure 11(b), the displacement field exhibits marked spatial heterogeneity. The upper-left corner of the slope crest and the shoulder region form the main displacement concentration zone, where the maximum displacement reaches 1.01 mm. Relatively large displacement bands, ranging from 0.50 mm to 0.84 mm, are also distributed along the exposed slope face and the middle-to-lower bank-toe region, whereas the deeper and inner parts of the slope show much smaller displacements, approximately 0.17 mm. These results indicate that, under the strong effect of thawing, the bank surface—particularly the slope shoulder and slope face—develops a clear tendency for outward deformation toward the free face.

Figure 12 presents the relationship between the factor of safety and maximum displacement of the bank slope under different conditions. Under the original natural condition without freeze-thaw or seepage effects, the bank slope exhibited the highest stability margin, with a factor of safety of approximately 2.8. After one freeze-thaw cycle, deterioration of the soil structure resulted in a pronounced decrease in the factor of safety to approximately 2.05. When seepage was considered, the factor of safety further decreased to approximately 1.4 because pore-water pressure generated by groundwater reduced the effective stress within the soil, bringing the slope closer to a critical stability state.

When freeze-thaw and seepage effects were considered simultaneously, the factor of safety decreased markedly to approximately 1.1, well below the adopted stability control threshold of 1.3 [33,34]. Comparison of the four conditions indicates that slope stability is jointly controlled by freeze-thaw and seepage effects, and that their combined action produces a stronger reduction in stability than either process considered individually.

From a mechanistic perspective, freeze-thaw action may promote the development of microcracks within the soil, providing preferential pathways for subsequent seepage. Seepage, in turn, can weaken the soil structure and intensify freeze-thaw-induced deterioration. The interaction between these processes may therefore progressively reduce the stability of the bank slope and increase its susceptibility to overall sliding failure during the spring thaw period.

5. Conclusions

This study investigated the embankments along the Inner Mongolia reach of the Yellow River, which are strongly affected by seasonal freeze-thaw processes. Basic physical tests and multivariable mechanical tests were conducted on sampled soils to examine the effects of freeze-thaw cycling, water content, and related factors on soil mechanical behavior. A thermo-hydro-mechanical (THM) coupled model for riverbank stability was subsequently developed, validated, and applied to the stability analysis of a typical riverbank cross-section. The main conclusions are as follows:

(1) Compared with water content and the number of freeze-thaw cycles, dry density was the dominant factor controlling soil shear strength. The soil exhibited its highest shear strength near the optimum water content. Increasing the number of freeze-thaw cycles caused pronounced degradation of soil cohesion, with the most substantial reduction occurring during the first three cycles; the rate of degradation decreased thereafter. By contrast, the internal friction angle was much less sensitive to freeze-thaw cycling.

(2) The developed THM coupled model satisfactorily reproduced the evolution of temperature, moisture, and displacement within the soil during freezing and thawing. When applied to the typical riverbank cross-section, the model captured the continuous evolution of the internal stress–strain response and displacement throughout the freeze-thaw period, identified the bank toe and near-surface slope regions as the most unfavorable locations for stability, and quantified variations in the factor of safety under different stages and influencing conditions.