{'fraction_non_alphanumeric': 0.033185840707964605, 'fraction_numerical': 0.01991150442477876, 'mean_word_length': 4.7786931818181815, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Characteristic analysis of the Nayong rock avalanche based on the seismic signal

A rock avalanche that destroyed 23 houses and killed 35 people occurred on August 28, 2017 in Nayong, southwest China. A discrete element model, MatDEM, was combined with the dynamic parameters of seismic signal inversion to determine the kinematic behavior of rock avalanches. The best fitting parameters were identified by comparing the velocity evolution process of numerical simulation with that of seismic signal inversion. The dynamic process obtained by modeling was compared with the frequency distribution spectrum of the nearest seismometer, with results revealing that the dynamic process agrees with the parameters inverted from the seismic signals. The simulation results further indicate that the movement process lasted for nearly 40 s, with a maximum speed of 40 m/s. The selected models and parameters can help provide a more accurate explanation of the dynamic processes of similar rock avalanches and are significant to hazard prediction in karst areas.


Introduction
Rock avalanches are gravitational instabilities that are highly destructive due to their bursting nature, long runout, high mobility, and entrainment capacity. Apart from conducting field observations and experiments, many studies have proposed dynamic models and numerical methods for predicting the post-failure behavior and movement process of long runout rock avalanches [1][2] [3] . However, obtaining the accurate parameters needed for these numerical models is quite difficult. Owing to the size effect, the experimental results of rock samples may be inconsistent with the actual material mechanical parameters during large-scale landslide movement, thus resulting in inaccurate simulation results. The validities of these models and approaches have been verified through comparisons with real events. Nevertheless, few models are considered effective due to the lack of direct time-dependent observation evidence regarding the movement processes of rock avalanches.
The seismic signals recorded by surrounding seismometers provide a possible approach to analyzing the movement processes and energy variations of rock avalanches. The seismic signals generated by landslides contain two distinct components reflecting different dynamic processes. On the one hand, the low-frequency component is produced by the loading and unloading cycle of solid soil under the action of body acceleration, and this can be used to invert the force history, the trajectory, and the mass of very large landslides [4][5] [6] . On the other hand, the high-frequency component is generated by the friction and collision between rocks and rock bed; it is used mostly for smaller events detected on regional networks [7][8] [9] . Some studies have recently focused on the application of seismic wave inversion to the landslide dynamics in numerical simulations [10] [11] . Using the dynamic parameters obtained by seismic signal inversion can help constrain the numerical simulation.
The current work analyzes seismic signals and provides more constraints for the dynamic process of the Nayong avalanche. The time-frequency distribution spectrum, which reflects the movement stage of the rock avalanches, can be obtained by Hilbert-Huang transform, and detailed dynamic characteristics can be determined through the force-time function. Afterwards, the runout behavior of the Nayong avalanche is simulated by the 3D discrete element model, MatDEM, and the simulated results are verified by comparing it with seismic and field investigation data. We expect that the method of combining seismic signals and the simulation model can help us gain a better understanding of the possible mechanisms of rock avalanches and conduct predictive simulations of similar events in bare karst areas.

The Nayong Rock Avalanche
The Nayong rock avalanche occurred in Southwest China, which comprises approximately 62 × 10 4 km 2 of bare karst topography [12] . In the past few years, several catastrophic rock avalanches, especially those that are triggered by mining activities and rainfall, occurred in the karst areas of China, bringing massive damages to the surrounding infrastructure and facilities [13] [14] . Unlike other rock avalanches, the Nayong rock avalanche showed a clear sign of rock collapse, and the local residents even recorded the whole process from rock collapse to debris flow with their smartphones [15] . Nevertheless, due to its size and runout distance far beyond the prediction, 35 people died and 23 houses were destroyed ( Figure 1). 1-1' crosssection line. [9] The study region has an elevation ranging from 2175 m a.s.l. at the mountain summit to 1800 m a.s.l at the Pusa Village, with a slope angle ranging from 10° to 55° (Figure 2). Three normal faults (F1 to F3, respectively, in Figure 2) are presented in this region. F1 and F2 evolve at the source area of the landslide; F1 and F2 have dip angles ranging from 63°-70° and 70°-75°, respectively. F3 is found 400-500 m northeast of the study area with a dip angle range of 75°-80°. Moreover, three major joints are identified at 345°/85°(J1), 84°/85°(J2), and 300°/85°(J3), respectively, with a bedding plane of 280°/8°, thus resulting in the crushing of the rock block.    Figure 3. Longitudinal profile of the Nayong rock avalanche [16] .

Seismic data
Six three-component broadband seismometers around Nayong County recorded the seismic signals generated from the rock avalanche ( Figure 4). The information details about the surrounding seismometers are shown in Table 1. According to the records of the closest station, ZJW station, with a distance of 5.59 km, the seismic signal started at 10:21:48 local time. Then, it increased to peak motion at 10:22:23, within a strong signal lasting for 35 s, before gradually fading into background noise after 10:24:00. The low-frequency component of the seismic signal is generated by the loading and unloading cycle of the solid Earth via the bulk acceleration and deceleration of the landslide mass; meanwhile, the high-frequency component is generated by the friction and collision between the rocks and the rock bed [4] . In this study, we used the Hilbert-Huang transform to obtain the time-frequency distribution spectrum of the seismic signal recorded by ZJW station ( Figure 5). The time-frequency distribution spectrum shows that the low-frequency signal was looming at 10:22:22, and amplitude was becoming more intense during 10:22:28.3-10:22:44.7. From 10:22:44.7, the amplitude gradually decreased until it finally disappeared in the background noise after 10:22:61.2. The high-frequency component (which was greater than 0.5 Hz) started to appear at 25 s, and the amplitude soon reached the peak. The highfrequency signal appeared suddenly at about 25s, and the large amplitude was maintained at 22-44 s, then started to decrease at 44.7 s. It further decreased after 61.2 s, and finally disappeared in the background at 80 s.

Method
The earthquake produced by landslides or rock avalanches can be represented by a single-force mechanism [17] . Therefore, the solid Earth can be considered as a slope, and the sliding mass can be a constant mass body sliding above the slope [4] (Figure 6).
For the sliding body, the force is exerted by the Earth's crust, and another force ( ) acts on the Earth's crust. Both are equal in magnitude and opposite in direction. Then, the ( ) is equal to where ( ) is equivalent to the single-force source for generating the long-period seismic signals of the rock avalanche. Thus, ( ) is called force-time function this study. ( ) can be determined by inverting the low-frequency seismic data [6] . The avalanche source is assumed to be a stationary point source. The impulse response set between the source and each station pair is known as Green's function. Hence, the seismic displacement records ( ) can be expressed by the convolution of forcetime function and Green's function [18] . Green's function calculation adopts a one-dimensional (1D) velocity model with the frequency range of 0.01-0.1 Hz.
For low-frequency signals, Green's functions ( ) can be calculated by the wave number integration method. Owing to the insensitivity of low-frequency signals to small-scale heterogeneity, the 1D generalized Earth model was used in this study, which is from Crust 1.0. The displacements of slide mass records are equal to = ( 1 cos + 2 sin ) + 3 , (3a) = ( 1 cos + 2 sin ) + 3 , (3b) = ( 1 sin − 2 cos ) , (3c) where is the azimuth angle of stations; 1 , 2 , and 3 are the N, E, and Z direction components of the force-time function, respectively; and are the vertical displacement components of generated by downward vertical and horizontal forces, respectively; and are the radial displacement components generated by the vertical and radial horizontal forces, respectively; and is the lateral displacement component generated by the lateral horizontal force.
To make the models smoother, the force-time function can be inverted by a damped least-squares approach [6] expressed as = ( * * + 2 ) −1 * , where * is the convolution matrix of Green's function, is the identity matrix, and is the damping coefficient. Once the force-time function is obtained, we can calculate the time-varying velocity of the sliding mass v(t). This can be done by integrating the force and displacement ( ) twice as follows: In this case, the results of the field investigation can be used to determine the mass of the sliding body. As shown in Figure 7, we can see that starting from 23.5 s, the speed of the sliding mass increased rapidly, reaching the maximum of 31. 8 m/s at 42.5 s. After that, its speed gradually decreased and almost reached zero at 75.8 s. The maximum movement distance of mass is 774 m.

The Discrete Element Model
The 3D discrete element model (MatDEM [19] ) was applied to simulate the kinematic behavior of sliding mass. This process was undertaken for investigating the detailed mechanisms of the Nayong avalanche. The motion of elements follows the Newtonian equation of motion, and the elements come into contact through breakable, linear elastic spring in normal and tangential directions. As shown in Figure 8, the spring force is equal to the stiffness ( n , s ) multiplied by the corresponding displacement ( n , ) expressed as n = n n , (7a) S = s s . (7b) The element bond will break, and the tensile force will cease to exist when the intergranular displacement exceeds the breaking displacement . However, when two particles return to the compressive contact, the repulsion still acts between them. The maximum shear forces ( S max ) allowed by the Coulomb criterion are expressed below. n max = n b intact bond (8a) S max = S0 − p n intact bond (8c) n = n n n < 0, broken bond (8d) S max = − p n broken bond (8e) S = S max (if S > S max ) (8f) In the equations above, S0 is the shear resistance between elements, and p is the intergranular friction coefficient. An artificial viscosity that can damped the rebound energy of the particles on boundary, is added to the model to avoid excessive accumulation of kinetic energy in a closed simulation system. This is expressed as where is the artificial viscosity, and ′ is the particle velocity. The time-stepping iterative algorithm was developed in the discrete element model to investigate and observed the dynamic evolution of elements [20] . The time step should be much smaller than the natural vibration period of the system so that the model's elastic properties can be accurately simulated. In this study, the time step is set to 0.02 times of the natural vibration period of the elements system.
MatDEM adopts the GPU matrix calculation to support the dynamic simulation of millions of discrete elements. Thus, discrete elements can also be used to construct the entrainable basement layer. In this study, we simulate the gravitational deposition of the discrete elements using the initial model. The discrete elements with a mean radius of 5 m have a certain initial velocity in a rectangular simulation box. These elements collide with one another under gravity, after which they deposit in a random position under artificial compression. The deposited elements are shaped according to the digital elevation model. To save computation, the elements in the lowest layer are defined as wall elements, which do not participate in the dynamic simulation process. The microparameters of discrete element model cannot directly use the macro-properties obtained from the laboratory due to the randomness of element accumulation and the size effect [21] . Liu et al. proposed a formula to convert laboratory data into MatDEM parameters [22] . Other discrete elements are divided into four layers based on the geological properties; in such a classification, the source area is divided separately due to the severe shattering ( Figure 9). Each layer has different mechanical properties (Table 2), which are obtained from the laboratory results converted by Liu's formula.

Simulation results
According to previous studies, the friction coefficient has a tremendous influence on the movement of landslides in the numerical model [23] . Consequently, the current study adopts different intergranular friction coefficients of source area elements (0.4, 0.6, and 0.8) to provide quantitative estimates of the initial conditions. Figure 10 shows  This phenomenon can be attributed to the kinetic energy of elements being dissipated by the collision between particles or between particles and the basement. These collisions cause some elements to spread further; hence, more elements are entrained into the movement. More elements rushed out of the main deposition body in the scenario of low-friction coefficient. Moreover, the displacement of elements at the front edge of the landslide mass is smaller than that at the center of the deposition mass. This indicates that the elements at the front edge of the landslide mass are not all from the source areas, and there are many elements entrained from the basement. These entrainment phenomena are also more intense in low-friction scenarios.

Best scenario
We can hardly distinguish which simulation result is most consistent with the actual event from Figure  10. Fortunately, we can better judge the numerical simulation based on the dynamic parameters from seismic signal inversion [10] . Figure 11 shows the velocity-distance pattern obtained from the numerical simulation of three scenarios and seismic signal inversion. The velocities of the simulation results are the average velocities of elements at the front edge of the sliding mass, whereas the velocity inversed from the seismic signal is the center of the sliding mass. We can see that the velocity of each scenario increases rapidly before the horizontal distance reaches 100 m. Furthermore, except for the data inversed by the seismic signal, the other scenarios reach the maximum speed when the horizontal distance reaches 200 m, i.e., before the 20 m-high hillslope. With the same displacement, velocities of sliding mass in different scenarios are just opposite to the order of the intergranular friction coefficient. Moreover, the friction coefficient only changes the absolute value of the sliding mass' velocity without influencing its changing law. This indicates that the changing law of the Nayong avalanche's velocity is mainly affected by the terrain, and the main obstacle to the sliding movement is the 20 m-high hillslope.   9 We can judge which scenario can best reflect the dynamic process of the seismic signal by calculating the area difference between the three curves and the seismic wave inversion curve via integration. Table 3 presents the results of the numerical model. As can be seen, the intergranular friction coefficient of source area elements equal to 0.6 is the best scenario that is closest to the real dynamic process of the landslide.

Dynamic process
The evolution of the Nayong rock avalanche simulated by the best MatDEM scenario is shown in Figure 12. The simulation starts from the overall rock collapse, namely, 10:22:23.5 on the day of the disaster. The primary movement process lasted for about 40 s, after which the sliding body only expands laterally. In Figure 9, the elements in the source area begin to rush, and the front edge of the landslide mass reaches the toe of the hillslope at 12.   Figure 12 shows the velocity evolution of the sliding mass. As can be seen, between 0 and 12.6s, the speed of the avalanche increases rapidly. Then, the speed of most elements begins to decline after 12.6 s. At 18.9 s, only the front edge of the sliding mass has visible speed. After 31.5 s, the speed of most elements is 0, and a small number of elements in the source area continue to fall. The maximum speed of the element, occurring at 12.6 s, is about 40 m/s. At this time, the landslide mass is passing through the 20 m-high hillslope.
For a comprehensive understanding of the dynamic process, we compare the time-frequency distribution spectra to the simulations that best fit the observations and the UAV video ( Figure 13). Both Fan [15] and Zhu [16] analyzed the failure process of the Nayong avalanche through UAV video, while the current paper studies the rock fall in the source area after the overall collapse. The magnitude of element velocity and the number of moving elements can indicate the intensity of energy  Figure 13, the changes in velocity and sliding volume are in good agreement with the seismic signal. At 6.3 s (10:22:29.8), the materials from the source reached high speed, colliding with one another violently, and affecting the basement layer. At this point, all the high-and lowfrequency signals become stronger. Between 6.3-21.5 s (10:22:29.8-10:22:45), a large amount of mass is entrained into the high-speed movement, resulting in both the high-and low-frequency signals maintaining a strong amplitude in this period. After 21.5 s (10:22:45), most of the materials stopped moving, and only the front end of the sliding body are still spreading forward, thus weakening the low-frequency signal. Meanwhile, the high-frequency signal remained strong due to the constant highspeed falling rocks from the source area colliding with accumulated materials. These rocks affected the basement severely, thus occasionally strengthening the low-frequency signal strong during this period. At 38.5 s (10:23:02), the avalanche movement stopped but the rockfall from the source area continued. Until 56.7s (10:23:20), the rockfall in the source area stopped, and the seismic signal faded into the background. Figure 13. Comparison of the seismogram and time-frequency distribution spectra of the ZJW station with the simulated velocity evolution of the best scenario. Snapshots from UAV video are added.

Conclusion
The interpretation of the seismic signals is used to calibrate the parameters of landslide simulation using the 3D discrete element model, MatDEM. An analysis of the intensity and time history of the seismic signals shows that the seismic signals generated by the Nayong avalanche have obvious longperiod components. The time-varying velocity and distance of the Nayong avalanche can be obtained via integration based on the force history. Combined with the results of seismic signal inversion, the best MatDEM simulation results are selected from the three scenarios. The simulation results indicate that the movement lasted for approximately 40 s at a maximum speed of 40 m/s. Furthermore, the dynamic process agrees with the parameters obtained from the seismic signal inversion. In summary, due to the uncertainty and suddenness of occurrence time, many large landslides have no direct evidence of their dynamic process. Conventionally, the parameter calibration of numerical simulation can only be based on the back analysis of field accumulation. For rock avalanches,  11 although the dynamic process from seismic signal inversion is affected by various factors, such as rockfall from the source area, it is still an effective method to constrain the discrete element modeling. Finally, using seismic signals can improve the modeling reliability and accuracy of the results, which is useful for predicting and assessing further landslide hazards.