Simulating Crack Development and Failure Characteristic of Toppling Rock Slope under Seismic Loading on Lancang River in China

Toppling rock slopes, induced by rapid and continuous downcutting of Lancang River, are widely distributed in the mountainous area of southwest China. To investigate the instability mechanism of 1# toppling rock slope of Huangdeng Hydropower Station under seismic loading, particle ﬂ ow code (PFC) is applied to simulate the dynamic response and failure mode. The study considers the particle characteristics of displacement, velocity, energy, and cracks. According to numerical results, the potential failure mechanism of toppling rock slope is identi ﬁ ed: multisliding surfaces form at the interfaces between the highly and moderately toppled rock mass and between the highly/moderately and weak toppled-crept rock mass; intersecting faults cut rock mass at the toe, leading to shear-toppling deformation; tension cracks develop, penetrate, and coalesce in the weak toppled-crept rock mass, resulting in tension-toppling-bending deformation. During the 2 to 5s of strong seismic intensity, crack increases sharply and energy of particles ﬂ uctuates greatly. The impacts of the amplitude of seismic loading and loading method in PFC are investigated. This study will provide a practically useful reference for seismic design of rock slopes.


Introduction
Toppling failure, as one of the most serious and hazardous landslide instability modes of rock slopes and underground caverns, usually occurs in mountainous and canyon areas, which is responsible for property loss and threats the safety of human life [1][2][3]. The characteristics of toppling failure are pervasive cracks and discontinuities developing parallel to the slope surface [4]. Landslides are commonly triggered by rainfalls, water level fluctuations, and earthquakes [5,6].
For earthquake-induced toppling failure and landslides, shaking-table tests and numerical simulations are effective approaches to study the characteristics of seismic wave propagation and dynamical response. Shaking-table tests can intuitively demonstrate toppling deformations and failure yet require more resources. Song et al. [7,8] studied the mechanism of instability of rock slopes with toppling discontinuous antidip joints under seismic loadings, with special attention focused on the relationship between local damage and landslides. However, the influence of joints on wave propagation and failure mode was not considered. Che et al. [9] and Song et al. [10] studied the influence of wave propagation on the stability of bedding rock slopes and other slopes containing nonpersistent joints through a combined approach using shaking-table test and finite element method (FEM). Discontinuity of wave at joint surfaces due to reflection and transmission was observed, but failure mode of slope was not studied. Other studies have investigated this aspect of wave propagation [11,12].
Discrete element method (DEM), including block discrete element (UDEC, 3DEC, and DDA) and particle discrete element (PFC2D, PFC3D), is more appropriate for simulating large-scale deformation of landslide than FEM and has been gradually applied to engineering projects. Currently, numerical simulations of toppling failure are usually carried out by block discrete element method, in which a set of discontinuous joints can be defined easily [13][14][15][16][17][18][19], but instability is induced mainly through weakening strength parameters of joints; it is difficult to implement such real condition as seismic loading. On the other hand, particle discrete element can be applied to better simulate landslides under real seismic loading condition [20,21]. Wang et al. [22] studied the characteristics of movement and deposits of Hongshiyan landslide under Ludian earthquake. Zhou et al. [23,24] investigated the formation process of the barrier lake formed by Yangjiagou landslide caused by Wenchuan earthquake. Tang et al. [25] and Chen et al. [26] simulated the dynamic deformation of Tsaoling landslide induced by Chi-Chi earthquake. To study the dynamic fragmentation process and the mechanism of transport of earthquake-induced landslides, Gao et al. [27] simulated the initiation, propagation, and coalescence of cracks inside rock mass and the process of dynamic disintegration of rock mass with the increase of sliding distance of the landslide in Tangjia Vally induced by Lushan earthquake in 2013. These studies focused on the characteristics of movement of sliding     bodies by implementing seismic loading in PFC2D/3D through changing the velocities of the bottom boundary. However, dynamic response and deformation inside rock slopes and characteristics of wave propagation were neglected. Furthermore, few studies implemented the equivalent nodal force of seismic loading and viscous boundary in particle discrete element. There is no publication yet on studying toppling deformation using the particle discrete element method.
In this paper, the #1 toppling rock slope of Huangdeng Hydropower Station is investigated. Based on the geological investigation, the toppling rock slope model is established using the particle flow method. The dynamic process of failure and the mechanism of instability under seismic loading are identified based on the development of cracks. The characteristics of dynamic response and wave propagation are investigated from two aspects: displacement response and acceleration response. The influences of amplitude of seismic loading and loading method in PFC2D are studied.
2. Characteristics of the #1 Toppling Rock Slope 2.1. Overview of the #1 Toppling Rock Slope. The #1 toppling rock mass slope is located on the right bank of Lancang River, about 700 m upstream of Huangdeng Hydropower Station in Yunnan Province, China [28,29]. The toppling rock slope is 1480-1830 m in elevation, 400-500 m in width, and 700 × 10 4 -800 × 10 4 m 3 in volume ( Figure 1). Detailed engineering geology conditions and formation process under effect of gravity and downcutting of the Lancang River can be found in other researches [30,31].
Currently, the toppling rock slope is basically in the stable state, with a potential instability mode of progressive failure in terms of displacement monitoring data. The toppling rock slope has complex deformation characteristics including shear sliding deformation in the fault zone at the toe and toppling deformation in the middle and rear parts of the rock slope.

Engineering Geological
Characteristics. The main exposed strata in the dam site area are Upper Triassic Xiaodingxi   and Quaternary avalanche deposit (Q col-dl ). As shown in Figure 2, the strata T3xd 4 -T3xd 5 mainly contain metamorphic basalt and metamorphic tuff, and T3xd 6 -T3xd 7 contain metamorphic volcanic breccia and metamorphic volcanic conglomerate. The tj4, tp1, and tp3 are thicker and more crushed metamorphic tuff interlayers, although tj4 develops along extrusion surface and more cleave develops well in tp1 and tp3. The t2 and t5 are also metamorphic tuff interlayers but have relatively good integrity.

Numerical Model and Seismic Loading
3.1. Numerical Model of the #1 Toppling Rock Slope. Necessary simplifications are made in the slope model based on field investigation and geology information analysis, where different zones are distinguished through the toppling extent of rock mass ( Figure 3). The main toppled rock mass can be classified into six zones according to six factors: dip angle, maximum width of the opening crack, unloading intensity, weathering intensity, and speed of the elastic wave [32]. The six zones include normal rock mass, highly toppled rock mass, moderately toppled rock mass, weakly toppled-crept rock mass, weakly toppled rock mass, and quaternary deluvium/quaternary alluvial deposits. Development of faults and joints leads to discontinuities and has a negative influence on the stability of slope.
A numerical model of A-A′ profile is established, as shown in Figure 4. The length and height are 654 m and 558 m, respectively, consistent with the actual sizes. In this study, the toppling rock slope is simulated by 34507 particles and the structures of faults and joints are modeled by fractures in PFC2D.

Contact Model and Calibration of Microparameters in PFC2D.
In PFC2D, rock mass can be modelled as an assembly of particles connected by the parallel bonded contact model (PBM). The parallel bonded model has been widely employed to study crack development and failure process of rock [33,34]. The nonpersistent joints are modeled by the smooth joint contact model (SJM), which can simulate the mechanical behavior of joints and are not affected by the contact direction of particles [35,36].
In general, the microparameters of particle and contact in PFC2D need to be calibrated to simulate the macrocharacteristics of a specific rock mass. However, there are no clear relationships between microparameters and macroparameters. The trial-and-error method has to be applied to link microparameters to the macroparameters obtained from uniaxial compressive tests. Based on the mechanical parameters of different stratum zones (shown in Table 1), including uniaxial compressive strength (UCS), Young's modulus (E), and Poisson's ratio (ν), a series of numerical tests are conducted and microparameters are obtained accordingly (shown in Table 2). Microparameters of PBM through intact rock are calibrated first, and the parameters of joints and faults are then obtained while keeping the PBM parameters unchanged. The calibration process of microparameters of the normal and jointed rock mass in bedrock zone are shown in Figures 5, satisfying UCS, E, and v obtained in the laboratory. Note that the dip of the joint is 81°, equal to that in the bedrock zone. The main microparameters influencing the mechanical characteristics of jointed rock mass are friction coefficient and dip angle in SJM.

Defining Monitoring Points.
In order to study the dynamic response, including displacement, velocity, and   5 Lithosphere acceleration, and the characteristics of wave propagation, 23 particle monitoring points are defined inside the slope. There are eight monitoring points located in alluvial/deluvium, eight points on the toppling rock surface, and other points are located in the horizontal and vertical directions inside the slope (Figure 6).

Boundary Conditions.
To simulate the infinite medium during process of dynamic analysis, the left and right boundaries should be set as viscous boundaries, which are capable of absorbing input seismic wave energy and eliminating wave reflection [37][38][39]. In PFC2D, the boundary force can be applied on the bottom boundary particles. The   6 Lithosphere relationship between boundary particle force and particle velocity is where R is the particle radius, ρ is the particle density, _ u is the particle velocity, and C is the wave speed of the rock.
According to the characteristics of wave propagation, C can be divided into two components (C p and C s ) through the following equations: Considering the transmission effect of the boundaries, the amplitude of the input wave should be doubled to prevent the amplitude from being halved after energy absorption. The equivalent force on the boundary is given by where _ UðtÞ is the velocity of the input seismic wave.

Seismic Input.
Seismic loading is applied on the bottom particles as particle contact forces through Equation (3), in which the velocity is determined by integrating the accelerations recorded at the Longtoushan Seismic Station during the Ludian earthquake (Figure 7), which can reflect regional geological conditions of southwest China. In this study, the amplitudes of the E-W and vertical components are adjusted to half of the original values (474.25 cm/s and 251.9 cm/s) to make seismic loading, and the impact of amplitude is in

Numerical Results and Analysis
4.1. Overall Instability of the Toppling Rock Slope. Figure 8 shows the displacement contour of the toppling rock slope, where the purple line represents tensile crack and lightyellow line represents shear crack. The alluvial/deluvium zones fail at 2 s and start to slide. Tension cracks develop preferentially at the toe and rear parts, while shear cracks concentrate mainly in the relatively denser joint area zone. As times go on, intersecting faults cut the rock mass, making the toe of slope unstable. At 3 s, a first-order sliding surface occurs in the highly toppled rock mass. At this time, only a

Characteristics of Dynamic
Response. The displacement of particles at different locations inside the slope is shown in Figure 9. The sliding mass in deluvium/alluvial travels farther than other parts since free slope surface has less resistance, as shown in Figure 9(a). The middle part of the deluvium is easy to slide for longer distances due to its own failure and being pushed by the upper sliding body above it. The largest travelling displacement is about 60 m. For the monitoring points inside the toppled rock mass along the surface (Figure 9(b)), points 3 and 4 of the highly toppled rock mass travel longer distances, implying that this part fails. Figures 9(c) and 9(d) show the characteristics of displacement in the horizontal and vertical directions, respectively. It is observed that the rock mass is relatively stable. Figure 10 shows the velocity of particles at different locations. Inside the toppling rock slope, the velocity timehistories are similar in both xand y-directions. As can be seen from Figures 10(c) and 10(d), these particles gradually return to a steady state at the end of seismic input.

Characteristics of Crack Development and Energy
Evolution. It is an effective tool to study the failure mechanism of a slope through the process of crack development. As shown in Figure 11, cracks increase sharply during 2-5 s, corresponding to the strongest motion of the seismic. Figure 8 at 10 s, the total number of cracks is 4214 with tension cracks (up to 3926) playing a dominant role, indicating that the failure mode of the slope under seismic is tension failure.

Referring to
The energy evolution of particles is monitored. During the same period of 2-5 s, energy fluctuates remarkably. In Figure 12, damping energy and slip friction energy consume most of the energy produced by seismic loading and gravity potential. There is very little change of kinetic energy and strain energy (elastic strain energy and parallel bond strain energy) in the steady state after 5 s.
A rose diagram is used to describe the distribution of the inclination angles of cracks (Figure 13(a)). It is seen that the inclination angles are basically distributed uniformly in 0°-180°. To better study the instability of the toppling rock slope, cracks in deluvium/alluvial are removed and the results are presented in Figure 13(b). It is observed that there are more cracks with inclination angles in 120°-180°, which are related to the slope angle and stratigraphic tendency.

Potential Failure Mode of the Toppled Rock Mass.
To better demonstrate the characteristics of local deformation, the region in the red rectangle shown in Figure 14 is magnified. Based on previous discussions and Figures 15 and 16, the failure mode of the middle highly toppled rock mass is multisliding surface traction landslide. At the toe, intersecting faults cut through the rock mass and toppling deformation occurs. Note that there are several large tension cracks parallel to the slope surface developing at the rear part of the slope. Under the influence of joints and tensile cracks, two rock blocks exhibit bend-toppling deformation as shown in the displacement contour plots. Figure 16 shows the evolution of velocity field and wave propagation in the toppled rock mass. At 3 s in the initial stage, velocity increases the fastest, leading to the damage of the highly toppled rock mass. It is obvious that wave propagation is discontinuous when passing through structural planes due to wave scattering. At the end of the seismic loading, the front part of the highly toppled rock mass slides, and the rear part starts to slide when the front part stops. This can be seen as the influence of traction.  Figure 14. When the PGA is small, deluvium/alluvial is unstable and slides, and only a few cracks develop. When the PGA is increased to 70% of the original value, large cracks penetrate joint surfaces at the bottom of the model. The front and middle parts are crushed, and toppling deformations, produced by tension cracks in the rear part, occur in several rock blocks.

The Impact of Input Mode of Seismic Loading.
For dynamic analysis in PFC2D, the dynamic input mode can be selected, but the final results may be different to some extent. Results of three loading modes are presented in Figure 17. In the previous sections, the loading is applied as forces on the bottom boundary particles, which are determined using equation (2). When the loading is applied as velocity obtained by integrating acceleration on the bottom boundary particles, crack development is not obvious at the interface between different strata. Meanwhile, the inclination dip of cracks is parallel to the slope surface. On the other hand, applying loading as velocity on the bottom     12 Lithosphere boundary rigid wall has been widely used in dynamic simulation; however, no significant failure deformation can be observed, and more shear cracks appear along joints and faults.

Conclusions
The potential failure mode of the toppling rock slope under seismic loading is studied using the particle flow method.
The following conclusions can be drawn from this study: (1) The toppling rock slope has unique and complex failure-deformation characteristics under seismic loading. The deluvium/alluvial is prone to instability and sliding. Multisliding surfaces traction landslide occurs in the toppled rock mass. A third-order sliding surface develops in the interface of the highly/moderately toppled rock mass zone and the weakly toppled-crept rock mass zone, followed by the bifurcation behavior of sliding surface. At the toe, intersecting faults cut through the rock mass, leading to shear-toppling deformation. At the rear part, large tension cracks nearly parallel to the slope surface develop, resulting in combined tension-bendingtoppling deformation (2) Cracks develop and energy changes drastically during 2 -5 s, when the earthquake has the strongest motion. Cracks are dominantly tension cracks, and shear cracks develop only along joints and faults. The energy dissipates mainly through damping and friction (3) The impact of seismic intensity is studied. Seismic loading with large intensity (in terms of PGA) leads to more cracks and crushed rocks but has no obvious influence on the failure mode (4) Different loading methods in PFC2D can make a difference in the numerical results. By comparison, loading by applying velocity on the bottom particles or on the bottom boundary rigid wall may cause instability along the interface of the highly and weakly toppled rock mass, but this is not consistent with the actual situation

Data Availability
All the data supporting the results of this study have been presented in the paper.

Conflicts of Interest
The authors declare that they have no conflicts of interest.