Three-Dimensional Numerical Investigation on the Seepage Field and Stability of Soil Slope Subjected to Snowmelt Infiltration

Cutting slope failures occur frequently along the high-speed railways in Northeast China during the construction due to snowmelt infiltration. This study addresses this issue by applying a three-dimensional numerical model. The influence of the depth of accumulated snow (ds), daily temperature variation (ΔT), and freeze-thaw (F-T) cycles on the seepage field and stability of cutting slopes is discussed. The results demonstrate that water seepage due to snowmelt infiltration primarily extends through the ground surface by about 10 m. The deep-seated instability is likely to occur under a prolonged and highly accumulated infiltration, while shallow failure is associated with intense, short-duration snowmelt infiltration. The maximum degree of saturation (Sr) and pore-water pressure (PWP) values are observed at the slope toe. Increasing ds and ΔT increase the Sr and PWP due to snowmelt infiltration and thereby decreases cutting slope stability. Compared to the ds and ΔT, the F-T cycle is more likely to cause slope failure. In addition, the F-T cycle also induces the reduction of soil strength and the crack propagation. Overall, the conducted study provided useful help toward the process of safer design for cutting slope along the high-speed railway in seasonally cold regions.


Introduction
Snowmelt infiltration affects the moisture content of soils in seasonally cold regions, and extensive infiltration due to continuous snowmelt over long periods of thawing may trigger landslides [1][2][3][4][5][6][7][8] (Figure 1, from a study on a snowmelt estimation method for road management), which can result in significant loss of life and damage to public and private property. Therefore, the influence of snowmelt infiltration on the soil slope stability in seasonally cold regions cannot be ignored.
This issue is particularly significant for ensuring the safe operation and maintenance of the high-speed railways constructed in Northeast China, which is a typical seasonally cold region [9][10][11]. While some studies have reported on the occurrence of cutting slope failures due to snowmelt infiltration during the construction of high-speed railways [12], the majority of studies have been restricted to investigating the stability of soil slopes subjected to rainfall [13][14][15][16][17] and those within permafrost regions [7,18]. Moreover, the limited efforts focused on exploring the stability of soil slopes under the influence of snowmelt infiltration have relied typically on two-dimensional (2D) numerical analyses [3,5,7]. However, the slope stability in seasonally cold regions has been demonstrated to be affected by freeze-thaw (F-T) cycles as well [10]. For example, cracking is a common natural phenomenon observed within clayey soils in seasonally cold regions as a result of multiple F-T cycles and significantly influences the mechanical and hydraulic behaviors of these soils [17,19]. The effects of F-T cycles on crack propagation have been further demonstrated by computerized tomography (CT) images of clayey soil in a cutting slope along a high-speed railway, where multiple F-T cycles are observed to induce an increase in the concentration of micro-voids in the soil microstructure [20]. Nonetheless, few studies have investigated the effects of soil cracking on slope stability. Moreover, the effects of F-T cycles generate complicated slope problems that necessitate an analysis of three-dimensional (3D) geometric effects. Therefore, 3D numerical analysis of soil slope stability under the combined effects of snowmelt infiltration and F-T cycling has significant potential for directing efforts toward mitigating these effects in seasonally cold regions.
The present study addresses the above discussed issues by applying a 3D finite element (FE) model to investigate the cutting slope stability of a typical section along the Jilin-Tumen-Hunchun high-speed railway subjected to snowmelt infiltration. We first validate the numerical analysis methodology by comparing the pore-water pressure (PWP), the degree of saturation (Sr), the soil slope stability, and the slip surface predicted by the model with those obtained in [21]. Then, we investigate the influences of the depth of accumulated snow (ds), daily temperature variation (ΔT), and F-T cycling on the seepage field (i.e., PWP and Sr) and the stability of cutting slopes arising from snowmelt infiltration. Finally, some important insights for cutting slope along the high-speed railway subjected to snowmelt infiltration are presented in conclusion.

Three-Dimensional Numerical Modeling
The geometry of the cutting slope model and the FE mesh used are presented in Figure 2. The profile consists of two soil layers, including a silty clay layer and a sandstone Moreover, the limited efforts focused on exploring the stability of soil slopes under the influence of snowmelt infiltration have relied typically on two-dimensional (2D) numerical analyses [3,5,7]. However, the slope stability in seasonally cold regions has been demonstrated to be affected by freeze-thaw (F-T) cycles as well [10]. For example, cracking is a common natural phenomenon observed within clayey soils in seasonally cold regions as a result of multiple F-T cycles and significantly influences the mechanical and hydraulic behaviors of these soils [17,19]. The effects of F-T cycles on crack propagation have been further demonstrated by computerized tomography (CT) images of clayey soil in a cutting slope along a high-speed railway, where multiple F-T cycles are observed to induce an increase in the concentration of micro-voids in the soil microstructure [20]. Nonetheless, few studies have investigated the effects of soil cracking on slope stability. Moreover, the effects of F-T cycles generate complicated slope problems that necessitate an analysis of three-dimensional (3D) geometric effects. Therefore, 3D numerical analysis of soil slope stability under the combined effects of snowmelt infiltration and F-T cycling has significant potential for directing efforts toward mitigating these effects in seasonally cold regions.
The present study addresses the above discussed issues by applying a 3D finite element (FE) model to investigate the cutting slope stability of a typical section along the Jilin-Tumen-Hunchun high-speed railway subjected to snowmelt infiltration. We first validate the numerical analysis methodology by comparing the pore-water pressure (PWP), the degree of saturation (S r ), the soil slope stability, and the slip surface predicted by the model with those obtained in [21]. Then, we investigate the influences of the depth of accumulated snow (d s ), daily temperature variation (∆T), and F-T cycling on the seepage field (i.e., PWP and S r ) and the stability of cutting slopes arising from snowmelt infiltration. Finally, some important insights for cutting slope along the high-speed railway subjected to snowmelt infiltration are presented in conclusion.

Three-Dimensional Numerical Modeling
The geometry of the cutting slope model and the FE mesh used are presented in Figure 2. The profile consists of two soil layers, including a silty clay layer and a sandstone layer, and a retaining structure composed of a sheet pile wall system with a pile length of 12 m. The initial ground-water table (GWT) is at depth ranging from about 24.8 m (10 m from the model bottom) to 22.8 m (12 m from the model bottom) within the slope profile, as illustrated in Figure 2a. A vertical view of the soil slope-sheet pile wall system and the detailed dimensions of the pile and sheet are presented in Figure 2b,c, respectively. A horizontal view of the model is presented in Figure 2d  direction. The pile and sheet were both assumed to be linearly elastic for simplicity [22,23]. The 8-node brick elements were employed to simulate the sheet pile walls. In addition, 8-node brick, trilinear displacement, and trilinear pore pressure elements were used to model the soil domain [24,25]. layer, and a retaining structure composed of a sheet pile wall system with a pile length of 12 m. The initial ground-water table (GWT) is at depth ranging from about 24.8 m (10 m from the model bottom) to 22.8 m (12 m from the model bottom) within the slope profile, as illustrated in Figure 2a. A vertical view of the soil slope-sheet pile wall system and the detailed dimensions of the pile and sheet are presented in Figure 2b, c, respectively. A horizontal view of the model is presented in Figure 2d and the model width is 5 m in the y direction. The pile and sheet were both assumed to be linearly elastic for simplicity [22,23]. The 8-node brick elements were employed to simulate the sheet pile walls. In addition, 8-node brick, trilinear displacement, and trilinear pore pressure elements were used to model the soil domain [24,25]. The soil water characteristic curves (SWCCs) of the silty clay and sandstone are presented in Figure 3a, whereas the corresponding relative hydraulic conductivities, given as k w /k s , where k w and k s are the unsaturated and saturated coefficients of permeability, respectively, are presented in Figure 3b. The SWCCs of the soils were obtained using the Fredlund and Xing (1994) [26] equation with a correction factor C (ψ) = 1, as recommended by Leong and Rahardjo (1997a) [27]. The values of k s used for silty clay and sandstone were set to 1.73 × 10 −6 and 10 −4 , respectively, and the ratios k w /k s were obtained using the indirect procedure presented by Leong and Rahardjo (1997b) [28]. The SWCC parameters (θ s , a, m, and n) and the parameters (p) used for calculating k w /k s for silty clay and sandstone are listed in Table 1.
The soil water characteristic curves (SWCCs) of the silty clay and sandstone are presented in Figure 3a, whereas the corresponding relative hydraulic conductivities, given as kw/ks, where kw and ks are the unsaturated and saturated coefficients of permeability, respectively, are presented in Figure 3b. The SWCCs of the soils were obtained using the Fredlund and Xing (1994) [26] equation with a correction factor C (ψ) = 1, as recommended by Leong and Rahardjo (1997a) [27]. The values of ks used for silty clay and sandstone were set to 1.73 × 10 −6 and 10 −4 , respectively, and the ratios kw/ks were obtained using the indirect procedure presented by Leong and Rahardjo (1997b) [28]. The SWCC parameters (θs, a, m, and n) and the parameters (p) used for calculating kw/ks for silty clay and sandstone are listed in Table 1.
Frictional contact was considered for the interfaces between the sheet pile walls and the soil layers [23,29], and the friction coefficients μ were 0.3 and 0.52 for the interfaces involving the silty clay and sandstone layers, respectively. Besides, stiff contact was used in the normal direction for the pile elements and sheet elements, whereas the friction coefficient in the tangential direction was assumed to be 0.7 [23].      Frictional contact was considered for the interfaces between the sheet pile walls and the soil layers [23,29], and the friction coefficients µ were 0.3 and 0.52 for the interfaces involving the silty clay and sandstone layers, respectively. Besides, stiff contact was used in the normal direction for the pile elements and sheet elements, whereas the friction coefficient in the tangential direction was assumed to be 0.7 [23].
Based on the previous work of Qi and Vanapalli (2015) [30], both the PWP at the x = 0 m and x = 71.15 m faces are assumed to increase linearly with respect to increasing depth below the GWT. Water at the lateral boundaries (i.e., y = 0 and y = 5 m) as well as at the bottom boundary of the computational domain is assumed to have a zero flux. In addition, all nodes at the lateral boundaries (i.e., y = 0 and y = 5 m) are constrained by displacement degrees of freedom (DOF) in the y direction, and the horizontal displacement DOF is constrained at the x = 0 m and x = 71.15 m faces. The displacement DOF of nodes at the z = 0 m face are fully fixed in all directions [23,31,32].

Snowmelt Infiltration
The climatic conditions used in the numerical model were based on meteorological monitoring data (i.e., monthly average snowfall and cumulative snowfall) obtained for Yanji, China, from 2010 to 2020, which are presented in Figure 4. We note that the accumulative snowfall depth d s are close to 1 m. We also note that d s can vary substantially. For example, the region of Yanji in Jilin Province experienced unusually high snowfall on 18 January 2016, and the maximum local snowfall was greater than 1 m. a (kPa) 28 10 p 4.7 4 m 1 1 n 1 1

Snowmelt Infiltration
The climatic conditions used in the numerical model were based on meteorological monitoring data (i.e., monthly average snowfall and cumulative snowfall) obtained for Yanji, China, from 2010 to 2020, which are presented in Figure 4. We note that the accumulative snowfall depth ds are close to 1 m. We also note that ds can vary substantially. For example, the region of Yanji in Jilin Province experienced unusually high snowfall on 18 January 2016, and the maximum local snowfall was greater than 1 m.  The process of snowmelt infiltration is illustrated in Figure 5a [33]. This infiltration process is known to be affected by many factors, such as soil temperature, soil freezing depth, the water content of the soil before winter, and the depth of the snow cover [34]. Because of the complexity of the snowmelt infiltration, therefore, the present work uses the snow water equivalent (SWE) to simplify the analysis of the process [35,36]. Here, the SWE is the depth of water that would cover the ground if the snow cover was in a liquid state and is used for light-, normal-, and heavy-snow years. In this study, the SWE is equivalent to daily rainfall. In addition, the effect of potential root biomass or incorporated residue is not considered in this study.
According to the monthly average snowfall and cumulative snowfalls in Yanji, China, the standard simulations apply an SWE rainfall intensity (Ir) of 3 mm/h, which is equivalent to ds = 0.72 m, to the top surface of the model slope for a period of 48 h, as illustrated in Figure 5b. This represents the standard value of Ir used in all simulations unless otherwise specified. The SWE rainfall intensity on the slope surface is Ir·cosα, where α is the slope angle. Here, Ir is less than the permeability of the silty clay layer. Therefore, all snowmelt infiltrated into the soil, and no run-off occurred in the simulations employing standard conditions. Moreover, the seepage field and stability of the unsaturated soil slope were investigated from 0 h to 96 h. The process of snowmelt infiltration is illustrated in Figure 5a [33]. This infiltration process is known to be affected by many factors, such as soil temperature, soil freezing depth, the water content of the soil before winter, and the depth of the snow cover [34]. Because of the complexity of the snowmelt infiltration, therefore, the present work uses the snow water equivalent (SWE) to simplify the analysis of the process [35,36]. Here, the SWE is the depth of water that would cover the ground if the snow cover was in a liquid state and is used for light-, normal-, and heavy-snow years. In this study, the SWE is equivalent to daily rainfall. In addition, the effect of potential root biomass or incorporated residue is not considered in this study.

FE Method with Shear Strength Reduction and Soil Model
The present work adopts the FE method with shear strength reduction technique because it has been demonstrated to be a reliable and robust approach for assessing the safety factors of slopes and locating the corresponding critical slip surfaces [37]. The Mohr-Coulomb model was employed to simulate the soil behaviors. Here, the shear strength (τ) of an unsaturated soil is expressed based on the Bishop effective stress principle outlined by Vanapalli et al. [38] as follows: where c′ is the effective cohesion, is the net total stress, ua is the pore-air pressure, uw  According to the monthly average snowfall and cumulative snowfalls in Yanji, China, the standard simulations apply an SWE rainfall intensity (I r ) of 3 mm/h, which is equivalent to d s = 0.72 m, to the top surface of the model slope for a period of 48 h, as illustrated in Figure 5b. This represents the standard value of I r used in all simulations unless otherwise specified. The SWE rainfall intensity on the slope surface is I r ·cosα, where α is the slope angle. Here, I r is less than the permeability of the silty clay layer. Therefore, all snowmelt infiltrated into the soil, and no run-off occurred in the simulations employing standard conditions. Moreover, the seepage field and stability of the unsaturated soil slope were investigated from 0 h to 96 h.

FE Method with Shear Strength Reduction and Soil Model
The present work adopts the FE method with shear strength reduction technique because it has been demonstrated to be a reliable and robust approach for assessing the safety factors of slopes and locating the corresponding critical slip surfaces [37]. The Mohr-Coulomb model was employed to simulate the soil behaviors. Here, the shear strength (τ) of an unsaturated soil is expressed based on the Bishop effective stress principle outlined by Vanapalli et al. [38] as follows: where c is the effective cohesion, σ n is the net total stress, u a is the pore-air pressure, u w is the pore-water pressure, ϕ is the effective angle of internal friction, and χ is a function of S r , which varies between 0 (completely dry condition) and 1 (fully saturated condition).
Here, χ can be expressed in terms of the volumetric water content of the soil (θ w ) as where θ s is the saturated volumetric water content, and θ r is the residual volumetric water content. The calculated shear strengths of the silty clay and sandstone layers used in the numerical analyses are displayed in Table 1.

Validation of the Numerical Analysis Methodology
The hydraulic and mechanical boundary conditions of the soil slope prevailing at the Pohang site used in the work of Oh and Lu [21] are presented in Figure 6. In addition, the values used for the geotechnical and hydromechanical properties of the Pohang soil are listed in Table 2.
hematic view of (a) the process of snowmelt infiltration (after [33]) (b) and the equivalent rainfall infiltration.

FE Method with Shear Strength Reduction and Soil Model
The present work adopts the FE method with shear strength reduction technique because it has been demonstrated to be a reliable and robust approach for assessing the safety factors of slopes and locating the corresponding critical slip surfaces [37]. The Mohr-Coulomb model was employed to simulate the soil behaviors. Here, the shear strength (τ) of an unsaturated soil is expressed based on the Bishop effective stress principle outlined by Vanapalli et al. [38] as follows: where c′ is the effective cohesion, is the net total stress, ua is the pore-air pressure, uw is the pore-water pressure, φ′ is the effective angle of internal friction, and χ is a function of Sr, which varies between 0 (completely dry condition) and 1 (fully saturated condition). Here, χ can be expressed in terms of the volumetric water content of the soil (θw) as where θs is the saturated volumetric water content, and θr is the residual volumetric water content. The calculated shear strengths of the silty clay and sandstone layers used in the numerical analyses are displayed in Table 1.

Validation of the Numerical Analysis Methodology
The hydraulic and mechanical boundary conditions of the soil slope prevailing at the Pohang site used in the work of Oh and Lu [21] are presented in Figure 6. In addition, the values used for the geotechnical and hydromechanical properties of the Pohang soil are listed in Table 2.    [21].

Type Parameter Pohang Soil
Strength The PWP and S r results obtained by the proposed numerical methodology for the Pohang site are presented Figure 7a,b, respectively. In addition, the potential slide surface of the soil slope obtained by the proposed methodology is compared with that obtained by Oh and Lu (2015) [21] in Figure 7c. Here, the three typical points denoted in Figure 7 as a, b, and c, respectively, represent the upper, middle, and lower regions along the failure surface.  The PWP and Sr results obtained by the proposed numerical methodology for the Pohang site are presented Figure 7a, b, respectively. In addition, the potential slide surface of the soil slope obtained by the proposed methodology is compared with that obtained by Oh and Lu (2015) [21] in Figure 7c. Here, the three typical points denoted in Figure 7 as a, b, and c, respectively, represent the upper, middle, and lower regions along the failure surface.  The values of PWP, S r , the factor of safety (FOS), and the slip surface in Oh and Lu [21] have been validated by filed evidences. Therefore, these results can be used to validate the correctness of numerical analysis methodology. We note that the PWP values obtained by the proposed methodology on 1 April 2011 are −85 kPa at point a, −55 kPa at point b, and −35 kPa at point c. We further note that the values of S r obtained at point c are always greater and approach a saturated state earlier than those observed at points a and b. The numerically determined values of S r on 1 April 2011 are 89%, 92.5%, and 95% at points a, b, and c, respectively. In comparison, the PWP values reported by Oh and Lu [21] were −80 kPa at point a, −60 kPa at point b, and −40 kPa at point c, whereas the S r values were 91.5%, 92.5%, and 94% at points a, b, and c, respectively. Accordingly, we can conclude that the results obtained by the proposed numerical methodology for the PWP and S r are in good agreement with the results presented in [21]. We further note from Figure 7c that the shape and position of the failure surface obtained by the proposed methodology compare satisfactorily with the results presented in [21]. Moreover, the factor of safety (FOS) of the soil slope at the Pohang site was calculated in [21] to be 2.8. Meanwhile, the FOS value predicted by the proposed methodology was 2.78, which again demonstrates that the results obtained by the proposed numerical methodology are in good agreement with the results presented in [21].
Therefore, this numerical analysis methodology can be used to investigate the seepage field and stability of soil slope subjected to snowmelt infiltration.

Snowmelt Infiltration Analysis
Before proceeding with an analysis of snowmelt infiltration, we must first evaluate the slope model in an initial state without snowmelt infiltration. The S r values and the PWP contours in the initial state are presented in Figure 8. Differences in the hydraulic properties of the silty clay and sandstone layers (i.e., the SWCCs and k s values) result in significant changes in S r at the interface between the two layers. Here, the value of S r at the upper surface layer is about 55%, and it increases along the slope surface (Figure 8a). The PWP displays a linear relation with the depth (Figure 8b). The FOS obtained for the cutting slope by the FE method with shear strength reduction technique was 2.81. The values of PWP, Sr, the factor of safety (FOS), and the slip surface in Oh and Lu [21] have been validated by filed evidences. Therefore, these results can be used to validate the correctness of numerical analysis methodology. We note that the PWP values obtained by the proposed methodology on 1 April 2011 are −85 kPa at point a, −55 kPa at point b, and −35 kPa at point c. We further note that the values of Sr obtained at point c are always greater and approach a saturated state earlier than those observed at points a and b. The numerically determined values of Sr on 1 April 2011 are 89%, 92.5%, and 95% at points a, b, and c, respectively. In comparison, the PWP values reported by Oh and Lu [21] were −80 kPa at point a, −60 kPa at point b, and −40 kPa at point c, whereas the Sr values were 91.5%, 92.5%, and 94% at points a, b, and c, respectively. Accordingly, we can conclude that the results obtained by the proposed numerical methodology for the PWP and Sr are in good agreement with the results presented in [21]. We further note from Figure 7c that the shape and position of the failure surface obtained by the proposed methodology compare satisfactorily with the results presented in [21]. Moreover, the factor of safety (FOS) of the soil slope at the Pohang site was calculated in [21] to be 2.8. Meanwhile, the FOS value predicted by the proposed methodology was 2.78, which again demonstrates that the results obtained by the proposed numerical methodology are in good agreement with the results presented in [21].
Therefore, this numerical analysis methodology can be used to investigate the seepage field and stability of soil slope subjected to snowmelt infiltration.

Snowmelt Infiltration Analysis
Before proceeding with an analysis of snowmelt infiltration, we must first evaluate the slope model in an initial state without snowmelt infiltration. The Sr values and the PWP contours in the initial state are presented in Figure 8. Differences in the hydraulic properties of the silty clay and sandstone layers (i.e., the SWCCs and ks values) result in significant changes in Sr at the interface between the two layers. Here, the value of Sr at the upper surface layer is about 55%, and it increases along the slope surface (Figure 8a). The PWP displays a linear relation with the depth (Figure 8b). The FOS obtained for the cutting slope by the FE method with shear strength reduction technique was 2.81.  The S r values and PWP contours obtained for the cutting slope under snowmelt infiltration (I r = 3 mm/h) are presented in Figures 9 and 10, respectively. The results in Figure 9a indicate that the value of S r at the upper surface layer is greater than 70% when snowmelt infiltration is terminated at time t = 48 h. As shown in Figure 9b, snowmelt infiltrates into the soil slope at t = 72 h, and the pressure head of the slope decreases, resulting in a decrease in S r . Furthermore, Figure 9c indicates that the continuous infiltration of snowmelt at t = 96 h results in further decrease in S r and PWP at the top and surface of the slope, and the seepage field gradually approaches the field observed in the initial state.
The Sr values and PWP contours obtained for the cutting slope under snowmelt infiltration (Ir = 3 mm/h) are presented in Figures 9 and 10, respectively. The results in Figure  9a indicate that the value of Sr at the upper surface layer is greater than 70% when snowmelt infiltration is terminated at time t = 48 h. As shown in Figure 9b, snowmelt infiltrates into the soil slope at t = 72 h, and the pressure head of the slope decreases, resulting in a decrease in Sr. Furthermore, Figure 9c indicates that the continuous infiltration of snowmelt at t = 96 h results in further decrease in Sr and PWP at the top and surface of the slope, and the seepage field gradually approaches the field observed in the initial state.
The above results indicate that due to the PWP diffusion process, the deep-seated landslide is likely to occur under a prolonged and highly accumulated infiltration, while the shallow slope failure is associated with intense, short-duration snowmelt infiltration.
In addition, the two locations denoted as L1 and L2 shown in Figure 11 were selected to illustrate the variations in the Sr and PWP in the inner region of the slope and the slope surface, respectively. Variations in Sr and PWP along the two lines under snowmelt infiltration are discussed in detail in later sections along with the computed FOS values.   The above results indicate that due to the PWP diffusion process, the deep-seated landslide is likely to occur under a prolonged and highly accumulated infiltration, while the shallow slope failure is associated with intense, short-duration snowmelt infiltration.
In addition, the two locations denoted as L1 and L2 shown in Figure 11 were selected to illustrate the variations in the S r and PWP in the inner region of the slope and the slope surface, respectively. Variations in S r and PWP along the two lines under snowmelt infiltration are discussed in detail in later sections along with the computed FOS values. (c)

Effect of Accumulated Snow Depth on Soil Slope Stability
The effect of d s on the seepage field and slope stability was investigated by applying d s values of 0.48 m, 0.72 m, 0.96 m, and 1.2 m to the numerical model. Here, the snowfall accumulated above the slope during the winter season will begin melting once the air temperature is greater 0 • C, and the snowmelt will infiltrate into the soil or runoff above the slope.
The obtained values of S r and PWP at lines L1 and L2 (t = 48 h) are plotted in Figure 12 with respect to d s . Obviously, the shapes of the S r and the PWP profiles along line L1 are quite different from those obtained in the initial state without snowmelt infiltration. As shown in Figure 12a

Effect of Accumulated Snow Depth on Soil Slope Stability
The effect of ds on the seepage field and slope stability was investigated by applying ds values of 0.48 m, 0.72 m, 0.96 m, and 1.2 m to the numerical model. Here, the snowfall accumulated above the slope during the winter season will begin melting once the air temperature is greater 0 °C, and the snowmelt will infiltrate into the soil or runoff above the slope.
The obtained values of Sr and PWP at lines L1 and L2 (t = 48 h) are plotted in Figure  12 with respect to ds. Obviously, the shapes of the Sr and the PWP profiles along line L1 are quite different from those obtained in the initial state without snowmelt infiltration. As shown in Figure 12a  We also note that the values of Sr and PWP initially decreased along the line L1 and then increased with further along the L1. The results in Figure 12c, d indicate that Sr and PWP increase along the slope surface, and their maximum values are observed at the slope toe. This may be a reasonable explanation as to why snowmelt-induced soil slope failures usually initiate from the slope toe. Compared with the values of Sr observed at the top and We also note that the values of S r and PWP initially decreased along the line L1 and then increased with further along the L1. The results in Figure 12c,d indicate that S r and PWP increase along the slope surface, and their maximum values are observed at the slope toe. This may be a reasonable explanation as to why snowmelt-induced soil slope failures usually initiate from the slope toe. Compared with the values of S r observed at the top and slope toe to that in the initial state, the corresponding values of S r obtained at d s = 1.2 m increase by about 24.3% and 17.4%, respectively. We note from Figure 12c,d that the S r and PWP values observed at a horizontal distance of approximately 15.6 m from the slope crest are located within a "stable region" and that the S r and PWP almost remain the same (about 3.2 m in x direction) for the initial state. However, no such "stable region" is observed with increasing d s .
The FOS values obtained for the cutting slope are plotted in Figure 13 with respect to d s . Here, the FOS of the cutting slope decreases significantly with increasing d s , which indicates that the soil slope stability is negatively affected by increasing snowmelt infiltration. In addition, the FOS obtained at the moment at which snowmelt infiltration is terminated (t = 48 h) is significantly less than that obtained at t = 96 h because the seepage field gradually tends toward its initial state with increasing time.
Water 2021, 13, x FOR PEER REVIEW (about 3.2 m in x direction) for the initial state. However, no such "stable region served with increasing ds.
The FOS values obtained for the cutting slope are plotted in Figure 13 with re ds. Here, the FOS of the cutting slope decreases significantly with increasing ds, w dicates that the soil slope stability is negatively affected by increasing snowmelt tion. In addition, the FOS obtained at the moment at which snowmelt infiltration nated (t = 48 h) is significantly less than that obtained at t = 96 h because the seepa gradually tends toward its initial state with increasing time.

Effect of Daily Temperature Variation on Soil Slope Stability
Daily temperature variation (ΔT) is a key factor affecting the distribution o content and the stability of soil slopes during spring thaw periods due to its e snowmelt and the resulting surface infiltration [5]. The degree-day factor (mm·°C is typically used to convert the degree-days to snowmelt or ice melt (i.e., the d snowmelt induced by a temperature increase of 1 °C). The value of the degree-da varies with respect to the melt period due to variations in the properties of the ice [39]. According to [36], an average value of the degree-day factor for snow in th of interest over the period considered in this study was 4.8 mm·°C −1 ·day −1 .
The values of Sr and PWP obtained along lines L1 and L2 are plotted in F with respect to ΔT. We note that increasing ΔT from 2 °C to 10 °C yields an incre and PWP along both lines L1 and L2. The results in Figure 14a, b indicate that th of ΔT on the seepage field along line L1 at depths greater than 10 m is not obviou ever, snowmelt infiltration is observed to affect the seepage field in the top 10 m th of the soil layer seriously along line L1. The values of Sr observed at the top and s at ΔT = 10 °C are about 10.8% and 7.3% greater than the corresponding values o in the initial state. However, these differences in the seepage field are small at ΔT

Effect of Daily Temperature Variation on Soil Slope Stability
Daily temperature variation (∆T) is a key factor affecting the distribution of water content and the stability of soil slopes during spring thaw periods due to its effect on snowmelt and the resulting surface infiltration [5]. The degree-day factor (mm· • C −1 ·day −1 ) is typically used to convert the degree-days to snowmelt or ice melt (i.e., the depth of snowmelt induced by a temperature increase of 1 • C). The value of the degree-day factor varies with respect to the melt period due to variations in the properties of the snow or ice [39]. According to [36], an average value of the degree-day factor for snow in the region of interest over the period considered in this study was 4.8 mm· • C −1 ·day −1 .
The values of S r and PWP obtained along lines L1 and L2 are plotted in Figure 14 with respect to ∆T. We note that increasing ∆T from 2 • C to 10 • C yields an increase in S r and PWP along both lines L1 and L2. The results in Figure 14a,b indicate that the effect of ∆T on the seepage field along line L1 at depths greater than 10 m is not obvious. However, snowmelt infiltration is observed to affect the seepage field in the top 10 m thickness of the soil layer seriously along line L1. The values of S r observed at the top and slope toe at ∆T = 10 • C are about 10.8% and 7.3% greater than the corresponding values observed in the initial state. However, these differences in the seepage field are small at ∆T = 2 • C. The FOS profile of the cutting slope is plotted in Figure 15 with respective ΔT. Here, the FOS decreases dramatically with increasing ΔT. This is because snowmelt infiltration increases with increasing ΔT, which increases the water content in the slope and thereby decreases slope stability. Again, we note that the FOS observed at the termination of snowmelt infiltration (t = 48 h) is significantly less than that observed at t = 96 h. However, we further note that the effect of ΔT on the FOS is relatively small compared with the corresponding effect of ds.  The FOS profile of the cutting slope is plotted in Figure 15 with respective ∆T. Here, the FOS decreases dramatically with increasing ∆T. This is because snowmelt infiltration increases with increasing ∆T, which increases the water content in the slope and thereby decreases slope stability. Again, we note that the FOS observed at the termination of snowmelt infiltration (t = 48 h) is significantly less than that observed at t = 96 h. However, we further note that the effect of ∆T on the FOS is relatively small compared with the corresponding effect of d s . The FOS profile of the cutting slope is plotted in Figure 15 with respective ΔT. Here the FOS decreases dramatically with increasing ΔT. This is because snowmelt infiltration increases with increasing ΔT, which increases the water content in the slope and thereby decreases slope stability. Again, we note that the FOS observed at the termination of snow melt infiltration (t = 48 h) is significantly less than that observed at t = 96 h. However, we further note that the effect of ΔT on the FOS is relatively small compared with the corre sponding effect of ds.

Effect of Freeze-Thaw Cycles on Soil Slope Stability
A schematic illustrating the process of crack propagation in a soil slope under F-T cycle is presented in Figure 16. The Figure 16a is the CT images of silty clay after various number of F-T cycles (N FT ) (typical Section 2 of soil samples is chosen to depict the microstructure variations, and a detailed introduction was provided by   [10]). It is obvious that as the N FT increases, the micro voids become clearly visible, as shown in the image in Figure 16a. After the F-T cycles, the soil aggregates separate and disrupt the interlocking of soil grains, thus leading to the damage in the microstructure of the soils, and more micro voids will generate. The increase in micro voids in the microstructure due to cyclic F-T deterioration results in generation of cracks at the surface of silty clay, which is increased as the N FT increases (Figure 16b). Due to the crack's generation in the surface of soils, more snowmelt water will infiltrate into the soil slope; therefore, referring to Zhang et al. [17], different SWE rainfall intensity (I r ) values of 3.2 mm/h, 4.0 mm/h, 4.8 mm/h, and 5.2 mm/h are employed for the first, third, seventh, and eleventh F-T cycles, respectively (Figure 16c). In addition, the analysis results presented in [20] have demonstrated that only the top 2.4-m thickness of the surface layer is substantially affected by crack propagation during multiple F-T cycles. Therefore, the present work considered crack propagation only within this depth zone (Figure 17), whereas the effective shear strengths of the silty clay layer below the 2.4-m level were maintained throughout.

Effect of Freeze-Thaw Cycles on Soil Slope Stability
A schematic illustrating the process of crack propagation in a soil slope under F-T cycle is presented in Figure 16. The Figure 16a is the CT images of silty clay after various number of F-T cycles (NFT) (typical Section 2 of soil samples is chosen to depict the microstructure variations, and a detailed introduction was provided by   [10]). It is obvious that as the NFT increases, the micro voids become clearly visible, as shown in the image in Figure 16a. After the F-T cycles, the soil aggregates separate and disrupt the interlocking of soil grains, thus leading to the damage in the microstructure of the soils, and more micro voids will generate. The increase in micro voids in the microstructure due to cyclic F-T deterioration results in generation of cracks at the surface of silty clay, which is increased as the NFT increases (Figure 16b). Due to the crack's generation in the surface of soils, more snowmelt water will infiltrate into the soil slope; therefore, referring to Zhang et al. [17], different SWE rainfall intensity (Ir) values of 3.2 mm/h, 4.0 mm/h, 4.8 mm/h, and 5.2 mm/h are employed for the first, third, seventh, and eleventh F-T cycles, respectively (Figure 16c). In addition, the analysis results presented in [20] have demonstrated that only the top 2.4-m thickness of the surface layer is substantially affected by crack propagation during multiple F-T cycles. Therefore, the present work considered crack propagation only within this depth zone (Figure 17), whereas the effective shear strengths of the silty clay layer below the 2.4-m level were maintained throughout. Figure 16. Schematic view of the crack propagation in soil slope: (a) CT images of silty clay after various freeze-thaw (F-T) cycles [10]; (b) crack propagation in the surface of soil slope (after [17]); and (c) SWE rainfall intensity (Ir) for soil slope exposed to different F-T cycle.  [17]); and (c) SWE rainfall intensity (I r ) for soil slope exposed to different F-T cycle.
The effective cohesion (c ) and effective internal friction angle (ϕ ) obtained from the laboratory tests are plotted in Figure 18 with respect to the number of F-T cycles (N FT ). We note that both c and ϕ decrease with increasing N FT , indicating that the shear strength of the soil decreases with increasing N FT . The effective cohesion (c′) and effective internal friction angle (φ′) obtained from the laboratory tests are plotted in Figure 18 with respect to the number of F-T cycles (NFT). We note that both c′ and φ′ decrease with increasing NFT, indicating that the shear strength of the soil decreases with increasing NFT. The values of Sr and PWP observed along lines L1 and L2 are plotted in Figure 19 with respect to NFT. The results indicate that Sr and PWP increase significantly with increasing NFT. The values of Sr observed at the top and slope toe at the eleventh F-T cycle are about 23.7% and 17.4% greater than the corresponding values observed in the initial state. Clearly, NFT has a significant effect on fluctuations in the water content, which influences the stability of the cutting slope, as demonstrated by the plot of the FOS with respect to NFT in Figure 20. Compared to the FOS under ds and ΔT, a lower FOS was observed under the effect of F-T cycle. In addition, Figure 19c, d demonstrate that Sr and PWP are significantly increased within the dotted box region of the figures. Such behavior may be due to the application of an Ir value that is close to the value of ks for the silty clay layer, which would induce runoff along the slope surface.
The results obtained here demonstrate that, among the various parameters considered, F-T cycling has a particularly significant impact on fluctuations in the soil water content and in turn heavily influences the stability of cutting slope in seasonally cold regions. Therefore, it is important to consider the F-T cycling in estimating the soil water distributions for the proper pre-judgement of slope stability.  The effective cohesion (c′) and effective internal friction angle (φ′) obtained f laboratory tests are plotted in Figure 18 with respect to the number of F-T cycles (N note that both c′ and φ′ decrease with increasing NFT, indicating that the shear str the soil decreases with increasing NFT.  Figure 20. Compared to the FOS under ds and ΔT, a lower FOS was observe the effect of F-T cycle. In addition, Figure 19c, d demonstrate that Sr and PWP are cantly increased within the dotted box region of the figures. Such behavior may b the application of an Ir value that is close to the value of ks for the silty clay laye would induce runoff along the slope surface.
The results obtained here demonstrate that, among the various parameters ered, F-T cycling has a particularly significant impact on fluctuations in the so content and in turn heavily influences the stability of cutting slope in seasonally gions. Therefore, it is important to consider the F-T cycling in estimating the so distributions for the proper pre-judgement of slope stability. The values of S r and PWP observed along lines L1 and L2 are plotted in Figure 19 with respect to N FT . The results indicate that S r and PWP increase significantly with increasing N FT . The values of S r observed at the top and slope toe at the eleventh F-T cycle are about 23.7% and 17.4% greater than the corresponding values observed in the initial state. Clearly, N FT has a significant effect on fluctuations in the water content, which influences the stability of the cutting slope, as demonstrated by the plot of the FOS with respect to N FT in Figure 20. Compared to the FOS under d s and ∆T, a lower FOS was observed under the effect of F-T cycle. In addition, Figure 19c,d demonstrate that S r and PWP are significantly increased within the dotted box region of the figures. Such behavior may be due to the application of an I r value that is close to the value of k s for the silty clay layer, which would induce runoff along the slope surface.   The results obtained here demonstrate that, among the various parameters considered, F-T cycling has a particularly significant impact on fluctuations in the soil water content and in turn heavily influences the stability of cutting slope in seasonally cold regions. Therefore, it is important to consider the F-T cycling in estimating the soil water distributions for the proper pre-judgement of slope stability.

Conclusions
The present work addressed the currently limited investigation of the impact of snowmelt infiltration on the seepage field and stability of cutting slopes. A 3D FE model was employed to explore the influence of the depth of accumulated snow, daily temperature variation, and the number of F-T cycles on the seepage field and stability. The numerical analysis methodology was validated by comparisons with the analysis results obtained by Oh and Lu [21]. The primary conclusions from the parameter study can be given as follows: (1) Due to the PWP diffusion process, the deep-seated instability is likely to occur under a prolonged and highly accumulated infiltration, while the shallow slope failure is associated with intense, short-duration snowmelt infiltration.
(2) For the given soil slope, snowmelt infiltration primarily affects the seepage field within a region about 10 m below the ground surface. Moreover, the maximum S r and PWP values are observed at the slope toe; this may be the reason why snowmelt-induced soil slope failures usually initiate from the slope toe.
(3) The water content distributions inside the soil slope and along the slope surface are seriously affected by snowmelt infiltration. Furthermore, the depth of accumulated snow, daily temperature variation, and number of F-T cycles have a negative effect on the slope stability. Besides, as these factors increase, the S r and PWP increase considerably. Therefore, these factors should be considered in landslide susceptibility prediction.
(4) The F-T cycle has a substantial impact on the seepage field and shear strength of soils, which in turn affects the soil slope stability. F-T cycle also induces the crack propagation in cutting slope. Compared to the depth of accumulated snow and daily temperature variation, the influence of F-T cycle is more likely to cause slope failure.