# Landslide Accumulation Ice-Snow Melting for Thermo-Hydromechanical Coupling and Numerical Simulation

Authors: Tongqiang Xiong; Jianlin Li; Lehua Wang; Huafeng Deng; Xiaoliang Xu
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

In this paper, we discussed the phase change coupling algorithm of accumulation bank slope under the action of ice-snow melting. We described the effect of temperature gradient on the water migration of soil. We simplified the stress balance, continuity, and energy equation in the coupled model. We discussed the variation law of temperature, seepage, and stress (deformation) field under different conditions of ice-snow melting on the bank slope of the accumulation body. Based on the three-field coupling energy balance equation of ice-snow melting with phase change, the simplified algorithm of three-field coupling is obtained. The simplified algorithm is applied to the coupling model of ice-snow thawing on indoor accumulation bank slope. We established a practical numerical model for the coupling analysis of temperature, seepage, and stress field. We established the coupled control differential equation of three fields. We investigated the three-dimensional numerical simulation of stress, displacement, plastic deformation, and other indicators. The results show that the numerical simulation results are in good agreement with the monitoring results. It is expected that the research results can more truly simulate the actual characteristics of ice and snow melting water on the bank slope of the Three Gorges Reservoir and provide reference for the prevention and prediction of extreme snow and ice disasters in the Three Gorges Reservoir area.


---

## Body

## 1. Introduction

A rapid huge landslide occurred in Yi Gong River in Bomê, Tibet, on April 9, 2000 [1]. The Yi Gong landslide is caused by the melting of snow and ice on Xuefeng Mountains, resulting in the formation of debris flow, which was a huge disaster [2]. Geological disasters caused by melting ice and snow are common. The volume of the Alps is larger than that of the Swiss Alps (108 m3). There are dozens of extremely large slopes [3]. Most of them are caused by rock mass loosening and water saturation caused by the melting of ice and snow during the retreat of Quaternary glaciers [4]. A large number of landslides occurred from February to March 2006 due to the melting of ice and snow [5–8]. The extreme ice-snow disaster climate rarity appears in the Three Gorges Reservoir area, and the probability is very small. The geological disaster caused by it has not been paid enough attention. There is a great difference with the alpine region, in the topography, geological environment, and climate environment. Under extreme snow and ice disaster conditions, the law of landslide disaster triggered by snow melting is similar to that of heavy rainfall [9, 10]. But it is not exactly the same. The time and scale of landslide induced by the two methods are also different [11–13]. The increase of its porosity makes the snow water fully infiltrate and soften the landslide mass and sliding surface. Compared with heavy rainfall [14], the ice-snow melting deformation is shallow, local, and slow [15]. The infiltration process of ice-snow melting is essentially a process of infiltration of water movement in the unsaturated area. The infiltration boundary of ice-snow melting is also a flow boundary [16]. This flow is variable. In the calculation process, it is necessary to adjust the infiltration flow according to the change of water content so as to realize the numerical simulation of ice-snow melting infiltration [17]. At this time, the snow cover only changes temperature and does not change phase [18]. Each time step is based on energy balance and water balance to calculate ice-snow melting, snow depth, and density which need to be recalculated [19–21]. The “coupled cycle” algorithm has the effect of three cycles coupling to temperature, seepage and stress field. It is equivalent to “heat flow-thermal-flow force” in series, and the two fields in each field are paralleled in the substep mode [22, 23]. Some scholars have studied the advance and review on frozen soil mechanics [24, 25]. They are calibration of snow parameters in SWAT: Comparison of three approaches in the Upper Adige River [26, 27]. We study models for polythermal ice sheets [28]. The melting process of ice sheet is the main body research on ice-snow melting [29–31]. They studysee page and stability analysis of slope under the condition of snowmelt infiltration [32]. The change of snow cover is reconstructed, whether the condition of multilayer snow cover is still satisfied [33, 34]. We verify enthalpy schemes for polythermal glaciers and ice sheet dimensional model [35–39]. They research data assimilation and prognostic whole ice sheet modeling with the variation derived, higher order, open source, and fully parallel ice sheet [40–45].By comparison of extreme ice-snow in the Three Gorges Reservoir area, the depth, thermal, and freeze-thaw processes are not the same. The rare ice-snow climate occurred in the Three Gorges Reservoir area is rare, and the geological disaster caused by it has not been paid enough attention. In addition, there are great differences in topography, geological, and weather environment and so on. In the middle-low latitudes of the Three Gorges Reservoir area, the theoretical, experimental, and numerical studies on the types of geological disasters and their influencing mechanisms under extreme ice and snow disasters are basically blank. In this paper, we discussed the phase change coupling algorithm of accumulation bank slope under the action of ice-snow melting. We described the effect of temperature gradient on the water migration of soil. We simplified the stress balance, continuity, and energy equation in the coupled model. We discussed the variation law of temperature, seepage, and stress (deformation) field under different conditions of ice-snow melting on the bank slope of accumulation body. Based on the three-field coupling energy balance equation of ice-snow melting with phase change, the simplified algorithm of three-field coupling is obtained. The simplified algorithm is applied to the coupling model of ice-snow thawing on the indoor accumulation bank slope.


## 2. Analysis of Ice-Snow Melting Infiltration

Under extreme snow and ice disaster conditions, the law of landslide disaster triggered by snow melting is similar to that of heavy rainfall. But it is not exactly the same. The time and scale of landslide induced by the two methods are also different. Firstly, the amount of snow melting cannot reach the level of heavy rainfall. It will not generate a large amount of water and rapid runoff, and no pressure infiltration and surface erosion will be formed. Secondly, because of the temperature difference between day and night, ice-snow melts slowly. This leads to slow infiltration. However, under the action of repeated freezing and thawing in day and night, the soil mass of the bank slope is under the action of repeated freezing and thawing. The increase of its porosity makes the snow water fully infiltrate and soften the landslide mass and sliding surface. Compared with heavy rainfall, the deformation of slope is shallow, local, and slow.The infiltration process of ice-snow melting is essentially a process of infiltration water movement in the unsaturated area. The infiltration boundary of ice-snow melting is also a flow boundary. The flow boundary is variable. In the calculation process, it is necessary to adjust the infiltration flow according to the change of water content. So, we realized the numerical simulation of ice-snow melting infiltration.


### 2.1. Phase Change Equation of Seepage

This paper is based on the bank slope model test under the combined action of indoor reservoir water and ice-snow melting. The interaction of ice-snow on the surface of the bank slope lasted for 15 days. The basic melting is over. Based on the interaction of multiple working conditions (gravity action of 7T ice-snow load and the ice-snow melting acted alone in the early stage and combined action of reservoir water and ice-snow melting), an unsaturated seepage calculation model is established. The distribution function of seepage field of soil is obtained as follows:(1)∂∂xkx∂H∂x+∂∂yky∂H∂y+Q=mwγw∂H∂t,where mw=∂θ/∂ua−uw, H=uw/γw+y, H is the total head, y is the elevation, kx and ky are x, y directional permeability coefficient, Q is boundary input flow, ua–uw is matric suction, γw is the gravity of water, t is time, and mw is the change rate of water content with pore water pressure.The soil water characteristic curve (SWCC) is divided intoN equal parts along the volume water content axis, and the matric suction at the midpoint of each segment was calculated as follows:(2)KθikskscAd∑j=1n2j+1−2iua−uw,where k(θi) is the permeability coefficient corresponding to the volume water content of the segment point; j is from 1 to m count; ua−uwj is j point matric suction value; and Ad is surface tension and water density relevant adjustment factors.


### 2.2. Definite Solution of Partial Differential Equation for Unsaturated Seepage

This paper takes water potential as a dependent variable. The water potential models of free infiltration and pressure infiltration are derived. The semianalytical and semi-iterative algorithm is used to solve the partial differential equation of unsaturated seepage:(3)Ch∂z∂t=−∂∂hKh∂z/∂h−∂Kh∂h,where h is the negative pressure head; Z is the infiltration depth; K is the infiltration coefficient; and C is the water melting rate. The nonlinear equation of unsaturated permeability coefficient K is the negative pressure head.The free infiltration stage of ice-snow melting generally occurs in the early stage of ice-snow melting. The definite solution conditions of free infiltration stage are as follows:(4)Ch∂z∂t=−∂∂hKh∂z/∂h−∂Kh∂h,h=h0,when,t=0,z≥0,kh1+∂h∂z=R,when,t>0,z=0,h=h0,when,t>0,z⟶∞,where R is the ice-snow melting strength.At this time, the pressure head of the surface changes with time, as the samez is the depth below the surface. It is assumed that the initial negative pressure is large. Unsaturated permeability coefficient K(h) is very small. Then, by formula (4),(5)Kh+Kh∂zih,t/∂h=R.In the light of formula (5) from h0 (t) to h(t) integral, we use α representing integral variable, and it can be sorted out as follows:(6)Z1h,t=∫h0thtKαR−Kαdα.In the light of formula (6) knowable, to get the result of the last iteration Z1 (h, t), h0 (t) has to be solved, that is the variation law of the surface negative pressure head with time. We derive formula (6).(7)−∂z1h,t∂t=Kh0tR−Kh0th0′t.We research a nonlinear equation of unsaturated infiltration with pressure head as an unknown number. In this paper, the average value is replaced by the real valueC(h) of melting water degree, and it approximately linearized as follows:(8)C¯=1hi+1−hi∫hihi+1Chdh.Then, we substitute formulas (7) and (8) into (4):(9)Kh0tR−Kh0th0′t=−1C¯∂∂hKh∂zi+1h,t/∂h+Kh.We are solving the integral of equation (9), integral limit from h0 to h. The initial pressure head h0 is assumed to be large. The permeability coefficient tends to 0, as the same K(h0)≈0:(10)Kh0tR−Kh0th0′th−h0=−1C¯∂∂hKh∂zi+1h,t/∂h+Kh.Leth equal to h0(t) be substituted into the above formula. By using the boundary conditions, the variation law of negative pressure head of surface soil with time in the free infiltration stage of ice-snow melting is obtained:(11)h0′t=−1C¯Rh0t−h0R−Kh0tKh0t.We study the one-dimensional unsaturated vertical infiltration model of ice-snow melting. Its basic ideas are as follows: thei times iteration result is Zi (h, t), derivative can get ∂zih,t/∂h, one time integral can get ∂zi+1h,t/∂h, one more integration can get the i + 1 times iteration result Zi+1 (h, t). Continuous iteration, until Zi+1(h, t) obtained by the two iterations is less than the allowable iteration error. According to the above equation:(12)Ch∂zih,t∂t=−∂∂hKh∂zi+1h,t/∂h−∂Kh∂h.In the light of formula (12) integral from h0 (t) to h(t), the initial pressure head h0 is assumed to be large. The permeability coefficient K (h0) tends to 0. It can be sorted out:(13)∂zi+1h,t∂h=−KhKh+∫h0hCh∂zih,t/∂tdh.In the light of formula (12) integral, the integral limit is from h0 to h, In order to avoid confusion of integral limits, In the first integration, the variable symbol h is changed into α, and in the second integration, the variable symbol is changed to β, The results of the known I iterations Zi (h, t) are obtained. The i + 1 iteration result is obtained Zi+1 (h, t).(14)zi+1h,t=∫hh0KβKβ+∫h0βCα∂ziα,t/∂tdαdβ.


### 2.3. Energy Conservation Equation of Ice-Snow Melting Phase Change

The temperature of snow cover on bank slope is usually less than 0°C. When the net energyQ of all soil mass in the surface layer of the bank slope is greater than the temperature of the bottom layer of snow cover, Qs, as same Q > Qs, the underlying snow began to melt. If the amount of melted snow is greater than the water holding capacity of the underlying snow, The melting snow began to flow. When the net energy of surface soil is 0Q < Qs, At this time, the snow cover only changes temperature and does not change phase; Each time step is based on energy balance and water balance to calculate ice-snow melting, snow depth and density need to be recalculated [17].The change of snow cover is reconstructed, Whether the condition of multilayer snow cover is still satisfied. Then the thickness of multilayer snow cover and snow water equivalent are recalculated. The three-field coupled energy balance equation with phase transition is as follows:(15)C3¯∂T∂t+C1¯vs¯∇T+C2¯vw¯∇T=∇λ¯∇T−β¯T∂trε∂t−nw∂pw∂TT∇vw¯+12βθε+Lρi∂nwi∂t,where C¯1, C¯2, and C¯3 are the thermal conductivity and melting coefficient of each component of the soil. ns, ni, and nw are the soil particles, the nonmelting ice-snow, and water volume contents; ρs, ρi, and ρw are the density of soil particles, nonmelting ice-snow, and water; nwi is the volume content of melted snow water; vw¯ is the migration rate of melted snow water; λs, λi, and λw are the heat conduction coefficients of soil particles, the snow water, and nonmelting ice-snow, respectively; Ci, Cs, and Cw are soil particles, thawing ice, and snow water.


## 2.1. Phase Change Equation of Seepage

This paper is based on the bank slope model test under the combined action of indoor reservoir water and ice-snow melting. The interaction of ice-snow on the surface of the bank slope lasted for 15 days. The basic melting is over. Based on the interaction of multiple working conditions (gravity action of 7T ice-snow load and the ice-snow melting acted alone in the early stage and combined action of reservoir water and ice-snow melting), an unsaturated seepage calculation model is established. The distribution function of seepage field of soil is obtained as follows:(1)∂∂xkx∂H∂x+∂∂yky∂H∂y+Q=mwγw∂H∂t,where mw=∂θ/∂ua−uw, H=uw/γw+y, H is the total head, y is the elevation, kx and ky are x, y directional permeability coefficient, Q is boundary input flow, ua–uw is matric suction, γw is the gravity of water, t is time, and mw is the change rate of water content with pore water pressure.The soil water characteristic curve (SWCC) is divided intoN equal parts along the volume water content axis, and the matric suction at the midpoint of each segment was calculated as follows:(2)KθikskscAd∑j=1n2j+1−2iua−uw,where k(θi) is the permeability coefficient corresponding to the volume water content of the segment point; j is from 1 to m count; ua−uwj is j point matric suction value; and Ad is surface tension and water density relevant adjustment factors.


## 2.2. Definite Solution of Partial Differential Equation for Unsaturated Seepage

This paper takes water potential as a dependent variable. The water potential models of free infiltration and pressure infiltration are derived. The semianalytical and semi-iterative algorithm is used to solve the partial differential equation of unsaturated seepage:(3)Ch∂z∂t=−∂∂hKh∂z/∂h−∂Kh∂h,where h is the negative pressure head; Z is the infiltration depth; K is the infiltration coefficient; and C is the water melting rate. The nonlinear equation of unsaturated permeability coefficient K is the negative pressure head.The free infiltration stage of ice-snow melting generally occurs in the early stage of ice-snow melting. The definite solution conditions of free infiltration stage are as follows:(4)Ch∂z∂t=−∂∂hKh∂z/∂h−∂Kh∂h,h=h0,when,t=0,z≥0,kh1+∂h∂z=R,when,t>0,z=0,h=h0,when,t>0,z⟶∞,where R is the ice-snow melting strength.At this time, the pressure head of the surface changes with time, as the samez is the depth below the surface. It is assumed that the initial negative pressure is large. Unsaturated permeability coefficient K(h) is very small. Then, by formula (4),(5)Kh+Kh∂zih,t/∂h=R.In the light of formula (5) from h0 (t) to h(t) integral, we use α representing integral variable, and it can be sorted out as follows:(6)Z1h,t=∫h0thtKαR−Kαdα.In the light of formula (6) knowable, to get the result of the last iteration Z1 (h, t), h0 (t) has to be solved, that is the variation law of the surface negative pressure head with time. We derive formula (6).(7)−∂z1h,t∂t=Kh0tR−Kh0th0′t.We research a nonlinear equation of unsaturated infiltration with pressure head as an unknown number. In this paper, the average value is replaced by the real valueC(h) of melting water degree, and it approximately linearized as follows:(8)C¯=1hi+1−hi∫hihi+1Chdh.Then, we substitute formulas (7) and (8) into (4):(9)Kh0tR−Kh0th0′t=−1C¯∂∂hKh∂zi+1h,t/∂h+Kh.We are solving the integral of equation (9), integral limit from h0 to h. The initial pressure head h0 is assumed to be large. The permeability coefficient tends to 0, as the same K(h0)≈0:(10)Kh0tR−Kh0th0′th−h0=−1C¯∂∂hKh∂zi+1h,t/∂h+Kh.Leth equal to h0(t) be substituted into the above formula. By using the boundary conditions, the variation law of negative pressure head of surface soil with time in the free infiltration stage of ice-snow melting is obtained:(11)h0′t=−1C¯Rh0t−h0R−Kh0tKh0t.We study the one-dimensional unsaturated vertical infiltration model of ice-snow melting. Its basic ideas are as follows: thei times iteration result is Zi (h, t), derivative can get ∂zih,t/∂h, one time integral can get ∂zi+1h,t/∂h, one more integration can get the i + 1 times iteration result Zi+1 (h, t). Continuous iteration, until Zi+1(h, t) obtained by the two iterations is less than the allowable iteration error. According to the above equation:(12)Ch∂zih,t∂t=−∂∂hKh∂zi+1h,t/∂h−∂Kh∂h.In the light of formula (12) integral from h0 (t) to h(t), the initial pressure head h0 is assumed to be large. The permeability coefficient K (h0) tends to 0. It can be sorted out:(13)∂zi+1h,t∂h=−KhKh+∫h0hCh∂zih,t/∂tdh.In the light of formula (12) integral, the integral limit is from h0 to h, In order to avoid confusion of integral limits, In the first integration, the variable symbol h is changed into α, and in the second integration, the variable symbol is changed to β, The results of the known I iterations Zi (h, t) are obtained. The i + 1 iteration result is obtained Zi+1 (h, t).(14)zi+1h,t=∫hh0KβKβ+∫h0βCα∂ziα,t/∂tdαdβ.


## 2.3. Energy Conservation Equation of Ice-Snow Melting Phase Change

The temperature of snow cover on bank slope is usually less than 0°C. When the net energyQ of all soil mass in the surface layer of the bank slope is greater than the temperature of the bottom layer of snow cover, Qs, as same Q > Qs, the underlying snow began to melt. If the amount of melted snow is greater than the water holding capacity of the underlying snow, The melting snow began to flow. When the net energy of surface soil is 0Q < Qs, At this time, the snow cover only changes temperature and does not change phase; Each time step is based on energy balance and water balance to calculate ice-snow melting, snow depth and density need to be recalculated [17].The change of snow cover is reconstructed, Whether the condition of multilayer snow cover is still satisfied. Then the thickness of multilayer snow cover and snow water equivalent are recalculated. The three-field coupled energy balance equation with phase transition is as follows:(15)C3¯∂T∂t+C1¯vs¯∇T+C2¯vw¯∇T=∇λ¯∇T−β¯T∂trε∂t−nw∂pw∂TT∇vw¯+12βθε+Lρi∂nwi∂t,where C¯1, C¯2, and C¯3 are the thermal conductivity and melting coefficient of each component of the soil. ns, ni, and nw are the soil particles, the nonmelting ice-snow, and water volume contents; ρs, ρi, and ρw are the density of soil particles, nonmelting ice-snow, and water; nwi is the volume content of melted snow water; vw¯ is the migration rate of melted snow water; λs, λi, and λw are the heat conduction coefficients of soil particles, the snow water, and nonmelting ice-snow, respectively; Ci, Cs, and Cw are soil particles, thawing ice, and snow water.


## 3. Coupled Mathematical Model of Soil with Melt Water

The THM coupling is temperature seepage stress coupling or thermal fluid solid coupling. Meltwater soil is a kind of complex soil mass encountered in reservoir bank slope engineering under extreme ice-snow disaster environment. It is a discontinuous medium composed of snow water with normal distribution and soil eroded by seepage due to temperature rise. The coupling of temperature field, seepage field, and stress field is a relatively complex problem. It is mainly manifested in the various components of geological environment of bank slope soil, that is, seepage field, stress field, and temperature field change with time and space. The coupling between the above components is a complex dynamic process. The traditional discrete medium hypothesis (neglecting the permeability of soil pore system) and the continuous medium hypothesis (averaging the melted snow water seepage into the soil seepage) are often difficult to reflect the coupling effect of temperature field seepage field stress field. To some extent, the three-field coupling model can simulate the seepage, heat transfer, and deformation of bank soil. According to the basic and physical property equation [17], the coupled control equations of temperature field, seepage field, and stress field can be derived. In recent years, many coupling models are proposed, such as simple heat conduction model, water flow model, water heat coupling model, and heat fluid solid coupling model. Considering the properties of rock mass, there are homogeneous model, fracture medium model, and so on. By comparison, the more comprehensive model is the thermal fluid solid coupling mathematical model which is proposed. We proposed and improved the coupling finite element method of temperature field, seepage field, and stress field of fractured rock mass.


### 3.1. Couple Theory of Seepage Field and Temperature Field

When one-dimensional heat transfer is considered, moisture transfer and wetting peak in soil change. There are two development directions of temperature field variation, that is, the heat conduction phenomenon of soil particles themselves and the heat carried by soil particles during water migration.The variation of one-dimensional temperature field is as follows:(16)qx=CwρwνxT−λx∂T∂x,where qx, λx, and νx are the latent heat flux density, thermal conductivity, and penetration velocity in X direction; Cw is the amount of heat fusion of water; ρw is the density of water; and T is the function of soil temperature.On the basis of equation (16), the three-dimensional temperature field equation is derived by adding temperature heat flux Qv,(17)∂∂xλx∂T∂x+∂∂xλy∂T∂y+∂∂xλz∂T∂z−Cwρw∂νxT∂x+∂νyT∂y+∂νzT∂z+qv=Cρ∂T∂t.Equation (17) expresses the interaction between seepage field and temperature field. With the increase of environmental temperature, the heat transfer and migration of soil mass and the redistribution of temperature field occur. A theoretical solution of seepage (or temperature) based on the boundary of seepage (or temperature) and satisfying the differential equation of motion is established.


### 3.2. Coupling Effect of Temperature Field and Seepage Field on Stress Field

The linear momentum balance equation can be derived from the principle of linear momentum balance:(18)∇⋅σij+ρfi=0ρDVsDt,where σij is the component of the total stress tensor; ρ is the equivalent density of bank slope medium (soil and groundwater), kg/m3; and fi is the volume force component of bank slope medium. There is the effective stress formula as follows:(19)σij=σij′+αpδij.Substituting (19) into (18), the static equilibrium equation is expressed by the effective stress:(20)∇⋅σij′+αpδijj+ρfi=0.From the general constitutive equation of thermoelasticity,(21)σij′=Deεij−εijT=Deεij−βsTs−Ts0δij,where De is the elastic matrix of bank slope soil; Ts and Ts0 are the temperature of bank slope and reference temperature; and βs is the linear thermal expansion coefficient of bank slope.The thermoelastic equilibrium equation can be obtained by substituting (21) into (20):(22)∇⋅Deεij−βsTs−Ts0δij+αpδij+ρfi=0.For the case of isotropic thermal expansion and mechanical properties of bank slope skeleton, the effective stress-strain relationship of bank slope considering the effect of temperature change can be written as follows:(23)σij′=2Gεij+λδijεkk−βTs−Ts0δij,where equation (23) is the component of effective stress tensor; εij is the component of the strain tensor; λ is the elastic constant; and β = (2G + 3λ) is the thermal stress coefficient.The isotropic thermoelastic equilibrium equation can be obtained by substituting equation (23) into (22):(24)∇⋅2Gεij+λδijεkk−βTs−Ts0δij+αpδij+ρfi=0.The incremental form of the constitutive equation is as follows:(25)dσij′=Depdεij−βsdTsδij,where Dep is the elastoplastic matrix of bank slope material.By substituting equation (24) into (23), an incremental thermoelastic-plastic equilibrium equation can be obtained:(26)∇⋅Depdεij−βsdTsδij+αpδij+ρfi=0.Equations (24)–(26) are the governing equations of deformation field of bank slope soil under the joint action of temperature field and seepage field. The deformation states of isotropic thermoelasticity are corresponding respectively. The last two items on the left of equations (24)–(26), respectively, reflect the influence of thermal stress and pore water pressure on the static equilibrium equation of bank slope. Therefore, the three static equilibrium equations under different deformation states fully embody and describe the coupling effect of temperature field and seepage field on stress field.


### 3.3. Seepage Field Model Coupled with Temperature Field and Stress Field

The influence of water head gradient and temperature gradient is also considered. The seepage flow flux (seepage velocity) consists of two parts: one part is the flow movement caused by head gradient, and the other part is the flow movement caused by temperature gradient. In the case of one-dimensional seepage, we study the expression of seepage velocity:(27)vx=−KdHdx−DTdTdx,where vx is the seepage velocity along the X direction of one-dimensional coordinate axis; H is the seepage head; T is the temperature; DT is the diffusion coefficient of water flow under the action of temperature difference; and K is the permeability coefficient and it is also a function of stress.(28)K=K0Te−ασ−γH,where K is the permeability coefficient of pressure σ; K0 is the permeability coefficient of pressure, and it is also a function of temperature; γH is the osmotic hydrostatic pressure; and α is a constant, which is determined by experiment. Equations (27) and (28) are substituted into the continuity equation of one-dimensional seepage. The following formula can be obtained:(29)ddxK0Te−ασ−γHdHdx+DTd2Tdx2=0.Equation (29) is the basic one-dimensional seepage equation considering both temperature and stress. The extension to 3D is as follows:(30)ddxK0xTe−ασ−γHdHdx+ddyK0yTe−ασ−γHdHdy+ddzK0zTe−ασ−γHdHdz+DTd2Tdx2+d2Tdy2+d2Tdz2=0.Equation (30) is the governing equation of three-dimensional seepage field coupled with temperature and stress, combined with certain definite solution conditions, and a three-dimensional seepage field model is formed.


### 3.4. Temperature Field Model Coupled with Seepage Field and Stress Field

Considering the change of soil thermodynamic properties caused by bank slope deformation, the thermal convection of seepage flow, and the heat exchange with the soil mass of bank slope, the governing equation of temperature field coupled with seepage and stress is as follows:(31)∂∂xλxσ∂T∂x+∂∂yλyσ∂T∂y+∂∂zλzσ∂T∂z−cwσρwσ∂vxHT∂x+∂vyHT∂y+∂vzHT∂z+QT=cwσρwσ∂T∂t,where c, ρ, and λ are the specific heat, density, and thermal conductivity of bank slope soil; cwandρw are the specific heat and density of water, which are functions of stress; vx, vy, and vz are the three components of seepage velocity, which are the function of water head; and QT is the source (sink) term of heat.Equation (31) combines with certain definite solution conditions to form a three-dimensional temperature field model.In the three-field coupling model shown in equations (26), (30), and (31), the relationship among seepage, stress, and temperature is “strong coupling.” Three-field iterative methods are used to solve the problem. In other words, the coupling relationship among seepage, stress, and temperature is regarded as “bridge,” the iterative calculation of seepage field, stress field, and temperature field is carried out, respectively (for example, firstly, the distribution of stress and temperature fields is assumed according to the definite solution conditions of stress and temperature fields; secondly, the stress field distribution is solved according to equation (30) from the distribution of seepage field and assumed temperature field; thirdly, the temperature field distribution is solved according to equation (31); so, we iterate over and over again), until convergence (that is, between two adjacent iteration cycles, the maximum errors of the distribution of seepage field, stress field and temperature field are all within the controllable accuracy). The solution of seepage field, stress field, and temperature field under the coupling effect of three fields can be obtained. It is more objective to reflect the essential attribute of the bank slope of ice-snow melting water accumulation body.


## 3.1. Couple Theory of Seepage Field and Temperature Field

When one-dimensional heat transfer is considered, moisture transfer and wetting peak in soil change. There are two development directions of temperature field variation, that is, the heat conduction phenomenon of soil particles themselves and the heat carried by soil particles during water migration.The variation of one-dimensional temperature field is as follows:(16)qx=CwρwνxT−λx∂T∂x,where qx, λx, and νx are the latent heat flux density, thermal conductivity, and penetration velocity in X direction; Cw is the amount of heat fusion of water; ρw is the density of water; and T is the function of soil temperature.On the basis of equation (16), the three-dimensional temperature field equation is derived by adding temperature heat flux Qv,(17)∂∂xλx∂T∂x+∂∂xλy∂T∂y+∂∂xλz∂T∂z−Cwρw∂νxT∂x+∂νyT∂y+∂νzT∂z+qv=Cρ∂T∂t.Equation (17) expresses the interaction between seepage field and temperature field. With the increase of environmental temperature, the heat transfer and migration of soil mass and the redistribution of temperature field occur. A theoretical solution of seepage (or temperature) based on the boundary of seepage (or temperature) and satisfying the differential equation of motion is established.


## 3.2. Coupling Effect of Temperature Field and Seepage Field on Stress Field

The linear momentum balance equation can be derived from the principle of linear momentum balance:(18)∇⋅σij+ρfi=0ρDVsDt,where σij is the component of the total stress tensor; ρ is the equivalent density of bank slope medium (soil and groundwater), kg/m3; and fi is the volume force component of bank slope medium. There is the effective stress formula as follows:(19)σij=σij′+αpδij.Substituting (19) into (18), the static equilibrium equation is expressed by the effective stress:(20)∇⋅σij′+αpδijj+ρfi=0.From the general constitutive equation of thermoelasticity,(21)σij′=Deεij−εijT=Deεij−βsTs−Ts0δij,where De is the elastic matrix of bank slope soil; Ts and Ts0 are the temperature of bank slope and reference temperature; and βs is the linear thermal expansion coefficient of bank slope.The thermoelastic equilibrium equation can be obtained by substituting (21) into (20):(22)∇⋅Deεij−βsTs−Ts0δij+αpδij+ρfi=0.For the case of isotropic thermal expansion and mechanical properties of bank slope skeleton, the effective stress-strain relationship of bank slope considering the effect of temperature change can be written as follows:(23)σij′=2Gεij+λδijεkk−βTs−Ts0δij,where equation (23) is the component of effective stress tensor; εij is the component of the strain tensor; λ is the elastic constant; and β = (2G + 3λ) is the thermal stress coefficient.The isotropic thermoelastic equilibrium equation can be obtained by substituting equation (23) into (22):(24)∇⋅2Gεij+λδijεkk−βTs−Ts0δij+αpδij+ρfi=0.The incremental form of the constitutive equation is as follows:(25)dσij′=Depdεij−βsdTsδij,where Dep is the elastoplastic matrix of bank slope material.By substituting equation (24) into (23), an incremental thermoelastic-plastic equilibrium equation can be obtained:(26)∇⋅Depdεij−βsdTsδij+αpδij+ρfi=0.Equations (24)–(26) are the governing equations of deformation field of bank slope soil under the joint action of temperature field and seepage field. The deformation states of isotropic thermoelasticity are corresponding respectively. The last two items on the left of equations (24)–(26), respectively, reflect the influence of thermal stress and pore water pressure on the static equilibrium equation of bank slope. Therefore, the three static equilibrium equations under different deformation states fully embody and describe the coupling effect of temperature field and seepage field on stress field.


## 3.3. Seepage Field Model Coupled with Temperature Field and Stress Field

The influence of water head gradient and temperature gradient is also considered. The seepage flow flux (seepage velocity) consists of two parts: one part is the flow movement caused by head gradient, and the other part is the flow movement caused by temperature gradient. In the case of one-dimensional seepage, we study the expression of seepage velocity:(27)vx=−KdHdx−DTdTdx,where vx is the seepage velocity along the X direction of one-dimensional coordinate axis; H is the seepage head; T is the temperature; DT is the diffusion coefficient of water flow under the action of temperature difference; and K is the permeability coefficient and it is also a function of stress.(28)K=K0Te−ασ−γH,where K is the permeability coefficient of pressure σ; K0 is the permeability coefficient of pressure, and it is also a function of temperature; γH is the osmotic hydrostatic pressure; and α is a constant, which is determined by experiment. Equations (27) and (28) are substituted into the continuity equation of one-dimensional seepage. The following formula can be obtained:(29)ddxK0Te−ασ−γHdHdx+DTd2Tdx2=0.Equation (29) is the basic one-dimensional seepage equation considering both temperature and stress. The extension to 3D is as follows:(30)ddxK0xTe−ασ−γHdHdx+ddyK0yTe−ασ−γHdHdy+ddzK0zTe−ασ−γHdHdz+DTd2Tdx2+d2Tdy2+d2Tdz2=0.Equation (30) is the governing equation of three-dimensional seepage field coupled with temperature and stress, combined with certain definite solution conditions, and a three-dimensional seepage field model is formed.


## 3.4. Temperature Field Model Coupled with Seepage Field and Stress Field

Considering the change of soil thermodynamic properties caused by bank slope deformation, the thermal convection of seepage flow, and the heat exchange with the soil mass of bank slope, the governing equation of temperature field coupled with seepage and stress is as follows:(31)∂∂xλxσ∂T∂x+∂∂yλyσ∂T∂y+∂∂zλzσ∂T∂z−cwσρwσ∂vxHT∂x+∂vyHT∂y+∂vzHT∂z+QT=cwσρwσ∂T∂t,where c, ρ, and λ are the specific heat, density, and thermal conductivity of bank slope soil; cwandρw are the specific heat and density of water, which are functions of stress; vx, vy, and vz are the three components of seepage velocity, which are the function of water head; and QT is the source (sink) term of heat.Equation (31) combines with certain definite solution conditions to form a three-dimensional temperature field model.In the three-field coupling model shown in equations (26), (30), and (31), the relationship among seepage, stress, and temperature is “strong coupling.” Three-field iterative methods are used to solve the problem. In other words, the coupling relationship among seepage, stress, and temperature is regarded as “bridge,” the iterative calculation of seepage field, stress field, and temperature field is carried out, respectively (for example, firstly, the distribution of stress and temperature fields is assumed according to the definite solution conditions of stress and temperature fields; secondly, the stress field distribution is solved according to equation (30) from the distribution of seepage field and assumed temperature field; thirdly, the temperature field distribution is solved according to equation (31); so, we iterate over and over again), until convergence (that is, between two adjacent iteration cycles, the maximum errors of the distribution of seepage field, stress field and temperature field are all within the controllable accuracy). The solution of seepage field, stress field, and temperature field under the coupling effect of three fields can be obtained. It is more objective to reflect the essential attribute of the bank slope of ice-snow melting water accumulation body.


## 4. Simplified Model of Ice-Snow Melting Coupling

In the process of melting ice and snow, when heat is always transferred from the high temperature zone to the low temperature zone, temperature rise of ice-water phase changes. In addition to the influence of external environment temperature, in the range of ice and snow melting temperature, it is also subject to the latent heat of ice-snow melting phase change. The system temperature is higher than 0°C in the process of ice-snow melting phase change. Until the latent heat is absorbed, the ice and snow melt into water.


### 4.1. Ice-Snow Melting Water Energy Change Stage

In the ice-snow melting stage of bank slope accumulation, it is assumed that the relationship of energy change in each stage is as follows:(i)
Warming stage before melting water phase change(ii)
Latent heat release before phase transition of ice-snow melting(iii)
Phase change and latent heat storage of ice-snow meltingIn other words, the heat of storing latent heat offsets the temperature rise of the system caused by the heat released from the external environment(iv)
The temperature rising stage after the phase change of ice-snow melting


### 4.2. Melting Water Phase Change of Ice and Snow

The basic assumptions of the deformation model of melt water soil include four parts:(1)
The linear elastic deformation or viscoelastic plastic deformation of the soil due to external force(2)
The volume expansion caused by the transformation of unfused ice and snow (including migration water) into water(3)
The deformation caused by the density change of soil particles and water caused by the sealed gas(4)
The driving force model of snow water is snow water seepage (erosion) modelThe flow chart of simplified algorithm of ice-snow melting phase transformation is shown in Figure1. The specific ideas are as follows.Figure 1
The simplified flow chart of ice-water phase change.The bottom temperature is higher than freezing temperature (T0) of snow cover. The model is in the state of “melting” detection of ice-snow melting phase change. At this time, the program traverses the unit temperature, when the temperature of a unit (T1) is higher than the freezing temperature (T0). Within Δt time, each step of calculation resets the unit temperature to the freezing temperature (T1 = T0) when the accumulated temperature difference is generated. ∑ΔT=∑T1−T0 corresponding accumulated heat ∑Q1=∑ΔT∗Cv. When it is not less than the latent heat of the unit, the phase transformation of ice-snow melting is completed. If the continuity of bidirectional phase transition is taken into account, the accumulated heat offsets each other (∑Q − ∑Q1 > 0), and the phase transformation of ice-snow melting is completed.


### 4.3. Coupled Cycle Simplification Algorithm

The temperature field, seepage field, and stress field change at the same time. The system is in equilibrium at all times. Therefore, in order to ensure the real-time synchronization of temperature field and seepage field, the time substep must be the same or multiple. Ensure that the total amount of time in each “coupling cycle” is equal. The corresponding stress field is also balanced in the process of synchronous change. The flow chart of coupling cycle algorithm is shown in Figure2. The core of the simplified phase change algorithm is that the accumulated heat corresponding to the unit temperature difference before and after the coupling calculation is not less than the latent heat of phase change. Before phase change detection and “coupling cycle,” temperature detection must be carried out. After the phase transformation, a temperature calculation is needed. It is used for phase change temperature update. Suppose 60 s is a large cycle, the “coupling cycle” time is 10 s, and the time substep is set 1 s (10 steps of calculation). That is, in the three-field coupling calculation, all three fields are in balance, and 60 s is a cycle calculation cycle.Figure 2
Chart of coupling cycle.The “coupled cycle” algorithm contains temperature field, seepage field, and stress field coupling in series and parallel connection, double play three times. It is equivalent to “thermal flow force” in series, and the two fields in each field are paralleled in substep mode.Firstly, the feasibility of the phase transformation algorithm and the three-field coupling algorithm of the melting water and soil is verified by the ice-snow melting. Secondly, the experiment was regressed to explore the different convective heat transfer coefficient, melting temperature, ambient temperature, porosity, and deformation law. Lastly, the seepage coefficient and erosion coefficient are introduced to deduce the deformation law under the environment of gradual melting and thawing settlement.


### 4.4. Numerical Analysis Method of Coupled Process

The stress balance, seepage, and phase change heat transfer of the model change with time (it is instantaneous effect), and the solution of the model is discretized in time domain. The energy change of convection is considered. It makes the coupling term of seepage and temperature asymmetric in the equivalent matrix of coupling model, considering the dissipation of thermal resistance of soil with water melting. Therefore, the effect of temperature on soil skeleton in the equilibrium equation is asymmetric to that of soil skeleton in energy equation. The interaction between seepage and stress is also asymmetric. The corresponding load vector is changed to make the equation symmetrical. In this paper, the cross-iteration method is used to solve the temperature field and stress field, respectively.


#### 4.4.1. Three-Field Coupling Relation

(i)
Influence of temperature field on the seepage fieldIt can be seen that the influence of temperature change on the driving force of water migration in seepage field includes two aspects. On the one hand, the change of temperature directly causes the change of water migration dynamics. On the other hand, the change of temperature first causes the change of the amount of ice and snow melt, and then the water migration power is changed.(ii)
Influence of seepage field on temperature fieldIn the influence of water migration on temperature field, the heat convection between water and soil skeleton caused by water migration is considered in the seepage equation (see the third item on the left of equation (15)).(iii)
Influence of stress field (deformation field) on temperature field.Changing temperature field caused by soil deformation of bank slope, it is mainly embodied by energy transformation. It includes the effect of thermal expansion and contraction and the thermal change caused by temperature change rate (see the second and fourth items on the right of equation (15).


#### 4.4.2. Numerical Analysis Process

The finite element program for solving the three-field coupling problem is very complicated. The finite element mesh and calculation time step must be determined according to the actual process, initial conditions, and boundary conditions. The flow chart for the transfer is shown in Figure3. The specific steps of this paper are as follows:(i)
Solving the Laplace equation according to boundary conditions of temperature field, the initial temperature field is obtained.(ii)
The calculation parameters are determined according to the initial temperature field (including water, heat, and force parameters). We judged the melting zone and the unfused zone.(iii)
According to the initial temperature field and boundary conditions, the strain energy is assumed to be 0. The temperature field of the next moment is obtained by solving the equation.(iv)
According to the water transport driving force model, calculate the water driving force. The water migration caused by water driving force is calculated, plus the amount of ice and snow melting in this time step. The additional deformation was caused by water migration and ice-water phase transition. The equivalent joint force is determined by the additional deformation.(v)
According to the boundary conditions, additional deformation of the equivalent nodal force and the modified calculation parameters occurred. We obtained the pore water pressure field. We are according to the total deformation field and the additional deformation field (Ua, Va, and εa) to solve the stress field in the soil.(vi)
According to this total deformation field, water component field and stress obtained field determined viscous dissipation energy. We coupled the term of thermal dissipation, heat transfer term, and latent heat of phase change. We obtained a new temperature field.(vii)
Repeat steps 4 to 6, loop through until stable.(viii)
The calculation is made step by step, until the time needed. We obtained the temperature field, water content field, and stress field of the whole process.Figure 3
The flow chart for the transfer.


## 4.1. Ice-Snow Melting Water Energy Change Stage

In the ice-snow melting stage of bank slope accumulation, it is assumed that the relationship of energy change in each stage is as follows:(i)
Warming stage before melting water phase change(ii)
Latent heat release before phase transition of ice-snow melting(iii)
Phase change and latent heat storage of ice-snow meltingIn other words, the heat of storing latent heat offsets the temperature rise of the system caused by the heat released from the external environment(iv)
The temperature rising stage after the phase change of ice-snow melting


## 4.2. Melting Water Phase Change of Ice and Snow

The basic assumptions of the deformation model of melt water soil include four parts:(1)
The linear elastic deformation or viscoelastic plastic deformation of the soil due to external force(2)
The volume expansion caused by the transformation of unfused ice and snow (including migration water) into water(3)
The deformation caused by the density change of soil particles and water caused by the sealed gas(4)
The driving force model of snow water is snow water seepage (erosion) modelThe flow chart of simplified algorithm of ice-snow melting phase transformation is shown in Figure1. The specific ideas are as follows.Figure 1
The simplified flow chart of ice-water phase change.The bottom temperature is higher than freezing temperature (T0) of snow cover. The model is in the state of “melting” detection of ice-snow melting phase change. At this time, the program traverses the unit temperature, when the temperature of a unit (T1) is higher than the freezing temperature (T0). Within Δt time, each step of calculation resets the unit temperature to the freezing temperature (T1 = T0) when the accumulated temperature difference is generated. ∑ΔT=∑T1−T0 corresponding accumulated heat ∑Q1=∑ΔT∗Cv. When it is not less than the latent heat of the unit, the phase transformation of ice-snow melting is completed. If the continuity of bidirectional phase transition is taken into account, the accumulated heat offsets each other (∑Q − ∑Q1 > 0), and the phase transformation of ice-snow melting is completed.


## 4.3. Coupled Cycle Simplification Algorithm

The temperature field, seepage field, and stress field change at the same time. The system is in equilibrium at all times. Therefore, in order to ensure the real-time synchronization of temperature field and seepage field, the time substep must be the same or multiple. Ensure that the total amount of time in each “coupling cycle” is equal. The corresponding stress field is also balanced in the process of synchronous change. The flow chart of coupling cycle algorithm is shown in Figure2. The core of the simplified phase change algorithm is that the accumulated heat corresponding to the unit temperature difference before and after the coupling calculation is not less than the latent heat of phase change. Before phase change detection and “coupling cycle,” temperature detection must be carried out. After the phase transformation, a temperature calculation is needed. It is used for phase change temperature update. Suppose 60 s is a large cycle, the “coupling cycle” time is 10 s, and the time substep is set 1 s (10 steps of calculation). That is, in the three-field coupling calculation, all three fields are in balance, and 60 s is a cycle calculation cycle.Figure 2
Chart of coupling cycle.The “coupled cycle” algorithm contains temperature field, seepage field, and stress field coupling in series and parallel connection, double play three times. It is equivalent to “thermal flow force” in series, and the two fields in each field are paralleled in substep mode.Firstly, the feasibility of the phase transformation algorithm and the three-field coupling algorithm of the melting water and soil is verified by the ice-snow melting. Secondly, the experiment was regressed to explore the different convective heat transfer coefficient, melting temperature, ambient temperature, porosity, and deformation law. Lastly, the seepage coefficient and erosion coefficient are introduced to deduce the deformation law under the environment of gradual melting and thawing settlement.


## 4.4. Numerical Analysis Method of Coupled Process

The stress balance, seepage, and phase change heat transfer of the model change with time (it is instantaneous effect), and the solution of the model is discretized in time domain. The energy change of convection is considered. It makes the coupling term of seepage and temperature asymmetric in the equivalent matrix of coupling model, considering the dissipation of thermal resistance of soil with water melting. Therefore, the effect of temperature on soil skeleton in the equilibrium equation is asymmetric to that of soil skeleton in energy equation. The interaction between seepage and stress is also asymmetric. The corresponding load vector is changed to make the equation symmetrical. In this paper, the cross-iteration method is used to solve the temperature field and stress field, respectively.


### 4.4.1. Three-Field Coupling Relation

(i)
Influence of temperature field on the seepage fieldIt can be seen that the influence of temperature change on the driving force of water migration in seepage field includes two aspects. On the one hand, the change of temperature directly causes the change of water migration dynamics. On the other hand, the change of temperature first causes the change of the amount of ice and snow melt, and then the water migration power is changed.(ii)
Influence of seepage field on temperature fieldIn the influence of water migration on temperature field, the heat convection between water and soil skeleton caused by water migration is considered in the seepage equation (see the third item on the left of equation (15)).(iii)
Influence of stress field (deformation field) on temperature field.Changing temperature field caused by soil deformation of bank slope, it is mainly embodied by energy transformation. It includes the effect of thermal expansion and contraction and the thermal change caused by temperature change rate (see the second and fourth items on the right of equation (15).


### 4.4.2. Numerical Analysis Process

The finite element program for solving the three-field coupling problem is very complicated. The finite element mesh and calculation time step must be determined according to the actual process, initial conditions, and boundary conditions. The flow chart for the transfer is shown in Figure3. The specific steps of this paper are as follows:(i)
Solving the Laplace equation according to boundary conditions of temperature field, the initial temperature field is obtained.(ii)
The calculation parameters are determined according to the initial temperature field (including water, heat, and force parameters). We judged the melting zone and the unfused zone.(iii)
According to the initial temperature field and boundary conditions, the strain energy is assumed to be 0. The temperature field of the next moment is obtained by solving the equation.(iv)
According to the water transport driving force model, calculate the water driving force. The water migration caused by water driving force is calculated, plus the amount of ice and snow melting in this time step. The additional deformation was caused by water migration and ice-water phase transition. The equivalent joint force is determined by the additional deformation.(v)
According to the boundary conditions, additional deformation of the equivalent nodal force and the modified calculation parameters occurred. We obtained the pore water pressure field. We are according to the total deformation field and the additional deformation field (Ua, Va, and εa) to solve the stress field in the soil.(vi)
According to this total deformation field, water component field and stress obtained field determined viscous dissipation energy. We coupled the term of thermal dissipation, heat transfer term, and latent heat of phase change. We obtained a new temperature field.(vii)
Repeat steps 4 to 6, loop through until stable.(viii)
The calculation is made step by step, until the time needed. We obtained the temperature field, water content field, and stress field of the whole process.Figure 3
The flow chart for the transfer.


## 4.4.1. Three-Field Coupling Relation

(i)
Influence of temperature field on the seepage fieldIt can be seen that the influence of temperature change on the driving force of water migration in seepage field includes two aspects. On the one hand, the change of temperature directly causes the change of water migration dynamics. On the other hand, the change of temperature first causes the change of the amount of ice and snow melt, and then the water migration power is changed.(ii)
Influence of seepage field on temperature fieldIn the influence of water migration on temperature field, the heat convection between water and soil skeleton caused by water migration is considered in the seepage equation (see the third item on the left of equation (15)).(iii)
Influence of stress field (deformation field) on temperature field.Changing temperature field caused by soil deformation of bank slope, it is mainly embodied by energy transformation. It includes the effect of thermal expansion and contraction and the thermal change caused by temperature change rate (see the second and fourth items on the right of equation (15).


## 4.4.2. Numerical Analysis Process

The finite element program for solving the three-field coupling problem is very complicated. The finite element mesh and calculation time step must be determined according to the actual process, initial conditions, and boundary conditions. The flow chart for the transfer is shown in Figure3. The specific steps of this paper are as follows:(i)
Solving the Laplace equation according to boundary conditions of temperature field, the initial temperature field is obtained.(ii)
The calculation parameters are determined according to the initial temperature field (including water, heat, and force parameters). We judged the melting zone and the unfused zone.(iii)
According to the initial temperature field and boundary conditions, the strain energy is assumed to be 0. The temperature field of the next moment is obtained by solving the equation.(iv)
According to the water transport driving force model, calculate the water driving force. The water migration caused by water driving force is calculated, plus the amount of ice and snow melting in this time step. The additional deformation was caused by water migration and ice-water phase transition. The equivalent joint force is determined by the additional deformation.(v)
According to the boundary conditions, additional deformation of the equivalent nodal force and the modified calculation parameters occurred. We obtained the pore water pressure field. We are according to the total deformation field and the additional deformation field (Ua, Va, and εa) to solve the stress field in the soil.(vi)
According to this total deformation field, water component field and stress obtained field determined viscous dissipation energy. We coupled the term of thermal dissipation, heat transfer term, and latent heat of phase change. We obtained a new temperature field.(vii)
Repeat steps 4 to 6, loop through until stable.(viii)
The calculation is made step by step, until the time needed. We obtained the temperature field, water content field, and stress field of the whole process.Figure 3
The flow chart for the transfer.


## 5. Three-Dimensional Finite Element Numerical Calculation of Landslide

### 5.1. Calculation Model, Boundary Conditions, and Calculation Parameters

The physical and mechanical parameters of rock and soil are determined. The physical and mechanical parameters of soil mainly include the permeability coefficient (k), the heavy density (γ), Poisson’s ratio (μ), the cohesion (c), the internal friction angle (φ), and so on. We determine main physical and mechanical parameters. The range of values is shown in Table 1.Table 1
Basic property indexes of soils.
SymbolMeaningQuantityγ (kg/m3)Heavy density21.4k (m/s)Permeability coefficient2.3 ∗ 10−6C (kPa)Cohesion17.2φ (°)Internal friction angle18.6Based on the physical and mechanical parameters of rock and soil, test data and recommended parameters are determined. The soil water parameters of temperature field and seepage field are obtained as follows in Table2.Table 2
The thermal-hydraulic parameter of soil and water.
SymbolMeaningQuantityCv/(J·kg−1·°C−1)Specific heat16λ/(W·m−1·°C−1)Thermal conductivity2.70kh/(cm·s−1)Permeability coefficient1.0 ∗ 10−7β/(C−1)Thawing settlement3.0 ∗ 10−6Boundary conditions are introduced as follows. The theory of saturated unsteady seepage and the coupling stress analyzed the water-thermal coupling of landslide. The conditions of seepage and stress calculation are determined in the calculation conditions.The initial seepage boundary conditions are as follows: the groundwater level in landslide often changes with the change of external factors (such as ice-snow melting water and reservoir water level). The initial groundwater level is determined by using the exploration data provided in the geological model study report. The pore pressure distribution of initial seepage field is determined by slope.The surrounding boundary of the model is calculated according to the permeability coefficient of different rock and soil on the boundary. The corresponding seepage velocity is used to simulate the permeable condition of the boundary.We studied the boundary conditions of unsteady seepage. The ice-snow melting caused infiltration and runoff generation on the slope. The boundary conditions are determined by the melting water strength and the infiltration rate. When the strength of ice-snow melting is less than the infiltration rate of bank slope soil, all the melting snow get into the slope. When the strength of ice-snow melting is greater than the infiltration rate of the bank slope, the bank slope surface occurs runoff. The research work should grasp the essence of the process. The seepage velocity is determined by the permeability coefficient and the melting strength. The boundary conditions of slope infiltration simulation can be determined by the above relationship. The bottom edge of the numerical calculation is impermeable boundary. The boundary around the no contact slope and reservoir water region is the free permeable boundary. The scope of 3D numerical calculation model is selected according to the geological conditions and landform characteristics of landslide. The calculation domain includes sliding mass, sliding zone, and bedrock. The whole computational domain is divided into 44890 hexahedral elements. There are 50864 nodes in total. The faults in the calculation model are simulated by solid elements. The 3D calculation model and mesh are shown in Figure4.Figure 4
The model computing grid.


### 5.2. Computing Method

Based on Section2.2, we have researched on the definite solution of partial differential equation of unsaturated seepage. We derived the negative pressure models of free infiltration and pressure infiltration of ice and snow-melted water. We have carried out the calculation and analysis of the variation law of the negative pressure in the soil during the infiltration process. We have developed a program for iterative calculation. We have realized the simultaneous solution of free infiltration and pressure infiltration. We have replaced the traditional boundary conditions and conversion methods. The calculation process is shown in Figure 5.Figure 5
Flow chart off ice-snow melting infiltration iterative models.We calculated the unsaturated infiltration parameters of ice-snow melting and the negative pressure change rate of the top soil based on the iterative mode. The program has realized the simultaneous solution of the model, which accords with the physical infiltration process.The test of ice-snow is assumed to be multilayer snow, and the melting of the bottom snow can affect from the upper snow cover to the whole snow layer. Firstly, we estimate the melting process in the lower layer, there is no lateral outflow, and all the snow water seeps into the soil. Secondly, we estimate the upper snowmelt process. Lastly, the total amount of multilayer outflow is determined as the total discharge.


### 5.3. Calculation Results and Analysis

During the test, elastic-plastic deformation on the soil body indicates that the model test was not repeatable. This model test mainly considered the combined action test of ice-snow melting water and reservoir water at the level of 100 cm under the extreme ice-snow disaster in winter. The similarity ratio of the model test concerning this test platform is 1 : 175. The model test conditions are as follows.We research the combined action of reservoir water and ice-snow melting. There are continuous works for 15 days.(i)
Build the model bank slope and calibrate the buried sensor which is then debugged to meet the standard. Stable initial settlement of bank slope is achieved when the similar materials of the bank slope reach the initial volumetric water content of the test, i.e., 22%–24%.(ii)
Evenly pave the ice and snow with a thickness of 30 cm in the 50 m2 test bank slope soil.(iii)
The water rises at a speed of 25 cm/h from the bottom (0 cm) of the dam to the high water level at 100 cm. It is a continuous work for 4 hours.(iv)
It is fixed at 175 m water level.We analyzed the stability under the ice-snow melting and reservoir water combined action. We calculate the three-dimensional numerical model of landslide. Then, we research the deformation and stability of different reservoir water levels and ice-snow melting condition [46].(1)
Condition 1: reservoir impoundment from 0 m to 175 m water level and ice-snow loadThis condition is a continuous work for 4 hours. The calculation results are shown in Figures6 and 7. The results are summarized as follows:(i)
The stress calculation results and analysisThe compressive stresses have greatly changed in the front edge of landslide. The maximum principal stress ranges from −1.568 to 0.005 MPa in the landslide mass. The maximum value appears at the trailing edge of the sliding mass. The minimum principal stress ranges from −3.743 to −0.023 MPa. The minimum value appears at the posterior middle part of the sliding mass. The maximum principal stress ranges from −1.940 to 0.026 MPa in the sliding belt. The maximum value appears at the trailing edge of the sliding belt. The minimum principal stress ranges from −4.090 to −0.086 MPa. The minimum value appears at the posterior middle part of the sliding belt. Reservoir impoundment from 0 m to 175 m water level and ice-snow load changed the seepage field. The coupling effect of the seepage field and stress field leads to the change of stress field. This is an important reason for the change of stress field.(ii)
Displacement calculation results and analysisAfter the impoundment to 175 m, a large horizontal displacement occurred in the middle and rear part of the sliding mass, and the maximum horizontal displacement was 2.492 cm. The horizontal displacement of the rear part of the sliding mass ranges from 0.002 cm to 1.657 cm. The horizontal displacement of the middle sliding mass ranges from 0.001 cm to 1.300 cm. The horizontal displacement of the front sliding mass ranges from 0.103 cm to 1.330 cm. This is mainly due to the floating effect of reservoir water on the sliding body, which reduces the antisliding force of the sliding body and causes the horizontal displacement of the sliding body. The downward vertical displacement occurs at the leading edge of the sliding mass, and the range is from 0.058 cm to 0.651 cm. This is mainly due to the self-weight of the sliding mass. There are some abrupt displacements at the junction between the upper edge of the sliding mass and the sliding bed. This is due to the inhomogeneous characteristics of the material properties and topography of the sliding mass and sliding bed in this area.(iii)
Calculation results and analysis of the plastic zoneThe plastic failure zone appears in local landslide; however, the landslide is in a relatively stable state [47](iv)
Calculation results and analysis of pore pressureThe pore pressure distribution diagram is shown in Figure7(a). The water level below 175 is directly affected by the pore pressure. The maximum pore pressure at the leading edge of the sliding mass is about 0.656 MPa. Its value appears at the height of 118 m at the front edge of the sliding mass. The pore pressure at the back edge of the landslide is mainly controlled by the groundwater level, and the maximum pore pressure is 4.912 MPa.(2)
Condition 2: 175 m reservoir water level and ice-snow loadThis condition is a continuous work for 1 hour. The calculation results are shown in Figures8 and 9. The results are summarized as follows:(i)
The stress calculation results analysis [48].The distribution of tensile and compressive stresses has greatly changed in the front edge of landslide. The compressive stresses have greatly changed in the front edge of landslide. In the landslide mass, the maximum principal stress ranges from −2.594 to 0.382 MPa. The maximum value appears at the middle of the sliding mass. The minimum principal stress ranges from −6.078 to −0.028 MPa. The minimum value appears at the posterior part of the sliding mass. The maximum principal stress ranges from −3.072 to 0.026 MPa in the sliding belt. The maximum value appears at the front of the sliding belt. The minimum principal stress ranges from −6.253 to −0.086 MPa. The minimum value appears at the posterior of the sliding belt.(ii)
Displacement calculation results and analysisThe horizontal displacement of the positive of the sliding mass ranges from 0 to 2.205 cm. The horizontal displacement of the rear of the sliding mass is relatively small, with the magnitude from 0.485 to 0.515 cm. The horizontal displacement of the middle of the sliding mass is relatively small, with the magnitude from 0 to 0.245 cm. The horizontal displacement of the middle of the sliding mass is relatively small, with the magnitude from 0.068 to 2.165 cm. The horizontal displacement of condition 2 is larger than condition 1. The trailing edge of the sliding mass have produced downward vertical displacement, with the magnitude from 0.678 to 18.6 cm. The vertical displacement of condition 2 is larger than condition 1.This is mainly due to the softening of the slope soaked in reservoir water, which causes the sliding body to slide downward. At the same time, the settlement occurs under the action of self-weight and melting water of ice-snow [49].(iii)
Calculation results and analysis of the plastic zoneThe plastic failure zone appears in local landslide; however, the landslide is in a relatively stable state in this condition(iv)
Calculation results and analysis of pore pressureThe groundwater is accordingly increased in front of the landslide. The seepage field greatly changes the front edge of landslide. The pore water pressure changes within the fluctuation range of reservoir water at the front edge. The pore water pressure has not changed in the back edge of landslide. This is due to the relatively long distance between the back edge of the landslide and the water level change zone of the reservoir.The ice-snow melting infiltration on the slope seepage field is concentrated in the slope toe. That is, the area where the transient saturation zone is first formed. The maximum horizontal displacement is at the toe. The horizontal displacement decreases with the increase of horizontal depth. The vertical displacement reaches the maximum at the top surface. The vertical depth decreases with the increase of depth. It should be noted that the horizontal displacement at the toe of the slope is relatively dense. It means that obvious slip has occurred at this point. At the same time, because the transient saturation zone is connected with the groundwater level, the pore water pressure increases due to the infiltration of ice-snow melting. As a result, the effective stress decreases at the toe.(3)
Condition 3: 175 m reservoir water level and ice-snow meltingThis condition is a continuous work for 355 hours. The calculation results are shown in Figures10 and 11. The results are summarized as follows:(i)
The stress calculation results and analysisThe tensile and compressive stresses are greatly changed in the front edge of landslide. The maximum principal stress ranges from −2.594 to 0.382 MPa in the middle landslide mass. The minimum principal stress ranges from −4.555 to 0.012 MPa in the posterior landslide mass. The maximum principal stress ranges from −2.414 to 0.029 MPa in the sliding belt. The maximum value appears at the front of the sliding belt. The minimum principal stress ranges from −4.649 to 0.012 MPa. The minimum value appears at the posterior of the sliding belt [50].The 175 m reservoir water level and ice-snow melting changed the seepage field. The effect of seepage field coupling stress field have led to the change of temperature field. This is an important reason for the change of temperature field.(ii)
Displacement calculation results and analysisThe horizontal displacement of the positive of the sliding mass ranges from 0.123 to 0.905 cm. The horizontal displacement of the rear of the sliding mass is relatively smaller, with the magnitude from 0.083 to 0.5 cm. The horizontal displacement of the middle of the sliding mass is relatively more lager, with the magnitude from 0.315 to 0.645 cm. The horizontal displacement of the front of the sliding mass is relatively the most maximum, with the magnitude from 0.368 to 2.165 cm. The trailing edge of the sliding mass have produced downward vertical displacement, with the magnitude from 0 .078 to 43.6 cm. The middle of the sliding mass have produced maximum downward vertical displacement, with the magnitude from 3.078 to 45.6 cm. The vertical displacement of condition 3 is larger than condition 2. This is mainly due to the downward sliding of the sliding mass caused by the infiltration of ice-snow melting. At the same time, the settlement occurs under the self-weight and ice-snow melting combined action. The vertical displacement enlarged the sliding mass.(iii)
Calculation results and analysis of plastic zoneThe plastic failure zone appears in local landslide; however, the landslide is in a relatively stable state in this condition.(iv)
Calculation results and analysis of pore pressureThe groundwater level is accordingly increased in front of the landslide. It is directly affected by reservoir water infiltration and ice-snow melting water. The reservoir water has directly affected groundwater distribution in front of the landslide. The phreatic line of groundwater level slowly rises. This is because it takes certain times. The difference of pore water pressure has dissipated in front of the landslide. The pore water pressure has not changed in the back edge of landslide. This is due to the relatively long distance between the back edge of the landslide and the water level change zone of the reservoir.Figure 6
Distribution of stress field under condition 1. (a) Maximum principal stress field. (b) Minimum principal stress field.
(a)(b)Figure 7
Deformation distribution diagram under condition 1. (a) Pore pressure distribution. (b) Displacement distribution.
(a)(b)Figure 8
Distribution of stress field under condition 2. (a) Maximum principal stress field. (b) Minimum principal stress field.
(a)(b)Figure 9
Deformation distribution diagram under condition 2. (a) Pore pressure distribution. (b) Displacement distribution.
(a)(b)Figure 10
Distribution of stress field under condition 3. (a) Maximum principal stress field. (b) Minimum principal stress field.
(a)(b)Figure 11
Deformation distribution diagram under condition 3. (a) Pore pressure distribution. (b) Displacement distribution.
(a)(b)The accumulation bank slope is studied under three working conditions. The data of stress, pore pressure, and displacement are effectively counted. With the continuous melting of snow, the matric suction gradually decreases in the unsaturated area. The overall stability gradually decreases the slope. The horizontal displacement and vertical displacement gradually increase on the slope.The deformation of the leading edge is larger than that of the trailing edge in the horizontal direction. In the vertical direction, the posterior margin is dislocated, and the front edge is uplifted. The results show that certain snow melting rate and snow melting amount have great influence on the seepage field of loose accumulation bank slope.


### 5.4. Comparison of Results and Discussion

The change data of horizontal and vertical displacement collected by numerical simulation are effectively counted. Draw the relationship curve between slope displacement and snow melting time, as shown in Figure12. The horizontal displacement of the slope increases gradually from the back edge to the front edge and reaches the maximum near the toe of the slope, while the vertical displacement is just the opposite, which is consistent with the actual state of the slope.Figure 12
The cumulative displacement of monitoring points. (a) Cumulative horizontal displacement. (b) Cumulative vertical displacement.
(a)(b)In this simulation, the hydraulic gradients as the same as the velocity fields are homogeneous with the depth; therefore, homogeneous erosion rates will be generated, which lead to homogeneous fine content distributions within the soil specimens. The evaluation profiles of the eroded fines mass ratio in the specimens under different constant seepage velocities are presented in Figure13.Figure 13
Variations of hydraulic gradient and time with average internal erosion velocity.We studied the numerical simulation results of the three working conditions. We collect numerical simulation results. We draw the distribution curve of pore water pressure of the ice-snow melting (as shown in Figure14). With the continuous melting of ice-snow, the horizontal displacement gradually increases. The deformation of the leading edge is greater than that of the trailing edge. It shows that the vertical displacement gradually increases with the continuous snowmelt. It is shows a downward dislocation at the trailing edge and an uplift at the front edge. The infiltration increased with ice-snow melting, the horizontal displacement and vertical displacement gradually increased, while the area of obvious displacement gradually reduced.Figure 14
Pore water pressure distribution curve of ice-snow melting. (a) Before ice-snow melting. (b) After ice-snow melting.
(a)(b)Under the combined effect of reservoir water, melting ice, and snow, there is infiltration on the bank slope. The dynamic evolution characteristics of time-varying stability (10 days interval) are shown in Figure15. The stability coefficient of landslide has three stages:(i)
When the reservoir water level rises, the stability coefficient gradually increases the bank slope.(ii)
During the period of rapid snow melting, the stability coefficient slowly decreases the bank slope.(iii)
The temperature rising in the daytime is a continuous infiltration process of ice-snow melting, and the stability coefficient rapidly decreases. With the decrease of temperature at night, the continuous infiltration capacity of ice-snow melting decreases, and the stability coefficient of landslide gradually stabilizes.Figure 15
Evolution characteristics of stability coefficient of ice-snow melting.The stability of the bank slope caused by ice-snow melting is a cyclic and gradual weakening process, and the stability coefficient is also gradually reduced.The accumulation landslide studied the key control factors of the combined action of reservoir water and ice-snow melting as follows:(i)
We study the ice-snow melting strength index. The maximum daily snow melting amount of 23 mm can be defined as the precipitation threshold.(ii)
We study the snow water infiltration rate index. The evolution process of shallow landslide is greatly affected by the dynamic change of groundwater seepage field.(iii)
We are study the penetration rate of shear plastic zone of landslide. The criterion is 90%.The relationship curve is between stability coefficient and time, under different infiltration rates as shown in Figure16.Figure 16
The relationship curve between stability coefficient and time, under different infiltration rates.


## 5.1. Calculation Model, Boundary Conditions, and Calculation Parameters

The physical and mechanical parameters of rock and soil are determined. The physical and mechanical parameters of soil mainly include the permeability coefficient (k), the heavy density (γ), Poisson’s ratio (μ), the cohesion (c), the internal friction angle (φ), and so on. We determine main physical and mechanical parameters. The range of values is shown in Table 1.Table 1
Basic property indexes of soils.
SymbolMeaningQuantityγ (kg/m3)Heavy density21.4k (m/s)Permeability coefficient2.3 ∗ 10−6C (kPa)Cohesion17.2φ (°)Internal friction angle18.6Based on the physical and mechanical parameters of rock and soil, test data and recommended parameters are determined. The soil water parameters of temperature field and seepage field are obtained as follows in Table2.Table 2
The thermal-hydraulic parameter of soil and water.
SymbolMeaningQuantityCv/(J·kg−1·°C−1)Specific heat16λ/(W·m−1·°C−1)Thermal conductivity2.70kh/(cm·s−1)Permeability coefficient1.0 ∗ 10−7β/(C−1)Thawing settlement3.0 ∗ 10−6Boundary conditions are introduced as follows. The theory of saturated unsteady seepage and the coupling stress analyzed the water-thermal coupling of landslide. The conditions of seepage and stress calculation are determined in the calculation conditions.The initial seepage boundary conditions are as follows: the groundwater level in landslide often changes with the change of external factors (such as ice-snow melting water and reservoir water level). The initial groundwater level is determined by using the exploration data provided in the geological model study report. The pore pressure distribution of initial seepage field is determined by slope.The surrounding boundary of the model is calculated according to the permeability coefficient of different rock and soil on the boundary. The corresponding seepage velocity is used to simulate the permeable condition of the boundary.We studied the boundary conditions of unsteady seepage. The ice-snow melting caused infiltration and runoff generation on the slope. The boundary conditions are determined by the melting water strength and the infiltration rate. When the strength of ice-snow melting is less than the infiltration rate of bank slope soil, all the melting snow get into the slope. When the strength of ice-snow melting is greater than the infiltration rate of the bank slope, the bank slope surface occurs runoff. The research work should grasp the essence of the process. The seepage velocity is determined by the permeability coefficient and the melting strength. The boundary conditions of slope infiltration simulation can be determined by the above relationship. The bottom edge of the numerical calculation is impermeable boundary. The boundary around the no contact slope and reservoir water region is the free permeable boundary. The scope of 3D numerical calculation model is selected according to the geological conditions and landform characteristics of landslide. The calculation domain includes sliding mass, sliding zone, and bedrock. The whole computational domain is divided into 44890 hexahedral elements. There are 50864 nodes in total. The faults in the calculation model are simulated by solid elements. The 3D calculation model and mesh are shown in Figure4.Figure 4
The model computing grid.


## 5.2. Computing Method

Based on Section2.2, we have researched on the definite solution of partial differential equation of unsaturated seepage. We derived the negative pressure models of free infiltration and pressure infiltration of ice and snow-melted water. We have carried out the calculation and analysis of the variation law of the negative pressure in the soil during the infiltration process. We have developed a program for iterative calculation. We have realized the simultaneous solution of free infiltration and pressure infiltration. We have replaced the traditional boundary conditions and conversion methods. The calculation process is shown in Figure 5.Figure 5
Flow chart off ice-snow melting infiltration iterative models.We calculated the unsaturated infiltration parameters of ice-snow melting and the negative pressure change rate of the top soil based on the iterative mode. The program has realized the simultaneous solution of the model, which accords with the physical infiltration process.The test of ice-snow is assumed to be multilayer snow, and the melting of the bottom snow can affect from the upper snow cover to the whole snow layer. Firstly, we estimate the melting process in the lower layer, there is no lateral outflow, and all the snow water seeps into the soil. Secondly, we estimate the upper snowmelt process. Lastly, the total amount of multilayer outflow is determined as the total discharge.


## 5.3. Calculation Results and Analysis

During the test, elastic-plastic deformation on the soil body indicates that the model test was not repeatable. This model test mainly considered the combined action test of ice-snow melting water and reservoir water at the level of 100 cm under the extreme ice-snow disaster in winter. The similarity ratio of the model test concerning this test platform is 1 : 175. The model test conditions are as follows.We research the combined action of reservoir water and ice-snow melting. There are continuous works for 15 days.(i)
Build the model bank slope and calibrate the buried sensor which is then debugged to meet the standard. Stable initial settlement of bank slope is achieved when the similar materials of the bank slope reach the initial volumetric water content of the test, i.e., 22%–24%.(ii)
Evenly pave the ice and snow with a thickness of 30 cm in the 50 m2 test bank slope soil.(iii)
The water rises at a speed of 25 cm/h from the bottom (0 cm) of the dam to the high water level at 100 cm. It is a continuous work for 4 hours.(iv)
It is fixed at 175 m water level.We analyzed the stability under the ice-snow melting and reservoir water combined action. We calculate the three-dimensional numerical model of landslide. Then, we research the deformation and stability of different reservoir water levels and ice-snow melting condition [46].(1)
Condition 1: reservoir impoundment from 0 m to 175 m water level and ice-snow loadThis condition is a continuous work for 4 hours. The calculation results are shown in Figures6 and 7. The results are summarized as follows:(i)
The stress calculation results and analysisThe compressive stresses have greatly changed in the front edge of landslide. The maximum principal stress ranges from −1.568 to 0.005 MPa in the landslide mass. The maximum value appears at the trailing edge of the sliding mass. The minimum principal stress ranges from −3.743 to −0.023 MPa. The minimum value appears at the posterior middle part of the sliding mass. The maximum principal stress ranges from −1.940 to 0.026 MPa in the sliding belt. The maximum value appears at the trailing edge of the sliding belt. The minimum principal stress ranges from −4.090 to −0.086 MPa. The minimum value appears at the posterior middle part of the sliding belt. Reservoir impoundment from 0 m to 175 m water level and ice-snow load changed the seepage field. The coupling effect of the seepage field and stress field leads to the change of stress field. This is an important reason for the change of stress field.(ii)
Displacement calculation results and analysisAfter the impoundment to 175 m, a large horizontal displacement occurred in the middle and rear part of the sliding mass, and the maximum horizontal displacement was 2.492 cm. The horizontal displacement of the rear part of the sliding mass ranges from 0.002 cm to 1.657 cm. The horizontal displacement of the middle sliding mass ranges from 0.001 cm to 1.300 cm. The horizontal displacement of the front sliding mass ranges from 0.103 cm to 1.330 cm. This is mainly due to the floating effect of reservoir water on the sliding body, which reduces the antisliding force of the sliding body and causes the horizontal displacement of the sliding body. The downward vertical displacement occurs at the leading edge of the sliding mass, and the range is from 0.058 cm to 0.651 cm. This is mainly due to the self-weight of the sliding mass. There are some abrupt displacements at the junction between the upper edge of the sliding mass and the sliding bed. This is due to the inhomogeneous characteristics of the material properties and topography of the sliding mass and sliding bed in this area.(iii)
Calculation results and analysis of the plastic zoneThe plastic failure zone appears in local landslide; however, the landslide is in a relatively stable state [47](iv)
Calculation results and analysis of pore pressureThe pore pressure distribution diagram is shown in Figure7(a). The water level below 175 is directly affected by the pore pressure. The maximum pore pressure at the leading edge of the sliding mass is about 0.656 MPa. Its value appears at the height of 118 m at the front edge of the sliding mass. The pore pressure at the back edge of the landslide is mainly controlled by the groundwater level, and the maximum pore pressure is 4.912 MPa.(2)
Condition 2: 175 m reservoir water level and ice-snow loadThis condition is a continuous work for 1 hour. The calculation results are shown in Figures8 and 9. The results are summarized as follows:(i)
The stress calculation results analysis [48].The distribution of tensile and compressive stresses has greatly changed in the front edge of landslide. The compressive stresses have greatly changed in the front edge of landslide. In the landslide mass, the maximum principal stress ranges from −2.594 to 0.382 MPa. The maximum value appears at the middle of the sliding mass. The minimum principal stress ranges from −6.078 to −0.028 MPa. The minimum value appears at the posterior part of the sliding mass. The maximum principal stress ranges from −3.072 to 0.026 MPa in the sliding belt. The maximum value appears at the front of the sliding belt. The minimum principal stress ranges from −6.253 to −0.086 MPa. The minimum value appears at the posterior of the sliding belt.(ii)
Displacement calculation results and analysisThe horizontal displacement of the positive of the sliding mass ranges from 0 to 2.205 cm. The horizontal displacement of the rear of the sliding mass is relatively small, with the magnitude from 0.485 to 0.515 cm. The horizontal displacement of the middle of the sliding mass is relatively small, with the magnitude from 0 to 0.245 cm. The horizontal displacement of the middle of the sliding mass is relatively small, with the magnitude from 0.068 to 2.165 cm. The horizontal displacement of condition 2 is larger than condition 1. The trailing edge of the sliding mass have produced downward vertical displacement, with the magnitude from 0.678 to 18.6 cm. The vertical displacement of condition 2 is larger than condition 1.This is mainly due to the softening of the slope soaked in reservoir water, which causes the sliding body to slide downward. At the same time, the settlement occurs under the action of self-weight and melting water of ice-snow [49].(iii)
Calculation results and analysis of the plastic zoneThe plastic failure zone appears in local landslide; however, the landslide is in a relatively stable state in this condition(iv)
Calculation results and analysis of pore pressureThe groundwater is accordingly increased in front of the landslide. The seepage field greatly changes the front edge of landslide. The pore water pressure changes within the fluctuation range of reservoir water at the front edge. The pore water pressure has not changed in the back edge of landslide. This is due to the relatively long distance between the back edge of the landslide and the water level change zone of the reservoir.The ice-snow melting infiltration on the slope seepage field is concentrated in the slope toe. That is, the area where the transient saturation zone is first formed. The maximum horizontal displacement is at the toe. The horizontal displacement decreases with the increase of horizontal depth. The vertical displacement reaches the maximum at the top surface. The vertical depth decreases with the increase of depth. It should be noted that the horizontal displacement at the toe of the slope is relatively dense. It means that obvious slip has occurred at this point. At the same time, because the transient saturation zone is connected with the groundwater level, the pore water pressure increases due to the infiltration of ice-snow melting. As a result, the effective stress decreases at the toe.(3)
Condition 3: 175 m reservoir water level and ice-snow meltingThis condition is a continuous work for 355 hours. The calculation results are shown in Figures10 and 11. The results are summarized as follows:(i)
The stress calculation results and analysisThe tensile and compressive stresses are greatly changed in the front edge of landslide. The maximum principal stress ranges from −2.594 to 0.382 MPa in the middle landslide mass. The minimum principal stress ranges from −4.555 to 0.012 MPa in the posterior landslide mass. The maximum principal stress ranges from −2.414 to 0.029 MPa in the sliding belt. The maximum value appears at the front of the sliding belt. The minimum principal stress ranges from −4.649 to 0.012 MPa. The minimum value appears at the posterior of the sliding belt [50].The 175 m reservoir water level and ice-snow melting changed the seepage field. The effect of seepage field coupling stress field have led to the change of temperature field. This is an important reason for the change of temperature field.(ii)
Displacement calculation results and analysisThe horizontal displacement of the positive of the sliding mass ranges from 0.123 to 0.905 cm. The horizontal displacement of the rear of the sliding mass is relatively smaller, with the magnitude from 0.083 to 0.5 cm. The horizontal displacement of the middle of the sliding mass is relatively more lager, with the magnitude from 0.315 to 0.645 cm. The horizontal displacement of the front of the sliding mass is relatively the most maximum, with the magnitude from 0.368 to 2.165 cm. The trailing edge of the sliding mass have produced downward vertical displacement, with the magnitude from 0 .078 to 43.6 cm. The middle of the sliding mass have produced maximum downward vertical displacement, with the magnitude from 3.078 to 45.6 cm. The vertical displacement of condition 3 is larger than condition 2. This is mainly due to the downward sliding of the sliding mass caused by the infiltration of ice-snow melting. At the same time, the settlement occurs under the self-weight and ice-snow melting combined action. The vertical displacement enlarged the sliding mass.(iii)
Calculation results and analysis of plastic zoneThe plastic failure zone appears in local landslide; however, the landslide is in a relatively stable state in this condition.(iv)
Calculation results and analysis of pore pressureThe groundwater level is accordingly increased in front of the landslide. It is directly affected by reservoir water infiltration and ice-snow melting water. The reservoir water has directly affected groundwater distribution in front of the landslide. The phreatic line of groundwater level slowly rises. This is because it takes certain times. The difference of pore water pressure has dissipated in front of the landslide. The pore water pressure has not changed in the back edge of landslide. This is due to the relatively long distance between the back edge of the landslide and the water level change zone of the reservoir.Figure 6
Distribution of stress field under condition 1. (a) Maximum principal stress field. (b) Minimum principal stress field.
(a)(b)Figure 7
Deformation distribution diagram under condition 1. (a) Pore pressure distribution. (b) Displacement distribution.
(a)(b)Figure 8
Distribution of stress field under condition 2. (a) Maximum principal stress field. (b) Minimum principal stress field.
(a)(b)Figure 9
Deformation distribution diagram under condition 2. (a) Pore pressure distribution. (b) Displacement distribution.
(a)(b)Figure 10
Distribution of stress field under condition 3. (a) Maximum principal stress field. (b) Minimum principal stress field.
(a)(b)Figure 11
Deformation distribution diagram under condition 3. (a) Pore pressure distribution. (b) Displacement distribution.
(a)(b)The accumulation bank slope is studied under three working conditions. The data of stress, pore pressure, and displacement are effectively counted. With the continuous melting of snow, the matric suction gradually decreases in the unsaturated area. The overall stability gradually decreases the slope. The horizontal displacement and vertical displacement gradually increase on the slope.The deformation of the leading edge is larger than that of the trailing edge in the horizontal direction. In the vertical direction, the posterior margin is dislocated, and the front edge is uplifted. The results show that certain snow melting rate and snow melting amount have great influence on the seepage field of loose accumulation bank slope.


## 5.4. Comparison of Results and Discussion

The change data of horizontal and vertical displacement collected by numerical simulation are effectively counted. Draw the relationship curve between slope displacement and snow melting time, as shown in Figure12. The horizontal displacement of the slope increases gradually from the back edge to the front edge and reaches the maximum near the toe of the slope, while the vertical displacement is just the opposite, which is consistent with the actual state of the slope.Figure 12
The cumulative displacement of monitoring points. (a) Cumulative horizontal displacement. (b) Cumulative vertical displacement.
(a)(b)In this simulation, the hydraulic gradients as the same as the velocity fields are homogeneous with the depth; therefore, homogeneous erosion rates will be generated, which lead to homogeneous fine content distributions within the soil specimens. The evaluation profiles of the eroded fines mass ratio in the specimens under different constant seepage velocities are presented in Figure13.Figure 13
Variations of hydraulic gradient and time with average internal erosion velocity.We studied the numerical simulation results of the three working conditions. We collect numerical simulation results. We draw the distribution curve of pore water pressure of the ice-snow melting (as shown in Figure14). With the continuous melting of ice-snow, the horizontal displacement gradually increases. The deformation of the leading edge is greater than that of the trailing edge. It shows that the vertical displacement gradually increases with the continuous snowmelt. It is shows a downward dislocation at the trailing edge and an uplift at the front edge. The infiltration increased with ice-snow melting, the horizontal displacement and vertical displacement gradually increased, while the area of obvious displacement gradually reduced.Figure 14
Pore water pressure distribution curve of ice-snow melting. (a) Before ice-snow melting. (b) After ice-snow melting.
(a)(b)Under the combined effect of reservoir water, melting ice, and snow, there is infiltration on the bank slope. The dynamic evolution characteristics of time-varying stability (10 days interval) are shown in Figure15. The stability coefficient of landslide has three stages:(i)
When the reservoir water level rises, the stability coefficient gradually increases the bank slope.(ii)
During the period of rapid snow melting, the stability coefficient slowly decreases the bank slope.(iii)
The temperature rising in the daytime is a continuous infiltration process of ice-snow melting, and the stability coefficient rapidly decreases. With the decrease of temperature at night, the continuous infiltration capacity of ice-snow melting decreases, and the stability coefficient of landslide gradually stabilizes.Figure 15
Evolution characteristics of stability coefficient of ice-snow melting.The stability of the bank slope caused by ice-snow melting is a cyclic and gradual weakening process, and the stability coefficient is also gradually reduced.The accumulation landslide studied the key control factors of the combined action of reservoir water and ice-snow melting as follows:(i)
We study the ice-snow melting strength index. The maximum daily snow melting amount of 23 mm can be defined as the precipitation threshold.(ii)
We study the snow water infiltration rate index. The evolution process of shallow landslide is greatly affected by the dynamic change of groundwater seepage field.(iii)
We are study the penetration rate of shear plastic zone of landslide. The criterion is 90%.The relationship curve is between stability coefficient and time, under different infiltration rates as shown in Figure16.Figure 16
The relationship curve between stability coefficient and time, under different infiltration rates.


## 6. Conclusions

(1)
We established practical coupled numerical model of stress field, seepage field, and temperature field of saturated soil with ice-snow melting. We constructed the viscoelastic constitutive relation of soil skeleton affected by temperature. The deformation of thawed soil is divided into ice-snow load and temperature gradient. We established the mass conservation equation of soil skeleton under the combined action of ice-snow load and temperature field.(2)
We established the three-dimensional evolution model of homogeneous loose accumulation slope under extreme snow and ice disaster conditions. We calculated the evolution characteristics of homogeneous loose accumulation bank slope under different rainfall and snow equivalent. We constructed the relationship between landslide stability coefficient and time, and it is under the condition of different ice-snow melting infiltration.(3)
We analyze the numerical simulation results, statistical effective data and drawing bank slope deformation curve. With the continuous ice-snow melting, the matric suction gradually decreases in the unsaturated zone. We studied the coupling effect of water and heat in the process of ice-snow thawing in three kinds of working conditions. We verified the correctness and rationality of the program.



---

*Source: 6664213-2021-03-19.xml*