There occurred a lot of catastrophic landslides in Southwest China karst areas in recent years. This study numerically simulated the failure process of Pusa landslide by distinct element method (UDEC). The Voronoi diagram algorithm was used to discretize the rock blocks in the upper part of the mining area, and the erosion in sliding source area caused by heavy rainfall was simulated by increasing the block density and reducing the joint strength. Results obtained using UDEC matched the failure process recorded by UAV, and the subsidence of monitoring points on the slope surface was well coincident with the InSAR results. The DEM results are as follows: firstly, under the action of underground coal mining, the overlying strata of the mountain deform toward the goaf and the deformation of the stratum increase significantly with the caving of the roof stratum; and the deep karst fissures in the overlying stratum have an important influence on the deformation of the mountain; the upper strata form a subsidence zone along the fissures, expanding the range of tensile fissures. Secondly, due to the excavation and unloading of the rock masses at the leading edge of the sliding area, the extrusion deformation occurs to the outside of the slope, resulting in the critical instability of sliding area. Finally, with the strength decreasing of rock masses and structural planes by heavy rainfall, the stability of the mountain continues to decrease until the sliding area collapses.

1. Introduction

Mining-induced landslide is a research that has been concerned for more than 100 years [1, 2]. With the continuous exploitation of coal resources in mountainous terrain, catastrophic landslides induced by underground coal mining occur frequently, resulting in a large number of economic losses and casualties [310]. The failure process and mechanisms of mining-induced landslides are becoming an increasingly important research [1113].

Although there were several methods for predicting rock movement induced by coal mining, such as influence function methods [14, 15] and physical test [16, 17], there is no doubt that numerical method can reveal the failure process more comprehensively [4, 18]. Finite element method (FEM) has been widely utilized in engineering, but it is more suitable for analyzing small deformation problems rather than large deformation problems like landslides [19]. On the other hand, discontinuous deformation analysis (DDA) [20] and distinct element method (DEM) [21] can simulate the large deformation of blocks. Do and Wu [13] and Xiong et al. [4] studied a mining-induced rock avalanche at Nattai North, Australia, using DDA and DEM methods, respectively, which verify that DDA and DEM can be applied to clarify the failure process of a mining-induced landslide. UDEC is a mature DEM software, verified by Cui et al. [22] that it can reveal the progressive failure process of Daguangbao landslide induced by Wenchuan earthquake. Meanwhile, compared with the DDA software, the Voronoi diagram algorithm has been successfully applied in UDEC for the stability analysis of different slope failure types [2325]. With the Voronoi diagram algorithm, the large rock block can be discretized into many small parts, which promotes a more accurate study of mining fracture propagation and sliding surface formation [26].

Furthermore, landslide is a progressive failure process, affected by several environmental factors. A large number of landslides indicate that the impact of heavy rainfall should not be ignored [18, 27, 28]; there are still a few studies about the impact of rainfall on landslide based on DEM. In addition, the unique geological structure of landslide area and the existence of karst fracture zone will reduce the sliding resistance of source area, which further increases the complexity of numerical method analysis [2932]. Moreover, fault is also a frequent geological structure of rock slope [33]. Because of the broken rock masses near the fault, the hydraulic characteristics of underground rock masses (permeability, porosity, etc.) [34, 35] are strongly affected. When there exist faults near the underground mining activities, fault in some cases has a significant impact on landslides, resulting in some irregular deformation responses. The Pusa landslide induced by underground coal mining in Guizhou Province, China, including heavy rainfall, fault, and complex karst geology, is a perfect candidate to study mining-induced landslide in karst area with heavy rainfall [36]. Although Zhu et al. [37] used DAN3D to study the runout behavior of Pusa landslide, there remains a lack of slope failure mechanism analysis induced by coal mining.

In summary, many researchers have used numerical methods and model tests to study the failure process of mining slopes. In addition to the slope factors such as weathering action, underground mining, karstification, rock mass structure, and lithological, it is also strongly affected by heavy rainfall factors, which has a stronger influence on karst slope. To our knowledge, there are no reports on the mining-induced landslide considering karst fissures, faults, and heavy rainfall in one model. It is necessary to comprehensively consider the above factors to study the instability mechanism of mining slope. In this paper, with the help of the Voronoi diagram algorithm, the UDEC software is used to simulate the Pusa landslide under mining operation. In this numerical model, the ground fissure characteristics before the landslide are considered, and the distribution of the mining areas is investigated. The slope erosion caused by heavy rainfall is simulated by increasing the density of the sliding source area and reducing the joint strength. The research was to reveal the failure mechanism of Pusa landslide, so that it can provide a reference for the evaluation of the long-term stability on mining-induced slopes.

2. Data and Method

The mega landslide in Pusa, Guizhou Province, China, is selected as the research object in this paper. Firstly, a survey involving the basic geological conditions was performed to understand the mechanism of landslide. Secondly, a detailed field investigation was conducted on the karst characteristics of the landslide and the distribution of the mined area, which further clarified the failure mechanism of the landslide. Finally, a discrete element numerical model was utilized to study the damage and contributing factors. In addition, the deformation results of the model were verified by the InSAR data provided by Chen et al. [38].

3. Engineering Geological Conditions

3.1. Geological Setting

Pusa landslide occurred on August 28, 2017, in a village of Zhangjiawan Town, Nayong County, Guizhou Province, China (N26°3804.55, E105°2656.14) [36]. More than 500000 m3 rock masses fell down from the mountain ridge, causing 26 deaths. The ridge of the Pusa landslide located at 2170 m a.s.l, and the foot of the mountain located at 1800 m a.s.l. The study area has a subtropical monsoon climate, with an average annual temperature of 14°C and annual precipitation slightly higher than 1200 mm.

As shown in Figure 1, the exposed formations in the study area include: (1) the Longtan Formation of the Upper Permian (P2l), consisting of silty mudstone, argillaceous siltstone, mudstone, and coal seams; (2) the Changxing–Dalong Formation of the Upper Permian (P2c + d), consisting of limestone, argillaceous limestone, and sandy shale; (3) the Yelang Formation of the Lower Triassic (T1y), mainly composed of limestone and marlstone, with a small amount of sandstone and mudstone.

The occurrence of rock layer is 180°∠8°, and the direction angle of the sliding direction is 320°. Based on field investigation, there are four structural planes in the overlying rock masses of slope body. Their occurrences are (1) 110°∠87°, (2) 326°∠86°, (3) 140°∠88°, and (4) 225°∠88°, respectively. According to the investigation by Fan et al. [36] and Chen et al. [38], there existed many long and deep ground fissures at the trailing edge of the mountain before 2015.

There are two faults F1 and F2 in the mining area, which have some effects on the strata and coal seams in the mining area: (1) fault F1 is close to the outside of mountain, belongs to a normal fault, strikes northeast, leans south-north, dips 63°~70°, and is about 1320 m long; (3) fault F2 is close to the inside of mountain, belongs to a reverse fault, strikes northeast, leans south-north, dips 70°~75°, and is about 709 m long.

There exist 6 mineable coal seams in the Longtan Formation, the attitude of which is consistent with that of the strata. From top to bottom, they are No. 6, No. 10, No. 14, No. 16, No. 18, and No. 20 coal seams, respectively. The inclination angle of the coal seams is generally 7-12°. The average thickness of coal seam is about 1.6 m, and the maximum total thickness is 9.6 m.

According to the field survey results, the Pusa landslide is mainly divided into three areas: the source area (zone A in Figure 2), the track-abrasion area (zone B in Figure 2), and the accumulation area (zone C in Figure 2).

The average height of the caving rock mass in the source area (zone A) is about 85 m, the width is about 145 m, the average thickness is about 40 m, and the total volume is about . The upper part of the rear wall is a steep medium-thick-thick limestone with a thickness of about 20 m. Since the lower rock mass has collapsed, two concave cavities have been formed in the lower part of the limestone. There are scratches on the middle and lower rock walls, with an inclination of 305°. The main rock mass near the shear outlet is siltstone and marl. Although it has been affected by multiple groups of structural planes, it still maintains a large block structure. The maximum size of the rock mass in the front and middle of the collapse is .

The track-abrasion area is mainly distributed between the elevations of 1920 m and 2030 m. When the unstable rock masses suddenly collapse, the rock scrapes the original loose deposits on the lower slope surface with huge potential and kinetic energy. The scraping action made the steep wall surface in this area retreat about 1.5 m. The track-abrasion area is about 180 m wide and 80 m high, and the track-abrasion area is about .

The deposition area was a relatively open gentle slope area before collapse. The middle part of the deposition area is distributed with a lot of rubble, mainly composed of siltstone and marl. The average length of deposition area is about 575 m, and the width is about 360 m. It is estimated that the thickness of the deposition area is about 1.5 m, and the amount of shovel scraped in the deposition area is about . After the collapse, a series of tensile fissures with a strike of 35° were produced behind the mountain (zone II in Figure 2). The largest tensile fissure is about 180 m long and 12 m deep, with an opening degree of about 34 m.

As shown in Figure 3, karst features are developed in the limestone layers at the top of the slope, including dissolution cavities, dissolution pipes, and dissolution fissures on the steep cliff walls. A series of dissolution phenomena can be clearly seen in the limestone on the ridge, including the dissolution of the rock out of the cavity, a series of small depressions distributed on the slopes, and the existence of limestone blocks in the accumulation. On the shoulder of the slope (about 2110 m above sea level), four groups of deep and large karst fissures are developed downward (shown in Figure 1), with a depth of 50-80 m and a width of 2-3 m.

3.2. Triggering Factors of Landslide
3.2.1. Coal Mining Activities

The coal mines in the study area mainly adopt inclined shaft development, the strike longwall retreating mining method is used to mine the mine, and the roof is treated by all caving method. From 2007 to December 2010, No. 16 coal seam was mainly under mining. From 2011 to 2012, the coal mine operation was suspended until No. 14 coal seam was started to mining in November 2013. After the mining of No. 14 coal seam was completed in September 2016, the mining of No. 10 coal seam began and lasted until the Pusa landslide occurred in August 2017. Meanwhile, the width of No. 10 and No. 14 mining areas is 90 m (in 1-1 cross section), and No. 14 mining area lies directly below No. 10.

3.2.2. Heavy Rainfall

Nayong County has a subtropical monsoon climate, with an annual average temperature of 13.6°C and an annual average rainfall of 1243.9 mm. As shown in Figure 5, the rainfall concentration period is from May to August. From June 21, 2017, to July 21, 2017, the rainfall in the study area increased sharply by 374 mm [36], which further deteriorated the stability of the relatively broken slope affected by goaf, karst fissure, and weathering. Rainwater flows into the slope with joints and fissures at the top, which accelerates the collapse and sliding of the slope from two aspects: (1) resulting in the increase of rock and soil bulk density and (2) promoting the deterioration of mechanical parameters of slope fissures and joints, such as tangential stiffness, tensile strength, cohesion, and internal friction angle [28].

4. Numerical Model

The two-dimensional distinct element code UDEC is widely utilized to study the failure mechanisms of rock slope [39]. In this code, blocks are separated by contacts that are considered analogous to discontinuities in rock masses; the rotation of rock blocks can also be simulated in this software.

4.1. Parameters of Rock Masses and Joints

In this model, the rock masses adopt the classic elastoplastic theory Mohr-Coulomb criterion, and the contacts adopt the Mohr-Coulomb criterion with residual strength. When the contact breaks, the residual strength parameters are used to replace the original mechanical parameters of contact. The simulation of the fracture process of the complete rock masses can restore the failure process of the slope more accurately, so the Voronoi grid division method [25] was used to discretize the block above the mining area, as shown in Figure 6. According to the study by Cui et al. [40], the parameters of rock masses with the Voronoi grid are listed in Table 1. The parameters of rock masses without the Voronoi grid are listed in Table 2. The parameters contacts, including interfaces and joints, can be found in Table 3.

4.2. DEM Model

The DEM model is established based on the cross section in Figure 1. The coal seams below No. 14 coal seam were mined before 2011; No. 10 and No. 14 coal seams were under mining after November 2013. However, the deep and large fractures behind the mountain were all produced before 2013. This means that the deep fissures behind the mountain are caused by the coal mining below No. 14 coal seam. Four deep fissures have been considered in the numerical model in order to be equivalent to the influence of the coal mining below No. 14 coal seam, so this numerical model before mining is in the actual situation after 2013. Except for No. 10 and No. 14 coal seams, the other coal seams have been mined out for several years before the landslide, so only the mining activities of No. 10 and No. 14 coal seams were analyzed in this model. These four fissures are 50~80 m deep and 1.5 ~ 2 m wide. The No. 14 coal seam is mined first and then the No. 10 coal seam. Each coal seam is set with 9 calculation steps, and each calculation step is 10 m. After the two seams are mined, rainfall conditions are applied to the sliding area. Monitoring points P1 and P2 (as shown in Figures 1 and 4) were monitored by satellites from September 2016 to August 2017 [41].

In this model, the stratum above No. 14 coal seam is discretized by the Voronoi grid, which will be more conducive to study the stratum movement and the crack expansion above the mining goaf. The influence of karst on the slope is reflected by setting four deep ground fissures and reduction of the parameters in weathering zone. In the field investigation, it was found that those four ground fissures (in Figure 1) were filled with many fillings, mainly including broken rock masses and dissolution weathering products. Therefore, the formation of ground fissures is mainly achieved by reducing the block parameters to the filling parameters. The geometry of the numerical model and boundary conditions can be found in Figure 6. It should be noted that point A and point B are, respectively, in the sliding source area and noncollapsing area, and the sliding surface is located between point A and point B.

5. Numerical Results

5.1. Characteristics of Rock Movement

The displacement cloud diagram and velocity vector diagram after coal mining and heavy rainfall are shown in Figures 7 and 8, respectively.

As shown in Figure 7(a), after the No. 14 coal seam was mined, the roof strata caved, and the subsidence basin of the overburden stratum expanded upward in a conical shape. The maximum deformation is roof caving, nearly 1.34 m. When the No. 10 coal seam was mined (Figure 7(b)), the deformation of overburden stratum was significantly increasing, forming a step-like bending subsidence zone, and continuously sliding along the deep ground fissures. Moreover, the front edge of slope was most affected. The rock masses at the slope toe were extruding out of the slope, which has an important impact on the stability of the mountain. The mining of No. 10 coal seam increases the differential deformation between the potential sliding area and stable area. However, the maximum displacement is about 3.95 m, approximately the same to the thickness of No. 14 and No. 10 coal seam, which indicates that the slope is still stable at this stage. After the heavy rainfall (Figure 7(c)), the differential settlement of the rock masses on both sides of the deep ground fissures increased sharply. The maximum displacement of sliding body is more than 6.8 m, locating at the lower position of the sliding surface.

No. 10 and No. 14 coal seams are located about 250 m directly below the top of the mountain. After the mining of the No. 14 coal seam, roof materials collapse into the mining panel; the upper rock layer bends and deforms, resulting in a maximum displacement of 1.34 m. The moving direction of the rock mass in the sliding source area is parallel to the free surface of the slope, without extruding outwards (shown in Figure 8(a)). However, the rock masses between the sliding source area and the fault are almost horizontally squeezed out. When the No. 10 coal seam was mined, the above phenomenon was obviously aggravated, with the maximum moving displacement of 3.96 m, and rock avalanches appeared in part of the rock mass in the sliding area (shown in Figure 8(b)). After the heavy rainfall, the slope was in a critical state of instability, and the rock mass displacement in the slid source area increased sharply and squeezed out (shown in Figure 8(c)).

Point A and point B in Figure 6 are selected to analyze the sliding era. Point A is located on the potential sliding body, point B is located behind the potential sliding body, and there is a deep ground fissures between these two points. Therefore, the deformation process of the potential sliding body can be quantitatively described by comparing the relative displacements of point A and point B.

It can be seen from Figure 9 that the horizontal and vertical displacements of point A are both larger than those of point B during coal seam mining and heavy rainfall. The vertical displacements of two monitoring points are always greater than the horizontal displacements at each stage. After the No. 14 coal seam is mined, the horizontal displacement and vertical displacement of point A are 1.6 times and 1.5 times of point B, respectively, the displacements of points A and B are not much different, and the potential sliding surface has not yet formed. When the No. 10 coal seam was mined, the horizontal displacement and vertical displacement of point A are 2.1 times and 1.9 times of point B, respectively, indicating that the potential sliding surface had been formed at this time. With the action of heavy rainfall, the ratio of horizontal and vertical displacements between A and B becomes 4.2 and 3.5, respectively, which means that the sliding body has begun to slide and the slope body has become unstable.

5.2. Characteristics of Stress Evolution

Figure 10 shows the distribution of maximum shear stress in the mining karst mountain. When No. 14 coal seam is mined (Figure 10(a)), the maximum shear stress at both ends of the void area is above 4.5 MPa. There is a large-scale stress redistribution in the karst mountain above the goaf. The roof strata broke due to the lack of support, causing the gravity of overburden rock masses to transfer to both sides of goaf. The pressure relief zone appears above the goaf, and the pressurized zone is formed on both sides of goaf; the pressurized zone extends up to the potential sliding surface.

After the No. 10 coal seam was mined (Figure 10(b)), the maximum shear stress increased to 4.97 MPa, and the height of the pressure relief zone also increased with the thickness of the mined coal seam. This indicates that when the thickness of the mined coal seam is not large, the rock layer above the goaf forms a composite beam structure to bear the overlying load. When the thickness of the mined coal seam increased, the roof strata collapse in a large area. The composite beam structure broke due to the large span, lying on the rock masses of the roof caving.

After heavy rainfall (Figure 10(c)), due to the strength reduction of the mechanical parameters of the overlying rock masses and joints, shear stress concentration appears on the potential sliding surface, which means that the sliding surface was slipping.

5.3. Comparison with Preevent Deformation Monitored by InSAR

The Pusa landslide, in Guizhou, China, occurred on 28 August 2017, and the preevent deformation has been recovered by Chen et al. [41] with multisensor synthetic aperture radar (SAR) imagery. The time series of satellite monitoring data is from September 2016 to August 2017; during which, only No. 10 coal seam was under mining (as shown in Figure 4).

To analyze the temporal evolution of the landslide deformation, deformation time series at two points (P1 and P2 in Figures 1 and 4) have been extracted. It should be noted that the monitored deformation is in the line of sight (LOS) direction, and the vertical deformation component can be retrieved as follows [42]: where is the vertical deformation component; and are, respectively, the deformation components along the descending and ascending LOS directions; and is radar incidence angle. According to research of Li et al. [43], the radar incidence angle is selected as 23°.

Due to the lack of the mining schedule of coal seam, the comparison between the InSAR results and the numerical results is selected after heavy rainfall. The comparisons are listed in Table 4 (positive values represent uplift, and negative values represent sinking down). The InSAR monitoring point P1 is located below the source area, and P2 is located behind the source area. It is clear that the DEM model of the Pusa landslide has a well prediction on the uplift of P1, but the error of prediction on point P2 is 31%.

6. Discussion on Failure Mechanism

The evolution process of rock stratum fissure caused by underground coal mining and heavy rainfall is shown in Figure 11. To represent the phenomenon of rock fissures more clearly, the deformation of the model is amplified by 6 times, that is, the “deformed factor” in Figure 11. The overlying rock stratum after coal seam mining can be regarded as beams or slabs. With the continuous mining activities, the span of goaf is getting wider, causing the increase of the bending moment at the midspan and the shear force at both ends of the support. When the stress at the most unfavorable position of the rock strata reaches load limit, the rock layer will break or even collapse. After the No. 14 coal seam was mined (Figure 11(a)), the roof strata collapsed, and the overlying rock stratum above the roof appeared a bending subsidence zone. There occurred local rock collapse in front of the slope body, forming a potential sliding surface. When the No. 10 coal seam is mined (Figure 11(b)), the roof caving zone keeps expanding, and the overburden rock stratum above the roof strata continues bending and sinking. There appeared rock layer dislocation near the fault, tensile fissures behind the slope, and compression-shear fractures in front of the slope, forming a potential sliding slope body. After the heavy rainfall (Figure 11(c)), the deformation of the potential sliding body increases sharply. The rock mass at the shoulder of the slope slides down along the deep ground fissure, the rock mass in the middle of the slope collapses and extrudes, and the slope is already unstable and damaged.

In general, based on the numerical simulation analysis of the Pusa landslide, the failure process of the Pusa landslide can be summarized as follows: (1)Mining operation before 2017 aggravated the deformation of the mountain and the propagation of fissures, resulting in the formation of deep and large fissures on the surface of the slope, which is also more conducive to water infiltration and dissolution development(2)Due to the superimposed mining disturbance of multiple coal seams, the rock mass at slope toe was squeezed out, the rock layer on the left side of the goaf was pulled apart, and a rock bridge was formed on the slope. The heavy rainfall made part of the ground fissures filled with water; the rock mass was gradually saturated, which increased the weight of sliding area and reduced the strength of the joints. The slope slipped down along the deep fissures and deformed, and the lower rock bridge was crushed and extruded. The slope has an overall dumping destruction

In addition, the study mountain is located on the Yunnan-Guizhou Plateau, and long-term sunshine and temperature differences will accelerate the fragmentation of the rock mass. This is also a factor that cannot be included in the numerical simulations in this paper.

7. Conclusions

In this paper, combined with the field investigation and numerical results, the movement characteristics of the overburden rock stratum and the evolution law of cracks in the Pusa landslide were studied. The failure process and failure mode of karst mountain collapse under the action of underground mining and heavy rainfall are proposed. The following main conclusions are drawn: (1)In the process of coal seam mining, the slope undergoes three stages: coal seam roof caving, overburden cantilever tension fracture, and overburden overall dumping and settlement. After the heavy rainfall, the source area began to slide, and the slope collapsed as a whole body. The failure mode of the Pusa landslide controlled by deep and large karst fissures can be summarized as slip and extrusion of the leading edge of the mountain, upward transmission of tensile fractures, collapse of the leading edge, and overall dumping(2)After the coal seam is mined, it is mainly manifested as loading on both sides of goaf and unloading above the goaf. With the continuous advancement of the mining space, the overburden rock layers break and deform, compacting the rock masses in the caving zone, also transferring the goaf clearance space to the sliding area and deep karst fissures(3)The long-term underground coal seam mining action, karst dissolution, and weathering action make the mountain stay in a critical instability state. The reduction of the mechanical strength of rock masses and structural planes caused by heavy rainfall before the collapse is the direct trigger factor of Pusa landslide

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


This work was funded by the National Natural Science Foundation of China (Grant No. 52074042) and the National Key R&D Program of China (Grant No. 2018YFC1504802).