Deformation Behavior of Mining beneath Flat and Sloping Terrains in Mountainous Areas

Slope structures and surface terrains are two signi ﬁ cant factors a ﬀ ecting the deformation behavior of mining slopes in mountainous areas. This research is aimed at investigating the deformation characteristics of a mining slope wielding Particle Flow Code (PFC), with 9 di ﬀ erent mining con ﬁ gurations (i.e., horizontal distance from extracted panel center to slope shoulder, D = − 200m, − 150 m, − 100m, − 50 m, 0 m, 50 m, 100 m, 150 m, and 200m). A representative slope in Faer Town, Liupanshui City, Guizhou Province, China, was selected, which was characterized by soft and hard interbedded rock strata. The results indicated that the overlying rock mass tended to move towards the sloping surface with mining beneath sloping terrain, which brought an asymmetrical subsidence funnel, and formed a wider relative disturbance range on the slope surface. With the vertical subsidence increasing additionally, the stability of the overall slope deteriorated. A safe mining range should be proposed based on evaluating the time-dependent deformation behavior at the slope shoulder and the overall slope stability.


Introduction
Longwall mining is one of the most generally adopted underground mining methods, particularly in mining areas with relatively uniform and thick coal beds [1]. Surface subsidence is the leading form of mining-induced geological hazards which has caused various adverse effects to the environment [2,3].
Analyzing the mechanism of ground movement and estimating their magnitudes and geometries have long been the main concerns in risk management of mining operations. A reliable prediction of ground subsidence caused by the mining operations remains a great challenge [4]. An extracted panel formed, the deformation of the overlying rock mass depends on many factors, i.e., bedding structure, thickness, strength, discontinuous geometric, mechanical characteristics of the panel, thickness of the mined coal seam, and width and length of the extracted panel [5][6][7][8]. In addition to the properties of the coal seams and the overlying rock masses, in situ stress conditions, groundwater conditions, terrain gradient, mining method, process of extraction, and distribution of pillars may all add complexity to the ground subsidence estimation [4,9,10]. The ground deformation induced by the mining operation is thus a multifactor coupling problem in temporal and spatial scales.
The prediction of ground subsidence can be performed by various methods, such as numerical simulation, physical modelling, influence function method, empirical approach, and analytical technique [10][11][12][13][14]. In addition, with the advancement in geographic information system (GIS) and remote sensing (RS) technologies, the spatial database can be constructed to analyze the shape and magnitude of the subsidence [15,16]. The capability of GPS network application in measuring ground horizontal displacement has made great contributions to monitoring and early warning systems for mining-induced geological hazards [17]. Application of geophysical methods enables geologists to develop a comprehensive understanding of fracture evolution in the overburden rock masses induced by shallow mining activities [18,19].
Many researchers performed numerical and physical modelling to estimate the subsidence of complex ground profiles. The physical modelling normally has difficulties dealing with an in situ stress state of rock mass (i.e., effect of gravitational force), which can be only simulated by geotechnical centrifuge [20,21]. Performing a large geotechnical centrifuge test, however, can be extremely costly. The numerical simulation has been widely used to analyze mining slopes with complex geometries and simulate discontinuous and nonlinear mechanical behavior of rock masses [4,9,[22][23][24][25][26][27]. The Discrete Element Method (DEM) is an efficient tool for analyzing instability of jointed rock slope, as proven by numerous successful cases [28][29][30].
There were in total 356 incidents of geological hazards (refer to Table 1) reported in the mountainous area in Liupanshui City, Guizhou Province, southwest of China. By the end of 2016, most of the hazards occurred on the P 3 l and T 1 f strata. A typical mining slope named Jianshanying slope in Faer Town, Liupanshui City, was selected. Wielding Particle Flow Code (PFC) to analyze deformation behavior with underground mining operations.

Mechanisms of Mining-Induced Slope Instability
Numerous previous studies reported that slopes subjected to underground mining are prone to caving and landslides [31][32][33][34]. There were 2 main models of mining-induced slope failures: (1) progressive rock falls and caving failures, which cause nearly vertical cliffs [25,35], and (2) rock masses extruded from the slope toe causing holistic instability [36]. With underground extraction performed beyond the slope shoulder, most instabilities are initiated [37,38]. According to Salmi et al. [4], the surface topography has a considerable impact on the mechanisms of mininginduced slope instability. Mining in hilly and mountainous terrains usually increases the risk of slope failure. In addition, mining beneath sloping terrains imposes an additional threat on ground subsidence, which usually occurs near a valley [35]. The coal seam extracted beneath flat terrains, the rock masses above the extracted panel were fragmented and caved into the panel immediately, and the fragmented rocks tend to fill the void forming a goaf. As a result of ascending step-loading imposed by the upper caving block, the stiffness of rock mass increases gradually [39]. The overburden strata remain intact and bend towards the extracted panel [40]. Owing to the expansibility of rock, the rock masses falling into the panel are subjected to lateral forces from the virgin strata, which rise gradually with the increase of depth and reach the maximum value at the coal pillar [41]. The boundary conditions on both sides of the extracted panel are identical, and hence, the magnitude and direction of lateral forces are completely symmetrical. Under the circumstances, there are different horizontal deformations in the rock masses, and it causes a symmetrical subsidence funnel on the surface [42]. In an extraction performed near or beneath a cliff, the lateral forces induced above the goaf (directed from the plateau towards the valley) are not counterbalanced by an identical force in the opposite direction. The rock mass near the valley has a greater horizontal displacement, and hence, an asymmetric subsidence funnel occurs on the slope surface.
The magnitude and shape of the surface subsidence, which is induced by mining operations under the condition of flat terrains, have been studied extensively. The localized deformation and overall instability constitute much more uncertainties attributed to the complex combination of terrains and structures in the sloping terrains. Studies on ground subsidence induced by mining activities in sloping terrains are still very limited. Several previous case studies of large-scale slope failures induced by mining include the Zhangjiawan collapse and Madaling landslide in Guizhou, China [36], and one of the largest contemporary landslides and mass movements reported at Nattai North, Australia [43], have drawn the global attention on the mass movement caused by mining activities and provoke the present study to be carried out.
This research is rooted in the exploration of mininginduced subsidence rules in P 3 l and T 1 f with Liupanshui City, Guizhou Province, China, as the typical. To be detailed, this research innovatively proposed 9 mining configurations for expounding the deformation behavior of mining beneath flat and sloping terrains by wielding Particle Flow Code, which is applied to the soft and hard interbedded and jointed slopes in the mountainous area. Furthermore, the time-dependent deformation was measured at the slope shoulder, as the transition part of the flat and sloping terrains, to propose a safe mining range, which was meaningful in the risk management of mining operations.

Case Study
3.1. Model Development. Liupanshui City in Guizhou Province, China, is known for its proven coal resources and reserves. The city which is known as the "Southwest Coal Sea" has developed multistage coal seams on the P 3 l stratum.
In 2015, there were more than 800 mines in the city, while more than 400 landslide and ground subsidence incidents have been reported, which were mainly caused by improperly planned mining activities. In the present study, a mining slope named Jianshanying slope in Faer Town, Liupanshui City, was selected as the case study. The specific geographical coordinates of the study area are E104°44′11″ and N26°18 ′ 20 ″ (Figures 1 and 2).
In the western part of the Guizhou plateau, a low mountainous terrain was formed because of the tectonic erosion. Typically, both steep and gentle structures were formed in the mining slope. Most of the coal-bearing strata are located in the flat terrains, while the interbedded sandstone and mudstone are mainly located in the sloping terrains. The surface terrain of the Jianshanying slope was reasonably generalized to simplify the subsequent numerical modelling pro-cesses. Three sets of dominant joints were considered in each strata, and one set has the same tendency as the strata. The general stratifications of the slope are presented in Figure 3. Coal and mudstone formed the relatively weak strata in the slope, however, the effect of the stratified structure was not considered in them.
Mining slopes are typically prone to time-dependent failures in the form of ground subsidence and slope sliding [44]. After going through a long process of mining, 6 coal layers had been mined out beneath the Jianshanying mining slope forming a total of 13 mining panels. These mining activities had caused severe impact on the stability of the slope. To simplify the analysis, the present study only focused on the impact of mining with the first coal layer on the slope insta-  Figure 4 and Table 2). The width of each extracted panel along the strata dip direction was within the range of 150~250 m, while the thickness of the coal seam was ranging from 2 to 4 m in  3.2. Particle Flow Code. Rock masses are discontinuous medium, and hence, the use of the DEM is justifiable [4]. PFC (Particle Flow Code), a popular program based on the DEM, is widely used to simulate the macroscopic characteristics of rock-soil masses. The soil/rock aggregates are modelled as either rigid disks (2D) or spheres (3D), and they are

Geofluids
connected by specific contact models as an equivalent model of rock-soil mass [45]. The PFC adopts the time step iterative calculation method ( Figure 5). Newton's second law and the law of force and displacement are repeatedly applied in the calculation for updating the motion state of units in real-time, and the contact force and torque between the updated units are further determined by the force-displacement relationship [46]. The law of force-displacement reflects the contact relationship between particles, also the relationship between the contact force and relative motion. In the PFC model, the contact force ball-ball and ball-wall can be divided into normal force and transverse force (Equation (1)). The particles move and rotate under the action of unbalanced forces and unbalanced torques (Equations (2) and (3) (8)). The definitions of the model parameters are summarized in Table 3: _ ω 3.3. Rock Mass Parameters. In particle flow simulation, the macroscopic mechanical behaviors of the rock and soil masses are governed by the microscopic mechanical properties of particles, nevertheless, there is a highly nonlinear relationship between them. Typically, the transformations of the macroscopic and microscopic parameters are carried out by means of biaxial compression tests [47]. The most common set of siltstone and pelitic siltstone was selected for parameter calibration to avoid the discreteness of rock samples. The stress-strain curve obtained from the PFC simulation under the condition of no confining pressure was reasonably consistent with that of the laboratory ( Figure 6). Both siltstone and pelitic siltstone showed significant brittle failure characteristics. The initial balance was carried out after gravity loading in the process of engineering scale simulation, and the increment of displacement and velocity during the process was cleared, and hence, the consistency of the stress-strain curve in the compaction stage was superfluous. The elastic modulus (E) and unconfined compressive strength (UCS) of siltstone and pelitic siltstone obtained from the laboratory and PFC are shown in Table 4; moreover, the fitting degree of their magnitude values is a measure index of calibration. Both mudstone and coal retrieved in the field had great discreteness with mechanical properties, which obstructed calibration by the PFC test. This research attempted to bring empirical values to the parameters of coal and mudstone, and the full mining model simulation was used in comparison with the actual situation. The parameter inversion method was used to adjust the rock block parameters as empirical.
The mechanical parameters of rock masses are generally smaller than those of intact rock with laboratory scale because of the size effect and discontinuity of rock masses [48]. Practically, the effects of bedding plane and dominant joints are often considered in simulation, and the equivalent jointed rock masses technology is applied. The smooth-joint model was chosen over the flat-joint model, which is poor in simulating the plane dilation mechanics, to reflect the constitutive relation. Furthermore, the smooth-joint model enabled the joint properties to a limited range on both sides of them, and a random joint model was formed to verify and correspond to the characteristics of the slope on-site. For this purpose, based on the laboratory mechanical experimental results, numerical simulation calibration testing, and the equivalent rock masses technique, the full mining model (the six-coal-seam mining model) simulation was used in comparison with the actual situation, which ensured that the simulation outputs were reasonably consistent with the actual field deformation (i.e., occurrence of deposition at the slope toe, presence of tensile cracks in the middle of the slope, and subsidence at the edge of slope), which was acquired through the UAV survey ( Figure 7). In this paper, parameters related to rock and soil masses materials were obtained and adjusted by parameter inversion. The calibrated microscopic mechanical parameters of the rock masses, which were adopted for the PFC simulation of the Jianshanying mining slope, are summarized in Table 5.  [49]. Obviously, tensile cracks were observed beneath the slope shoulder with D = −100 m. The panel was extracted directly below the slope shoulder (D = 0 m); the tensile cracks had developed in front of the slope shoulder and extended to the sloping terrain (Figure 8(e)). These results implied that there was an "aggregation" on the extension of tensile cracks beneath the slope shoulder.
According to Salmi et al. [50], neglecting the effect of stratum bedding in simulating mining-induced subsidence in flat terrains would yield a wider but shallower subsidence trough as compared with the field conditions. Therefore, the strata bedding surface and joints should be carefully modelled to improve the simulation outputs ( Figure 3).
The sum of upside and downside crack angles (γ + β) was wielded to characterize the relative disturbance range of overburden rock masses. The term "relative" was used to indicate that the thickness of rock strata, which is above the extracted panel, was not taken into consideration. A low value of "γ + β " indicated a large relative disturbance range, and vice versa (refer to Figure 9). The extracted panel is partially or completely located beneath the sloping terrains; the relative disturbance ranges for the cases were greater than those beneath the flat terrains. With the upside crack angle decreasing, the subsidence trough of the latter was wider than that of the former (see in Figure 8). Moreover, the propagations of tensile cracks at the extracted panel boundary were almost parallel with all mining configurations. The crushing of the coal pillar resulted in an increment in the distance between the position of boundary tensile cracks and the center of the extracted panel and hence altered the crack angle ( Figure 10). The center of the extracted panel positioned in front of the slope shoulder; the otherness between upside and downside crack angles increased. The upside crack angle reached the minimum value of 64°, and the summation angle of "γ + β" reached the minimum value of 151°with D = 0 m; moreover, the rock masses above the extracted panel were disturbed to the greatest.

Evaluation of Horizontal Displacement.
Coal mining causes significant vertical deformations. For materials which are characterized by low compactness and high expansibility in the subsidence area, the lateral compression of the strata surrounding the extracted panel would increase and cause an expansion to the sloping terrains. Subsidence immediately causes lateral deformation with the constraining forces of surrounding rock mass. A lower confining pressure makes the effect of lateral deformation more prominent [41]. Therefore, symmetrical and high constraining forces make the lateral deformation inappreciable [51]. In sloping terrains, the overburden rock masses produce relatively low lateral constraining forces, which are insufficient to offset the dilatational forces of rock masses caving into the panel. As a result, the disturbed rock mass would displace towards the sloping terrains ( Figure 10). The extraction panels are located at different positions; Figure 11 shows the lateral deformations of rock masses. The lateral deformation of overburden rock mass above the extracted panel was not symmetrical, with flat-sloping terrains as simulated. The lateral deformation beneath the flat terrain side was lower than that of the sloping side. The extraction panel was located close to the sloping terrain; the lateral deformation was intensified and caused an outcrop towards the sloping surface.
Bedding planes provided a suitable path for lateral movement of strata in both flat and sloping terrains. A bedding plane, with a low bonding strength, provided less resistance in the direction of the overburden material movement and hence caused the sliding between the layers. The "zigzag" horizontal displacement change zone can be seen in Figures 11 and 12.
Interestingly, lateral deformations towards the slope inner part were observed in the mudstone layer with D < 100, with    6 Geofluids the maximum value reaching 0.98 m (Figure 12(b)). However, there was no similar phenomenon in the mudstone directly above the extracted panel ( Figure 11). With the in situ stress releasing, the materials at the slope shoulder poorly cemented were further loosed. Holding a more complete and dense layered structure, the siltstone and pelitic siltstone were subjected to stick-slip resistances along the bedding plane in the lateral motion. Conversely, the resultant force, including the gravitational force, redistributed stress, and the cementing force between the materials, leading to the deformation of mudstone. A potential through slip plane appeared on the slope with D = −100 m (Figure 11(c)), which was initiated from the inner boundary of the extracted panel, extended upward to the goaf and the thin mudstone layer, and subsequently spread from the outer boundary of the extracted panel to the toe of the slope, causing the overall instability of slope. A horizontal displacement was observed of 0.4 m at the toe.

Evaluation of Surface
Subsidence. The maximum surface subsidence is consistently located above the inner part of the extracted panel under various mining configurations in the countertilt slope. Furthermore, the sloping surface with thinner overburden materials has a larger maximum subsidence area, as indicated by the cases of D ≥ 0 m (Figure 13(a)). These results proved that the slope has reached the "sufficient mining conditions" with D ≥ 0 m. For obtaining the increment in surface subsidence beneath the sloping terrains, the maximum subsidence (W 0 ) of flat terrains was brought to the present research, which was referred to Equation (9) proposed by Zou [52] under the "sufficient mining conditions":

Geofluids
where q is the subsidence coefficient under the "sufficient mining conditions," m is the thickness of the mining coal seam, and α is the dip angle of the coal seam.
Wielding the lithology comprehensive evaluation index ðpÞ to characterize the degree of influence of lithology on surface subsidence [52]: where h i is the thickness of overburden rock strata and Q i is the lithologic classification index of overburden rock mass. The value of Q i ranges from 0 to 1 for the first mining slope according to the hardness of the lithology. The Q i values for the coal seam, mudstone, pelitic siltstone, and siltstone in this research were set at 0.9, 0.8, 0.4, and 0.05, respectively. Based on the comprehensive evaluation index of lithology ðpÞ, the subsidence coefficient ðqÞ of flat terrains under the "sufficient mining conditions" can be acquired by the following: Khanal et al. [53] suggested that the ratio of subsidence to thickness of overburden rock mass (S/T) could be positively correlated with the width to depth ratio of the mine (W/D 1 ). In this research, the thickness (T) and width (W) were fixed at 6 m and 200 m, respectively. The findings from the mining  8 Geofluids configurations of D < 100 m showed consistency with that of Khanal et al. [53], but an opposite trend was observed with D > 100 m ( Figure 13).

Evolution of Deformation at Slope
Shoulder. The deformation behavior of the slope shoulder, a transition from flat terrains to sloping terrains, has been focused on in this research. The monitoring data of M13 was selected to characterize the deformation of rock mass at the slope shoulder.

Discussion
This research investigated the mechanical mechanisms of mining activities in flat and sloping terrains, with the considerations of deformation of overburden rock masses and propagation of tensile cracks. A model was established based on the typical mining-induced slope structure of a case study in southwest China. To the authors' knowledge, this is the first research reported on the use of PFC for analyzing the mining slope deformation behavior, with the extracted panel located beneath various complex terrains (i.e., flat terrain, slope shoulder, and sloping terrain). The time-dependent deformation characteristics of the rock mass at the slope shoulder were studied in detail. The sum of upside and downside crack angles was proposed to characterize the relative disturbance range of overburden rock mass, and reasonable and safe configurations of single-layered mining operations in mountainous areas were put forward.
It should be noted that the findings from the present research were handicapped by several limitations, such as the width of the extracted panel along the inclined strata was remained constant at 200 m, and the interval of two adjacent extracted panels in all mining configurations was kept at 50 m. In addition, the microstrength parameters of coal and pelitic siltstone were not derived from the uniaxial compression simulation by the PFC. Owing to these limitations, the functional relationship between surface subsidence and geo-logical and geotechnical factors cannot be fully revealed in this research. These issues can be solved if the following future improvements are taken: (i) setting the extracted panel width as an independent variable and shortening the interval, (ii) increasing the number of samples for coal and pelitic siltstone and acquiring the mechanical parameters by wielding the support vector machine (SVM) coupled with the microstrength parameters by PFC, and (iii) setting the number of configurations of mining in sloping terrains to be 5 to 10 times of independent variables, for finding the regression relationship between the maximum subsidence and geological and geotechnical factors.   13 Geofluids comparison with the actual situation, which ensured the rationality of microstrength parameters of rock masses and joints. The following conclusions can be drawn for the study of deformation behavior of mining slopes in mountainous areas with gentle anti-incline overburden rock strata:

Conclusion
(i) Mining in mountainous areas usually meets with the risk of slope instability. The overlying rock masses tended to move towards the sloping surface with mining beneath sloping terrain, which brought an asymmetrical subsidence funnel, and formed a wider relative disturbance range on the slope surface. In particular, the rock masses above the extracted panel were disturbed to the greatest with D = 0 m (ii) The constraining forces of overburden rock masses towards the valley decrease; mining beneath sloping terrains usually acquire larger subsidence (up to 3 m) and additional horizontal displacement (up to 1.4 m) than flat terrains. The "zigzag" horizontal displacement change zone formed with the control of bedding planes (iii) It is important to remain the center position of the extracted panel behind the slope shoulder, beyond which the deformation rate and ultimate value of the rock mass at the slope shoulder would increase drastically. In addition, the boundary of the extracted panel should also be avoided to cross over the slope shoulder (i.e., D ≤ −100 m) to prevent an overall slope instability

Data Availability
The simulation codes and slope information used to support the findings of this study are available from the corresponding author upon request.