Numerical Analysis of the Effects of Crack Characteristics on the Stress and Deformation of Unsaturated Soil Slopes

Abstract

Cracks induced by evaporation or rainfall have a great influence on the stability of unsaturated soil slopes, which can lead to landslides during the rainfall process. In order to study the effect of crack characteristics on the evolution of stress and deformation of unsaturated soil slopes, a series of numerical analyses under different conditions were performed using a coupled elastoplastic finite element program that we developed for unsaturated soil. When carrying out the numerical analyses, the effective stress for unsaturated soil proposed by Bishop and an elastoplastic double-hardening constitutive model for the soil skeleton were employed. The varying parameters, including the crack location, the discharge speed, evaporation rate, infiltration rate, and tensile strength, were investigated to study the coupling process of pore water pressure and deformation in the process of evaporation and rainfall infiltration. The numerical results showed that the minimum pore water pressure of the soil slope at the end of evaporation/rainfall decreased gradually and the crack width increased gradually as the crack set closer to the slope; the larger the discharge speed of pore air, the greater the crack width. With the increase in the evaporation rate, the pore water pressure of the soil slope reduced and the crack initiated earlier and became wider. As the infiltration rate increased, the pore water pressure of the soil slope and the crack width increased, but the decreasing duration became shorter. The change of tensile strength had little effect on the pore water pressure, but the development of the crack width changed with evaporation and rainfall infiltration.

1. Introduction

Slope instability is related to many factors, such as geological structure, geotechnical properties, climatic conditions, topography, surface water effects, groundwater activities, human engineering activities, and earthquakes. Among them, rainfall infiltration is one of the most important conditions that cause slope instability. The occurrence and development of slope failures, such as collapse, landslide, and debris flow, are mostly controlled by factors such as rainfall infiltration [1]. During the process of evaporation and rainfall infiltration, the development and closure of slope cracks in unsaturated soils may be induced, which will have a great influence on the evolution of stress and deformation of the unsaturated soil slope.

Due to the complexity of field tests on unsaturated soil slopes, there is scant research in this area. Blatz et al. [2] initiated a project including a field investigation program, a laboratory testing program, and advanced numerical modeling to identify the cause of two shallow slope failures. A residual soil slope in Singapore was equipped with pore water pressure, water content, and rainfall measuring devices, where studies were carried out under natural and simulated rainfalls [3]. Timpong et al. [4] described the development of a new in-flight ground water table control system in a centrifuge, which can be used to investigate the mechanism of slope failure induced by ground water table changes. Ng et al. [5] investigated the effects of pole transpiration on rainfall-induced slope hydrology through centrifuge model tests by using real branch cuttings. Jeong et al. [6] investigated rainfall-induced landslides on partially saturated soil slopes using the 2011 Umyeonsan landslides at the center of Seoul, Korea and then carried out an integrated analysis of rainfall-induced landslides through laboratory tests, field tests, and numerical analysis. Experimental works and numerical analysis were conducted to determine the critical condition of the slope stability due to the evolution of shear strength parameters [7]. Ismail et al. [8] described a study in which simulated rainfall events were used with a two-dimensional soil column to study the response of unsaturated soil behavior based on different slope angles.

In order to deal with complex initial and boundary conditions, numerical methods have been applied more often in stress and deformation analyses of unsaturated soil slopes [9,10,11,12,13,14,15]. To investigate the influence of various rainfall events and initial ground conditions on transient seepage and, hence, slope stability, a parametric study was carried out by Ng et al. [16]. Numerical models were used to study how infiltration into a slope varied with respect to rainfall intensity and how this infiltration affected the stability of the slope [17]. A safety factor was calculated based on the smoothed stress field obtained from finite element analysis, and an optimization technique was used to search for a critical slip surface [18]. Shen [19] used the simplified consolidation theory of unsaturated soil to simulate the development of pore pressure and deformation of a canal slope excavated in unsaturated expansive soil in Zaoyang City during artificial rainfalls. Based on the theory of unsaturated soil, the seepage and stability of a slope under rainfall infiltration were studied with the finite element method [20]. A multiphase coupled elasto-viscoplastic finite element analysis formulation, based on the theory of porous media, was used to describe the rainfall infiltration process into a one-dimensional soil column [21]. Zhan et al. [22] developed an analytical solution for simulating rainfall infiltration into an infinite unsaturated soil slope based on the general partial differential equation for water flow through unsaturated soils. Liu et al. [23] used a coupled elastoplastic finite element analysis based on simplified consolidation theory for unsaturated soils to investigate the coupling processes of water infiltration and deformation. A numerical analysis of the relationship between three rainfall patterns and the anisotropic ratios was designed to investigate the effects of the anisotropic ratio on the stability of slopes using the reliability index approach [24]. A novel approach developed by Arairo et al. [25] to predict the behavior of unsaturated soils was discussed to investigate the effect of rainfall events on the stability of soil slopes. With a suction-stress-based effective stress representation, Vahedifard et al. [26] performed a stability analysis of unsaturated engineered and natural slopes effectively in the same manner as the classical limit equilibrium (LE) methodologies. Kim and Jeong [27] presented a numerical investigation to study the hydromechanical response of a shallow landslide in unsaturated slopes subjected to rainfall infiltration using a coupled model. Gofar and Rahardjo [28] presented results of saturated and unsaturated stability analyses of typical residual slopes subjected to rainfall infiltration, which corresponded to a 50 year rainfall return period.

Some researchers [29,30] have studied crack-containing unsaturated soil slopes. Finite element simulation was used to analyze the influences of cracks’ position, depth, and seepage characteristics on a slope’s rainfall infiltration [31]. Considering the cracking of expansive soils, the seepage characteristics of expansive soil slope were studied under the circumstance of rainfall infiltration [32]. Based on the generalized consolidation theory of unsaturated soils, coupled analysis of deformation of the soil skeleton and the movement of pore water and pore air in stiff fissured clay and swelling soil were performed for a selected slope [33]. Liu et al. [34] proposed a new analytical model for describing the progressive slow movement of natural slopes based on the relationship between velocity and the viscoplastic strain rate. Wang [35] derived the formulae of secondary crack spacing and secondary trend crack spacing after stress analysis. Li [36] investigated crack development, characterized crack geometrical parameters under natural atmosphere conditions by field tests, and then studied the permeability tensor and REV (representative element volume) for saturated soils containing random crack networks through numerical simulation. The dynamic development of cracks in expansive soil during drying and wetting has been measured in the laboratory to study the hydraulic properties of cracks by Cao et al. [37].

From the abovementioned work, we can see that though some in situ experiments and numerical analyses have been carried out on the coupling process of stress and deformation of unsaturated soil slopes, few studies have been done on the influence of crack characteristics on stress and deformation, which can affect the stability of unsaturated soil slopes and was investigated here. In this study, the effects of crack characteristics, including the location of the crack, discharge speed, evaporation rate, infiltration rate, and tensile strength, on the stress and deformation of unsaturated soil slopes were explored numerically by using a coupled finite element method developed by the authors.

2. Materials and Methods

Based on the simplified consolidation theory of unsaturated soil proposed by Shen [38], the coupling elastoplastic finite element analysis method is used to analyze the seepage–deformation coupling process of unsaturated soils here. By using the reduced suction and elastoplastic constitutive equations, the simplified consolidation theory of unsaturated soils is incorporated into the finite element program in Fortran. And the development of pore water pressure and deformation under evaporation and rainfall infiltration conditions can be calculated through this numerical method.

2.1. Simplified Consolidation Theory for Unsaturated Soils

The detailed calculation principle can be found in the works of Shen [38] and Liu and Xing [39].

2.1.1. Effective Stress Formula

Bishop’s formula was used to represent the principle of effective stress for unsaturated soils:

$${\sigma }^{\prime }=\sigma -u_a+\chi \left(u_a-u_w\right)$$

(1)

where ${\sigma }^{\prime }$ is the effective stress, $\sigma$ is the total stress, $u_a$ is the pore air pressure, $u_w$ is the pore water pressure, the term $s=u_a-u_w$ is called the matrix suction, and $\chi$ is called the coefficient of reduced suction. The coefficient of reduced suction varies with the matrix suction in the manner proposed here as follows:

$$\chi ={(\frac{s}{s_e})}^{-m_2}$$

(2)

where $s_e$ is the air entry value, and $m_2$ is a material constant.

2.1.2. Pore Air Pressure Assumption

We defined

$$n_a=\left[1-\left(1-c_h\right)S_r\right]n$$

(3)

as the ratio of the pore air of a unit soil element to signify the air content in the pores of a soil element within a unit volume, where $n$ is the porosity, $c_h$ is the Henry coefficient of solubility, and $S_r$ is the degree of saturation of the soil element. Based on Boyle’s law, ${\rho }_a={\rho }_{a0}(1+u_a/p_a)$, and the air pressure in the pores can be obtained:

$$u_a=\left(\frac{n_{a0}}{n_a}-1\right)p_a,$$

(4)

where $n_{a0}=\left[1-\left(1-c_h\right)S_r\right]n_0$ is the initial ratio of the pore air of a unit soil element, and $p_a$ is the atmospheric pressure. Under the conditions that the air in the pores is discharged partially, it was assumed that the mass of air discharged per unit time is $\Delta q_a$, and we defined the discharge speed of pore air as follows:

$$\xi =\frac{\Delta q_a}{{\rho }_a\Delta n_a}\left(1-\xi \right)\Delta n_a.$$

(5)

So, we can obtain the expression for the increment of pore air pressure as follows:

$$\Delta u_a=-\frac{p_a+u_a}{n_a}\left(1-\xi \right)\Delta n_a.$$

(6)

If $\xi$ is constant, by integrating Equation (6), we have

$$u_a=\left[{\left(\frac{n_{a0}}{n_a}\right)}^{\left(1-\xi \right)}-1\right]p_a.$$

(7)

When the air cannot discharge, $\xi =0$, and Equation (7) can be reduced to Equation (4); when $\xi =1$, $u_a=0$, which corresponds to the conditions that air is discharged completely.

2.1.3. Governing Equations

If the effect of temperature is not considered and the flow of dissolved air in the pore water and vapor in the pore air is ignored, the consolidation equations for unsaturated soils are as follows:

• Equilibrium equations:

  $$\left[L\right]\left\{\Delta \sigma \right\}+\left\{\Delta b\right\}=0;$$

  (8)
• Continuous equations of pore water:

  $$\frac{\partial }{\partial t}\left(S_{r}n\right)=\mathrm{div}\left[k_w\mathrm{grad}\left(\frac{u_w}{{\rho }_{w}g}+z\right)\right];$$

  (9)
• Continuous equations of pore air:

  $$\frac{\partial }{\partial t}\left[{\rho }_a\left(1-S_r\right)n+{\rho }_a{c}_h{S}_{r}n\right]=\mathrm{div}\left[{\rho }_a{k}_a\mathrm{grad}\left(\frac{u_a}{g}\right)\right];$$

  (10)
• The relationship of effective stress displacement:

  $$\left\{\Delta {\sigma }^{\prime }\right\}=\left[D\right]{\left[L\right]}^T\left\{\Delta U\right\};$$

  (11)
• The relationship of the saturation matrix suction:

  $$S_r=f_r\left(s\right);$$

  (12)
• The coefficient of permeability of the pore water:

  $$k_w=f_w\left(s\right);$$

  (13)
• The coefficient of permeability of pore air:

  $$k_a=f_a\left(s\right)$$

  (14)

  where $\left[L\right]=\left[\begin{array}{cccccc}\frac{\partial }{\partial x}& 0& 0& \frac{\partial }{\partial y}& 0& \frac{\partial }{\partial z}\\ 0& \frac{\partial }{\partial y}& 0& \frac{\partial }{\partial x}& \frac{\partial }{\partial z}& 0\\ 0& 0& \frac{\partial }{\partial z}& 0& \frac{\partial }{\partial y}& \frac{\partial }{\partial x}\end{array}\right]$, $\left\{\Delta \sigma \right\}$ is the increment of total stress, $\left\{\Delta b\right\}$ represents the increments of effective stress, $n$ is the porosity, ${\rho }_w$ is the density of pore water, ${\rho }_a$ is the density of air in the pore, $g$ is the gravitational acceleration, $\left\{\Delta {\sigma }^{\prime }\right\}$ is the increment of effective stress, $\left[D\right]$ is the matrix of stress–strain, and $\left\{\Delta U\right\}$ is the increment of displacement.

2.1.4. Simplified Consolidation Equations for Unsaturated Soils in 2D

In consideration of $\Delta {\epsilon }_v=-\Delta n,$ $\Delta n_a=\frac{\partial n_a}{\partial S_r}\frac{\partial S_r}{\partial s}\left(\Delta u_a-\Delta u_w\right)-\frac{\partial n_a}{\partial n}\Delta {\epsilon }_v,$ the increment of total stress from Equation (1) can be described as follows:

$$\Delta \sigma =\Delta {\sigma }^{\prime }+A_1\Delta u_w+A_2\Delta {\epsilon }_v$$

(15)

where $A_1=\frac{\chi +s\frac{\partial \chi }{\partial s}+P\frac{\partial n_a}{\partial S_r}\frac{\partial S_r}{\partial s}}{1+P\frac{\partial n_a}{\partial S_r}\frac{\partial S_r}{\partial s}}$ and $A_2=\frac{\left(\chi +s\frac{\partial \chi }{\partial s}-1\right)P\frac{\partial n_a}{\partial n}}{1+P\frac{\partial n_a}{\partial S_r}\frac{\partial S_r}{\partial s}},$ $P=\frac{\left(1-\xi \right)\left(p_a+u_a\right)}{n_a}, \frac{\partial n_a}{\partial S_r}=-\left(1-c_h\right)n, b \frac{\partial n_a}{\partial n}=-\left(1-c_h\right)S_r$. Substituting the above equations into Equation (8), we have:

$$\left(d_{11}+A_2\right)\frac{{\partial }^2\Delta u_x}{\partial x^2}+\left(d_{14}+d_{41}\right)\frac{{\partial }^2\Delta u_x}{\partial x\partial z}+d_{44}\frac{{\partial }^2\Delta u_x}{\partial z^2}+d_{14}\frac{{\partial }^2\Delta u_z}{\partial x^2}+\left(d_{12}+d_{44}+A_2\right)\frac{{\partial }^2\Delta u_z}{\partial x\partial z}+d_{42}\frac{{\partial }^2\Delta u_z}{\partial x\partial z}-A_1\frac{\partial u}{\partial x}=\Delta F_x,$$

(16)

$$d_{41}\frac{{\partial }^2\Delta u_x}{\partial x^2}+\left(d_{21}+d_{44}+A_2\right)\frac{{\partial }^2\Delta u_x}{\partial x\partial z}+d_{24}\frac{{\partial }^2\Delta u_x}{\partial z^2}+d_{44}\frac{{\partial }^2\Delta u_z}{\partial x^2}+\left(d_{24}+d_{42}\right)\frac{{\partial }^2\Delta u_z}{\partial x\partial z}+\left(d_{22}+A_2\right)\frac{{\partial }^2\Delta u_z}{\partial z^2}-A_1\frac{\partial u_w}{\partial z}=\Delta F_z$$

(17)

where $\Delta u_x$ is the increment of the horizontal displacement, $\Delta u_z$ is the increment of the vertical displacement, $\Delta F_x$ are the increments of loads in the horizontal direction, $\Delta F_z$ are the increments of loads in the vertical direction, and $d_{11},$ $d_{12}$ … are the elements of the elastoplastic matrix of stress–strain. The continuous equation of pore water is formulated as follows:

$$\mu n\frac{\partial u_w}{\partial t}=-\frac{\partial }{\partial x}{k}_{wx}\frac{\partial h}{\partial x}-\frac{\partial }{\partial z}{k}_{wz}\frac{\partial h}{\partial z}+S_r\frac{\partial {\epsilon }_v}{\partial t}$$

(18)

where $h=u_w/{\rho }_{w}g+z,$ $k_{wx}$ is the coefficient of permeability in the horizontal direction and $k_{wz}$ is the coefficient of permeability in the vertical direction, $\mu =\partial S_r/\partial u_w.$

2.2. Constitutive Equations

2.2.1. Double-Hardening Model for Soil Skeleton

Regardless of the influence of temperature, the double-hardening elastoplastic model for saturated soils [39] was used to describe the mechanical features of the soil skeleton of unsaturated soils.

Set ${{\sigma }_m}^{\prime }=\frac{1}{3}{{\sigma }^{\prime }}_{kk},$ ${{\sigma }^{\prime }}_s=\sqrt{\frac{3}{2}{s^{\prime }}_{ij}{s^{\prime }}_{ij}},$ ${s^{\prime }}_{ij}={{\sigma }^{\prime }}_{ij}-{{\sigma }^{\prime }}_{kk}{\delta }_{ij},$ ${\epsilon }_s=\sqrt{\frac{2}{3}{e}_{ij}{e}_{ij}},$ $e_{ij}={\epsilon }_{ij}-\frac{1}{3}{\epsilon }_{kk}{\delta }_{ij}$. The yield function of the model is expressed as follows:

$$F\left({\sigma }^{\prime },{\epsilon }_v^p,{\epsilon }_s^p\right)=\frac{{{\sigma }^{\prime }}_m}{1-{\left[\frac{\eta }{\alpha \left({\epsilon }_s^p\right)}\right]}^m}-p\left({\epsilon }_v^p\right)$$

(19)

where $\eta ={{\sigma }^{\prime }}_s/{{\sigma }^{\prime }}_m$, and $m$ is the parameter of yield function; when $m=1.2,$ the shape of the yield surface is close to an ellipse. $p$ and $\alpha$ are the two hardening parameters, which evolve with plastic volumetric strain ${\epsilon }_v^p$ and plastic shear strain ${\epsilon }_s^p$, respectively, as follows:

$$p=p_0\exp \left(\frac{{\epsilon }_v^p}{c_c-c_e}\right),$$

(20)

$$\alpha ={\alpha }_m-\left({\alpha }_m-{\alpha }_0\right)\exp \left(\frac{{\epsilon }_s^p}{c_a}\right)$$

(21)

where ${\epsilon }_v^p={\epsilon }_{kk}^p,$ ${\epsilon }_s^p=\sqrt{\frac{2}{3}{e}_{ij}^p}{e}_{ij}^p,$ $e_{ij}^p={\epsilon }_{ij}^p-\frac{1}{3}{\epsilon }_{kk}^p{\delta }_{ij},$ $c_c$ is the slope of the compressional curve, $c_e$ is the slope of the rebound curve, and $p_0$ is the reference pressure when ${\epsilon }_v^p=0$. Equation (20) is in the same form as the hardening parameter of the original Cam-clay model. In Equation (21), ${\alpha }_m=\left(\sqrt[m]{1+m}\right)\sin \phi ,$ where $\phi$ is internal frictional angle, and ${\alpha }_0$ and $c_a$ are two other parameters that can be determined by an unloading triaxial compression test, in which the axial load is kept constant and the confining pressure is reduced gradually.

Assuming that the flow rule is associated, the plastic strain increment can be determined by the use of elastic-plastic theory as follows:

$$d{\epsilon }_p^{ij}=d\lambda \frac{\partial F}{\partial {{\sigma }^{\prime }}_{ij}}$$

(22)

or

$$d{\epsilon }_v^p=d\lambda \frac{\partial F}{\partial {{\sigma }^{\prime }}_m},$$

(23)

$$d{\epsilon }_s^p=\frac{3}{2}d\lambda \frac{\partial F}{\partial {{\sigma }^{\prime }}_s}$$

(24)

where $d\lambda$ is the plastic multiplier, which can be derived from the consistency conditions:

$$\frac{\partial F}{\partial {{\sigma }^{\prime }}_{ij}}d{{\sigma }^{\prime }}_{ij}+\frac{\partial F}{\partial {\epsilon }_v^p}d{\epsilon }_v^p+\frac{\partial F}{\partial {\epsilon }_s^p}d{\epsilon }_s^p=0.$$

(25)

Substituting for the plastic volumetric strain increment $d{\epsilon }_v^p$ and the plastic shear strain increment $d{\epsilon }_s^p$ in Equations (23) and (24), $d\lambda$ is obtained as follows:

$$d\lambda =\frac{\frac{\partial F}{\partial {{\sigma }^{\prime }}_{ij}}d{{\sigma }^{\prime }}_{ij}}{H}$$

(26)

where the hardening modulus $H$ is

$$H=-\frac{3}{2}\frac{\partial F}{\partial {\epsilon }_s^p}\frac{\partial F}{\partial {{\sigma }^{\prime }}_s}-\frac{\partial F}{\partial {\epsilon }_v^p}\frac{\partial F}{\partial {{\sigma }^{\prime }}_m}=-\frac{3}{2}\frac{\partial F}{\partial \alpha }\frac{\partial \alpha }{\partial {\epsilon }_s^p}\frac{\partial F}{\partial {{\sigma }^{\prime }}_s}-\frac{\partial F}{\partial p}\frac{\partial p}{\partial {\epsilon }_v^p}\frac{\partial F}{\partial {{\sigma }^{\prime }}_m}.$$

(27)

2.2.2. Soil-Water Characteristic Curve

The soil-water characteristic curve is divided into two sections. When the value of matrix suctions is smaller than that of the air entry suction $s_e,$ the soil can be assumed to be quasi-saturated and the degree of saturation of quasi-saturated soil is assumed to be $S_{r1},$ so the degree of saturation can be expressed by using the Hilf formulation.

When the value of matrix suction s is greater than that of the air entry suction $s_e,$ the degree of saturation is computed as follows [40]:

$$S_r=S_{r0}+\left(S_{r1}-S_{r0}\right){\left(\frac{s}{s_e}\right)}^{-m_1}.$$

(28)

During the process of drying shrinkage or absorbing water, the parameters $S_{r0}$, $S_{r1},$ and $s_e$ may be different. By the derivation of Equation (28), we have

$$\mu =\left(S_{r1}-S_{r0}\right)\frac{m_1}{s}{\left(\frac{s}{s_e}\right)}^{-m_1}.$$

(29)

The coefficient of permeability of unsaturated soils is calculated as follows:

$$k_w=k_{ws}\exp \left(-c_k\frac{s-s_e}{p_a}\right)$$

(30)

where $k_{ws}$ is the coefficient of permeability of saturated soils, and $c_k$ is constant. When $s<s_e,$ $k_w=k_{ws}.$

2.3. Formulations of the Finite Element Equations

Isoperimetric elements were implemented with eight-node interpolating functions for the displacements and four-node interpolating functions for the pore water pressures, which can be expressed as follows:

$$u_x={{\displaystyle \sum }}_{i=1}^8{N}_i{u}_{xi},$$

(31)

$$u_z={{\displaystyle \sum }}_{i=1}^8{N}_i{u}_{zi},$$

(32)

$$u_w={{\displaystyle \sum }}_{i=1}^4\overline{N_i}{u}_{wi}$$

(33)

where $u_{xi},$ $u_{zi}$, and $u_{wi}$ are the nodal variables at nodal point $i.$ The weak forms of Equations (31)–(33) are discretized in space and solved by the finite element method as follows:

$${{\displaystyle \sum }}_{j=1}^{N_i}\left[k_{ij}^{11}\Delta u_{xj}+k_{ij}^{12}\Delta u_{zj}+k_{ij}^{13}\Delta {\overline{h}}_j\right]=\Delta F_{xi},$$

(34)

$${{\displaystyle \sum }}_{j=1}^{N_i}\left[k_{ij}^{21}\Delta u_{xj}+k_{ij}^{22}\Delta u_{zj}+k_{ij}^{23}\Delta h_j\right]=\Delta F_{zi},$$

(35)

$${{\displaystyle \sum }}_{j=1}^{N_i}\left[k_{ij}^{31}\Delta u_{xj}+k_{ij}^{32}\Delta u_{zj}+k_{ij}^{33}\left(h_{j0}+\beta \Delta h_j\right)+s_{ij}\Delta h_j\right]=\Delta Q_i$$

(36)

where $N_i$ is the total number of nodal points; $\Delta F_{xi},$ $\Delta F_{zi}$, and $\Delta Q_i$ are the load increments and flux increment at node $i,$ respectively; $h_{i0}$ is the initial value of the water head at node $i,$ and $\beta$ is an integral parameter ($\beta =2/3$ was used here). The coefficients in the upper equations are as follows:

$$k_{ij}^{11}={\displaystyle \int }\left[\left(d_{11}+A_2\right)\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+d_{44}\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial z}+d_{14}\left(\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial z}+\frac{\partial N_j}{\partial x}\frac{\partial N_i}{\partial z}\right)\right]dxdz,$$

$$k_{ij}^{12}={\displaystyle \int }\left[\left(d_{12}+A_2\right)\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial z}+d_{14}\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+d_{24}\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z}+d_{44}\frac{\partial N_j}{\partial x}\frac{\partial N_i}{\partial z}\right]dxdz,$$

$$k_{ij}^{21}={\displaystyle \int }\left[\left(d_{12}+A_2\right)\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial z}+d_{14}\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+d_{24}\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z}+d_{44}\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial z}\right]dxdz,$$

$$k_{ij}^{22}={\displaystyle \int }\left[\left(d_{22}+A_2\right)\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z}+d_{44}\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+d_{24}\left(\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z}+\frac{\partial N_j}{\partial z}\frac{\partial N_i}{\partial x}\right)\right]dxdz,$$

$$k_{ij}^{33}=-{\displaystyle \int }\left[k_x\frac{\partial N_i}{\partial x}\frac{\partial N_j}{\partial x}+k_z\frac{\partial N_i}{\partial z}\frac{\partial N_j}{\partial z}\right]dxdz,k_{ij}^{13}=-{\rho }_{w}g{\displaystyle \int }{A}_1\frac{\partial \overline{N_i}}{\partial r}\overline{N_j}dzdz,$$

$$k_{ij}^{23}=-{\rho }_{w}g{\displaystyle \int }{A}_1\frac{\partial \overline{N_i}}{\partial z}\overline{N_j}dxdz,k_{ij}^{31}=-{\rho }_{w}g{\displaystyle \int }{S}_r\frac{\partial \overline{N_j}}{\partial x}\overline{N_i}dxdz,k_{ij}^{32}=-{\rho }_{w}g{\displaystyle \int }{S}_r\frac{\partial \overline{N_j}}{\partial z}\overline{N_i}dxdz,$$

and

$$s_{ij}=-{\rho }_{w}g{\displaystyle \int }\left[c_s{N}_i{N}_j\right]dxdz$$

where $c_s=m_{f}n$ (saturated) or $c_s=\frac{\mu n}{S_r}$ (unsaturated).

2.4. Computational Model

The seepage–deformation processes under evaporation and rainfall infiltration conditions were simulated for an unsaturated soil slope which was 20 m in depth and 54 m in width using the coupled elastoplastic finite element program in Fortran. The computational mesh and boundary conditions are shown in Figure 1.

For the lateral surfaces, there were undrained and constrained boundaries. For the bottom surface, also the surface of the groundwater table, there was a constrained boundary; the pore water pressure was zero all the time with the surface of the groundwater table. For the upper surfaces composed of three surfaces, it drained freely, and evaporation/rainfall occurred on these surfaces.

When t = 0, the soil slope was saturated and in equilibrium with the weight stress state. The calculation simulated a long geological process without ground pressure and formed a soil slope through natural erosion or manual excavation. The excavation was carried out in eight steps, each of which was 1 billion h; in each row was dug a row of units, and the groundwater level was reduced. Evaporation occurred 1100 days after the completion of the excavation, which was carried out in 22 steps, each of which was 50 days. The next 300 days of rainfall were divided into 30 steps, and each step was 10 days. The computed parameters were as follows: $\gamma$ = 20 kN/m³, $k_0$ = 0.6 cm/s, $\nu$ = 0.3077, $k_{ws}$ = 0.001 cm/s, $c_k$ = 0.2, $c_c$ = 0.0332, $c_e$ = 0.0064, $n_0$ = 0.412, $S_{r1}$ = 0.7, $S_{r0}$ = 0.96, $m_1$ = $m_2$ = 0.1, $s_e$ = 5 kPa, ${\alpha }_m$ = 1.0, ${\alpha }_0$ = 0.75, $c_a$ = 0.05, $\xi$ = 0.6, ${\sigma }_b$ = 5 kPa, the evaporation rate was 0.3 mm/day, and the infiltration rate was 0.5 mm/day. The pore water pressure distribution and displacement distribution at the end of evaporation without cracks were obtained as shown in Figure 2. The pore water pressure distribution and displacement distribution at the end of rainfall are shown in Figure 3.

During the simulation, we observed that evaporation led to a decrease in the pore water pressure, which was most pronounced at the right side of the top of the slope. After the completion of the rainfall, the pore water pressure increased greatly, and the maximum displacement of the soil slope decreased. From the simulation results of evaporation and rainfall infiltration of this noncracked unsaturated soil slope, we found that the above method can be used to simulate the seepage–deformation coupling process of unsaturated soils.

Simulation of Crack Propagation

A dual-node technique was used to simulate crack propagation [41]. The point where a crack may show was set as two separate nodes stuck together. The coordinates of these two points were the same but belonged to neighboring elements. Before crack propagation, the two nodes had the same degree of freedom. Therefore, the two nodes were regarded as a single node. When the total horizontal stress satisfied the cracking condition (over the tensile strength), the two nodes separated and had their own degree of freedom. Once the cracked section became free face, evaporation took place on that surface, and horizontal stress of other nodes at the bottom of the two cracked nodes progressively increased to reach the tensile strength. Therefore, these bottom nodes separated from one another eventually. The crack propagation would be restrained if the crack moved near the groundwater table.

The soil slope with a crack located at a distance of 1400 cm from the left edge of the slope and the same calculated parameters as the noncracked slope were used as the baseline group, with variations of some parameters (listed in

Table 1. Parameters changed in subsequent analysis.

Name	Symbol	Unit
Crack location	\	cm
Discharge speed of pore air	$\xi$	\
Evaporation rate	\	mm/day
Infiltration rate	\	mm/day
Tensile strength	${\sigma }_b$	kPa

) for the subsequent analysis.

3. Results and Discussion

3.1. Influence of Crack Location

The pore water pressure distribution and displacement distribution at the end of evaporation of the baseline group were obtained as shown in Figure 4. The pore water pressure distribution and displacement distribution at the end of rainfall are shown in Figure 5.

Figure 6 shows the pore water pressure distribution and displacement distribution at the end of evaporation when the crack was located 600 cm from the left edge of the soil slope. The pore water pressure distribution and displacement distribution at the end of rainfall are shown in Figure 7.

From the above figures, it can be seen that after the completion of evaporation, the appearance of the crack led to the discontinuity of pore water pressure distribution, the pore water pressure decreased, and then the pore water pressure rose significantly after the completion of rainfall and changed the most at the right side of the crack. It also caused the reduction of the influence on the right side of the slope. Then, the slopes with cracks at 1000 and 1800 cm from the left edge of the soil slope were analyzed as well. We found that the minimum pore water pressure of the soil slope decreased gradually and the crack width increased gradually as the crack location moved to the right side of the slope.

The curve of the crack width over time with cracks located at different locations is shown in Figure 8.

Figure 8 shows that the position of the crack had no influence on the development trend of the crack width. The cracks all appeared after 100 days of evaporation. When the evaporation developed, the crack developed quickly and then slowly. As for rain, the crack width slowly decreased and then slowly increased, and the maximum crack width appeared when the rainfall ended. Further, with the rightward movement of the crack, the crack width gradually increased, because the more the crack position was on the right slope, the closer the crack was to the slope surface, and the smaller the constraints were on the soil body.

3.2. Influence of Discharge Speed of Pore Air

The numerical analysis of the unsaturated soil slopes with the discharge speeds of pore air of 0.2, 0.4, 0.8, and 1.0 respectively showed that there was not much difference between stress and deformation under each discharge speed, so they are not individually listed. Figure 9 shows the variation of crack width with time under different discharge speeds, and the crack-width–time curve under different discharge speeds after 800 days is shown in Figure 10.

As can be seen from Figure 9 and Figure 10, the greater the discharge speed, the greater the crack width, but the overall difference is not obvious. This is because the increase of the discharge speed caused the pore air to be expelled more quickly. The drier the soil was, the easier it was for the stress to reach the tensile strength, resulting in an increase in the crack width.

3.3. Influence of Evaporation Rate

Figure 11 shows the pore water pressure and displacement distribution at the end of evaporation with an evaporation rate of 0.1 mm/day. Figure 12 shows the pore water pressure and displacement distribution at the end of rainfall with an evaporation rate of 0.1 mm/day.

Comparing Figure 11 and Figure 12 with Figure 4 and Figure 5, it is obvious that the negative pore water pressure increased as the evaporation rate changed from 0.3 to 0.1 mm/day, and the total displacement of the slope decreased.

Then, we analyzed the condition of the evaporation rate being 0.5 mm/day. The variation of crack width with time under different evaporation rates is shown in Figure 13.

As can be seen in Figure 13, the greater the evaporation rate, the wider the crack, and the earlier the crack appeared. This is because when the evaporation rate was high, the pore water was reduced, so the negative pore water pressure decreased and the stress reached tensile stress quickly and finally made the crack width increase.

3.4. Influence of Infiltration Rate

For the numerical analysis of unsaturated soil slopes under different infiltration rate conditions, since the evaporation rate had not changed, the pore water pressure and the displacement distribution were the same as when the evaporation ended.

Figure 14 shows the pore water pressure distribution and displacement distribution at the end of rainfall with an infiltration rate of 0.1 mm/day.

Figure 15 shows the pore water pressure distribution and displacement distribution at the end of rainfall with an infiltration rate of 0.8 mm/day.

From the above figures, we can see that with the increase of the infiltration rate, the pore water pressure and displacement increased at the end of rainfall, and the change rules were similar.

The crack-width–time curve during the rainfall period is shown in Figure 16.

From Figure 16, we can see that the three lines show a decreasing trend at the initial phase but then increase till the end. As the infiltration rate increased, the crack width increased, the decreasing period became shorter, and the subsequent increasing rate increased. This is because at the beginning of the rain, the water infiltrated into the pore of the soils and slightly increased the tensile stress; but with the water infiltrating, the pore water stress increased, and a sustainable crack developed.

3.5. Influence of Tensile Strength

A numerical analysis was carried out on unsaturated soil slopes with tensile strengths of 1, 3, and 5 kPa (e.g., [42,43,44]).

Figure 17 shows the pore water pressure and displacement distributions at the end of evaporation under the condition of 1 kPa tensile strength, and Figure 18 shows the pore water pressure and displacement distributions at the end of rainfall under the condition of 1 kPa tensile strength.

Figure 17 and Figure 18 show that the larger the tensile strength, the smaller the pore water pressure. After the rainfall had completed, the pore water pressure around the crack increased greatly.

The crack-width–time curve under different tensile strength conditions is shown in Figure 19, and the crack-width–time curve under different tensile strength conditions after 800 days is shown in Figure 20.

From Figure 19 and Figure 20, we can see that when the evaporation developed, the crack width under 5 kPa of tensile stress was bigger than that under 3 kPa of tensile stress, and the crack width under 3 kPa was bigger than that under 1 kPa but almost the same. When it began raining, the crack width decreased as the tensile stress increased. This shows that increasing the tensile strength has a certain effect on reducing the crack width.

4. Conclusions

In this study, the finite element method was used to numerically analyze the development of cracks in unsaturated soil slopes, and the effect of crack characteristics on the stress and deformation of unsaturated soil slopes was studied. The main conclusions are as follows:

The finite element method considering the effective stress for unsaturated soil proposed by Bishop and an elastoplastic double-hardening constitutive model for the soil skeleton can be used to simulate the formation and propagation of cracks of unsaturated soil slopes.

The pore water pressure decreased and the crack width increased gradually with the right position of the crack. The larger the discharge speed, the greater the crack width. The negative pore water pressure increased as the evaporation rate increased, the total displacement of the slope decreased, and the crack widened and appeared earlier.

With the increase of the infiltration rate, the pore water pressure increased at the end of rainfall. The crack width at the rainfall stage firstly decreased and then increased, and as the infiltration rate increased, the decreasing stage became shorter and the crack width increased. The larger the tensile strength, the smaller the pore water pressure.