Deformation and Stability Characteristics of Layered Rock Slope Affected by Rainfall Based on Anisotropy of Strength and Hydraulic Conductivity

Abstract

The strength and hydraulic conductivity anisotropy of rock slopes have a great impact on the slope stability. This study took a layered rock slope in Pulang, Southwestern China as a case study. The strength conversion equations of the seriously weathered rock mass were proposed. Then, considering the anisotropy ratio and anisotropy angle (dip angle of bedding plane) of strength and hydraulic conductivity, the deformation and stability characteristics of rock slope were calculated and compared with field monitoring data. The results showed that the sensitivity analysis of strength and hydraulic conductivity anisotropy could successfully predict the occurrence time, horizontal displacement (HD), and the scope of the rock landslide. When the anisotropy ratio was 0.01 and the dip angle was 30°, the calculated HD and scope of the landslide were consistent with the field monitoring data, which verified the feasibility of the strength conversion equations. The maximum horizontal displacement (MHD) reached the maximum value at the dip angle of 30°, and the MHD reached the minimum value at the dip angle of 60°. When the dip angle was 30°, the overall factor of safety (FS) and the minimum factor of safety (MFS) of the rock slope were the smallest. By assuming that the layered rock slope was homogeneous, the HD and MHD would be underestimated and FS and MFS would be overestimated. The obtained results are likely to provide a theoretical basis for the prediction and monitoring of layered rock landslides.

1. Introduction

Rainfall is one of the main factors that trigger landslides [1]. In the past two decades, research on rainfall-induced soil landslides has been gaining popularity [2,3,4,5,6], but research on rock landslides is less reported [1,7,8,9,10,11,12]. According to available statistical data, more than 90% of rock slope failures are related to fluid hydraulic conduction [13]. The strength parameters of rock slopes such as Poisson’s ratio and elastic modulus also influence the stability of rock slopes [14]. Collins and Znidarcic [14] pointed out that the depths of slope failures were not only controlled by hydraulic characteristics, but also by strength characteristics; thus, both factors should be considered to assess the stability of slopes.

Anisotropy is everywhere but isotropy is rare [15]. Barton and Quadros [15] questioned whether the a priori assumption of homogeneous-isotropic-elastic behavior had any significant place in the scientific practice of realistic rock mechanics. Hence, the physical and mechanical behaviors of anisotropic rocks, such as sandstone [16,17] shale [18,19,20], slate [21,22], and phyllite [23,24], have become a hot topic over the last decades [25]. Seriously weathered carbonaceous slate, as a kind of layered rock, has an obvious bedding structure [26]. The layered structure of slate rock induces strength anisotropy. The rock strength in the direction parallel to the layers is considerably smaller than in other directions [27]. Based on the triaxial compression test of shale, Geng et al. [28] found out that the brittleness anisotropy of Longmaxi shale samples was dependent on the bedding angle. Ding et al. [25] conducted laboratory tests of anisotropic slate with various dip angles and found out that the uniaxial compression strength exhibited a typical U-type trend with the variance of dip angle, the elastic modulus first decreased and then increased, and the variation of Poisson’s ratio showed an opposite trend with the variation in elastic modulus. Based on uniaxial and triaxial compression tests, Gholami and Rasouli [22] found out that the strength of slate exhibited obvious U-shaped anisotropy with an increasing dip angle. Chen et al. [21] found out that the apparent elastic modulus increased remarkably with the dip angle, and the compressive strength at any confining pressure varied in a typical U-shaped trend. The above research mainly focuses on the anisotropy of layered rock mass, but the strength anisotropy of actual rock slope is less studied.

The deformation and hydraulic conduction of layered rock masses display anisotropic characteristics [29]. The ratio of the hydraulic conductivity in the parallel direction to the vertical direction of the bedding plane in the western foothills of Taiwan was 10 to 100 [1]. Dong et al. [11] found out that the anisotropy ratio and dip angle of the bedding plane had significant effects on the pore pressure distributions and slope stability; however, the above research is restricted to dip angles of 30, 0, and -30°, and the influence of anisotropy ratio was ignored. Yeh et al. [1] also found the heterogeneity in hydraulic conduction and strength for layered rock masses. It was considered that the effects of anisotropy of hydraulic conductivity and strength on groundwater flow and slope stability were mainly reflected in the direction of the bedding plane. However, the rock slope was seriously weathered, while the Poisson’s ratio and elastic modulus used for numerical simulation were from intact rocks, which would lead to inaccurate results. Yu et al. [3] and Xia et al. [9] studied the seepage, deformation, and stability of soil slopes based on the anisotropy of hydraulic conductivity; however, the influence of the anisotropy of elastic modulus and Poisson’s ratio on the slope deformation and stability was ignored.

In view of the shortcomings of previous studies, the theory of fluid–solid coupling and the unsaturated seepage are introduced. Based on rock uniaxial compression tests, the strength conversion equations of weathered rock are proposed. A seriously weathered layered rock slope in the Pulang region of Southwestern China was selected as the study objective. By considering different anisotropy ratios and anisotropy angles of hydraulic conductivity and strength, numerical simulations were carried out to analyze the stability and deformation characteristics of rock slope. Major parameters, including the horizontal displacement (HD), maximum horizontal displacement (MHD), the factors of safety (FS), the minimum factors of safety (MFS), and factors of safety of landslide (FSL) were proposed to characterize the slope deformation and stability pattern. Finally, by comparing the actual detection data of the rock landslide with the numerical simulation results, the feasibility of the strength conversion formula and the accuracy of the landslide prediction were verified. The research results could provide some references for monitoring and predicting rock landslides.

2. Materials and Methods

2.1. Control Differential Equation

To simulate the seepage processes of rock slopes, rainfall infiltration was regarded as the variable flow boundary of the unsaturated zone, and a saturated–unsaturated seepage theory was applied to perform the simulations. Under the principle of mass conservation, the saturated–unsaturated seepage control differential equation is as follows [3,9,30]:

$$\frac{\partial }{\partial x}(k_x\frac{\partial H}{\partial x})+\frac{\partial }{\partial y}(k_x\frac{\partial H}{\partial y})+Q=m_w{\gamma }_w\frac{\partial H}{\partial t}$$

(1)

where H is the total head, the unit is m; x and y are the coordinates in the direction of x and y; k_x is the hydraulic conductivity in the x-direction, the unit is m/s; k_y is the hydraulic conductivity in the y-direction, the unit is m/s; Q is the applied boundary flux, m_w is the slope of the storage curve, t is the time, the unit is s; γ_w is the unit weight of water, the unit is N/m³.

Applying the Galerkin method of weighted residual to the governing differential equation, the finite element seepage equation in two-dimension can be derived as [3]:

$$\tau {\displaystyle \underset{A}{\int }({\left[B\right]}^T\left[C\right]\left[B\right])dA\{H\}+}{{\displaystyle \underset{A}{\int }(\lambda \langle N\rangle }}^T\langle N\rangle )dA\{H\},t=q\tau {\displaystyle \underset{L}{\int }({\langle N\rangle }^T)}dL$$

(2)

where [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 m_wγ_w; 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]\left\{H\right\}+[M]\left\{H\right\},t=\left\{Q\right\}$$

(3)

where [K] is the element characteristic matrix; [M] is the element mass matrix; and {Q} is the element applied flux vector.

2.2. Theory of Fluid–Solid Coupling

For a two-dimensional space, the incremental stress–strain relationship can be written as [3]:

$$\text{{}\Delta \sigma \text{}}=[D]\text{{}\Delta \epsilon \}-[D]\{m_H\left\}\right\{u_a-u_w\}+\text{{}\Delta u_a\}$$

(4)

where ε is the normal strain; γ is the engineering shear strain; σ is the normal stress, the unit is kN; τ is the shear stress, the unit is kN; u_a is the pore-air pressure, the unit is kPa; u_w is the pore water pressure, the unit is kPa; [D] is the drained constitutive matrix; and {m_H} can be expressed as:

$${\left\{m_H\right\}}^T=\langle \begin{array}{cccc}\frac{1}{H^{\prime }}& \frac{1}{H^{\prime }}& \frac{1}{H^{\prime }}& 0\end{array}\rangle$$

(5)

where H’ is the unsaturated rock modulus for rock structure with respect to matric suction (u_a − u_w).

2.3. Establishment of the Factor of Safety for Unsaturated Rock

This study adopted a Morgenstern–Price (M–P) method based on the limit equilibrium theory to calculate the factors of safety. The M–P method [31] for the analysis of rock slope stability offered various user-specified inter-slice force functions and considered both the shear and normal inter-slice forces, which made the analysis more consistent with the actual engineering situation [32,33,34]. It was found that the method strictly satisfied the force balance and torque balance, and the calculation accuracy was found to be high.

Taking the rock slice as the research object, as shown in Figure 1, when the boundary conditions E₀ = 0 and E_n = 0, according to the balance of the force in the normal and tangential directions of the slip surface of the rock slice, it can be obtained:

$$S_i=(N_i\tan \phi +cb\sec \alpha )/Fs$$

(6)

$$N_i=(W_i+\lambda f_{i-1}{E}_{i-1}-\lambda f_i{E}_i+Q_i\cos {\omega }_i+J_i\sin {\beta }_i)\cos \alpha +E_i-E_{i-1}+Q_i\sin {\omega }_i-J_i\cos {\beta }_i)\sin {\alpha }_i$$

(7)

where N_i is the effective normal force on the slip surface of the rock slice and the unit is kN; S_i is the shear force on the slip surface of the rock slice and the unit is kN, and the safety factor Fs can be derived as:

$$Fs=\frac{{\displaystyle \sum _{i=1}^{n-1}(R_i^{}{\displaystyle \prod _{j=i}^{n-1}{\psi }_j)+R_n}}}{{\displaystyle \sum _{i=1}^{n-1}(T_i^{}{\displaystyle \prod _{j=i}^{n-1}{\psi }_j)+T_n}}}$$

(8)

where R_i is the resistance force, $R_i=[W_i\cos {\alpha }_i+Q_i\cos ({\omega }_i-{\alpha }_i)+J_i\sin ({\beta }_i-{\alpha }_i)]\tan {\phi }_i+c_i{b}_i\sec {\alpha }_i$ and the unit is kPa; T_i is the sliding force, $T_i=W_i\sin \alpha -Q_i\sin ({\omega }_i-{\alpha }_i)+J_i\cos ({\beta }_i-{\alpha }_i)$; ψ_i is the transfer coefficient, ${\psi }_{i-1}=\frac{(\sin {\alpha }_i-\lambda f_{i-1}\cos {\alpha }_i)\tan {\phi }_i+(\cos {\alpha }_i+\lambda f_{i-1}\sin {\alpha }_i)Fs}{(\sin {\alpha }_{i-1}-\lambda f_{i-1}\cos {\alpha }_{i-1})\tan {\phi }_{i-1}+(\cos {\alpha }_{i-1}+\lambda f_{i-1}\sin {\alpha }_{i-1})Fs}$; c_i is the effective cohesion for every rock slice and the unit is kPa; φ_i represents the shear strength angle for every rock slice, the unit is the degree. J_i is the seepage force of the rock slip which adopts a simplified seepage field treatment method, assuming that its position of action is located at the center of gravity of the rock below the underwater line, and the distance from the center of the slip surface of the rock strip is h_wi/2, and h_wi is the distance from the underwater line to the slip surface. E_i and E_i−1 are the horizontal effective forces between the two sides of the rock slice and the unit is kN; λf_iE_i and λf_i−1E_i−1 are the shear forces between the two sides of the rock strip, f_i is the function of the inter-strip force, and λ is the proportional coefficient. The distances between the action positions on both sides and the center of the slip surface of the rock strip are Z_i and Z_i−1, respectively; W_i is divided into two parts based on the underwater line: W_i = W_i1 + W_i2, where W_i is the gravity of the rock above the underwater line and W_i2 is floating weight of the rock below the underwater line. Q_i is an external force on the rock surface; ω_i is the angle between the direction of the external force and the normal direction; α_i is the angle of the slip surface.

2.4. Establishment of the Numerical Calculation Model

The layered rock slope in the Pulang region of Southwestern China is shown in Figure 2. In Figure 3, the stratum lithology of the slope is seriously weathered carbonaceous slate with a layered structure. This finite element model of the slope body is 260 m high and 540 m long. The left and right borders extend horizontally for 50 m, and the lower border extends downward for 50 m. The grid model is divided into 18,233 nodes and 18,196 elements. To study the characteristics of the deformation at different positions, three monitoring lines were set at x = 150 m (top of the slope), x = 250 m (middle of the slope), and x = 330 m (bottom of the slope). This rock slope experienced a landslide on 4 July 2019 and the scope of the landslide is shown in Figure 2a. Displacement monitoring stations were set between the bottom and the middle of the slope: They were located 1 m below the slope and were set at x = 345 m (BPK2-1), x = 370 m (BPK2-2), and x = 395 m (BPK2-3).

In Figure 3, the boundary conditions are as follows: AB and FD are the fixed water level boundaries of 175 and 90 m, respectively; CD indicates the boundary of the rainfall infiltration; BC is the free boundary with a limited flow; AF is an impervious boundary. Since the rock landslide occurred on 4 July 2019, to restore the rainfall boundary, the average daily rainfall from 1 March to 2 August 2019 (155 days) in the Pulang region was adopted as the simulated rainfall conditions, as detailed in Figure 4.

2.5. Determination of the Maximum Initial Matric Suction

Since the studied rock slope had a strong degree of weathered and the rock layers were so intensely fractured that the properties of a porous medium could be assumed, it had a certain degree of matric suction. The extent of the initial maximum matric suction may have varied, and the seepage states and stability levels of slopes may have also been quite different. Therefore, to make the initial matric suction more in line with the actual situation (refer to previous studies [9,35]), this study set the different maximum initial matric suctions as −50, −75, and −100 kPa, respectively, and the annual average rainfall of the Pulang area as 0.06 mm/h in the steady-state comparative analysis process. As detailed in Figure 5, the changes in the matric suction of the average annual rainfall were similar to the setting of the maximum initial matric suction of −75 kPa. Therefore, the initial maximum matric suction of all this study’s simulation experiments was set as −75 kPa, which was determined to be in line with the real situation in the study area.

2.6. Establishment of Unsaturated Permeability Coefficient

Due to the difficulties of measuring the unsaturated characteristic curves of rock masses, the rock mass was assumed as porous media or unsaturated soils, and the water retention curves of classic unsaturated porous media [36,37,38,39], and relative hydraulic conduction models [40] were applied to describe the unsaturated flow in fissures and weakly permeable rocks. Therefore, this study adopted the van Genuchten model to describe unsaturated features of rock masses. The van Genuchten model is as follows [37]:

$$\theta ={\theta }_r+\frac{{\theta }_s-{\theta }_r}{{[1+{(p/a)}^n]}^m},$$

(9)

$$k=k_w\frac{{[1-{\left|p/a\right|}^{n-1}{(1+{\left|p/a\right|}^n)}^{-m}]}^2}{{(1+{\left|p/a\right|}^n)}^{\frac{m}{2}}}.$$

(10)

where p is the suction and the unit is kPa; θ indicates the adjusted volumetric water content, the unit is m⁻¹; θ_r is the residual volumetric water content, the unit is m⁻¹; θ_s indicates the saturated volumetric water content, the unit is m⁻¹; a is the fitting parameter closely related to the air-entry value of the unsaturated rock mass and the unit is kPa; n and m (m = 1 − 1/n and n > 1) are fitting parameters that control the slope at the inflection point in the volumetric water content function [18]. k_w is the saturated hydraulic conductivity, unit is m/s; k is the adjusted hydraulic conductivity, the unit is m/s.

Parameters such as a, n, m, and θ_r are known to have major influences on the unsaturated flow of rock slope [9]. The slope in the present study is located on the west side of the slope in a previous study [9], and the difference of unsaturated parameters between the two rock slopes is negligible. Therefore, the parameters of the water retention curve of the present slope were referred from a previous study [9]. The unsaturated parameter values are shown in

Table 1. Parameters of the water retention curve.

Layer	Materials	Fitting Parameters					Hydraulic Conduction Coefficient
		A (kPa)	m	n	θs	θr	k (m/s)
I	Strongly weathered carbonaceous slate	10	0.33	1.5	0.242	0.001	8.08 × 10−5

[9]. Accordingly, the water retention curve of the rock mass is shown in Figure 6.

2.7. Definition of Anisotropy and the Calculation Conditions

The different dip angles of the bedding plane of the rock slopes may result in different hydraulic conductivity [9]. In this study, the hydraulic conductivity matrix [C] is expressed as follows:

$$\left[C\right]=\left[\begin{array}{cc}{C}_{11}& C_{12}\\ C_{21}& C_{22}\end{array}\right]$$

(11)

In Equation (11), C₁₁ = k_xcos²β + k_ysin²β, C₂₂ = k_xsin²β + k_ycos²β, and C₁₂ = C₂₁ = k_xsinβcosβ + k_ysinβcosβ. The k_x and k_y, as well as the dip angle β, were defined according to Figure 3. In the present study, k_x represents the horizontal hydraulic conductivity; k_y is the vertical hydraulic conductivity; β is the direction between k_x and x-axis. Therefore, if β = 0°, then [C] will be reduced to:

$$\left[C\right]=\left[\begin{array}{cc}{k}_x& 0\\ 0& k_y\end{array}\right]$$

(12)

Equation (12) was adopted in this study [4] with only the anisotropy ratio k_r = k_y/k_x considered. However, for layered rock masses, due to the existence of dip angles of the bedding plane, the conditions that the anisotropy angle β was not equal to 0 were also considered. Therefore, in this study, to restore a reasonable rock formation state, anisotropy angles ranging between 0 and 90° in the layered rock slope in the Pulang area were taken into consideration. As shown in Figure 3, the anisotropic angle of hydraulic conductivity was equivalent to the dip angle of the layered rock slope.

In addition, the dip angle of the rock slope also affected the strength parameters of the rock mass, such as elastic modulus and Poisson’s ratio. The rock was processed into rock specimens with dip angles of 0, 15, 30, 45, 60, 75, and 90°, as shown in Figure 7. The WAW-1000kN microcomputer-controlled electro-hydraulic Servo universal testing machine was used to adopt the uniaxial compression test.

In Figure 8, through the uniaxial compression experiment, the elastic modulus of the rock specimens with different dip angles is in the shape of “U”, reaching the maximum when the dip angle is 90°, and reaching the minimum when the dip angle is 30° [21]. The influence of dip angle variance on Poisson’s ratio was minimal [21].

The strength parameters of the rock mass were not only affected by the dip angle, but also by the integrity of the rock mass [42]. Since the stratum lithology of the rock slope in the present study belonged to seriously weathered carbonaceous slate, the strength parameters of the intact rock could not be directly applied to the numerical simulation calculation. According to damage mechanics [42,43], the rock masses with different weathered states had different elastic modulus. There was an assumption that with the change of the dip angle, the change trends of the elastic modulus of different weathered states were similar. Therefore, the conversion equation of elastic modulus between intact rocks and seriously weathered rock masses was expressed as follows:

$$\frac{E^d-E_{\min }^d}{E_{\max }^d-E_{\min }^d}=\frac{E-E_{\min }}{E_{\max }-E_{\min }}$$

(13)

From Equation (13), the following can be derived:

$$E^d=\frac{E-E_{\min }}{E_{\max }-E_{\min }}(E_{\max }^d-E_{\min }^d)+E_{\min }^d$$

(14)

The same assumption for Poisson’s ratio is expressed as follows:

$${\upsilon }^d=\frac{\upsilon -{\upsilon }_{\min }}{{\upsilon }_{\max }-{\upsilon }_{\min }}({\upsilon }_{\max }^d-{\upsilon }_{\min }^d)+{\upsilon }_{\min }^d$$

(15)

E and v are the elastic modulus and Poisson’s ratio of the intact rock; E^d and v^d are the elastic modulus and Poisson’s ratio of the seriously weathered rock mass; E_max, v_max, E_min, and v_min are the maximum and minimum elastic modulus and Poisson’s ratio of the intact rock; E^d_max, v^d_max, E^d_min, and v^d_min are the maximum and minimum elastic modulus and Poisson’s ratio of the seriously weathered rock mass, respectively. E^d_max, v^d_max, E^d_min, and v^d_min were measured by the uniaxial compression experiment. According to Equations (13)–(15), the strength parameters of the seriously weathered rock mass for different dip angles in the present study are shown in Figure 8.

Based on previous research [1,3,9,12,44], the anisotropy ratios of hydraulic conductivity were set as 0.002, 0.01, 0.02, 0.1, and 1, and the dip angles (anisotropy angle of hydraulic conductivity) were set as 0, 15, 30, 45, 60, 75, and 90°. The failure criterion of the rock layered slope simulation adopted the Mohr–Coulomb criterion. Correspondingly, the elastic modulus and Poisson’s ratio of the rock slope were also set according to Figure 8. The average values of the cohesion and friction angle for the whole rock slope were measured by the direct shear experiment on site, as detailed in

Table 2. Different calculation conditions.

Rock Types	Anisotropy Ratio kr = ky/kx	Anisotropic Angle/ Dip Angle β (°)	Elastic Modulus Ed (GPa)	Poisson Ratio vd	Unit Weight (kN/m3)	Cohesion (kPa)	Friction Angle (°)
seriously weathered carbonaceous slate	$\left[\begin{array}{c}0.002\\ 0.01\\ 0.02\\ 0.1\\ 1\end{array}\right]$	$\left[\begin{array}{c}0\\ 15\\ 30\\ 45\\ 60\\ 75\\ 90\end{array}\right]$	$\left[\begin{array}{c}2.902\\ 3.014\\ 2.645\\ 2.728\\ 2.953\\ 4.428\\ 5.561\end{array}\right]$	$\left[\begin{array}{c}0.3569\\ 0.3546\\ 0.3500\\ 0.3731\\ 0.3638\\ 0.3800\\ 0.3615\end{array}\right]$	22.4	93.6	33.3

.

3. Results

3.1. Effects of the Hydraulic Conductivity and Strength Anisotropy on the Deformation Characteristics

3.1.1. Analysis of the Horizontal Displacement of Monitoring Line

Based on the fluid–solid coupling theory, the layered rock slope produced the corresponding displacement under the action of the fluid. The seepage field exerted normal seepage pressure and tangential tension on the rock mass skeleton and affected the rock mass stress distribution. Under such stress conditions, the rock mass deformed accordingly. At the same time, the stress field affected the pores of the rock mass, and the hydraulic conductivity changed with the change of the pore size; thus, the stress field also affected the seepage field of the rock mass. According to the calculation conditions in Table 2, a total number of 70 numerical simulations was carried out, and 105 sections of horizontal displacement (HD) variation were obtained. The variations of the HD of layered rock slope with different β values and k_r = 0.01 and k_r = 0.1 are shown in Figure 9 to illustrate the influence of anisotropy direction β on deformation characteristics. Different k_r values with β = 0, β = 30, and β = 90° are also displayed in Figure 10 to show the impacts of anisotropy ratio k_r.

In Figure 9, the initial state is the distribution of the horizontal displacement (HD) of the rock slope before rainfall, and the remainder shows the distribution of the HD of day 155. The HDs of day 155 changed differently at different positions. At the top and middle of the slope, the HDs first increased and then decreased with the elevation. The maximum HDs at the top and middle of the slope were approximately 50 and 15 m from the surface, respectively. At the slope bottom, the HDs increased with elevation. The HDs at the bottom and middle of the slope were larger than the top of the slope. According to previous studies [3,9,30], the hydrodynamic pressure generated by rainfall had a large impact on the bottom and middle parts of the slope, and the landslide between the bottom and middle of the slope, which verified the results of the numerical simulation. Generally speaking, the hydraulic conductivity and strength anisotropy characteristics also had great impacts on the variation of HD. When the dip angle β reached 30°, the HD reached the maximum. When the HDs at different positions of the rock slope reached the minimum, the corresponding dip angles of the rock formation were different.

The trend of HD distribution in Figure 10 was similar in Figure 9. It is worth noting that when the dip angle β reached 30°, as the anisotropy ratio k_r increased, the HD decreased. When the k_r was smaller than 0.01, the maximum HD of the slope top appeared at 10 m from the surface. When the anisotropy angle was greater than 30° and k_r = 1, the HD was the smallest. Therefore, considering the layered rock slope as a homogeneous medium may underestimate its horizontal displacement.

In Figure 11 and Figure 12, the initial state is the distribution of the pore water pressure (PWP) of the rock slope before rainfall and the remainder shows the distribution of the PWP of day 155. With the elevation, the initial PWP at the slope top and slope middle first decreased and then remained stable, and the initial PWP at the slope bottom kept decreasing. When the rock slope was regarded as an isotropic medium, the PWP distribution trend of day 155 was similar to the initial PWP distribution trend. The corresponding HDs meanwhile reached the minimum. When the rock slope was regarded as a heterogeneous medium, the PWP distribution trend of day 155 was different from the initial PWP distribution trend, and this difference was mainly manifested in the shallow layer. At the top and middle of the slope, as the k_r decreased, the shallow PWP gradually increased, and at the slope bottom, as the k_r decreased, the shallow PWP gradually decreased. The distribution of shallow HD was also similar to shallow PWP. At the top, middle, and bottom of the slope, as the β increased, the shallow PWP gradually increased. However, the shallow HD first increased, then decreased, and reached the maximum dip angle of 30°. Since the HD was not only affected by water pressure but also by its strength parameters (elastic modulus and Poisson’s ratio), the variations of HD were not synchronous with PWP.

3.1.2. Analysis of the Horizontal Displacement of Monitoring Point

The variations of horizontal displacement for the monitoring point of the layered rock slope during rainfall are shown in Figure 13. BPK2-1, BPK 2-2, and BPK 2-3 were the field monitoring data obtained by the displacement monitoring station. As a landslide occurred just near the field station during the monitoring activity, the deformation behavior of the rock slope could be detected [45]. In correspondence to the rainfall event of 4 July 2019 (day 91), the monitoring point of the rock slope experienced a sudden increase in horizontal displacement. Through numerical simulations, it was found that the HD of the monitoring point was underestimated when the rock slope was considered isotropic, compared to anisotropy during the rainfall period. Under varied k_r and β, the HDs calculated by numerical simulations all increased abruptly on 4 July 2019, indicating that the rock slope model could accurately predict the occurrence time of the landslide. It is worth noting that when k_r = 0.01 and β = 30°, the calculated HD was consistent with the field monitoring data.

3.1.3. Analysis of the Maximum Horizontal Displacement

In order to study the deformation at different slope positions under varied k_r and β conditions more prominently, the maximum horizontal displacement (MHD) of each monitoring line was defined. Figure 14a–c shows that MHD at the different positions of the slope is greatly affected by k_r and β. When β = 30°, the MHD reached the maximum, and the elastic modulus and Poisson’s ratio of the rock slope were the smallest in this case. When the anisotropy angle was close to 60°, and the MHD reached its minimum. When the rock slope was considered to be isotropic (k_r = 1), the simulated MHD was almost negligible. Previous studies only took k_r or β into consideration, but ignored their collective effect.

Table 3. Differences between only considering k_r or β and considering both k_r and β.

MHD	Only Considering kr	Only Considering β	Considering Both kr and β
Slope top	23%	118%	2.25 × 105%
Slope middle	86%	127%	3.38 × 105%
Slope bottom	81%	123%	1.17 × 105%

shows the MHD of the differences between only considering k_r or β and considering both k_r and β compared with the isotropy condition.

As can be seen from Table 3, the values of MHD for the rock slope are quite different under the conditions when only considering k_r or β and when considering both k_r and β. Therefore, the k_r and β must be considered in the deformation analysis of rock slope, and MHD is more affected by β than k_r.

3.2. Effects of the Hydraulic Conductivity and Strength Anisotropy on the Stability of the Rock Slope

3.2.1. Analysis of the Factor of Safety

The variations of the factor of safety (FS) for the rock slope under varied anisotropy ratios and dip angles are shown in Figure 15. When β = 30°, the overall FS of the rock slope was the smallest. During the period of heavy rainfall (from 21 June to 26 July 2019), the overall FSs obviously decreased and the FS of the rock slope decreased with the increase of β. Xia et al. [9] conducted relevant simulations and the results were similar to current research results. This is because a larger dip angle of the bedding plane would lead to a smaller hydraulic conductivity in the horizontal direction and a larger hydraulic conductivity in the vertical direction, which would promote the rainwater infiltrating into the deep part of the slope body [3,9,30], increase the height of the groundwater level and pore water pressure, and finally reduce the effective stress and shear strength. When the rainfall was heavy and the rainfall time was short, the FS increased sharply. This is because the shallow layer of the rock slope belongs to the unsaturated zone, and the hydraulic conductivity is small. A large amount of rainwater rushing into the rock slope for a short time caused blockage, which led to the accumulation of rainwater on the surface of the rock slope. The growth rate of pore water pressure inside the slope was less than the growth rate of the rainfall water pressure on the slope surface, which formed a reverse pressure on the rock slope. This would make it more difficult for the slope to transfer from an unsaturated state to a saturated state. The negative pressure zone gradually increased, leading to a corresponding increase in the FS of the rock slope.

It can be seen from Figure 16 that during the period of heavy rainfall (from 21 June to 26 July 2019), when β = 0°, the overall FSs of the rock slope decreased with the increase of k_r. In the horizontal rock bedding plane, the hydraulic conductivity in the horizontal direction reached the maximum, but as k_r increased, the hydraulic conductivity of the vertical direction increased, and the wet front of rain could reach the deep part of the slope. During the period of heavy rainfall, when β = 90°, the smaller the anisotropy ratio k_r, the more notable the decrease of the FS was. In the vertical rock bedding plane, the hydraulic conductivity of the vertical direction reached the maximum, and as k_r decreased, the hydraulic conductivity of the horizontal direction decreased, and the rainwater infiltration capacity was much greater than the horizontal drainage capacity. When the anisotropy angle β was close to 30°, with the change of the anisotropy ratio k_r, the change law of FS was not notable. However, when rock slope was treated as an isotropic medium, the FS was overestimated.

3.2.2. Analysis of the Minimum Factor of Safety

In order to study the stability of rock slopes under varied k_r and β conditions more prominently, the minimum factor of safety (MFS) and the factor of safety of landslide (FSL) on 4 July 2019 were defined. Variations in the factor of safety under varied k_r and β are shown in Figure 17. When the k_r is between 0.002 and 0.02, the MFS and FSL reached the minimum and maximum at β = 30 and β = 0°, respectively. It is worth noting that when the β was greater than 75°, the MFS reached the minimum, and the reason is discussed in Section 3.2.1.

Table 4. Differences between only considering k_r or β and considering both k_r and β.

MHD	Only Considering kr	Only Considering β	Considering Both kr and β
FSL	7%	16%	16.4%
MFS	6%	17%	17.6%

shows the MFS and FSL of the differences between only considering k_r or β and considering both k_r and β compared with the isotropy condition.

As can be seen from Table 4, the values of MFS and FSL for rock slope are similar under the conditions when only considering β to considering both k_r and β. Therefore, we can consider only the effect of β but ignore k_r to simplify the calculation.

4. Discussion

The strength parameters of seriously weathered rock slopes like the elastic modulus and Poisson’s ratio were obtained by the conversion Equations (13)–(15). Through numerical simulations, the anisotropic rock slope model could predict the time of landslide occurrence. Figure 18 shows the cloud map of the horizontal displacement when the landslide occurred on 4 July 2019 (day 91). It could be found that when k_r = 0.01 and β = 30°, the simulated landslide scope almost coincided with the field landslide scope. When the rock slope was regarded as a homogeneous medium (k_r = 1), the distance between the simulated landslide scope and the field landslide scope was about 45 m, and the horizontal displacement was significantly small, which made it difficult to predict the occurrence of the landslide. According to the field survey, the field dip angle of the rock slope was 32°, as shown in Figure 19, which was similar to the dip angle of 30° in the simulation. Thus, in general, the landslide scope, horizontal displacement, dip angle, and landslide occurrence time of the field investigation were consistent with the results of the simulation. Therefore, the feasibility of the conversion Equations (11)–(14) of elastic modulus and Poisson’s ratio was also verified.

5. Conclusions

A layered rock slope in the Pulang Region, southwestern China, is selected as a case study. A numerical model of the landslide is constructed by considering different anisotropy ratios and anisotropy angles of hydraulic conductivity and strength. The research results were as follows:

(1)

The strength conversion equations of elastic modulus and Poisson’s ratio were feasible, and the rock slope model could accurately predict the occurrence time, horizontal displacement, and scope of the landslide;

(2)

The different anisotropy ratios and dip angles of the bedding plane were found to have major impacts on the deformation and stability of the layered rock slope;

(3)

The horizontal displacement (HD) and maximum horizontal displacement (MHD) were determined to characterize the deformation characteristics of the rock slope. Considering the layered rock slope as a homogeneous medium could underestimate its HD and MHD. When the dip angle was 30°, the MHD reached the maximum. When the anisotropy angle was close to 60°, and the MHD reached its minimum;

(4)

The factor of safety (FS), the minimum factor of safety (MFS), and the factor of safety of landslide (FSL) were determined to characterize the stability characteristics of the rock slope. When the dip angle was 30°, the FS, MFS, and FSL of the rock slope reached the minimum. However, when the rock slope was treated as an isotropic medium, the FS, MFS, and FSL were overestimated;

(5)

The changing law of the anisotropy ratio was not obvious, and it was difficult to verify by field data. This could be a focus of future research work. The obtained results are likely to provide a theoretical basis for the prediction and monitoring of layered rock landslide.