**1. Introduction**

Slope stability is an important engineering problem in the geotechnical field, for example, the slope stability of the excavation of a foundation pit [1], the regulation of a riparian slope [2], and the stability of a high slope under complex geological environment in a reservoir area [3], etc. Slope failures mainly include external and internal causes. The external causes include rainfall [4–6], reservoir water level fluctuations [7–9], earthquakes [10–12], human activities such as excavation or blasting [1,3–15], etc. The internal causes are mainly affected by the properties of the slope soil, which include soil types [16–18], unsaturated characteristics [19–21], soil strength [22–24], etc. Rainfall is the key factor triggering the landslide, which accounts for 51% of all the landslide disasters, according to some relevant investigations [25]. Rainfall usually occurs in rainy seasons, which is concentrated and has a long duration and leads to landslide easily. The main reasons for slope instability caused by rainfall are as follows: (1) Rainfall increases the groundwater level inside the slope, which reduces the effective stress and shear strength of the soil. (2) Rainfall increases the slide force of the slope, which aggravates the slope instability. The main consequence for slope instability caused by rainfall is reflected in these

two aspects: (1) Landslide causes damage to buildings in the area where the slope is located. (2) Landslide threatens the lives and properties of the residents around the disaster area. For example, on 16 September 2011, heavy rainfall (250 mm/day, cumulative 430 mm) in Nanjiang County, Sichuan Province, China, induced thousands of landslides, and on 29 April 2015, the groundwater level in Dangchuan, Heifangtai, Gansu Province, China, rose due to heavy rainfall, which caused a large-scale slope failure [26]. Therefore, it is important to grasp the law and influencing factors of slope instability caused by rainfall in order to correctly understand the mechanism of rainfall infiltration and to prevent and control landslide disasters.

Scholars have conducted large numbers of research studies on the rainfall infiltration mechanism. The research results mainly focused on theory, experiments, and numerical simulations. With respect to theoretical research, Green-Ampt Semi-Analytical Method [27] was the first method to describe the transient infiltration process of rain water in unsaturated soils. It is assumed that the wetting front moves down along the depth direction and the velocity remains unchanged. The volume of water content of the soil after the wetting front is θ0 (completely saturated), and the initial water content of the soil before the front is θ*i*. The rainfall intensity is always higher than the infiltration capacity of the soil, and the actual infiltration rate is equal to the infiltration capacity of the soil. Mein et al. [28] improved the Green-Ampt model and rainfall infiltration process was divided into two periods. The first stage is the free infiltration period, and the rainfall intensity in this stage is less than the soil infiltration capacity. The second stage is the ponding infiltration period, and the rainfall intensity is greater than the soil infiltration capacity. Chu et al. [29] divided the non-uniform rainfall process into several uniform periods, and calculated the infiltration process according to whether there will be ponding. Chen et al. [30] derived a uniform slope rainfall infiltration model based on the Green-Ampt model. For experimental research, Wu et al. [31] carried out the laboratory model test of landslide under artificial rainfall, and the influence of rain water infiltration on the slope failure was analyzed. Li et al. [32] studied the influence of rainfall on the internal mechanical response characteristics of slope based on fiber grating monitoring technology. Zhang et al. [33] systematically investigated the stability of Xiakou slope under rainfall based on the field monitoring data. In the aspect of numerical simulation, Hao et al. [34] simulated the variations of stability of a typical slope under rainfall based on the limit equilibrium method. Wang et al. [35] used XFEM (extended finite element method) to simulate the crack propagation process in a slope under heavy rain. Dou et al. [36] considered the spatial variability of hydraulic conductivity of slope soil based on the Monte Carlo method, and simulated the seepage characteristics and slope stability. However, most of the previous studies regarded the slope materials as isotropic materials. According to the SEM (scanning electron microscope) microcosmic study of Song et al. [37], the anisotropy of permeability coe fficient is caused by the flocculation microstructure of clay and other soils, meanwhile, soil permeability coe fficient anisotropy is greatly a ffected by dry density and freeze-thaw cycles [38]. Generally speaking, the hydraulic conductivity anisotropy ratio (*k*x/*<sup>k</sup>*y) can reach 2–10, and it is possible to reach 100 for clayey soil [39]. There exist joint cracks in rock slope [40,41], which lead to strong anisotropy of seepage characteristics. The coe fficient anisotropy not only has grea<sup>t</sup> influences on the transient seepage but also has an impact on the safety factor of the slope. According to Mahmood et al. [42], the di fference of slope safety factors between considering and not considering soil anisotropy will be about 40% [43]. But most of the previous studies ignored the seepage anisotropy of slope soil, and the research results of soil anisotropy of slope were few and incomplete. Yeh et al. [44] took into account the hydraulic conductivity anisotropy ratio of soil slope, and simulated the seepage characteristics and local safety factors, but ignored the hydraulic conductivity anisotropy direction. In fact, the horizontal permeability coe fficient *kx* and vertical permeability coe fficient *ky* coincide with the natural coordinate axis only in some special cases such as layered crushed earth dam or naturally deposited layered soil. In nature, there are more cases where the anisotropy principal direction does not coincide with the coordinate axis. The actual conditions cannot be accurately reflected only by considering the hydraulic conductivity anisotropy ratio.

In view of the shortcomings of previous studies, the mathematical definition of hydraulic conductivity anisotropy ratio and direction are firstly described in this paper. The SEEP/W and SLOPE/W modules in Geo-studio were utilized to carry out the numerical analysis of a homogeneous slope in Luogang District, Guangzhou City, China [45]. Geo-studio is a professional software suitable for analyzing the seepage and stability of soil slopes, and the numerical results were consistent with the experimental results and field investigations. For example, Jiang et al. [46] analyzed the seepage and stability of a cracked slope with SEEP/W and SLOPE/W. Duong T.T. [2] used the Geo-studio program to study the e ffects of soil hydraulic conductivity and rainfall intensity on riverbank stability. Muqdad Al-Juboori et al. [47–49] conducted the machine learning research based on SEEP/W. Therefore, this software was utilized in this paper to analyze the e ffect of hydraulic conductivity anisotropy. Two kinds of soil (clay and sand) and the hydraulic conductivity anisotropy ratio and direction were included. Then the volume water content, rainfall infiltration depth (RID), rising height of groundwater (RHG) and maximum water content of the surface (MWCS) of three typical sections of slope (top, middle, and toe) were analyzed in detail. Finally, the safety factors (SF) of the slope under di fferent conditions were evaluated. The research results provide some references for the understanding of the seepage anisotropy law and prevention of landslides.

#### **2. Methods and Theory**

#### *2.1. Theory of Unsaturated Seepage*

The SEEP/W module in Geo-studio was utilized to simulate the rainfall infiltration, and the seepage control equation in the SEEP/W was derived from the saturated and unsaturated Darcy's law [50], which can be expressed as

$$\frac{\partial}{\partial \mathbf{x}} (k\_x \frac{\partial H}{\partial \mathbf{x}}) + \frac{\partial}{\partial y} (k\_y \frac{\partial H}{\partial y}) + Q = m\_w \gamma\_w \frac{\partial H}{\partial t}.\tag{1}$$

In Equation (1), *x* and *y* are the coordinates in the direction of *x* and *y*, *kx* is the hydraulic conductivity in the x direction, *ky* is the hydraulic conductivity in the y direction, *H* is the total head, *Q* is the applied boundary flux, *t* is the time, *mw* is the slope of the storage curve, and γ*w* is the unit weight of water.

Applying the Galerkin method of weighed residual to the governing di fferential equation, the finite element for two-dimensional seepage equation can be derived as

$$\tau \int\_{A} \left( \left[ B \right]^{T} \left[ \mathbf{C} \right] \middle| B \right] dA \{ H \} + \tau \int\_{A} \left( \lambda \left< \mathbf{N} \right>^{T} \left< \mathbf{N} \right> \right) dA \{ H \}, t = q\tau \int\_{L} \left( \left< \mathbf{N} \right>^{T} \right) dL. \tag{2}$$

In Equation (2), *[B]* is the gradient matrix, *[C]* is the element hydraulic conductivity matrix, *[H]* is the vector of nodal heads, <*N*> is the vector of interpolating function, *q* is the unit flux across the edge of an element, τ is the thickness of an element, λ is the storage term for a transient seepage equal to *mw*γ*<sup>w</sup>*, A is a designation for summation over the area of an element, and *L* is a designation for summation over the edge of an element.

In an abbreviated form, the finite element seepage equation can be expressed as

$$[K]\{H\} + [M]\{H\}, t = \langle Q \rangle. \tag{3}$$

In Equation (3), [ *K*] is the element characteristic matrix, [ *M*] is the element mass matrix, and { *Q*} is the element applied flux vector.

#### *2.2. E*ff*ect of Negative Pore-Water Pressures*

In locations above the groundwater table, the pore-water pressure in the soil is negative relative to the pore-air pressure. This negative pore-water pressure is commonly referred to as the matric suction of the soil. Under negative pore-water pressure conditions the shear strength may not change at the same rate as for total and positive pore-water pressure changes. Therefore, a modified form of the Mohr–Coulomb equation must be used to describe the shear strength of an unsaturated soil (i.e., the soil with negative pore-water pressures). The shear strength equation is [51]

$$s = c' + \sigma\_{\text{tr}} \tan q' + (u\_d - u\_{\text{w}}) \tan q\_{\text{b}}.\tag{4}$$

In Equation (4), *s* is the unsaturated shear strength, *c*' is the cohesive strength, ϕ' is the frictional strength, ϕ*b* is an angle defining the increase in shear strength for an increase in suction, *u*a is the pore-air pressure, and *u*w is the pore-water pressure.

#### *2.3. Safety Factor for Unsaturated Soil*

SLOPE/W adopts the Morgenstern–Price method based on limit equilibrium theory to calculate the safety factor. The modified method strictly satisfies the force balance and torque balance, and the calculation accuracy is high. The expression is listed below:

$$F\_{\sf S} = \frac{\sum\_{i=1}^{n\_{\sf S}} \frac{c'\_i b\_i + (\mathcal{W}\_i + P\_i \cos \beta\_i - u\_b b\_i) \tan \phi'\_i + (u\_a - u\_{\sf W}) b\_i \tan \phi\_{\sf b}}{\left[1 + \left(\tan \varphi'\_i \tan \alpha\_i\right) / F\_{\sf s}\right] \cos \alpha\_i}}{\sum\_{i=1}^{n\_{\sf S}} W\_{\sf i} \sin \alpha\_i - r\_i P\_i} \,. \tag{5}$$

In Equation (5), *ci*' is the cohesive strength for every soil slice, *i* is the soil slice number, *Wi* is the weight of every soil slice, *Pi* is the water pressure, β*i* is the angle of the bottom of the soil slice, *bi* is the length of every soil slice, ϕ*i*' is the frictional strength for every soil slice, *ri* is the radius of the sliding arc, and *F*s is the safety factor.

#### **3. Numerical Model Framework**

#### *3.1. Numerical Model and Boundary Conditions*

The case study is a homogeneous slope in Luogang District, Guangzhou City, China [46]. The slope height is 16 m, which is divided into 2 grades with a width of 2 m. In order to reduce the influence of boundary conditions, the range was extended and the model was divided into 13,484 nodes and 13,279 units, which are shown in Figure 1.

**Figure 1.** Illustration of slope model.

To investigate the seepage characteristics at different positions, three sections were set to reflect the effect of hydraulic conductivity anisotropy, whose positions were *x* = 73 m (top of the slope), *x* = 61 m (middle of the slope), and *x* = 50 m (toe of the slope). The boundary conditions were as follows: AB and GH were the fixed water level boundaries of 9 m and 24 m, respectively. CDEF was the rainfall infiltration boundary. BC and GF were small flux boundaries. AH was the impervious boundary.

#### *3.2. Unsaturated Soil Properties*

The soil-water characteristic curves (SWCC) adopted the Fredlund and Xing model, which can be written as [52]

$$k\_w = k\_s \frac{\sum\_{i=j}^{N} \frac{\Theta(\varepsilon^y) - \Theta(\Psi')}{\varepsilon^{y\_i}} \Theta'(\varepsilon^{y\_i})}{\sum\_{i=1}^{N} \frac{\Theta(\varepsilon^y) - \Theta\_s}{\varepsilon^{y\_i}} \Theta'(\varepsilon^{y\_i})}.\tag{6}$$

In Equation (6), *kw* is the calculated conductivity for a specified water content or negative pore-water pressure, *ks* is the measured saturated conductivity, Θ*s* is the volumetric water content, e is the natural number 2.71828, *y* is a dummy variable of integration representing the logarithm of negative pore-water pressure, *i* is the interval between the range of *j* to *N*, *j* is the least negative pore-water pressure to be described by the final function, *N* is the maximum negative pore-water pressure to be described by the final function, Ψ is the suction corresponding to the jth interval, Θ' is the first derivative of the equation, and Θ can be described as

$$\Theta = \mathbb{C}(\mathbb{V}) \frac{\Theta\_s}{\left\{ \ln \left[ \mathbf{e} + \left( \frac{\mathbb{V}}{a} \right)^n \right] \right\}^m} . \tag{7}$$

In Equation (7), a is the air-entry value of the soil, *n* is a parameter that controls the slope at the inflection point in the volumetric water content function, *m* is a parameter that is related to the retention capacity, and *C(*Ψ*)* is a correcting function defined as

$$C(\Psi) = 1 - \frac{\ln(1 + \frac{\mathcal{V}}{C\_r})}{\ln(1 + \frac{1000000}{C\_r})}.\tag{8}$$

In Equation (8), *Cr* is a constant related to the matric suction corresponding to the retention capacity.

The sandy soil and clayey soil were selected for analysis, which represent the high and low permeability [53], as shown in Figure 2. The unsaturated parameter values are shown in Table 1 [54], and the SWCC curves are shown in Figure 3.

**Figure 2.** Typical appearance of sandy and clayey soil. (**a**) Clayey soil. (**b**) Sandy soil.

**Table 1.** Unsaturated parameter values.

**Figure 3.** SWCC curves. (**a**) Permeability coefficient function. (**b**) Volume water content function.

#### *3.3. Definition of Anisotropy and Calculation Conditions*

Previous studies mostly ignored the anisotropy ratio and direction. In fact, anisotropy widely exists in the soil. For the hydraulic conductivity matrix *[C]* in Equation (2), it can be expressed as

$$\begin{bmatrix} \mathbf{C} \end{bmatrix} = \begin{bmatrix} \mathbf{C}\_{11} & \mathbf{C}\_{12} \\ \mathbf{C}\_{21} & \mathbf{C}\_{22} \end{bmatrix} \tag{9}$$

In Equation (9), *C*11 = *kx*cos2α + *ky*sin2α, *C*22 = *kx*sin2α + *ky*cos2α, and *C*12 = *C*21 = *kx*sinαcosα + *ky*sin<sup>α</sup>cosα. The *kx*, *ky*, the anisotropy direction α can be defined according to Figure 2. The *kx* is the horizontal permeability coefficient, *ky* is the vertical permeability coefficient, and α is the direction between *ky* and *y* axis. When α = 0◦, [*C*] is reduced to

$$[\mathbb{C}] = \begin{bmatrix} k\_{\mathbf{x}} & 0 \\ 0 & k\_{\mathbf{y}} \end{bmatrix}. \tag{10}$$

Equation (10) was adopted in [44], with only considering the anisotropy ratio *kr* = *kx*/*ky*. However, the definition of anisotropy not only includes the ratio *kr* but also the direction α. Previous investigations ignored the anisotropy, especially the anisotropy direction α.

To completely discuss the anisotropy of sandy soil and clayey soil, including the anisotropy ratio *kr* and the anisotropy direction α, the calculation conditions are shown in Table 2, which includes the anisotropy ratio *kr* = 1, 10, 50, 100, and the anisotropy direction α = 0◦, 15◦, 30◦, 45◦, 60◦, 75◦, and 90◦. The range of *kr* and α were selected according to Gilbert el al. [39]. To reflect the influence of the heavy rainfall, the rainfall intensity is set to 10−<sup>6</sup> m/s, and the rainfall duration time is set to 120 h. Meanwhile, 120 h of rainfall stop was also considered.


**Table 2.** Different calculation conditions.

#### **4. Results and Discussions**

## *4.1. Initial Conditions*

Initial conditions are important for further numerical simulations. In order to determine the initial conditions more accurately in this paper, the maximum negative pore water pressure of −25 kPa, −50 kPa, and −75 kPa and the specified annual average rainfall infiltration were numerically simulated, and the pore pressure variation of sandy soil and clayey soil are shown in Figure 4. Under the annual average rainfall infiltration, the initial pore pressure of sand and clay was obviously different in that the initial pore pressure of sandy soil slope was slightly larger than clayey soil slope, but the distribution along the elevation was similar, which was reflected in that the initial pore water pressure firstly remained unchanged then gradually increased along the elevation. The maximum negative pore-water pressure was close to −50 kPa, which was selected as the initial condition of all the calculation conditions in this paper.

**Figure 4.** Initial pore pressure distribution. (**a**) Sandy soil. (**b**) Clayey soil.

#### *4.2. Variation of Volumetric Water Content*

According to the calculation conditions in Table 2, we carried out a total number of 54 numerical simulations, and obtained, in all, 162 sections of volumetric water content variation. For ease of reading, this section will carry out the discussion based on the classification of sandy soil and clayey soil slope. The variations of volumetric water content of different α values with *kr* = 10 and *kr* = 100 are shown in Figures 5 and 6 to illustrate the influence of anisotropy direction α on seepage characteristics, and different *kr* values with α = 0◦, 45◦, and 90◦ are also displayed in Figures 7 and 8 to show the impact of anisotropy ratio *kr*. The volumetric water content of the 120th h is only shown in Figures 5–8.

**Figure 5.** Variation of volumetric water content under different α values for sandy soil. (**a**) Top of the slope with *kr* = 10. (**b**) Middle of the slope with *kr* = 10. (**c**) Toe of the slope with *kr* = 10. (**d**) Top of the slope with *kr* = 100. (**e**) Middle of the slope with *kr* = 100. (**f**) Toe of the slope with *kr* = 100.

**Figure 6.** *Cont.*

**Figure 6.** Variation of volumetric water content under different *kr* values for sandy soil. (**a**) Top of the slope with α = 0◦. (**b**) Middle of the slope with α = 0◦. (**c**) Toe of the slope with α = 0◦. (**d**) Top of the slope with α = 45◦. (**e**) Middle of the slope with α = 45◦. (**f**) Toe of the slope with α = 45◦. (**g**) Top of the slope with α = 90◦. (**h**) Middle of the slope with α = 90◦. (**i**) Toe of the slope with α = 90◦.

**Figure 7.** Variation of volumetric water content under different α values for clayey soil. (**a**) Top of the slope with *kr* = 10. (**b**) Middle of the slope with *kr* = 10. (**c**) Toe of the slope with *kr* = 10. (**d**) Top of the slope with *kr* = 100. (**e**) Middle of the slope with *kr* = 100. (**f**) Toe of the slope with *kr* = 100.

**Figure 8.** Variation of volumetric water content under different *kr* values for clayey soil. (**a**) Top of the slope with α = 0◦. (**b**) Middle of the slope with α = 0◦. (**c**) Toe of the slope with α = 0◦. (**d**) Top of the slope with α = 45◦. (**e**) Middle of the slope with α = 45◦. (**f**) Toe of the slope with α = 45◦. (**g**) Top of the slope with α = 90◦. (**h**) Middle of the slope with α = 90◦. (**i**) Toe of the slope with α = 90◦.
