Investigating the Role of Earthquakes on the Stability of Dangerous Rock Masses and Rockfall Dynamics

In recent years, earthquake rockfalls have occurred frequently all over the world, resulting in heavy casualties and property losses. Unfortunately, research on rockfall dynamics upon earthquake events is still rare due to limited field and experimental data, and restricted to numerical simulations of only two dimensions. In order to primarily reveal the role of earthquakes on rockfall, this study focused on the stability of rock blocks on an inclined slope, and following rockfall dynamics by three-dimensional discontinuous deformation analysis (3-D DDA). First, earthquake input methods were discussed and implemented for the triangulated regular network (TRN) in 3-D DDA. The effectiveness was verified by comparison with analytical solution results of a single block on the inclined slope under seismic loads. Further, by discussing the variations of the boundary chart of failure modes, it indicated that the block was more prone to slide even with a large friction angle, became instability under seismic conditions. Moreover, a regular dodecahedron rock block was released on a stochastic roughness slope with two platforms through parallel realizations. The indices of the movement characteristics of the block, such as runout distance, lateral displacement range, and resting position, were investigated. The results showed that the maximum runout distance was not sensitive to seismic load, but the lateral displacement range was significantly sensitive to seismic load and increased appreciably. Through the 3D-DDA numerical simulations, both rock stability and rockfall behaviors under earthquake conditions could be better understood. Furthermore, it will be helpful to by analyzing trajectories and kinetic energies predict earthquake rockfall disasters and design reasonable protective countermeasures under earthquake scenarios.


INTRODUCTION
Dangerous rock refers to a rock mass cut into by multiply sets of joints, generally in an unstable state or a state of limited equilibrium under the action of gravity, an earthquake, rainfall, or other external forces (Chen et al., 2003). It is usually located on a high and steep slope or cliff, suddenly loses its stability, and forms a rockfall or even a rock avalanche. After one or more combinations of falling, rebound, jumping, rolling, or sliding, the falling rock masses move rapidly downward, finally coming to rest in relatively flat areas or being stopped by an obstacle. Dangerous rocks are difficult to identify, although the 3S technologies and AI algorithms that have been developed aid substantially (Huang et al., 2017a;Huang et al., 2021a;Huang et al., 2021b) due to the small volume and uncertain occurrences (Dorren, 2003). In addition, Due to the high-altitude and extremely energetic characteristics, rockfalls usually have destructive characteristics (Volkwein et al., 2011), which can pose a serious threat to roads, railways, buildings, and other structures by means of impact and burial (He et al., 2014;Yuan and Pei, 2014;Hou et al., 2010).
Among the external forces triggering instability, earthquake is one of the major causes of rockfall. In recent years, rockfall and collapse disasters caused by earthquakes have occurred frequently all over the world, resulting in heavy casualties and property losses. There were co-seismic and post-seismic disasters in the recent 2008 Ms 8.0 Wenchuan earthquake (Yin et al., 2009;Li and Kong, 2011;Su et al., 2012;Cheng and Su, 2014), 2013 Ms 7.0 Lushan earthquake (Chen et al., 2013a;Pei and Huang, 2013;Cheng and Zheng, 2014;, 2014 Ms 6.5 Ludian earthquake, and 2017 Ms 7.0 Jiuzhaigou earthquake (Cheng et al., 2018). The effect of an earthquake on rockfall is manifested as showing large scale, huge quantity, and tremendous destructiveness, occurring even with small slope angles under small magnitudes (Qiu et al., 2009). For example, there were 1884 new potential rockfalls after the Wenchuan earthquake, more than triple those that were revealed before 2008. Rockfalls comprise an even larger proportion of the secondary geohazards induced by the Lushan earthquake . Meanwhile, statistical data shows that there were 15 typical earthquakes (magnitude greater than 5.0) that occurred in Southwest China (Wang, 2019) within the last 10 years. More than 100 earthquakes with a magnitude greater than 3.0 occurred in Sichuan Province within 2020. Therefore, it is undoubtedly of great theoretical significance and engineering application value to study the failure mechanism of dangerous rock masses and the following rockfall dynamics under earthquakes. Marzorati et al. (2002) investigated the earthquake-induced collapsed boulders in Umbria and Marche, Italy in 1997, and studied the relationship between the frequency of rockfalls and the environment and seismic parameters by using mathematical statistics. Li et al. (Pei et al., 2011) investigated the extremely longrunout rock avalanches and rockfalls induced by the Wenchuan earthquake, and pointed out that under the action of seismic force, the collapsed rock mass was thrown out at a certain initial speed, which is different from the ones induced by gravity only. The former has high initial speed, large kinetic energy, and a wide influence range of rockfalls. Unfortunately, research on rockfall dynamics upon earthquake events was still rare involved due to the long recurrence period and limited field data on rockfalls during earthquakes.
Laboratory experiments were also rare or missing in the literature. Huang et al. (2019) obtained the rockfall data under earthquake by the shaking table test and proposed a prediction model based on an improved KNN algorithm.
Stability analysis and trajectory modeling using the discrete numerical methods (the distinct element methods, UDEC, PFC, and et al.; Discontinuous deformation analysis, DDA) have the advantage of representing the physical details of dangerous rock masses and following rockfalls, especially under the seismic scenario.
Based on UDEC, Teng et al. (2013) investigated the influence of seismic intensity on the stability of a rock slope with two groups of joints. An actual seismic recording was implemented. Gao (Gao, 2015) used UDEC to reproduce the collapse process of dangerous rock mass under different terrain, geology, and seismic input conditions. The results were summarized to describe the influence of seismic input on the collapse mode and time. Xue (Xue, 2016) used UDEC to establish a simplified model of twodimensional expressway high rock slope, discussed the failure mechanism of high rock slope under seismic load. The rockfall dynamics after the collapse were recorded and analyzed. In addition, the impact force to the retaining wall was calculated according to the law of conservation of momentum.
Hatzor and Feintuch (Hatzor and Feintuch, 2001) applied different forms of sine wave horizontal acceleration to an inclined slope in 2001 and studied the sliding characteristics of blocks on the slope, and proved that DDA can simulate earthquake rock collapses and rockfalls; in 2004 (Hatzor et al., 2004), a twodimensional DDA method was proposed to add time-varying acceleration as a physical force to each block. This method can accurately predict the failure mode of key blocks of jointed rock slopes. Kamai and Hatzor (Kamai and Hatzor, 2008) established a two-dimensional slider model, inputting the displacement time history to the bedrock, studied the dynamic characteristics of the bedrock underlying earthquakes, and obtained the numerical results which are in good agreement with the theoretical solution. Zhang (Zhang, 2011) focused on the collapses developed along a highway during the Wenchuan earthquake as a research object, used the 2-D DDA to study the failure mechanism of rockfalls and rock avalanches under strong earthquakes. Huang et al. (Huang et al., 2016) focused on that the friction coefficient should be gradually reduced in the numerical simulation of seismic rockfall. Based on this, an improved two-dimensional DDA method was proposed. Chen and Wei (Chen and Wei, 2017) analyzed the phenomenon of collapse and rockfalls caused by the failure of dangerous rock mass due to earthquake by 2-D DDA, studied the rockfall dynamics and put forward protection suggestions. Huang et al. (Huang et al., 2017b) studied the influence of seismic load on the rockfall dynamics after the collapse of dangerous rock mass using the 2-D DDA method. The results implied that the influence of vertical seismic load on collapsing blocks is even greater than that of horizontal seismic load. Zhou et al. (2021) simulated the rockfall dynamics of an actual potential rockfall site, and evidenced the special thrown out phenomenon under seismic conditions.
The previous investigation of earthquake rockfall was mostly limited to 2-D numerical simulation. The important lateral displacement could not be obtained. In addition, the stochastic random nature of rockfall was also ignored. Beyabanaki et al. (2009) compared the block displacement solution obtained by the 3-D DDA method with the analytical solution under dynamic load to verify the effectiveness of the three-dimensional DDA method. Chen et al. (2013b) compared the differences in the prediction of rockfall dynamics between 2-D and 3-D DDA simulations and proposed a TRN (triangulated regular network) in 3-D DDA, which could be feasibly used to demonstrate the random roughness of slope surface by disturbing the elevation of each gridpoint. Based on the above, this paper Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 824889 implemented a virtual depth to TRN and conducted acceleration wave on its cells to develop a new rockfall simulator based on 3-D DDA. The algorithms were first validated based on comparison to the theoretical solution of block stability on an inclined slope. The variation of failure mode was also used to demonstrate the influence of seismic loading on the stability of the dangerous rock masses. The role of an earthquake on rockfall dynamics was then discussed especially in terms of runout distance and lateral displacement.

Block Displacements and Deformations
For a 3-D DDA formulation, the linear system has twelve degrees of freedom for each individual block with an arbitrary polyhedral shape: where (u 0 , v 0 , w 0 ) are the rigid body translations of a specific point (x 0 , y 0 , z 0 ) of Block i, (r x , r y , r z ) are the rigid body rotations of Block i with the center of rotation at (x 0 , y 0 , z 0 ), and (ε x , ε y , ε z ), and (γ yz , γ zx , γ xy ) are the normal and shear strains on Block i, respectively. The displacement of an arbitrary point (x, y, z) in Block i is: where [T i (x, y, z)]is the formula for the displacement function of Block i, which is given by: (3)

Equations of Motion
The total potential energy is the summation of all of the potential energy sources from the block stiffness, the initial stress, the point loads, the body loads, the inertia forces, the constraint springs of fixed points or measured displacements, and the contact forces between the blocks. The equations of motion for a system of n blocks is derived by minimizing the total potential energy: where K ij and F i are 12 × 12 stiffness submatrices and 12 × 1 loading vectors, respectively. The submatrices K ii depend on the material properties of Block i, and K ij (i ≠ j) is defined by the contacts between Blocks i and j. For additional details, the interested reader is referred to Shi (Shi, 1988;Shi, 2001).

Triangulated Regular Network for the Terrain of a Slope
Kinematic analyses of rockslides only consider the kinematic behavior of the displaced rocks, and the slope merely serves as the boundary for these rocks to interact with. Moreover, the ground motion and deformation that are induced by the rockslide are neglected due to their small values; for example, the velocities induced by a 2.5 × 10 6 m 3 rockslide are on the order of 10 −6 m/s (Favreau et al., 2010). Therefore, it is not necessary to generate a slope as blocks. A simpler slope model for the terrain of a slope, such as a triangulated regular network (TRN), can play a similar role as the boundary of the blocky slope.
The terrain of the slope in the raster format is represented mathematically by a set of elevation points in the form of an m × n matrix (where m and n are the rows and columns of the matrix, respectively) of elevations. The faces in 3-D DDA are defined as plane polygons. Therefore, it is necessary to divide each cell into two triangles (Figure 1). For the special contact treatment and the corresponding validations, the interested reader is referred to Wang et al. (2017).

Basic Assumptions of Seismic Input and DDA Implementation
The DDA can consider seismic loading in three forms: acceleration, velocity, and displacement time histories. Loading of time-dependent acceleration is widely used by many researchers to conduct earthquake analysis (Hatzor and Feintuch, 2001;Sasaki et al., 2004;Favreau et al., 2010;Wu;Ning and Zhao, 2013). In practice, acceleration loading is usually used to the failure rockfall (method 1) or the slope (method 2) as volume forces.
When the rockfall mass travels over the slope, method 1 could generate the wrong displacement (Wu). In 2-D, method 1 is more suitable for the pseudo-static analysis for the stability problem of dangerous rock; method 2 has been validated to be more adaptable for rock avalanche and rockfall problems Zhou et al., 2021).
However, the TRN only performs the rigid boundary without density. It is difficult to generate volume forces by acceleration loadings. Hence, a virtual depth of 1 m was set for each TRN cell. In order to simulate the main features of 3-D rockfall, This paper quotes some basic assumptions proposed by Zhang et al. (2013): (1) Blocks are treated as 3D polyhedrons; (2) Rockfall masses could be divided into smaller blocks by preexisting joints; (3) Geometry of the slope model is presented by TRN with a virtual depth of 1 m; (4) Earthquake forces only act on the TRN cells.
(5) Density of TRN cells is assumed as 105 times to avoid the influences of gravity and the impact of falling blocks. (6) Gravity acceleration acting to TRN is set to zero to make updown free vibrations.
A TRN slope with virtual depth was then constructed based on the raster data shown in Figure 1A. The real seismic acceleration record ( Figure 2A) from MZQP stations during the Wenchuan earthquake was horizontally-vertically loaded to the TRN virtual blocks. Relative displacement time histories are also shown in Figure 2B. The 3-D DDA results are notably consistent with the theoretical integration results.

State Boundary
A block on an incline (Figure 3) has four different possible states: ① static stability, ② downslope sliding, ③ toppling, and ④ toppling and sliding simultaneously. Its state is controlled by the geometry of both the block and the inclined plane, and the frictional resistance of the interface between them, the three of which are defined by three angles as follows ( Figure 3): δ, the block aspect angle defined by the ratio of the block width and height; α, the inclination angle of the slope; and ϕ, the friction angle of the interface between the slope and the block.
Define the pseudo-static force F k 1 mg acting on the block center, and k 1 tan β. The state chart of the block using limit equilibrium analysis is given as follows: ① the critical angle between stability and sliding: ② the one between stability and toppling: ③ the one between sliding and sliding + toppling: ϕ δ ④ the one between toppling and sliding + toppling: tan α + β 3 cos 2 δ tan ϕ − 3 sin δ cos δ + tan ϕ 3 sin 2 δ − 3 sin δ cos δ tan ϕ + 1 , if δ defined by the shape of the block is given, the four boundaries can be obtained (Figure 4). The detailed derivation could be found in Yang (2018) (Yang, 2018).

Verification of the Boundaries With DDA
To directly compare with the theoretical solution results, method 1 was adopted here to validate the dynamic accuracy of 3-D DDA with TRN.
The numerical control parameters used in the DDA verification are provided in Table 1.
Then 3d-DDA stability analysis of a single block under the same horizontal seismic load was conducted. By taking ϕ 30°as an example, the state of the block when tan β 0.2 are listed in Table 2.

Influence of Seismic Load
Compared with the condition of only gravity, the total boundary figure remains similar, but there is a translation determined by both the size (β) and direction of the pseudo-static earthquake force as shown in Figure 5. In the figure, the black lines BA, BO, BG, and BH are the boundaries with gravity only, and the red lines CD, CE, CF and CH are the boundaries with horizontal pseudo-static force. When the boundary moves a unit to the left, the stable region decreases significantly, and the sliding region increases by the area of ABCD, indicating that the original stable state changed into sliding after adding seismic load, that is, the block was more prone to slide. The toppling region decreases by the area of BGFI, indicating that the original toppling state was transformed into a toppling + sliding state under seismic conditions. The original stable state was transformed into toppling after adding seismic load as shown in the CIOR region. The calculated BGFI area is greater than CIOE, indicating that the proportion of toppling was reduced. The area where sliding and toppling occur simultaneously increases the area of the CBGF part, indicating that the application of horizontal seismic force had a great influence, and the stability was significantly decreased.
Assuming that the seismic acceleration was 0.2 g, k 1 tan β 0.2 (β 11.31°), ϕ varied from 20°to 50°. The instability mode maps with only gravity and seismic load were drawn in Figure 6. Where the green shadow depicts the area of sliding increase ΔS S+ , the orange shadow describes the area of toppling decrease ΔS T− , and the blue shadow demonstrates the area of toppling increase ΔS T+ .
The change of the sliding part is listed in Table 3. It implies the area increased after the seismic load was applied, while the rate shows an increasing trend with the increase of friction angle.
The variation of the toppling area is demonstrated in Table 4. It indicates ΔS T− was always greater than ΔS T+ after applying seismic load. The difference was a constant, only related to acceleration value, manifested by a small rate of change.
The results indicate that under the action of seismic load, the block on the inclined slope was more prone to slide, losing its stability.

Model Establishment
The TRN grid size was 2 m. The height of the slope was 30 m with a slope ratio of 1:1. In order to quickly dissipate the kinetic energy Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 824889 of the falling rock under gravity, two 10 m wide platforms were set at 6-16 m and 22-32 m, respectively. After several trial experiments, the width of the slope was determined as 40 m while the total length of the model was 60 m (Figure 7). To represent the apparent randomness during the impactrebound, a stochastic roughness was applied to the slope surface. The roughness parameter R was defined based on the concept proposed by Romkens and Wang (Romkens and Wang, 1986). R was 1/4 to generate significant stochastic rockfall dynamics. The implement can be found in Chen (2018) (Chen, 2018). It was achieved by applying gridpoint elevation variations which obeyed normal distributions with expectation μ 0, standard deviation σ 200 mm.
The simulations were performed using a regular dodecahedron rock block (volume of 1 m 3 ) initially located 2.5 m (center) above the top of the slope. The numerical procedure began with releasing the rock to the slope, generating an initial velocity. The rock then bounced, rolled over the slope. It ended when the rock finally stopped.

Parameter Setting
The TRN cell with virtual depth was rigid here with only friction angles. According to the field experiment conducted by the Japan Road Association in 2000, a friction coefficient of 0.11-0.20 was adopted, which is equivalent to an average dynamic friction angle of 10°, and its standard deviation was set at about 0.05.
To achieve energy dissipation during impact, the postmodification proposed by Chen (Chen et al., 2013b) was adopted in this paper. The rebound velocity (V x , V y , V z ) (the velocity when contact between rock block and surface TRN ends) was DDA simulated and then modified to be not greater than R v times of incident velocity (V′ x , V′ y , V′ z ) (the velocity recorded when block-TRN contact occurs).
R v value was assigned due to the surface material of each TRN cell.
Considering that R v is obviously dependent on both the modulus and incident angle of v in the normal direction, the R v (scaled) proposed by Tetsuya et al. (2009) was utilized. The R v (scaled) was 0.95, the standard deviation was 0.05, and K was 12 m/s.
SF 1 The volume of falling rock was 1 m 3 . The density was 2.0 × 10 3 kg/m 3 , the Young's modulus was 1.0 × 10 4 MPa and the Poisson's ratio was 0.2.
The values of physical parameters of falling rock and control parameters for DDA simulations are listed in Table 5.

The Number of Simulations
Considering the uncertainty raised by the nature of such a stochastic model, a limited number of observations will amplify the random fluctuations in rockfall dynamics. To avoid this, the model must to be performed with a number of parallel realizations. The procedure under gravity conditions was repeated until the relative error of each observation variable (runout distance, velocity, jumping height, lateral displacement, and platform resistance no.)  Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 824889 8 was less than 1%. As the no. of realization increased, the relative error of lateral displacement decreased greatly, and the one of each other variable also decreased ( Figure 8).
The simulations were performed; each block was released 250 times and all realizations were drawn to extract for the distribution of each variable. The modeled maximum value was used to estimate the uncertainty representing the variability in the stochastic results. Although there may still be a potential undersampling, it found that the statistics data (maximum value and distribution) from 250 parallel realizations could present the stochastic rockfall dynamics.

Earthquake Input
The control parameters of the earthquake are acceleration and frequency using a sine wave. The orthogonal test simulation was carried out by control variables. The duration was 25 s. Condition 1 was the reference group without seismic load. The calculation conditions are shown in Table 6. It can be seen from Figure 9 and Figure 12 that by comparing with the case without seismic load, the maximum runout distance was not sensitive to seismic load, but the lateral displacement range increased largely. When the seismic frequency was 2 Hz, as acceleration increased, the maximum runout distance changed little, even showing a decrease by 1.59%; but the lateral displacement range increased significantly, by 38.92%. When  Frontiers in Earth Science | www.frontiersin.org January 2022 | Volume 9 | Article 824889 acceleration was 0.5 g, with the increase of seismic frequency, the maximum runout distance and lateral displacement range decreased by 4.42 and 17.70%, respectively. The results indicate the strong earthquake not only provided large initial velocities to falling rocks but also increased their kinetic energy during the travelling stage.

CONCLUSION
In this paper, we first reviewed the basic theory of 3-D DDA, including block displacements, deformations, and equations of motion. Then we discussed the TRN in 3-D DDA simulation for stochastic rockfall analysis and proposed a realization of seismic load input in terms of acceleration for it by introducing a virtual depth. On this basis, the failure mode chart boundary of a single block under seismic load was first analyzed by 3-D DDA and compared to an analytical solution to verify the accuracy. The results from 3-DDA with TRN showed a good agreement with theoretical ones. The influence of seismic loading on the stability of the block was further discussed. It implies that the block was more prone to sliding even if the friction angle was significantly large that it decreased its stability under seismic loading.
Finally, we used the 3-D DDA with scholastic TRN to simulate and analyze the rockfall dynamics under seismic waves. Compared with the case without seismic load, the maximum runout distance was not sensitive to seismic load, but the lateral displacement range was significantly sensitive to seismic load and increased appreciably. 2) Strong earthquakes not only provided a large initial velocity for rockfall but also increased the kinetic energy of rockfall in the traveling stage. Through the 3D-DDA numerical simulations, rockfall phenomena and behaviors could be better understood. It is helpful for the prediction of rockfall disasters and the reasonable design of protective countermeasures.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusion of this article will be made available by the authors, without undue reservation.