Dynamic triggering mechanism of the Pusa mining-induced landslide in Nayong County, Guizhou Province, China

Abstract Catastrophic landslides have occurred frequently in the karst areas of Southwestern China for these years, which makes it urgently necessary to reveal their triggering mechanisms. The universal distinct element code (UDEC) was used to simulate whole response of the Pusa landslide during its underlying coal-seam mining. The results are as follows. First, under the given pre-conditions, when coal seams beneath the original slope were mined from bottom to top, caving zone, fracture zone, and deformation zone overlying the goaf gradually expanded; and as the uppermost layer was mined, the deformation, which was apparently characterized by subsidence and clockwise rotation, overwhelmed the slope. Second, under action of the heavy rainfall, the deformed slope, which was affected by the goaf, behaved further progressive deformation characterized by subsidence and clockwise rotation as a whole. Third, with continuous action of near stoping blasting vibrations, the more-deformed slope, which had been influenced by the goaf and the heavy rainfall, behaved continuous progressive deformation, critical failure, shattering, ejection, debris flow, and accumulation. The deformation exhibited significant subsidence, clockwise rotation, and collapse. These responses were verified by the field reconnaissance. It is the first time to explore dynamic triggering mechanism of the landslide triggered by the vibrations.


Introduction
Disastrous mountain landslides have occurred frequently in the karst areas of Southwestern China for these years. The events have a significant impact on human activities and lead to many deaths and injuries (Yin et al. 2011;Jiao et al. 2013;Li et al. 2015;Li et al. 2016). Studies have shown that in addition to the effects of heavy rainfall and groundwater, earthquakes, infrastructure construction, and the frequent exploitation of mineral resources in the area can also tremendously affects the stability of a slope. In particular, the underground mining of coal and other mineral resources has a particularly significant impact on the formation of the landslides in the areas, which seriously threatens public safety and mining production. For instance, 74 people were killed in the landslide that hit Jiweishan Village, Wulong County, Chongqing Municipality on 5 June 2009; on 28 August 2017, the Pusa landslide in Zhangjiawan Town, Nayong County, Guizhou Province killed 35 people (Yin et al. 2010;Ma et al. 2018;Fan et al. 2019). Under the comprehensive influence of the above factors, the major landslides frequently occurring in this area possess distinct regional characteristics.
Based on detailed field investigations and analyses of typical landslides in this areas, controlled by the geological structure and karstification, the landforms in this areas are dominated by high mountains and canyons, which complicates local geological and ore-forming conditions. The orebody generally exhibits a shallow burial depth and has a large inclination occasionally, so a landslide effect triggered by the underground mining is more significant. Therefore, numerous studies have been carried out to reveal the formation mechanism of typical catastrophic landslides in this areas under the synergistic influence of various factors by using theoretical modeling, in-situ measurements, experiments with centrifuges and shaking tables, finite element, finite difference, and distinct element numerical simulations (Peng and Luo 1989;Shu and Bhattacharyya 1992;Feng et al. 2014;Tang et al. 2015;Deng et al. 2016;Xu et al. 2016;Yang et al. 2016;Zhao et al. 2016;Tao et al. 2018;Zhu et al. 2019). Moreover, work in the UK, Australia, the Czech Republic, Brazil, Turkey, Thailand, and India have also applied similar technical methods to study landslides caused by underground mining in order to reveal their formation mechanism and main controlling factors and to provide technical support for the prevention and control of landslides (Jones et al. 1991;Donnelly et al. 2001;Zahiri et al. 2006;Ergi_Nal et al. 2008;Singh et al. 2008;Fuenkajorn and Archeeploha 2010;Marschalko et al. 2012;Yilmaz and Marschalko 2012;Unlu et al. 2013;Lana 2014;Salmi et al. 2017a;Salmi et al. 2017b;Wardle and McNabb 1992). Finally, 2 main models of mininginduced landslides are summarized, i.e. (1) progressive rock falls and caving failures, which cause nearly vertical cliffs, and (2) rock masses extruded from a slope toe causing a holistic instability Zhao et al. 2016;Wardle and McNabb 1992). Furthermore, mechanisms of mining-induced landslides and ground subsidences under different geological and mining settings were fully illustrated and generalized by detailed case studies (Salmi et al. 2017a;Zhao et al. 2016;Zhao et al. 2021).
Further analysis revealed that the landslides in the areas were not only influenced by some factors such as the micro-geomorphology, the lithological association of the slope, the rock mass structure, the weathering, the goaf in the underlying coal seams, and karstification, but it were also severely affected by other factors such as driving vibrations, stoping blasting vibrations, heavy rainfall, and natural earthquakes. Unfortunately, there are currently few studies on the latter factors. Therefore, it is important to study the dynamic triggering mechanism and the main factors controlling the landslides induced by underground mining in the karst areas.
In this study, by simulating the Pusa landslide in Zhangjiawan Town, Nayong County, Guizhou Province, the distinct element numerical method was used to reveal firstly the dynamic triggering mechanism of this landslide under the actions of various factors such as the goaf, heavy rainfall, and stoping blasting vibrations, especially. The results of this study are important for improving our understanding of the triggering mechanism and identifying the main controlling factors of the landslides in these areas.

Methods
Firstly, a desk study involving a basic geological settings, a preliminary understanding of the landslide, was fully performed. Secondly, a detailed site investigation was conducted, which further made the formation mechanism understandable. Finally, a numerical modeling was employed to identify the failure and the contributing factors. Moreover, the modeling results were verified by that of the site work.

Landslide features and geological conditions
The landslide occurred near the Pusa Village, Zhangjiawan Town, Nayong County, Bijie City, Guizhou Province. This landslide destroyed most of the residential area (Figures 1, 2). In addition, 26 people were killed, 9 were missing (later confirmed to be killed), and 8 were injured. Therefore, this is defined as a catastrophic mountain landslide. In addition, this landslide had a maximum length of about 610 m and an average length of about 575 m along the main sliding direction (about 310 ). It was 360-380 m wide in the middle, and the accumulation body was about 4 m thick on average with a volume of about 82.3 Â 10 4 m 3 . Based on the volume, this was a medium-sized landslide. According to the field investigation and the relevant research results (Fan et al. 2019), the strata (from top to bottom) in the original slope and their main lithologies are as follows: (1) the micrite, marlstone, silty mudstone, and argillaceous siltstone of the Lower Triassic Yelang Formation (T 1 y); (2) the marlstone and argillaceous siltstone, where several thin coal seams are interbedded, of the Upper Permian Changxing-Dalong Formation (P 3 c þ d); and (3) the silty mudstone, carbonaceous mudstone, argillaceous siltstone, and coal seams of the Longtan Formation (P 3 l).
In terms of the geological structure, the landslide was located in the southwestern section of the northwest wing of the Santang syncline in Zhijin County, where the strata form a monocline, and the occurrences of the coal seams are consistent with that of the strata. Affected by faults (fault F1 and F2 intersect the strata of the original slope), the occurrences of the strata vary significantly, with dip directions of 138-187 and dip angles ranging from 7 to 10 . Figure 3 shows geological settings of the area. Figure 4 shows the geological profile of the landslide along the main sliding direction.

Factors influencing the landslide
In this study, the factors influencing this landside are divided into geological factors and triggering factors, which are discussed in the following sections. In fact, the geological factors are preconditions for the landslide, which make the original slope slide probably. They deteriorated stability of the original slope gradually through a longtime influence. However, the triggering factors degraded the slope stability quickly to a critical failure.

Micro-landform
According to the field investigation, as shown in Figure 5, the landslide originated from the steep cliff at the top of the original slope. The 210-m high and steep free surface, which was nearly vertical, provided a favorable micro-landform and precondition for the later rock sliding. The adjacent cliff which is located to the southwest    of the landslide is 260-300 m high, as shown in Figures 2 and 5. Moreover, the adjacent cliff is unstable because of some deep-seated karst fractures which developed inside.

Lithological setting
Generally, the lithological association of the original slope is characterized by a hard upper part and a soft lower part from top to bottom (Figure 4), that is, the micrite and marl ( Figure 6) in the upper and middle parts of the slope are hard, while the part from the foot of the slope to the coal-bearing strata are composed of relatively soft silty mudstone (Figure 7), carbonaceous mudstone, argillaceous siltstone, and coal seams. Such a lithological association is very likely to result in differential weathering ( Figure 8) within the slope, which reduces the stability of the slope. This was also a precondition for the landslide. A detailed lithological distribution of the original slope is illustrated in Figure 4. The rock strata of the Lower Triassic Yelang Formation (T 1 y) are hard and the rock strata of the Upper Permian Longtan Formation (P 3 l) are relatively soft.

Structure and tectonics of the slope
Due to the synergistic influence of physical weathering, karstification, tectonic stress, and mining, structural planes universally developed in the original Pusa slope.
1. Figure 4 shows the rock-stratum interfaces with dips up-slope. Such an anti-dip rock-stratum interfaces possess a better water-holding capacity than down-dip interfaces, and groundwater flowed into the slope along the anti-dip interfaces. During this process, due to the gradual expansion of the pores in the  rock-stratum interfaces caused by karstification, stratification of the limestone in the slope became more significant (Figure 9), which constantly reduced the overall stability of the slope. 2. According to the detailed site geological investigation, nearly vertical joints, with a maximum length of about 30-50 m, universally developed in the rock strata of the original Pusa slope (Figure 9), most of which crossed 2-3 rock strata, with a density of about 2-3 joints/10 m. Under the synergistic influence of physical weathering, karstification, and mining stress (including stoping blasting vibrations), these joints were prone to further development, which made the rock mass structure of the original slope (especially in the shoulder of the slope) evolve from medium-thick strata to fragments, thereby further deteriorating the stability of the original slope. 3. Based on the in-situ geological investigation, deep karst cracks generally developed in the study area ( Figure 10 hydrostatic pressure expansion effect after water filling; and (c.) the tangential drag effect of the seepage flow in a larger range after water filling. 4. The original Pusa slope was crossed by 2 reverse faults at the foot of the slope (fault F 1 and F 2 in Figures 3, 4), and the strata on both sides of the faults were staggered by the high and steep fault plane (dip angle of about 60 ). Consequently, the mining of the coal seams beneath the slope was affected to some extent. Moreover, since the faults exhibited compressed deformation, it is inferred that the compressive stress field of the faults exerted an impact on the adjacent slope to a certain extent, which aggravated the development of the structural planes of the slope and weakened its stability.

Weathering
According to the detailed field investigation and the analysis of the related literature, weathering effect on the structure of the original slope was mainly manifested in the following 2 aspects.
1. Physical weathering (Figures 8, 11): thermal expansion and cold shrinkage of the slope's rock mass and freezing and thawing of fissure flow in the structural plane in the slope caused by a temperature change led to further development (starting from scratch or expansion of scale) of the structural planes in the slope, which reduced the stability of the original slope. 2. Karst chemical weathering: this action refers to how the weak acid flow formed by dissolution of CO 2 in the atmosphere dissolves CaCO 3 in the rock mass on both sides of the structural planes as it passes through the planes. This action further developed the structural planes in the original slope and even formed penetrating controlling structural planes ( Figure 12), thereby greatly reduced the stability of the original slope.

Triggering factors
Based on the on-site investigation and analysis of the relevant literature, main triggering factors of this landslide include goaf, heavy rainfall, and blasting vibrations during the underground mining. data shows that from 1995 to 2012, main coal seams mined in the area were C20 and C16, while from 2013 to 2017, C14 and C10 coal seams were mined were, and C6 was speculated to be partially mined. Figure 4 shows the goaf formed by the mining of the coal seams. The coal seams were mined by using inclined shaft development and long-wall mining method, blasting driving and mining were adopted and the roof was managed by using a collapse method. After occurrence of the disastrous landslide, nearby underground mining was obliged to be stopped because of its obvious influence to the village and the slope ( Figure 14). According to the mining pressure theory (Qian 2003), after underground coal mining, the overlying strata, which are divided into the caving zone, the fracture zone, and the deformation zone from bottom to top, are deformed due to redistribution of the strata stress. If these strata touch the ground surface, mining collapse (also known as mining subsidence) will occur, and ground cracks will form at the edge of the collapse area.
For the original Pusa slope, if this deformation zone develops into the whole or part of the slope, the stability of the slope will deteriorate to some extent, thereby partially deforming (damaging) the slope or even destroying it and causing sliding.

Heavy rainfall
As this area entered the rainy season in June 2017, the accumulated rainfall surged from 147.9 mm on June 21 to 522.2 mm on July 21. This further deteriorated the stability of the original slope, which had been broken by the goaf, karst cracks, and weathering. That is, (1) the rain water entered the cracks, spreading all over the bedding and throughout the joints, which increased the amount of water in the rock mass of the slope, thereby enlarging the volume weight of the slope and increasing the sliding force.
(2) Similar to the hydrodynamic effect of the water flow partially accumulated in the large, deep karst cracks, the water flow partially aggregated in the slope caused the hydrostatic pressure expansion and the tangential drag effect of the seepage in the cracks increasing. This reduced the effective stress and the stability of part of the slope and caused the third deterioration effect, that is, (3) declines of mechanical properties of the rock mass composing the slope and deteriorations of mechanical parameters of the slope cracks, e.g. normal & shear stiffness, tensile strength, cohesion, and internal friction angle.

Stoping blasting vibrations
As mentioned above, during the shaft & roadway driving and coal stoping, blasting was adopted as a feasible method to facilitate the work. However, as shown in Figure  4, the minimum distance between the main coal seam mined (C6) and the foot of the slope was only about 50 m, and that between coal seams C10 and C14 and the foot of the slope was about 80 m. Under such conditions, the vibration waves formed during the mining blasting significantly affected the stability of the slope. This was specifically reflected in the periodic tension-compression and shear effect caused when the vibration waves passed through the slope. This was manifested as horizontal and vertical acceleration, which significantly deteriorated the stability of the partial-broken slope. According to interviews with some survived local residents of the Pusa village during our post-accident field investigation, the village and adjacent slopes always underwent obvious vibrations accompanied with strong underground booms induced by the stoping blastings in recent years.  Figure 15 shows local magnifications of this numerical model, from which it can be seen that four large, deep structural planes developed in the slope shoulder, with depths ranging from 50 m to 80 m. C6, C10, C14, C16, C18, and C20 are the 6 main coal seams beneath the slope.
To reveal the dynamic response characteristics and mechanism of this slope under the synergistic influence of the goaf, the large and deep karst cracks, the heavy rainfall, and the mining blasting vibrations, 3 columns of monitoring points were set up in the shoulder of the slope model ( Figure 16) to record the variations of the displacement, acceleration, and velocity of the slope's shoulder under the integrated action of the above-mentioned influencing factors.

Constitutive model and boundary conditions
Based on the actual failure characteristics of the Pusa landslide, the distinct element numerical simulation was used to reveal the dynamic failure process of the discontinuous and large deformation, i.e. the progressive deformation and dynamic run-out response, of the slope under the synergistic influence of the goaf, the large and deep karst cracks, the heavy rainfall, and the mining vibrations. Thus, for the elements in the numerical model of the slope (meshes in the middle and the upper part of the numerical model shown in Figure 15), the relative displacement between each two elements was much larger than the elastic-plastic deformation of a single element. Therefore, for the part of the slope overlying the goaf, which would experience a large relative displacement, the rock mass was regarded to be composed of discrete rigid blocks (the region around point B in Figure 15) in the numerical model. In other words, a rigid constitutive model was adopted to reveal the relative displacements between the elements while ignoring their own small elastic-plastic deformations. However, the elements constituting the strata beneath all the goaf were characterized by a relative displacement that was significantly smaller than that of the elements making up the slope. Thus, this part was regarded as an elastic-plastic body, and a plastic constitutive model (the region around point A in Figure 15) and the Mohr-Coulomb yield criterion were adopted.
As far as the boundary conditions were concerned, the left and right sides of the numerical model were constrained by a zero horizontal displacement, while the bottom of the model was controlled by a zero horizontal and vertical displacements. In the subsequent simulation of the dynamic response of the slope triggered by the mining blasting vibrations, inner viscous and outer free field boundaries were applied on both lateral sides of the numerical model in order to prevent the vibration waves propagating outward from being reflected back into the interior of the model; that is, necessary energy divergence was allowed. Meanwhile, a viscous boundary was also used at the bottom for inputting dynamic stresses transformed from the velocity records.

Mechanical parameters
To facilitate calculations of the numerical model, the lithology of the slope was generalized, that is, the slope was mainly composed of micrite, marl, argillaceous siltstone, silty mudstone, coal seams, karst fissure filling, and a weathering zone at the slope shoulder (the area near point B in Figure 15, or the potential landslide mass). The structural planes were divided into the bedding planes, each representative lithological joints, and the fault planes. Table 1 and 2 show the physical and mechanical parameters of the typical rock masses in the slope and the mechanical parameters of the structural planes, respectively. The above parameters were assigned by the Jregion and Group keywords built in the distinct element software. Additional, dip angles of the control joints and fault planes were also used for the parameters assigning.
Selection of the parameters was conducted as follows.
(1) Since the elements composing the upper and the middle parts of the model are rigid bodies, all of the parameters in Table 1, except for the rock densities, are for reference only. (2) To avoid re-cracking of the fault planes in the slope during the simulation, the adopted values of mechanical parameters of the planes are relatively large. (3) To avoid excessive overlapping of the numerical elements during the dynamic simulation caused by the blasting vibrations, some type's joints have large normal stiffnesses. (4) Because fillings in the large, deep karst cracks were almost disintegrated before those cracks occurred above the goaf, their mechanical parameters were close to 0, and the deterioration of the parameters under the subsequent action of heavy rainfall was not considered. Due to weak water-stagnant effects of the karst cracks and the inside fillings in the slope shoulder, increment of the density of the fillings during the heavy rainfall was not taken into account either. (5) Since weak water-stagnant capacity of the weathering zone, only density increasing of the rock mass and mechanical parameters deteriorating of the cracks were considered under action of the heavy rainfall.  However, hydrostatic pressure was not considered. (6) Joint mechanical parameters of underlying strata adjacent to rock-bedding planes were used as the parameters of the planes.

Deformation of overlying strata and slope under action of the goaf
According to the mining data, the main underlying coal seams C20 and C16 were mined in this area from 1995 to 2012, while from 2013 to 2017, before the slope slided, the main underlying coal seams C14 and C10 were mined. However, the coal mining in this area was seriously out of the boundaries, and therefore, it is speculated that coal seam C6 had been partially mined before the slope failure. The red rectangular areas in Figure 15 are the goafs formed by the mining in the different coal seams, including C20, C16, C14, C10, and C6. Therefore, in this numerical simulation, the mining sequence in the underlying coal seams was set as C20!C16!C14!C10!C6. Figures 17-19 show the deformation characteristics of the overlying strata and different parts of the slope, which are described below.

Deformation Characteristics of the Overlying Strata
As the coal seams beneath the original Pusa slope were successively mined from bottom to top, i.e. C20!C16!C14!C10!C6, the overlying strata were damaged due to aredistribution of the strata stress. Moreover, as more coal seams were mined, deformation of the slope also gradually increased (Figure 17a-e). When coal seam C6 was mined in the end, obvious caving zone, fracture zone, and deformation zone were observed in the strata overlying the coal seams (Figure 17e and zone A, B, and C in Figure 18), and bed-separation zones were also obviously developed in the above-mentioned 3 zones (Figure 18). When the coal seam C6 was mined, the deformation zone reached the foot of the slope (zone C in Figure 17e), and at this time, the maximum vertical displacement according to records of the monitoring points in the slope shoulder reached À0.883 m (point A in Figure 19, the negative value denotes downward direction of the displacement, which is opposite to the positive direction of the Y axis denoting upward direction of the displacement).  The different degrees of failures in the overlying strata and the slope deformation fully show that the mining of the underground coal seams directly affected the surface slope.

Deformation Characteristics of the Slope
After the mining of coal seam C6, no obvious damage was observed in the original Pusa slope. However, the deformation accumulated to a great extent.
From the vertical displacement (Y-Dis cluster) curves including R, M, and L clusters corresponding to stage V in Figure 19, it can be seen that the maximum vertical displacement of the right side in the slope shoulder (R-cluster curves in Figure 19, monitoring points 12-17 in Figure 16) reached À0.883 m (point A in Figure 19), the maximum vertical displacement of the middle in the slope shoulder (M-cluster curves in Figure 19, monitoring points 5-11 in Figure 16) reached À0.638 m (point B in Figure 19), and the maximum vertical displacement of the left side in the slope shoulder (L-cluster curves in Figure 19, monitoring points 1-4 in Figure 16) reached À0.380 m (point C in Figure 19), which indicated that the vertical displacement of the right side in the slope shoulder was largest. Figure 19 shows the vertical displacement curves including the R, M, and L clusters corresponding to the simulation stages I (C20 mining) !II (C16 mining) !III (C14 mining) !IV (C10 mining) !V (C6 mining), from which it can be seen that as the number of coal seams mined increases, the vertical displacements on the left, middle, and right side of the slope shoulder gradually increases (the red arrow in Figure 17 also shows this variation trend), revealing the cumulative effect of the layer-by-layer mining of the coal seams on the vertical displacement of the slope shoulder. These vertical displacement curves reveal the severe subsidence deformation in the original slope after the mining of the 5 underlying coal seams.
From the horizontal displacement curves composing the X-Dis cluster of the slope shoulder corresponding to stages I!II!III!IV!V in Figure 19, it can be seen that as the number of the coal seams mined increases, the absolute value of the horizontal displacement of each monitoring point in the slope shoulder gradually increases. Specifically, the horizontal displacement at the top of the slope remains positive, which means the horizontal displacement is towards the positive X-axis. Nevertheless, the middle and lower parts of the slope shoulder exhibit negative horizontal displacements along the X-axis (the magenta arrows in Figure 17), which indicates that the overall deformation of the slope shoulder during this process is characterized by apparent clockwise rotation as well as subsidence deformation (Figure 17).
In a word, its macroscopic deformation is characterized by subsidence deformation accompanied by a clockwise rotation.

Deformation of the slope under action of the heavy rainfall
Based on the related meteorological data, after it entered the rainy season in June 2017, the accumulated rainfall surged from 147.9 mm on June 21 to 522.2 mm on July 21, which further deteriorated the condition of the slope with subsidence and rotation deformation that had been influenced by the goaf, karst cracks, and weathering.
It should be noted that since the top and middle parts of the slope were broken by structural planes such as the large, deep karst cracks and the dense joints, the hydrostatic pressure in the slope was insignificant during the heavy rainfall, i.e. the groundwater flow could quickly overflow by seeping along the slope and through the slope cracks. At this time, the heavy rainfall slightly increased the dead weight of the slope (the increase or even saturation of the groundwater contained in the rock mass) and decreased the mechanical parameters of the slope's joints (the joints experienced groundwater seepage). Therefore, instead of applying a fluid-solid coupling calculation in the distinct element numerical simulation, an equivalent method of increasing the volume weight of the slope (potential slide mass or weathering zone in the slope shoulder) and reducing the mechanical parameters of its joints was used to simulate the effect of the heavy rainfall to the slope deformation. Table 1 and 2 show the physical and mechanical parameters of the typical strata and structural planes in the slope, respectively, the simulation results of which are shown in Figure 17f and stage VI in Figure 19. It can be seen that under the influence of the heavy rainfall, the horizontal and vertical displacements in the slope shoulder further increased. However, the variation trend remained as before, that is, the absolute value (negative value) of the vertical displacement of the slope shoulder continued to increase, and the absolute value of the horizontal displacement (positive at the top of the slope, negative in the middle and lower parts of the slope shoulder) continued to increase, continuing the original trend of subsidence and rotation deformation of the slope shoulder.

Dynamic response of the slope triggered by the mining blasting
According to the related mining data, in the area, blasting stoping and blasting driving prevailed. Even worse, the minimum distance between the upmost mineable coal seam and the foot of the slope was only about 50 m. During the mining, the vibration wave produced by the mining blasting imposed periodic tension, compression, and shear influences on the slope. This exerted a more serious dynamic impact to the slope which had been significantly deformed by the mutual action of the goaf, the large and deep karst cracks, and the heavy rainfall. Under such conditions, the subsequent distinct element numerical simulation was carried out based on the above simulation results. Figure 20 shows the time-history records of site-measured horizontal and vertical velocities of the blasting vibration waves. A single blasting lasted for 3 s, and the blasting was conducted every 6 s. The peak horizontal velocities were 0.1651 cm/s and À0.1254 cm/s, and the peak vertical velocities were 0.1501 cm/s and À0.1140 cm/s. To determine response characteristics and dynamic triggering mechanism of the slope (especially the weathering zone in the shoulder, i.e. the potential sliding mass), 10 blasting events with a total duration of 60 s were input into the numerical model after the action of the goaf and the heavy rainfall.
It should be noted that when the dynamic response analysis of the slope during the subsequent blasting was conducted in the distinct element numerical model, the mechanical parameters of the rock mass were not updated (e.g. the dynamic volume and shear modulus should be used), which to a certain extent strengthened the damage effect of the blasting vibrations to the slope (because the dynamic mechanical parameters of the rock mass are often higher than the static mechanical parameters). Nonetheless, this trend was offset by the decreased blasting frequency. An equivalent blasting interval input into this model was 6 s, however the blasting interval in actual mining was much less, often many times simultaneously. Moreover, this numerical simulation focused on determining whether or not the slope would meet its critical failure under the action of the blasting, rather than the process or details, in order to identify main controlling factors, while the actual number of blasting events required to form a landslide needs to be analyzed further. Figure 21 shows the distinct element numerical simulation results under such working conditions. As can be seen from the figure, under the action of repeated stoping blasting, the dynamic response characteristics of the slope developed from macroscopic deformation to macroscopic failure, and the macroscopic collapse was also characterized by significant subsidence and clockwise rotation. In addition, under the continuous actions of horizontal and vertical velocities, the slope also exhibited the macroscopic characteristics of shattering, ejection, debris flow, and accumulation. A detailed dynamic response of the slope is illustrated as follows.
(1) When the slope was subjected to the first cycle of stoping blasting vibrations, i.e. the first 6 s of blasting, during which the vibration waves acted on the slope (Figures 21a-f and points a-b in Figure 20), the slope experienced overall collapse and a further crushing process along the large, deep karst cracks. The process could be refined as follows: when the vibration lasted for 0-4.60 s, the slope was progressively deformed under the forces of horizontal and vertical vibrations produced by the blasting and reached the state of critical failure at 4.60 s. At this moment, the original slope formed a larger sliding surface along the large, deep karst crack (the area enclosed by the yellow broken line in Figure 21e). When the vibration lasted for 6 s, the slope, which had already developed an overall sliding, began to break up further and collapse more severely (Figure 21f). During this process, the displacement of the landslide had the following characteristics: During 0-4.60 s, the slide mass displaced downward along the large, deep karst crack; while from 4.60 s onward, apart from the overall displacement downward, there was an obvious horizontal displacement at the top of the original slope. In other words, the overall displacement of the slide mass was characterized by apparent subsidence, bottom collapse, and clockwise rotation.
(2) When the slope experienced the second cycle of stoping blasting vibrations, that was, when the vibration waves acted on the slope from 6 s to 12 s (Figures 21f-l), the slide mass went through further shattering and ejection. The overall displacement of the slope was still characterized by subsidence and clockwise rotation. It should be noted out that during this process, since the slide mass was already separated, it presented obvious dynamic response characteristics of ejection under the action of the horizontal and vertical vibrations produced by the vibration waves (Figures 21g-j), that was, part of the broken masses were thrown up into the air. This was different from the overall dynamic response characteristics of the slope during the first cycle of the stoping blasting vibrations. The first cycle of the vibrations mainly led to progressive deformation and critical failure of the slope, and the slope was still an intact one, and therefore, the slope was not partially thrown up into the air. However, the slope experienced consistent horizontal and vertical vibration intensities throughout the first and second cycles. As the second cycle of the stoping blasting vibrations went on, the slide mass entered response stages of debris flow and accumulation (Figures 21k-l), and the overall displacement was still characterized by subsidence and clockwise rotation.
In summary, after 2 cycles of stoping blasting vibrations with a total duration of 12 s, the Pusa slope underwent progressive deformation, critical failure, shattering, ejection, debris flow, and accumulation process. During this process, the overall displacement of the landslide was characterized by significant subsidence, collapse, and clockwise rotation. The simulation indicates that the stoping blasting is a main controlling factor that triggered the landslide.
Additionally, we consent for identifying information/images in the paper in an online open-access publication.

Discussion
According to interviews with some survived local residents of the Pusa village during our post-accident field investigation, the village and adjacent slopes always underwent obvious vibrations accompanied with strong underground booms formed by the underground blasting in recent years. The vibrations of the slopes were more obvious during rainy seasons and some minor rock falls occurred continuously. In fact, based on records of the underground mining in underlain coal seams near the initial slope, stoping blasting was continued to be adopted even if its influence to the village and the adjacent slopes was destructive. Even worse, with decreasing depth of the underground mining and decreasing distance from the mining face to the original slopes, the vibrations, which was induced by the stoping blastings, of the slope became more and more obviously. Therefore, the strong underground booms sensed on the ground surface were formed by the stoping blastings, which also induced the vibrations of the village and the adjacent slopes. Furthermore, from 21 June 2017 to 21 July 2017, some times heavy rainfalls occurred, which made the vibrations of the initial slopes more seriously. Five weeks later, on August 28, the Pusa landslide happened eventually. Until its occurrence, the underground mining was obliged to be stopped. As a result, as far as the key triggering factor was concerned, in the light of the field investigation, the vibrations induced by the continuous stoping blastings was inferred to make the landslide eventually. It is in accordance with the result of the numerical modeling.
In addition, the following points should be noted.
1. Since the elements were mainly set as rectangles cut by orthogonal joints in the distinct element numerical model, the characteristics of the long-distance debris flow of the landslide were not fully revealed. We plan to modify the elements at the top of the slope into Voronoi polygons in order to determine the longdistance debris flow effect and to make the simulation results more in line with the actual failure process of the slope. 2. In this simulation, the slope experienced overall collapse and sliding after undergoing only 2 cycles of blasting vibrations, which is not completely consistent with the frequent blasting actions in actual mining. In actual mining, the blasting positions were usually far away from the slope; however, the vibration waves input into this simulation was close to the slope because of its purpose of exploring influence of the stoping blastings near the village and the initial slope, and therefore, the simulation result was qualitative rather than quantitative. Specifically, blasting vibrations are determined to be the key controlling factor of the dynamic response mechanism of the slope, while the corresponding quantitative analysis remains to be conducted further. However, it is the first time to explore dynamic formation mechanism of the landslide triggered by the local underground coal mining, especially by the vibrations induced by the continuous stoping blastings.

Conclusions
Based on the modeling results and related analyses, the following conclusions are drawn. First of all, the micro-landform, the lithological association of the slope, the structure and tectonics of the rock mass, and the weathering are pre-conditions that caused this landslide. When the coal seams beneath the original Pusa slope were mined from bottom to top, the caving zone, fracture zone, and deformation zone in the strata overlying the goaf gradually expanded. As the uppermost coal seam (C6) was mined, the deformation zone in the strata overlying the goaf overwhelmed the slope. During this period, the slope was in a phase of cumulative deformation, which was apparently characterized by subsidence and clockwise rotation. Under action of the heavy rainfall, the slope, which had been affected by the goaf, entered the further cumulative deformation phase, the macroscopic deformation of which was also characterized by subsidence and clockwise rotation. Finally, when the stoping blasting vibrations (the triggering factor) lasted for 2 cycles, the deformed slope, which had been impacted by the goaf and the heavy rainfall, entered the dynamic response stages of continuous progressive deformation, critical failure, shattering, ejection, debris flow, and accumulation. During this process, the macroscopic displacement of the slide mass exhibited significant subsidence, collapse, and clockwise rotation. The continuous vibrations caused by the underground blastings acted as the key controlling factor that triggered the landslide.

Data availability statement
The data that support the findings of this study are available from the corresponding author, [Fangpeng Cui], upon reasonable request.