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

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.


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, Water 2020, 12 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.

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]: 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 3 . 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]: (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: where [K] is the element characteristic matrix; [M] is the element mass matrix; and {Q} is the element applied flux vector.

Theory of Fluid-Solid Coupling
For a two-dimensional space, the incremental stress-strain relationship can be written as [3]: 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: where H' is the unsaturated rock modulus for rock structure with respect to matric suction (u a − u w ).

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 = 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: 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 ; 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 i E i and λf i−1 E 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. Taking the rock slice as the research object, as shown in Figure 1, when the boundary conditions E0 = 0 and En = 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: where Ni is the effective normal force on the slip surface of the rock slice and the unit is kN; Si 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: where Ri is the resistance force, and the unit is kPa; Ti is the sliding force, ; ψi is the transfer coefficient, ; ci 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. Ji 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 hwi/2, and hwi is the distance from the underwater line to the slip surface. Ei and Ei-1 are the horizontal effective forces between the two sides of the rock slice and the unit is kN; λfiEi and λfi-1Ei-1 are the shear forces between the two sides of the rock strip, fi 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 Zi and Zi-1, respectively; Wi is divided into two parts based on the underwater line: where Wi is the gravity of the rock above the underwater line and Wi2 is floating weight of the rock below the underwater line. Qi 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.

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

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.  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. 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.

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.

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 Water 2020, 12, 3056 7 of 24 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.

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]: where P is the suction and the unit is kPa; θ indicates the adjusted volumetric water content, the unit is m −1 ; θr is the residual volumetric water content, the unit is m −1 ; θs indicates the saturated volumetric water content, the unit is m -1 ; 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]. kw is the saturated hydraulic conductivity, unit is m/s; k is the adjusted hydraulic conductivity, the unit is m/s.

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]: where p is the suction and the unit is kPa; θ indicates the adjusted volumetric water content, the unit is m −1 ; θ r is the residual volumetric water content, the unit is m −1 ; θ s indicates the saturated volumetric water content, the unit is m −1 ; 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 [9]. Accordingly, the water retention curve of the rock mass is shown in Figure 6. 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 [9]. Accordingly, the water retention curve of the rock mass is shown in Figure 6.

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: [ ] 11 12 In Equation (11), C11 = kxcos 2 β + kysin 2 β, C22 = kxsin 2 β + kycos 2 β, and C12 = C21 = kxsinβcosβ + kysinβcosβ. The kx and ky, as well as the dip angle β, were defined according to Figure 3. In the present study, kx represents the horizontal hydraulic conductivity; ky is the vertical hydraulic conductivity; β is the direction between kx and x-axis. Therefore, if β = 0°, then [C] will be reduced to: Equation (12) was adopted in this study [4] with only the anisotropy ratio kr = ky/kx 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

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: In Equation (11), C 11 = k x cos 2 β + k y sin 2 β, C 22 = k x sin 2 β + k y cos 2 β, and C 12 = C 21 = k x sinβcosβ + k y sinβ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: 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.
Water 2020, 12, 3056 9 of 24 angles of 0, 15, 30, 45, 60, 75, and 90°, as shown in Figure 7. The WAW-1000kN microcomputercontrolled 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].  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].
angles of 0, 15, 30, 45, 60, 75, and 90°, as shown in Figure 7. The WAW-1000kN microcomputercontrolled 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: From Equation (13), the following can be derived: The same assumption for Poisson's ratio is expressed as follows:  (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 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: From Equation (13), the following can be derived: The same assumption for Poisson's ratio is expressed as follows:  (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.

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 .
Water 2020, 12, 3056 11 of 24 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.

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 kr = 0.01 and kr = 0.1 are shown in Figure 9 to illustrate the influence of anisotropy direction β on deformation characteristics. Different kr values with β = 0, β = 30, and β = 90° are also displayed in Figure 10 to show the impacts of anisotropy ratio kr.  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 kr increased, the HD decreased. When the kr 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 kr = 1, the HD was the smallest. Therefore, considering the layered rock slope as a homogeneous medium may underestimate its horizontal displacement.
In Figures 11 and 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 kr decreased, the shallow PWP gradually increased, and at the slope bottom, as the kr 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 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 Figures 11 and 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.
Water 2020, 12,3056 14 of 24 also by its strength parameters (elastic modulus and Poisson's ratio), the variations of HD were not synchronous with PWP.

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 kr 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 kr = 0.01 and β = 30°, the calculated HD was consistent with the field monitoring data.

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.

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 kr 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 kr = 0.01 and β = 30°, the calculated HD was consistent with the field monitoring data.

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 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 kr or β and when considering both kr and β. Therefore, the kr and β must be considered in the deformation analysis of rock slope, and MHD is more affected by β than kr.

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.  Table 3. Differences between only considering k r or β and considering both k r and β. 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.

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.

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 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.

Analysis of the Minimum Factor of Safety
In order to study the stability of rock slopes under varied kr 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 kr and β are shown in Figure 17. When the kr 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 shows the MFS and FSL of the differences between only considering kr or β and considering both kr 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 kr and β. Therefore, we can consider only the effect of β but ignore kr to simplify the calculation.

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 kr = 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 (kr = 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.  Table 4. Differences between only considering k r or β and considering both k r and β. 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.

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.

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;

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;

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.