Effects of Three-Directional Seismic Wave on Dynamic Response and Failure Behavior of High-Steep Rock Slide

The landslide triggered by earthquakes can cause severe infrastructure losses or even fatalities. The high-steep rock slide is the most common type of landslide in the earthquake area. In an earthquake, the ground moves randomly in all directions, two horizontal directions (East-West (EW) direction, North-South (NS) direction) and one vertical direction (Up-Down (UD) direction). Even though extensive studies have been carried out on the earthquake-triggered landslide, the effects of each single seismic wave and the three-directional seismic waves are not considered. This study aims to evaluate the effects of different types of the seismic waves on the dynamic response and failure behavior of the high-steep rock slide. To investigate the effects of each single seismic wave and three-directional seismic wave, this study presents a numerical model with four types of seismic waves, e.g., East-West (EW) direction, North-South (NS) direction, Up-Down (UD) direction, and three-directional wave (EW_NS_UD). The numerical results revealed that the types of the seismic waves have significantly different effects on the dynamic process, failure behavior, run-out distance, velocity, and deposition of the high-steep rock slide.


Introduction
The earthquake-triggered landslide represents a dangerous natural hazard causing significant damage to property and infrastructure and personal casualties [1,2]. There have been many reports on the earthquake-triggered landslides in the past few years [3][4][5]. For example, on 8 October 2005, the Kashmir earthquake triggered a large number of landslides in the northern Pakistan. The most disastrous landslide destroyed three villages and caused nearly 1000 deaths [6]. A total of 60,104 landslides induced by the 2008 Wenchuan earthquake (Ms = 8.0) directly killed over 20,000 people [7]. Less than five years later, the 2017 Jiuzhaigou earthquake shocked the Sichuan province, China, and induced as many as 179 landslides [8,9]. In spite of their significant damage to property and human casualties, the cause and effect of earthquake-triggered landslides is not well understood [10,11].
The process of earthquake-triggered landslides is complex, which is also called the dynamic response of the landslide [12]. When the seismic wave reaches the ground surface, the dynamic behavior of the seismic wave will be affected and modified by the soil.
Thus, it appears possible to reduce the seismic response by intervening and modifying the characteristics of the soil [13][14][15]. The failure behavior of the earthquake-triggered landslide is also closely related to its dynamic response under the seismic waves. The earliest study on the effect of the seismic waves can be traced back to 1965, when Newmark proposed a method permitting a rapid estimate of the dams and embankments displacement in an earthquake [16]. After that, many studies focusing on the effects of seismic waves on the earthquake-triggered landslide dynamic response and failure behavior have been published [17][18][19][20][21]. For example, Chuang et al. developed a nowcasting model for estimating the probability of landslides induced by earthquakes [22]. Zou et al. analyzed fault (TZF), the Xueshan fault (XSF), and the Minjiang fault (MJF), the north-east and the north-west corners are the Songpan-Ganzi orogenic zone, the Animaqing orogenic zone (anticline), and the Malkang thrust-detachment rock slice (syncline). The middle part of the Jiuzhaigou area is the Motianling nappe in the west orogenic zone. The Minjiang fault (MJF) is located in the west side of the epicenter. The overall trend is north-south, and the section is inclined to the west. It shows upthrust and left-lateral strike-slip action, indicating that the Minjiang fault (MJF) is not the seismogenic fault of the earthquake. However, it has a certain restrictive effect on the aftershocks, surface, and the western boundary of the earthquake landslide distribution of the earthquake. The Tazang fault (TZF) is located in the north side of the epicenter. The overall trend is northwest-southeast, and the fracture dip angle is 50 • -60 • . It also had a specific limit on the aftershocks of the earthquake, surface deformation, and the northern boundary of the distribution of the earthquake area. Therefore, it can be inferred that the seismogenic fault of this earthquake was the Huya fault (HYF).
which caused disasters in eight counties in the Sichuan and Gansu provinces. The earthquake caused 25 deaths, injured 525 people, caused six people to become lost, and damaged 73,671 houses.
The Jiuzhaigou area is located at the junction of the Danba-Wenchuan, Songpan-Ganzi and the Motianling in the West Qinling tectonic zone ( Figure 1). Bounded by the Tazang fault (TZF), the Xueshan fault (XSF), and the Minjiang fault (MJF), the north-east and the north-west corners are the Songpan-Ganzi orogenic zone, the Animaqing orogenic zone (anticline), and the Malkang thrust-detachment rock slice (syncline). The middle part of the Jiuzhaigou area is the Motianling nappe in the west orogenic zone. The Minjiang fault (MJF) is located in the west side of the epicenter. The overall trend is northsouth, and the section is inclined to the west. It shows upthrust and left-lateral strike-slip action, indicating that the Minjiang fault (MJF) is not the seismogenic fault of the earthquake. However, it has a certain restrictive effect on the aftershocks, surface, and the western boundary of the earthquake landslide distribution of the earthquake. The Tazang fault (TZF) is located in the north side of the epicenter. The overall trend is northwest-southeast, and the fracture dip angle is 50°-60°. It also had a specific limit on the aftershocks of the earthquake, surface deformation, and the northern boundary of the distribution of the earthquake area. Therefore, it can be inferred that the seismogenic fault of this earthquake was the Huya fault (HYF).

The High-Steep Rock Slide
The high-steep slope is one of the most common typical landforms in mountainous areas, which refers to the slope with an approximately vertical angle. The high-steep slope is generally composed of horizontal sedimentary rocks. Due to its steep terrain characteristic, rock slides often develop on the surface of the high-steep slope. The high-steep rock slide is one of the most common type of the landslide in the Jiuzhaigou area. For example,

The High-Steep Rock Slide
The high-steep slope is one of the most common typical landforms in mountainous areas, which refers to the slope with an approximately vertical angle. The high-steep slope is generally composed of horizontal sedimentary rocks. Due to its steep terrain characteristic, rock slides often develop on the surface of the high-steep slope. The high-steep rock slide is one of the most common type of the landslide in the Jiuzhaigou area. For example, Figure 2a shows the perilous rock mass at a high-steep slope on the right bank of the Swan-sea central highway (33.09 • N, 103.66 • E). Figure 2b shows the overall view of the rock slide at multiple points on the right bank of the Swan-sea central highway. Appl. Sci. 2022, 11, x FOR PEER REVIEW 4 of 25 Figure 2a shows the perilous rock mass at a high-steep slope on the right bank of the Swan-sea central highway (33.09° N, 103.66° E). Figure 2b shows the overall view of the rock slide at multiple points on the right bank of the Swan-sea central highway. The perilous rock mass on the surface of the high-steep slope is cut by vertical tensile joints or steeply dipping columnar joints [32]. It separates from the stable base rock and forms tall, narrow, upright, laminated rock columns. Finally, the rock columns are divided into blocks by two sets of nearly orthogonal structural planes, as shown in Figure 3. The strata is the general term for all stratified rocks. The high-steep slope is most prone to rock slide. Under the action of the earthquake, the perilous rock mass suddenly separates from the base rock. The stability of the high-steep rock slide is mainly affected by the overturning moment, as shown in Figure 4. When the overturning moment is greater than the antititling moment, a toppling failure occurs, and then it will continue to collapse, roll, fall, etc. The perilous rock mass on the surface of the high-steep slope is cut by vertical tensile joints or steeply dipping columnar joints [32]. It separates from the stable base rock and forms tall, narrow, upright, laminated rock columns. Finally, the rock columns are divided into blocks by two sets of nearly orthogonal structural planes, as shown in Figure 3. The strata is the general term for all stratified rocks. The high-steep slope is most prone to rock slide. Under the action of the earthquake, the perilous rock mass suddenly separates from the base rock. The stability of the high-steep rock slide is mainly affected by the overturning moment, as shown in Figure 4. When the overturning moment is greater than the anti-titling moment, a toppling failure occurs, and then it will continue to collapse, roll, fall, etc.

Overview
The process of the landslide failure involves the separation of the blocks and the reconstruction of the contact surface. Therefore, the numerical method requires that it allows the blocks can move freely. The discrete element method (DEM) is a numerical method especially suitable for stress analysis of jointed rock masses [33]. The key concept of DEM is based on Newton's law, which treats the research object as an assemblage of the rigid or deformable blocks, and the location of each block is updated during the entire deformation/motion process [34]. Therefore, the DEM is widely used in the process of the dynamic response and failure behavior of the high-steep rock slide. The DEM program called 3DEC is used for numerical modeling.
To investigate effects of the dynamic response and failure behavior of the high-steep rock slide under different seismic wavs, a typical high-steep rock slide numerical model is built in 3DEC. As illustrated in Figure 5, the research process is divided into three main parts: preprocessing of the numerical modeling (steps 1-2), numerical modeling (steps 3-6), and analysis of the numerical results (step 7). Specifically, the research steps are (1) filtering and baseline correction of the seismic waves; (2) calculation of the appropriate mesh size; (3) building of the high-steep rock slide; (4) determination of the numerical parameters; (5) inputting of the processed seismic waves; (6) setting of the monitors of the key blocks; and (7) analyzing of the numerical results.

Details of Dynamic Modeling Considerations in the Discrete Element Method
There are three issues should be considered when preparing a DEM model for a dynamic analysis, e.g., boundary conditions, local damping, and max unbalanced force. The boundary conditions are essential in 3DEC dynamic analysis. Figure 6 shows the setting of the slope boundary conditions, in which the viscous boundary is to set up normal and tangential viscous dampers on the boundary of the model. The normal and tangential viscous dampers can provide unbalanced force to absorb external incident waves. The unbalanced force will be applied to the viscous boundary, and the forces in the three coordinate directions are as follows:

Details of Dynamic Modeling Considerations in the Discrete Element Method
There are three issues should be considered when preparing a DEM model for a dynamic analysis, e.g., boundary conditions, local damping, and max unbalanced force. The boundary conditions are essential in 3DEC dynamic analysis. Figure 6 shows the setting of the slope boundary conditions, in which the viscous boundary is to set up normal and tangential viscous dampers on the boundary of the model. The normal and tangential viscous dampers can provide unbalanced force to absorb external incident waves. The unbalanced force will be applied to the viscous boundary, and the forces in the three coordinate directions are as follows: where, ρ is the mass density; C s and C p are the velocities of longitudinal and transverse waves; and A is the area affected by the viscous boundary grid. v m  (1) where, ρ is the mass density; Cs and Cp are the velocities of longitudinal and transverse waves; and A is the area affected by the viscous boundary grid. , , / , , and are the three direction velocity components of the lateral boundary nodes of the main grid/viscous boundary network, respectively. , , and are the forces on the viscous boundary grid. Therefore, the wave propagation will not be affected by the distortion caused by the boundary conditions. The local damping is used in the 3DEC dynamic analysis to reduce the unnecessary numerical vibrations block motion system. In general, the vibration amplitude of any vibration system will decrease due to air resistance, friction, etc. However, in the 3DEC dynamic analysis, the kinetic energy cannot be dissipated naturally. Therefore, the appropriate local damping is added to ensure the normal vibrations of the numerical model. It takes many repeated trials to obtain an appropriate value. In this study, the local damping ratio is set to 0.005.
The maximum unbalanced force is monitored to detect whether the system has reached a state of equilibrium. The ideal model balance means that the nodal force vector at the center of the block is exactly zero. However, the maximum unbalanced force in numerical analysis can never reach zero. Generally, it is considered that the maximum unbalanced force is sufficiently small compared with the representative force of the model. It represents the numerical modeling has reached a state of equilibrium. The maximum unbalanced force in this study is set to 0.00002.

Building of the High-Steep Rock Slide Numerical Model
As shown in Figure 2a To facilitate the study of the effect of three-directional seismic waves on the highsteep rock slide, the numerical model is simplified comparing with the actual situation. The local damping is used in the 3DEC dynamic analysis to reduce the unnecessary numerical vibrations block motion system. In general, the vibration amplitude of any vibration system will decrease due to air resistance, friction, etc. However, in the 3DEC dynamic analysis, the kinetic energy cannot be dissipated naturally. Therefore, the appropriate local damping is added to ensure the normal vibrations of the numerical model. It takes many repeated trials to obtain an appropriate value. In this study, the local damping ratio is set to 0.005.
The maximum unbalanced force is monitored to detect whether the system has reached a state of equilibrium. The ideal model balance means that the nodal force vector at the center of the block is exactly zero. However, the maximum unbalanced force in numerical analysis can never reach zero. Generally, it is considered that the maximum unbalanced force is sufficiently small compared with the representative force of the model. It represents the numerical modeling has reached a state of equilibrium. The maximum unbalanced force in this study is set to 0.00002.

Building of the High-Steep Rock Slide Numerical Model
As shown in Figure 2a To facilitate the study of the effect of three-directional seismic waves on the high-steep rock slide, the numerical model is simplified comparing with the actual situation. Figure 7a shows the numerical model according to Figure 2a. The numerical modeling adopts the left-handed coordinate system. In the discrete element analysis, multiple sets of horizontal and longitudinal structural planes are often used to form discrete blocks so that each block can move freely. One joint set can divide the perilous rock into planes, while two joint sets can separate the perilous into columns. Therefore, we need three joint sets to cut the rock into blocks. Three sets of discontinuities cut the perilous rock mass, creating an assemblage of blocks with a spacing of 4 m. Three joints set have dip direction/dip angles of 90 • ∠20 • , 0 • ∠90 • , and 90 • ∠90 • . Considering the lateral movement of the landslide, the size of the base rock is 200 m in width, which is much larger than the 20-m width perilous rock mass. Furthermore, seven monitors are set to explain the ground motion response of the failure process and the movement characteristics clearly, as shown in Figure 7b. Appl Figure 7a shows the numerical model according to Figure 2a. The numerical modeling adopts the left-handed coordinate system. In the discrete element analysis, multiple sets of horizontal and longitudinal structural planes are often used to form discrete blocks so that each block can move freely. One joint set can divide the perilous rock into planes, while two joint sets can separate the perilous into columns. Therefore, we need three joint sets to cut the rock into blocks. Three sets of discontinuities cut the perilous rock mass, creating an assemblage of blocks with a spacing of 4 m. Three joints set have dip direction/dip angles of 90° 20°, 0° 90°, and 90° 90°. Considering the lateral movement of the landslide, the size of the base rock is 200 m in width, which is much larger than the 20-m width perilous rock mass. Furthermore, seven monitors are set to explain the ground motion response of the failure process and the movement characteristics clearly, as shown in Figure 7b.  Table 1 illustrates the physical and mechanical parameters for the numerical analysis [26]. All blocks and joints follow the linearly elastic model and Coulomb-slip model, respectively. It is assumed that the numerical model is in a residual state because only the parameters after slope failure can be obtained.   Table 1 illustrates the physical and mechanical parameters for the numerical analysis [26]. All blocks and joints follow the linearly elastic model and Coulomb-slip model, respectively. It is assumed that the numerical model is in a residual state because only the parameters after slope failure can be obtained.   Figure 8 shows the original EW (East-West direction), NS (North-South direction), and UD (Up-Down direction) seismic waves used in the numerical model. These curves represent the acceleration time histories of the seismic waves. They are obtained from the 2017 Jiuzhaigou earthquake, corresponding to the X, Z, and Y directions of the dynamic load input during the simulation process. The geographic coordinates of the selected seismic wave station are 32.8° N 105.4° E. The dataset is provided by the China Earthquake Networks Center, National Earthquake Data Center (http://data.earthquake.cn (accessed on 20 March 2020)). The frequency component of the input waveform, and the wave velocity characteristics of the rock mass will affect the numerical accuracy of wave propagation. To describe the wave propagation in the model accurately, the study by Kuhlemeye and Lysmer shows that the zone size must be less than 1/8~1/10 of the wavelength, corresponding to the highest frequency of the input waveform [35]. On account of the maximum frequency of seismic waves has a more significant impact on the zone size, filtering seismic waves is required. The higher the maximum frequency is, the smaller the zone size is. The aim of filtering seismic wave is to filter out the high-frequency components of the original waveform and reduce the maximum frequency of small seismic waves, thereby increasing the minimum zone size for less calculation time.

Processing of the Input Seismic Wave Loadings
As an illustration, Figure 9a shows the seismic wave filtering of the EW seismic wave. The frequency corresponding to peak acceleration is mainly concentrated in the range of 10~15 Hz. Therefore, the part of the waveform greater than 20 Hz can be filtered by Wizard software. The maximum frequency of the filtered waveform is about 20 Hz, as shown in the red line in The frequency component of the input waveform, and the wave velocity characteristics of the rock mass will affect the numerical accuracy of wave propagation. To describe the wave propagation in the model accurately, the study by Kuhlemeye and Lysmer shows that the zone size must be less than 1/8~1/10 of the wavelength, corresponding to the highest frequency of the input waveform [35]. On account of the maximum frequency of seismic waves has a more significant impact on the zone size, filtering seismic waves is required. The higher the maximum frequency is, the smaller the zone size is. The aim of filtering seismic wave is to filter out the high-frequency components of the original waveform and reduce the maximum frequency of small seismic waves, thereby increasing the minimum zone size for less calculation time.
As an illustration, Figure 9a shows the seismic wave filtering of the EW seismic wave. The frequency corresponding to peak acceleration is mainly concentrated in the range of 10~15 Hz. Therefore, the part of the waveform greater than 20 Hz can be filtered by Wizard software. The maximum frequency of the filtered waveform is about 20 Hz, as shown in the red line in Figure 9a. The equation for calculating the shear wave velocity is C s = G/ρ, and the equation for calculating the longitudinal wave velocity is C p = K+4G/3 ρ , where G, K and ρ represent the model shear modulus, bulk modulus, and density, respectively. According to the input model parameters, the shear and longitudinal wave velocities of this model can be estimated to be 2435 m/s and 4551 m/s, respectively. According to the calculation equation of wavelength λ = C s / f (λ = C p / f ), where C s (C p ), f, and λ represent the shear (longitudinal) wave velocity, highest waveform frequency, and wavelength, respectively. And the wavelength of this seismic wave can be estimated to be from 121 m to 228 m. Therefore, the zone size is set to 5 m for the perilous rock mass and 15 m for the base rock, satisfying the calculation requirements. velocity and residual displacement will continue to appear at the end of the dynamic calculation. Baseline correction is usually conducted to eliminate the integration errors for time-domain seismic analysis by adding a low-frequency waveform to the original velocity time history. Figure 9b shows the baseline correction progress for EW seismic waves, and Figure 10 illustrates the final input of seismic waves of EW, NS, UD, and EW_NS_UD in 3DEC.  The velocity-time history within 45-55 s when the acceleration peak appears is intercepted as the input dynamic loading, and the bottom boundary of the model is the seismic application source. In the 3DEC seismic dynamic analysis, if the final input velocity is not zero, or the final displacement obtained by integrating the input velocity is not zero, the velocity and residual displacement will continue to appear at the end of the dynamic calculation. Baseline correction is usually conducted to eliminate the integration errors for time-domain seismic analysis by adding a low-frequency waveform to the original velocity time history. Figure 9b shows the baseline correction progress for EW seismic waves, and Figure 10 illustrates the final input of seismic waves of EW, NS, UD, and EW_NS_UD in 3DEC. Appl

Dynamic Response of the High-Steep Rock Slide
As an illustration, the velocity contour of the EW seismic wave propagation process is depicted in Figure 11. The numerical results show that the dynamic response of the slope does not occur at the same time. First, the velocity is generated at the model's bottom, as shown in Figure 11a, and then transmits from the bottom to the top. The propaga-

Dynamic Response of the High-Steep Rock Slide
As an illustration, the velocity contour of the EW seismic wave propagation process is depicted in Figure 11. The numerical results show that the dynamic response of the slope does not occur at the same time. First, the velocity is generated at the model's bottom, as shown in Figure 11a, and then transmits from the bottom to the top. The propagation process has a strip-like characteristic: with the growing altitude, the velocity increases.
When the velocity field reaches the slope surface, the surface velocity is larger than the internal velocity at the same altitude, as shown in Figure 11b,c. It can be inferred that the input seismic waves are affected by vertical and surface amplification of the slopes. The other three cases have the similar dynamic response phenomenon. Appl. Sci. 2022, 11, x FOR PEER REVIEW 13 of 25 tion process has a strip-like characteristic: with the growing altitude, the velocity increases. When the velocity field reaches the slope surface, the surface velocity is larger than the internal velocity at the same altitude, as shown in Figure 11b,c. It can be inferred that the input seismic waves are affected by vertical and surface amplification of the slopes. The other three cases have the similar dynamic response phenomenon.  Figure 11d shows the velocity field transmits to the top after 7000 steps and lasts about 0.13 s, equal to the calculation result. The dynamic response steps/time of the other three cases are NS: 8000/0.15 s, UD: 4000/0.07 s, and EW_NS_UD: 4000/0.07 s, also confirming the calculation results. The simulation results reveal that the seismic wave propagation time is affected by the type of seismic wave, and the propagation time of the EW_NS_UD seismic wave is consistent with the propagation time of the UD seismic single wave. We infer that the first velocity field reaching the top under EW_NS_UD seismic wave is the UD seismic wave, and the three seismic waves act independently. To verify this hypothesis, the velocity contour of UD wave, EW_NS_UD wave, and EW wave under 4000 steps are recorded, respectively. The velocity fields of monitor 1 and monitor 6 are selected for comparative analysis, as shown in Figure 12. Here, let Vn denote the velocity of monitor n. It can be seen that the velocity of monitor 1/6 under the EW wave is both 0. It shows that under the EW single wave, the velocity field has not yet reached the position of monitor 1 and monitor 6. The velocities of monitor 1/6 under the UD seismic wave are 0.12 and 0.12 m/s, respectively. The velocities of monitor 1/6 under the EW_NS_UD wave are 0.13 and 0.14 m/s, respectively. Those two results are almost consistent. That is to say, at 0.07 s, the velocity fields generated by the UD wave and the EW_NS_UD wave at monitor 1 and monitor 6 are the same. The numerical results satisfy the hypothesis.  Figure 11d shows the velocity field transmits to the top after 7000 steps and lasts about 0.13 s, equal to the calculation result. The dynamic response steps/time of the other three cases are NS: 8000/0.15 s, UD: 4000/0.07 s, and EW_NS_UD: 4000/0.07 s, also confirming the calculation results. The simulation results reveal that the seismic wave propagation time is affected by the type of seismic wave, and the propagation time of the EW_NS_UD seismic wave is consistent with the propagation time of the UD seismic single wave. We infer that the first velocity field reaching the top under EW_NS_UD seismic wave is the UD seismic wave, and the three seismic waves act independently. To verify this hypothesis, the velocity contour of UD wave, EW_NS_UD wave, and EW wave under 4000 steps are recorded, respectively. The velocity fields of monitor 1 and monitor 6 are selected for comparative analysis, as shown in Figure 12. Here, let V n denote the velocity of monitor n. It can be seen that the velocity of monitor 1/6 under the EW wave is both 0. It shows that under the EW single wave, the velocity field has not yet reached the position of monitor 1 and monitor 6. The velocities of monitor 1/6 under the UD seismic wave are 0.12 and 0.12 m/s, respectively. The velocities of monitor 1/6 under the EW_NS_UD wave are 0.13 and 0.14 m/s, respectively. Those two results are almost consistent. That is to say, at 0.07 s, the velocity fields generated by the UD wave and the EW_NS_UD wave at monitor 1 and monitor 6 are the same. The numerical results satisfy the hypothesis. Appl

Slope Failure Scenario of the High-Steep Rock Slide
The failure behaviors of high-steep rock slide under cases 1/4 and cases 2/3 are different, which will be described separately below. Figure 13a,b show the post-failure behaviors of landslides in cases 1/4. The high-steep rock slide can be described by four-stage scenario.
(1) Stage 1 (peeling): In the time of 0-4/0-3 s, the front edge rock peels and collapses. At this stage, the front edge rock slides vertically. The rock mass collapses after it crashes the base rock.
(2) Stage 2 (toppling): In the time of 5-14/4-12 s, the rear rock cracks propagate, and then, the rock dumps. The crack propagation occurs on the sliding surface in contact with the base rock, and then, the rock begins to rotate with the bottom of the rock as the pivot. At the critical time 1, two vertical cracks develop inside the dangerous rock. At the critical time 2, the upper parts of the rear two blocks dump at the 1/3 height of the blocks. At the critical time 3, the front block breaks at the same position.
(3) Stage 3 (collapsing): In the time of 15-37/14-32 s, the perilous rock ejects and disintegrates. Stage 3 has two sub-stages: at the first sub-stage (15-21/14-17 s), more than 1/3 height rock mass crashes the slope, forming collapse debris. Then, it rolls down and leaves the bottom rock remaining. At the second sub-stage (21-37/17-32 s), with the complete dumping of the upper block, the lower block is gradually out of balance and then collapses at the critical time 4, causing the secondary collapse of the perilous rock mass.
(4) Stage 4 (depositing): This stage is the deposition process, and the convergence time is 79/82 s. Figure 13c,d show the different failure behaviors of the rock slides in cases 2/3, which are not completely failing. The first stage is similar to that of cases 1/4, and both of the first stage times of NS/UD are 0-3 s. At the second stage, only the front crack propagates and dumps at the critical time 1 (NS/UD: 14/16 s). Compared to the dumping point with the cases 1/4, the cases 2/3 dumps at the bottom of rock mass, and the dumping phenomenon is not very apparent. Then, rock mass collides and accumulates until the system reaches equilibrium at 126/77 s, NS/UD.
Comparing the failure behavior of these four cases, it can be observed that EW and EW_NS_UD seismic waves have more vital damage to slope stability than NS and UD seismic waves. It can be concluded that the EW wave plays a dominant role in the failure behavior of landslides in three-directional seismic waves. Still, the existence of NS/UD seismic waves will accelerate the failure process of landslides.

Slope Failure Scenario of the High-Steep Rock Slide
The failure behaviors of high-steep rock slide under cases 1/4 and cases 2/3 are different, which will be described separately below. Figure 13a Comparing the failure behavior of these four cases, it can be observed that EW and EW_NS_UD seismic waves have more vital damage to slope stability than NS and UD seismic waves. It can be concluded that the EW wave plays a dominant role in the failure behavior of landslides in three-directional seismic waves. Still, the existence of NS/UD seismic waves will accelerate the failure process of landslides. Appl Figure 14 shows the velocity-time history curves in four cases, and Figure 15 shows the comparison of maximum velocities of the monitors and their appearance time in these four cases. Here, let V n denote the velocity of monitor n (n = 1-7). From the two horizontally arranged monitors 1/2/3 and monitors 4/7, we can observe that, in cases 2, 3, and 4, V 3 (dark blue line) is always greater than V 2 (red line), while it is approximately equal to V 2 in case 1. V 1 (black line) is the maximum value of the three monitors in case 1, while it is the minimum value in case 2. In conclusion, the velocity of the edge of the slope is more significant than that of the inside. V 7 (light blue line) is always 0 in cases 2/3. Therefore, only V 4 (green line) and V 7 in cases 1/4 are compared and analyzed. The result shows that V 7 is greater than V 4 in case 1, while V 7 and V 4 are approximately equal in case 4. We investigated that the velocity relationship between V 4 and V 7 was not apparent. The time for the max velocity appearance is that V 3 is earlier than V 2 , earlier than V 1 , and V 4 is earlier than V 7 .

Analysis of the Velocities of the Monitors
From the longitudinally arranged monitors 3/4/5/6 and monitors 1/7, it can be found that except for V 4 being more significant than V 3 in case 1, V 3 > V 4 > V 5 (purple line) in other cases. V 1 is more significant than V 7 in case 4, while it is the opposite in case 1. We conclude that the height of the landslide is positively correlated with the peak velocity. The velocity relationship between V 6 (yellow line) and the others is not apparent, which is just like V 4 and V 7 , so we analyzed that the velocity of the sliding surface is not related to others. The time for the max velocity appearance is that V 3 is earlier than V 4 , earlier than V 5 , earlier than V 6 , and V 1 is earlier than V 7 . Figure 16 shows the displacement time history curves in four cases. The monitors without displacement are not shown in the figure. The number of effective monitors is 7/6/5/4 for four cases, respectively, and the displacement trend is similar. Here, let D n denote the displacement of monitor n (n = 1-7). Appl Figure 14 shows the velocity-time history curves in four cases, and Figure 15 shows the comparison of maximum velocities of the monitors and their appearance time in these four cases. Here, let Vn denote the velocity of monitor n (n = 1-7). From the two horizontally arranged monitors 1/2/3 and monitors 4/7, we can observe that, in cases 2, 3, and 4, V3 (dark blue line) is always greater than V2 (red line), while it is approximately equal to V2 in case 1. V1 (black line) is the maximum value of the three monitors in case 1, while it is the minimum value in case 2. In conclusion, the velocity of the edge of the slope is more significant than that of the inside. V7 (light blue line) is always 0 in cases 2/3. Therefore, only V4 (green line) and V7 in cases 1/4 are compared and analyzed. The result shows that V7 is greater than V4 in case 1, while V7 and V4 are approximately equal in case 4. We investigated that the velocity relationship between V4 and V7 was not apparent. The time for the max velocity appearance is that V3 is earlier than V2, earlier than V1, and V4 is earlier than V7.

Analysis of the Velocities of the Monitors
From the longitudinally arranged monitors 3/4/5/6 and monitors 1/7, it can be found that except for V4 being more significant than V3 in case 1, V3 > V4 > V5 (purple line) in other cases. V1 is more significant than V7 in case 4, while it is the opposite in case 1. We conclude that the height of the landslide is positively correlated with the peak velocity. The velocity relationship between V6 (yellow line) and the others is not apparent, which is just like V4 and V7, so we analyzed that the velocity of the sliding surface is not related to others. The time for the max velocity appearance is that V3 is earlier than V4, earlier than V5, earlier than V6, and V1 is earlier than V7.  Figure 16 shows the displacement time history curves in four cases. The monitors without displacement are not shown in the figure. The number of effective monitors is 7/6/5/4 for four cases, respectively, and the displacement trend is similar. Here, let Dn denote the displacement of monitor n (n = 1-7).  Figure 16 shows the displacement time history curves in four cases. The monitors without displacement are not shown in the figure. The number of effective monitors is 7/6/5/4 for four cases, respectively, and the displacement trend is similar. Here, let Dn denote the displacement of monitor n (n = 1-7). As an illustration, Figure 17 shows the characteristics of case 1. From the two horizontally arranged monitors 1/2/3 ( Figure 17a) and monitors 4/7 (Figure 17b), we can observe that the value of displacement is D3 > D1 > D2 (D4 > D7), and the time for reaching the maximum displacement is that D3 is earlier than D2, which is earlier than D1 (D4 earlier than D7). The results consist with the previous section conclusion, "the velocity of the slope's edge is larger than that of the inside". The results show that (1) the rear rock has a more significant displacement than the front rock, and the displacement of the edge of the slope is more significant than that of the inside; and (2) the stable time is inversely related to the distance from the edge rock. The other three cases also have the similar conclusion. It is worth mentioning that in case 4, D1 is the largest, followed by D3 and D2.

Analysis of the Displacement of the Monitors
From the longitudinally arranged monitors 3/4/5/6 ( Figure 17c) and monitors 1/7 (Figure 17d), we can observe the displacement characteristics along the vertical direction. The order of displacement from large to small is D3 ≈ D4 > D5 > D6 (D1 > D7), and the stable time is D3 < D4 < D5 < D6 (D1 < D7). The results show that the displacement values are related to slope altitude under seismic loading, while stable time is inversely related, as is the case in the other three cases. As an illustration, Figure 17 shows the characteristics of case 1. From the two horizontally arranged monitors 1/2/3 ( Figure 17a) and monitors 4/7 (Figure 17b), we can observe that the value of displacement is D 3 > D 1 > D 2 (D 4 > D 7 ), and the time for reaching the maximum displacement is that D 3 is earlier than D 2 , which is earlier than D 1 (D 4 earlier than D 7 ). The results consist with the previous section conclusion, "the velocity of the slope's edge is larger than that of the inside". The results show that (1) the rear rock has a more significant displacement than the front rock, and the displacement of the edge of the slope is more significant than that of the inside; and (2) the stable time is inversely related to the distance from the edge rock. The other three cases also have the similar conclusion. It is worth mentioning that in case 4, D 1 is the largest, followed by D 3 and D 2 .
From the longitudinally arranged monitors 3/4/5/6 ( Figure 17c) and monitors 1/7 (Figure 17d), we can observe the displacement characteristics along the vertical direction. The order of displacement from large to small is D 3 ≈ D 4 > D 5 > D 6 (D 1 > D 7 ), and the stable time is D 3 < D 4 < D 5 < D 6 (D 1 < D 7 ). The results show that the displacement values are related to slope altitude under seismic loading, while stable time is inversely related, as is the case in the other three cases. Figure 18a shows no significant difference in the displacement value of D3 in four cases. The difference in displacement gradually increases as the altitude decreases. The order of D4, which is case 1 > 2 > 3 > 4, has the largest displacement difference in different cases. D5 is the same as D4, and the displacement order is case 2 > 3 > 1 ≈ 4. In conclusion, the response of different seismic waves at low altitudes is different, and the effects of seismic waves in high-altitude areas are equivalent. Figure 18a shows no significant difference in the displacement value of D 3 in four cases. The difference in displacement gradually increases as the altitude decreases. The order of D 4 , which is case 1 > 2 > 3 > 4, has the largest displacement difference in different cases. D 5 is the same as D 4 , and the displacement order is case 2 > 3 > 1 ≈ 4. In conclusion, the response of different seismic waves at low altitudes is different, and the effects of seismic waves in high-altitude areas are equivalent.  Figure 19 shows the rock slide deposition in four cases, and Table 2 presents the maximum displacement value (X direction and Z direction) and the main deposition area value. It shows that the deposition area of case 1 and case 4 are similar, which proves the view that the EW wave plays a dominant role in the three-directional seismic wave mentioned before. The farthest distance (X direction) of the landslide is in case 1, while in other cases, they are very close. It shows that the farthest landslide distance in the three-directional seismic wave is not necessarily more significant than others, for the possible reason that various waves may be superimposed or offset. In addition, the farthest distance (Z direction) in case 2 and the nearest is case 1. It shows that the displacement of the rock slide along the loading direction is the largest, while the vertical loading direction is the smallest.   Figure 19 shows the rock slide deposition in four cases, and Table 2 presents the maximum displacement value (X direction and Z direction) and the main deposition area value. It shows that the deposition area of case 1 and case 4 are similar, which proves the view that the EW wave plays a dominant role in the three-directional seismic wave mentioned before. The farthest distance (X direction) of the landslide is in case 1, while in other cases, they are very close. It shows that the farthest landslide distance in the threedirectional seismic wave is not necessarily more significant than others, for the possible reason that various waves may be superimposed or offset. In addition, the farthest distance (Z direction) in case 2 and the nearest is case 1. It shows that the displacement of the rock slide along the loading direction is the largest, while the vertical loading direction is the smallest.

Impact of the Local Damping Ratio
The local damping in this study is set to 0.005, which is obtained from many repeated trials. As mentioned in Section 3.2, local damping is used to reduce unnecessary numerical vibrations. If the local damping is too large, the rock slide will remain stable under all types of seismic waves. If the local damping is too small, the rock slide will float in the air, which is obviously unrealistic.
In the local damping ratio sensitivity analysis, a three-directional seismic wave is adopted. The value of damping ratio is conducted with 0.001, 0.005, and 0.01. The stability of the perilous rock mass increases with the growing value of local damping ratio ( Figure  20). When the value of the local damping is 0.001, the rock slide will float in the air, and the 3DEC computations will diverge (Figure 20a). When the value of the local damping is 0.01, the rock slide is not completely failing (Figure 20b).

Impact of the Local Damping Ratio
The local damping in this study is set to 0.005, which is obtained from many repeated trials. As mentioned in Section 3.2, local damping is used to reduce unnecessary numerical vibrations. If the local damping is too large, the rock slide will remain stable under all types of seismic waves. If the local damping is too small, the rock slide will float in the air, which is obviously unrealistic.
In the local damping ratio sensitivity analysis, a three-directional seismic wave is adopted. The value of damping ratio is conducted with 0.001, 0.005, and 0.01. The stability of the perilous rock mass increases with the growing value of local damping ratio ( Figure 20). When the value of the local damping is 0.001, the rock slide will float in the air, and the 3DEC computations will diverge (Figure 20a). When the value of the local damping is 0.01, the rock slide is not completely failing (Figure 20b). Appl

Impact of the Dip Angle of Joints
The dip angle of joints might influence the failure behavior of the rock slide. In the joint dip angle sensitivity analysis, three-directional seismic wave is adopted. The dip angle in this study is conducted with 0°, 20°, and 60°. The failure behavior of joint dip angle 20° is shown as Figure 13b. Figure 21 shows the critical failure behavior of joint dip angle 0° and 60°. The failure behavior of the rock slide changes from toppling movement to sliding movement with the increasing joint dip angle.

Discussions
This study evaluates the dynamic response of a typical type of landslides, that is, high-steep rock slide, under different types of seismic waves. Part of the analysis of the numerical results have been discussed in previous section, and they will be described as conclusions in next section as well. In addition, there are still some discussions worth of attention.
(1) Selection of the input seismic waves: At time of the Jiuzhaigou earthquake, China Earthquake Networks Center collected the seismic waves of 66 stations. The seismic waves selection criteria in this study are as follows: the station as close to the study area as possible. However, the recorded seismic waves are quite different in the total 66 stations, e.g., the amplitude. Therefore, it is worth discussing whether it is available to select seismic waves only by the value of distance.
(2) The validity of the typical high-steep rock slide numerical modeling: This study is to explore the regularity of the dynamic response of the typical highsteep rock slide under different seismic waves. Therefore, the numerical model is simplified based on the actual landslide. However, the conclusions investigated in this study are consistent with those obtained in other studies. For example, seismic waves have edge

Impact of the Dip Angle of Joints
The dip angle of joints might influence the failure behavior of the rock slide. In the joint dip angle sensitivity analysis, three-directional seismic wave is adopted. The dip angle in this study is conducted with 0 • , 20 • , and 60 • . The failure behavior of joint dip angle 20 • is shown as Figure 13b. Figure 21 shows the critical failure behavior of joint dip angle 0 • and 60 • . The failure behavior of the rock slide changes from toppling movement to sliding movement with the increasing joint dip angle.

Impact of the Dip Angle of Joints
The dip angle of joints might influence the failure behavior of the rock slide. In the joint dip angle sensitivity analysis, three-directional seismic wave is adopted. The dip angle in this study is conducted with 0°, 20°, and 60°. The failure behavior of joint dip angle 20° is shown as Figure 13b. Figure 21 shows the critical failure behavior of joint dip angle 0° and 60°. The failure behavior of the rock slide changes from toppling movement to sliding movement with the increasing joint dip angle.

Discussions
This study evaluates the dynamic response of a typical type of landslides, that is, high-steep rock slide, under different types of seismic waves. Part of the analysis of the numerical results have been discussed in previous section, and they will be described as conclusions in next section as well. In addition, there are still some discussions worth of attention.
(1) Selection of the input seismic waves: At time of the Jiuzhaigou earthquake, China Earthquake Networks Center collected the seismic waves of 66 stations. The seismic waves selection criteria in this study are as follows: the station as close to the study area as possible. However, the recorded seismic waves are quite different in the total 66 stations, e.g., the amplitude. Therefore, it is worth discussing whether it is available to select seismic waves only by the value of distance.
(2) The validity of the typical high-steep rock slide numerical modeling: This study is to explore the regularity of the dynamic response of the typical highsteep rock slide under different seismic waves. Therefore, the numerical model is simplified based on the actual landslide. However, the conclusions investigated in this study are consistent with those obtained in other studies. For example, seismic waves have edge

Discussions
This study evaluates the dynamic response of a typical type of landslides, that is, high-steep rock slide, under different types of seismic waves. Part of the analysis of the numerical results have been discussed in previous section, and they will be described as conclusions in next section as well. In addition, there are still some discussions worth of attention.
(1) Selection of the input seismic waves: At time of the Jiuzhaigou earthquake, China Earthquake Networks Center collected the seismic waves of 66 stations. The seismic waves selection criteria in this study are as follows: the station as close to the study area as possible. However, the recorded seismic waves are quite different in the total 66 stations, e.g., the amplitude. Therefore, it is worth discussing whether it is available to select seismic waves only by the value of distance.
(2) The validity of the typical high-steep rock slide numerical modeling: This study is to explore the regularity of the dynamic response of the typical high-steep rock slide under different seismic waves. Therefore, the numerical model is simplified based on the actual landslide. However, the conclusions investigated in this study are consistent with those obtained in other studies. For example, seismic waves have edge and altitude amplification effects [36], and the Up-Down seismic waves have little impact on the dynamic response of the landslide [29]. The available of the numerical results can be verified. The follow-up work will study a specific landslide to explore how to select seismic wave loading in 3DEC numerical modeling.
(3) Effects of the orientation of the slope: Firstly, it should be noted that the free-field boundary in 3DEC requires that the all five surfaces (except the top surface) of the numerical model should be parallel to the coordinate axis. Moreover, the seismic wave loading is limited to the coordinate axis; that is, it can only be loaded along X, Y, and Z directions. Owing to the left-handed coordinate system adopted in the numerical modeling, the Up-Down seismic wave is determined along the Y direction. The horizontal seismic waves (East-West and North-South seismic waves) are corresponding to the X or Z directions. However, the actual orientation of the slope is not exactly towards north, south, east, or west. It may lead to the inconsistency of the seismic wave loading in numerical modeling with the actual vibration of the slope. For example, this study assumed that the orientation of the slope is towards east, corresponding to the X direction. Therefore, the input X-directional seismic wave loading can be regarded as the East-West seismic wave. Once the actual slope is not towards east, the X direction is not actually the East-West seismic wave.
(4) The reason for the numerical simulation time longer than that of the actual situation: In general, the failure process of the rock slide occurs very fast and lasts for a few seconds. However, the convergence time of the numerical simulation is much longer than that of actual rock slide. The reason for the numerical simulation time longer than that of the actual situation might be the 3DEC calculation principle. The calculating unit of 3DEC is called "block". The overall system can be influenced by each single "block". Therefore, it may take a relatively longer time to ensure each "block" is in a stable state. However, in the general failure process of the landslide, the time of the failure process of the rock slide can be estimated roughly by observing the rock slide stops.

Conclusions
The effects of the four types of the seismic waves on the dynamic response and failure behavior of landslides were numerically investigated. The following conclusions were drawn: (1) The velocity contour shows that the velocity increases with the growing altitude, and the surface velocity is more significant than the internal velocity at the same altitude. In other words, the dynamic response of the slope is affected by vertical and surface amplification. By contrast to the velocity field of East-West, Up-Down, and three-directional seismic waves, it can be investigated that each seismic wave actually affects the landslide independently and will not interfere with each other.
(2) By comparing the failure behavior of four cases, the East-West and three-directional seismic wave have more vital damages to slope stability than North-South and Up-Down seismic wave. The East-West wave plays a dominant role in the failure behavior of landslides in the three-directional seismic wave. Still, the existence of North-South and Up-Down seismic wave will accelerate the failure process of landslides.
(3) The displacement/velocity of the slope's edge is more significant than that of its inside, and the height of the landslide is positively correlated with the displacement/velocity. As for the distinctive effects of different waves, the effects of seismic waves on the displacement of landslides are different at low altitude, while the effects of seismic waves are equivalent at high altitude.
(4) The farthest distance (X direction) of the rock slide is in East-West seismic wave case, while in other cases, they are very close. The reason the earthquake can trigger longdistance movement of landslide is that the slope can obtain energy from the continuous shaking of the ground surface and convert this energy into kinetic energy, which is reflected in the block movement. In addition, the displacement of the landslide along the dynamic loading direction is the largest, while the vertical loading direction is the smallest.
(5) The local damping ratio and the joints dip angle can affect the failure behavior of the rock slide. The stability of the perilous rock mass increases with the growing value of local damping ratio. The failure behavior of the rock slide changes from toppling movement to sliding movement with the increasing joint dip angle.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.