5.1 Numerical model
Based on the energy conservation law, the heat transfer equation of soil considering heat conduction and heat convection (Zhang et al. 2016) is as follows:
\({\rho }_{\text{s}}C\frac{\partial T}{\partial t}=\nabla \left(\lambda \nabla T\right)+{C}_{\text{w}}{\rho }_{\text{w}}\nabla [-K\nabla \left(\frac{p}{{\gamma }_{\text{w}}}\right)\bullet T]\)
|
(1)
|
where \(C\) and \(\lambda\) are equivalent specific heat and thermal conductivity of soil, and \({C}_{\text{w}}\) is the specific heat of the water. \({\rho }_{\text{s}}\) is the soil density, and \({\rho }_{\text{w}}\) is the water density. \({\gamma }_{\text{w}}\) is water bulk density, \(p\) is pore water pressure, and \(K\) is the soil permeability coefficient. In this study, the temperature range of the ice-water phase transition is assumed to be \({T}_{\text{b}}\le T\le {T}_{\text{p}}\), where \({T}_{\text{b}}\)=-1.5 ℃ and \({T}_{\text{p}}\)=-0.05 ℃. Latent heat of the phase transition has a great influence on heat transfer, which can be considered by the sensible heat capacity method. The expression is as follows:
\(C=\left\{\begin{array}{ll}{C}_{\text{u}}& T>{T}_{\text{p}}\\ {C}_{\text{f}}+\frac{{C}_{\text{u}}-{C}_{\text{f}}}{{T}_{\text{p}}-{T}_{\text{b}}}\left[T-{T}_{\text{b}}\right]+\frac{L}{1+{w}_{0}}\frac{\partial \left({w}_{0}-{w}_{\text{u}}\right)}{\partial T}& {T}_{\text{b}}\le T\le {T}_{\text{p}}\\ {C}_{\text{f}}& T<{T}_{\text{b}}\end{array}\right.\)
|
(2)
|
\(\lambda =\left\{\begin{array}{ll}{\lambda }_{\text{u}}& T>{T}_{\text{p}}\\ {\lambda }_{\text{f}}+\frac{{\lambda }_{\text{u}}-{\lambda }_{\text{f}}}{{T}_{\text{p}}-{T}_{\text{b}}}\left[T-{T}_{\text{b}}\right]& {T}_{\text{b}}\le T\le {T}_{\text{p}}\\ {\lambda }_{\text{f}}& T<{T}_{\text{b}}\end{array}\right.\)
|
(3)
|
Unfrozen water content is a function of temperature, and the expression is as follows:
\({w}_{\text{u}}=a\bullet {\left|T\right|}^{-b}\)
|
(4)
|
where a and b are test parameters.
Assuming that ice is incompressible, the total stress of frozen soil in a saturated state is composed of effective stress and pore water pressure. According to Darcy's law, the seepage equation is as follows:
\(\frac{\partial }{\partial t}\left({\rho }_{\text{w}}{\epsilon }_{\text{e}\text{l}}\right)+\nabla \bullet \left[{-\rho }_{\text{w}}\frac{k}{\mu }\left(\nabla p+g\nabla D\right)\right]={Q}_{\text{m}}\)
|
(5)
|
where \({\epsilon }_{\text{e}\text{l}}\) is the elastic volumetric strain of soil, \(g\) is the gravitational acceleration, \(k\) is the permeability. \(\mu\) is the viscosity coefficient of water, \(\nabla D\) is the relative height, and \({Q}_{\text{m}}\) is the source term.
When the soil temperature is higher than the freezing point, the permeability of the soil varies with the volumetric strain. When the soil temperature is below freezing point, the permeability coefficient is regarded as a function of soil temperature and is expressed as follows:
\(K=\left\{\begin{array}{ll}\frac{{K}_{-1}}{({T}_{\text{p}}-T{)}^{\beta }}& T<{T}_{\text{p}}\\ \frac{{K}_{\text{u}}}{1+{\epsilon }_{\text{e}\text{l}}}{\left(1+\frac{{\epsilon }_{\text{e}\text{l}}}{{n}_{0}}\right)}^{3}& T\ge {T}_{\text{p}}\end{array}\right.\)
|
(6)
|
where \({K}_{-1}\) is the permeability coefficient of soil at -1℃, \({K}_{\text{u}}\) is the water conductivity coefficient of thawed soil, \({n}_{0}\) is the initial porosity, and \(\beta\) = 2.7 is the test parameter.
Because the sliding rate of the slope is small, the deformation of the sliding body can be regarded as a quasi-static process. Shear stress \(\tau\) of soil is as follows:
\(\tau =c+{\sigma }_{\text{d}}tan\phi\)
|
(7)
|
where \(c\) is the cohesive force of soil, \({\sigma }_{\text{d}}\) is the normal stress of soil, and \(\phi\) is the internal friction angle of soil.
The damage constitutive equation of soil is as follows:
\({\sigma }_{\text{d}}=\left[1-d\left(\kappa \right)\right]E{\epsilon }_{\text{e}\text{l}}\)
|
(8)
|
\({\epsilon }_{\text{e}\text{l}}=\epsilon -{\epsilon }_{\text{t}\text{h}}-{\epsilon }_{\text{t}\text{r}}-{\epsilon }_{\text{h}\text{p}}-{\epsilon }_{0}\)
|
(9)
|
The strain displacement equation is as follows:
\(\epsilon =\frac{1}{2}\left[\nabla \mathbf{u}+{\left(\nabla \mathbf{u}\right)}^{\text{T}}\right]\)
|
(10)
|
where \(\epsilon\) is the total strain of soil. E is the elastic modulus of the soil in the undamaged state. \({\epsilon }_{\text{t}\text{h}}\) is the strain caused by thermal expansion, \({\epsilon }_{\text{t}\text{r}}\) is the strain caused by the ice-water phase transition, \({\epsilon }_{\text{h}\text{p}}\) is the strain caused by matrix potential, and \({\epsilon }_{0}\) is the initial strain. In the sliding state of the slope, these strains are small and negligible for the total strain (Liu and Yu 2011; Shan et al. 2022a). So \(\epsilon\) is approximately equal to \({\epsilon }_{\text{e}\text{l}}\). In addition, \(\kappa\) is a state variable, and \(d\left(\kappa \right)\) is a scalar damage variable, expressed as follows:
\(d\left(\kappa \right)=1-{e}^{-{\left(F/{F}_{0}\right)}^{\beta }}\)
|
(11)
|
where F is the nominal stress based on the Mohr-coulomb failure criterion and \(\beta\) and \({F}_{0}\) are unknown parameters of the Weibull distribution (Li et al. 2009).
Equations (1–11) constitute a hydro-thermal-mechanical coupling model of saturated frozen soil.
5.2 Soil parameters
According to laboratory test results, physical parameters of the soil layers are shown in Table 1, mechanical parameters in Table 2, and seepage parameters in Table 3. Other mechanical parameters are expressed as follows:
\(\begin{array}{ll}E=\left\{\begin{array}{l}{a}_{1}+{b}_{1}{\left|T\right|}^{0.6}\\ {a}_{1}\end{array}\right.& \begin{array}{l}T\le {T}_{p}\\ T>{T}_{p}\end{array}\\ v=\left\{\begin{array}{l}{a}_{2}+{b}_{2}\left|T\right|\\ {a}_{2}\end{array}\right.& \begin{array}{l}T\le {T}_{p}\\ T>{T}_{p}\end{array}\\ c=\left\{\begin{array}{l}{a}_{3}+{b}_{3}\left|T\right|\\ {a}_{3}\end{array}\right.& \begin{array}{l}T\le {T}_{p}\\ T>{T}_{p}\end{array}\\ \phi =\left\{\begin{array}{l}{a}_{4}+{b}_{4}\left|T\right|\\ {a}_{4}\end{array}\right.& \begin{array}{l}T\le {T}_{p}\\ T>{T}_{p}\end{array}\end{array}\)
|
(12)
|
Table 1
Physical parameters of the soil layers.
Soil layer
|
\({\lambda }_{\text{f}}\)
\(\text{W}/\left(\text{m}\bullet \text{K}\right)\)
|
\({\lambda }_{\text{u}}\)
\(\text{W}/\left(\text{m}\bullet \text{K}\right)\)
|
\({C}_{\text{f}}\)
\(\text{J}/\left({\text{m}}^{3}\bullet \text{K}\right)\)
|
\({C}_{\text{u}}\)
\(\text{J}/\left({\text{m}}^{3}\bullet \text{K}\right)\)
|
\({\rho }_{\text{s}}\)
\(\text{k}\text{g}/{\text{m}}^{3}\)
|
\({n}_{0}\)
|
\({w}_{0}\)
%
|
Embankment fill
|
1.980
|
1.920
|
\(1.625\times {10}^{6}\)
|
\(1.982\times {10}^{6}\)
|
1940
|
0.259
|
14.2
|
Cumulose soil
|
1.027
|
0.883
|
\(2.947\times {10}^{6}\)
|
\(3.436\times {10}^{6}\)
|
1802
|
0.496
|
36.7
|
Weathered sand
|
1.651
|
1.447
|
\(2.019\times {10}^{6}\)
|
\(3.175\times {10}^{6}\)
|
1821
|
0.483
|
24.6
|
Gravelly sand
|
1.870
|
1.692
|
\(2.215\times {10}^{6}\)
|
\(3.334\times {10}^{6}\)
|
1837
|
0.479
|
48.9
|
Siltstone
|
1.443
|
1.213
|
\(1.949\times {10}^{6}\)
|
\(2.771\times {10}^{6}\)
|
1845
|
0.447
|
26.3
|
Silty clay
Mudstone
|
1.229
1.973
|
1.105
1.836
|
\(2.011\times {10}^{6}\)
\(2.062\times {10}^{6}\)
|
\(2.968\times {10}^{6}\)
\(2.821\times {10}^{6}\)
|
1799
1811
|
0.491
0.313
|
31.1
29.9
|
Table 2 Mechanical parameters of the soil layers.
Soil layer
|
\({a}_{1}\left(\text{M}\text{P}\text{a}\right)\)
|
\({b}_{1}\)
|
\({a}_{2}\)
|
\({b}_{2}\)
|
\({a}_{3}\)(MPa)
|
\({b}_{3}\)
|
\({a}_{4}\)(o)
|
\({b}_{4}\)
|
\({F}_{0}\)
|
\(\beta\)
|
Embankment fill
|
48.43
|
27.26
|
0.35
|
-0.007
|
0.014
|
0.032
|
15.00
|
0.35
|
3.86
|
1.10
|
Cumulose soil
|
1.26
|
0.24
|
0.48
|
-0.012
|
0.068
|
0.034
|
9.45
|
0.56
|
2.28
|
1.14
|
Weathered sand
|
28.71
|
4.19
|
0.28
|
-0.002
|
0.016
|
0.005
|
26.26
|
0.61
|
3.91
|
1.12
|
Gravelly sand
|
24.37
|
4.68
|
0.25
|
-0.004
|
0.010
|
0.019
|
22.39
|
0.68
|
3.77
|
1.12
|
Siltstone
|
58.73
|
8.18
|
0.32
|
-0.002
|
0.011
|
0.006
|
32.16
|
0.54
|
3.17
|
1.11
|
Silty clay
Mudstone
|
46.92
87.64
|
24.17
22.76
|
0.32
0.25
|
-0.007
-0.001
|
0.097
0.142
|
0.025
0.058
|
18.49
26.74
|
0.45
0.59
|
3.29
3.43
|
1.14
1.11
|
Table 3 Seepage parameters of the soil layers.
Soil layer
|
\({K}_{-1}\left(\text{m}/\text{s}\right)\)
|
\({K}_{\text{u}}\left(\text{m}/\text{s}\right)\)
|
a
|
b
|
Embankment fill
|
\(3.76\times {10}^{-10}\)
|
\(4.32\times {10}^{-6}\)
|
0.082
|
0.28
|
Cumulose soil
|
\(7.71\times {10}^{-09}\)
|
\(2.35\times {10}^{-5}\)
|
0.242
|
0.36
|
Weathered sand
|
\(5.54\times {10}^{-11}\)
|
\(7.56\times {10}^{-8}\)
|
0.133
|
0.32
|
Gravelly sand
Siltstone
Silty clay
mudstone
|
\(3.92\times {10}^{-11}\)
\(9.27\times {10}^{-11}\)
\(3.66\times {10}^{-11}\)
\(1.45\times {10}^{-13}\)
|
\(7.28\times {10}^{-8}\)
\(9.82\times {10}^{-8}\)
\(6.19\times {10}^{-8}\)
\(4.18\times {10}^{-10}\)
|
0.256
0.114
0.172
0.153
|
0.31
0.38
0.36
0.34
|
5.3 Boundary and initial conditions
Temperature boundary conditions of slope can be expressed as:
\(T={T}_{0}+Asin\left(2\pi t+{\alpha }_{0}\right)+\varDelta Tt\)
|
(13)
|
where α0 is the initial phase angle determined by the model starting time (October 15, 2009), and its value is π. T is time. ΔT is the climate warming rate, and the value is 0.052 ℃/a (Shan et al. 2020). T0 is the mean annual temperature of the slope surface, with a value of 0.44 ℃. A is the annual temperature variation amplitude, with a value of 20.63 ℃. The two lateral sides of the geometric model (Fig. 9) are adiabatic and waterproof boundaries, and the horizontal displacement is restricted. The bottom boundary is impermeable, the geothermal flux is 0.03 W/m2, and the vertical and horizontal displacements are restricted. Except for asphalt pavement, pore water pressure at other positions of slope surface is 0. The entire slope is a free boundary and the displacement is not restricted. It should be noted that asphalt pavement is an impermeable boundary and its displacement is not restricted.
We had buried temperature sensors at the foot of the embankment, but the landslide destroyed them. Therefore, the initial hydrothermal field distribution was obtained by steady-state calculation using the slope hydrothermal and deformation boundary conditions on October 15, 2009. Then, the deformation, pore water pressure change, and effective stress distribution of the slope from October 15, 2009, to October 15, 2014, are calculated considering climate change.
It should be noted that the meltwater recharge area of the slope was determined by the high-density electrical diagram of the HH' section (Fig. 12). From July 19, 2010, to October 1, 2014, the pore water pressure monitoring results showed that meltwater recharge decreased year by year with the increase of time, which was reflected by the source term \({Q}_{\text{m}}\) in the seepage equation:
\({Q}_{\text{m}}=\left[0.032\left(1-0.2t\right)+0.011sin\left(2\pi t\right)+0.162\right]\times {10}^{-10}\)
|
(14)
|