In global mining, for a number of reasons, the number of mining operations considering an open pit (OP)–underground (UG) transition is significantly increasing and most of them (85%) are ore mines [1,2]. The application of such a mining technique in coal mines is uncommon due to difficult natural conditions, mainly the much lower mechanical parameters of rock mass compared to ore mines [3]. In the decision-making process, various factors should be taken under consideration. One of them is an assessment of the impact of UG mining on the stability of pit slopes. UG operations cause surface subsidence, deformations, and changes in hydrogeological conditions [4]. In OP–UG operation, the UG-induced subsidence can lead to slope deformation. Historical slope failure incidents, including the Palabora mine in South Africa [5], Ernest Henry mine in Australia [6], and Los Bronces mine in Chile [7], provide a clear warning regarding the impact of caving operations on an OP slope.

Limited studies were also conducted to prove the impact of longwall coal mining on an OP slope. Peng et al. [8] established a mined slope calculation model with the influence analysis of UG mining on slope, distortion, and stress distribution. Based on this model, the stability calculating method was analyzed and the judgment formulae were elicited, with the determination principle and judgment of minimum protecting coal pillar width of the slope drawn. Nguyen and Niedbalski [1] reported an attempt at evaluating the influence of longwall operations on the slope stability of a still-functioning OP coal mine. In this study, the slope face displacements were calculated, taking into account some considerations regarding UG operation, such as exploitation depth, direction, and the distance from the slope plane (face). Based on the outcomes, a protecting coal pillar was determined to avoid impact on the OP slope. Payne et al. [9] introduced a case study that has experienced a significant effect of longwall subsidence on the open cut highwall movement when the longwall operation approached their final position (stop line), close to the open cut highwall. Results from the monitoring (radar and laser scanners) found the highwall is displaced to magnitudes unlike those typically measured in open-cut mining and is in direct contrast to typical longwall subsidence behavior. Recommendations are made for mine and pit designs for future punch longwall layouts.

Generally, surface subsidence prediction methods can be divided into two categories: empirical methods and theoretical numerical methods. The accuracy and efficiency of surface subsidence prediction methods depend on location, knowledge of geologic and to pographic conditions, and the level of mining experience [10,11,12,13]. Empirical methods including graphical, profile function, and influence function methods are developed based on the mathematical assumptions between the surface movements and the mining geometry. The influence functions, such as Keinhorst's method [14], Bals’ method [15], and Knothe's method [16,17], are widely used due to their applicability to complex mine geometry and various types of mining situations and possibility of validation mathematically. Theoretical numerical methods are based on modeling the UG-induced surface subsidence mechanisms. Numerical methods enable the presentation of more complex relationships between strain and deformation based on the modeled properties of the rock mass. They can also incorporate the influence of varying surface topography, lithology, and large-scale structure [18]. The first numerical applications are presented in [19,20,21,22]. In recent years, with the development of computer science, numerical methods have become increasingly popular for solving the problems of geoengineering. They prove particularly useful when their results are used for comparison with the results of observations under natural conditions. Issues related to displacement and surface deformation predictions are also solved by applying numerical methods. The implementation of numerical methods for surface deformation prediction in coal mining using the longwall UG method was performed by many other researchers [23,24,25,26,27,28,29,30,31,32].

Currently, Vietnamese coal mining is a very important part of the national economy and hard coal is mined using both the OP method and the UG method. The decision to switch from OP mining to UG mining was taken for several reasons, including technological possibilities, economic and environmental conditions, and social aspects [33]. Coal production from the UG mines in Quang Ninh province is steadily increasing and is predicted to be approximately 75% of the total coal production in 2030 [34, 35]. Cao Son is an OP mine located in the Cam Pha coal basin, Quang Ninh. The Cao Son OP mine decided to change its operation system from OP to UG (Khe Cham II–IV). According to the latest operation plan, simultaneous operation started in 2022 and will be conducted until 2038 [33]. In order to evaluate the stability of the Cao Son slope, the surface and subsurface strain fields generated by UG operation at Khe Cham II–IV were investigated by using numerical modeling with an FDM (Finite Difference Method) code, FLAC2D [36]. In addition, the factor of safety as a slope stability measurement was performed for certain stages of OP slope with UG advance. Based on the obtained results, slope stability evaluations and practical recommendations for Cao Son and Khe Cham II–IV operation are presented.

The Cao Son mine is located in the Cam Pha coal basin in the Quang Ninh province, Vietnam. The OP is about 150 km east of the Vietnamese capital city, Hanoi (Figure 1).

The Cao Son OP mine decided to make the transition from OP to UG in 2015. The Khe Cham II–IV UG operation is planned to extract coal seams No. 9 (2022–2036) and No. 8 (2025–2040) using the longwall mining method. Simultaneous operation (Cao Son and Khe Cham II–IV) will be conducted until 2038, when the Cao Son OP mine ends its mining activities in coal seam No. 10 (Figure 2).

The Quaternary period rock mass of an average thickness of 20 m mainly consists of gravel, sand, and rolled fill. The mechanical properties of the Quaternary period rock mass are density 1800–2200 (kg/m^{3}), cohesion 25–130 (kPa), and friction angle 9°–31°. Below, the Triassic rock mass consists of hard rocks: sandstone (48%), mudstone (25%), conglomerate (16%), coal (anthracite) (8%), and claystone (3%) [33]. Thickness of the Triassic period rock mass ranges from 500 to 700 m. Average inclination angle of the rock layers is 10°–15°. The lithology and distribution of rock mass in the Cao Son and Khe Cham II–IV region is shown in Figure 3.

Tables 1 and 2 show the mechanical properties of rocks and faults in the studied region.

Mechanical properties of rock mass in the studied region [33].

^{t} |
^{3}) |
|||||
---|---|---|---|---|---|---|

Mudstone | 2.31 | 1.54 | 1.87 | 30 | 1.34 | 2620 |

Claystone | 2.23 | 1.34 | 2.05 | 26 | 1.14 | 2600 |

Anthracite | 2.17 | 1.36 | 2.14 | 27 | 1.22 | 1500 |

Conglomerate | 4.76 | 3.57 | 3.23 | 28 | 2.27 | 2510 |

Sandstone | 3.91 | 2.46 | 3.56 | 28 | 1.96 | 2600 |

Strength parameters of the fault in the studied region [33].

Value | 12–20 | 3.4–4.4 | 0 |

The FLAC2D program is based on the finite difference method, uses an explicit, time-marching method to solve the algebraic equations, but implicit, matrix-oriented solution schemes are more common in finite elements [36].

For the purposes of this piece of work, numerical calculations were carried out using the Coulomb–Mohr elastic–plastic model. The failure envelope for this model corresponds to a Mohr–Coulomb criterion (shear yield function) with tension cut-off (tensile yield function). The shear flow rule is nonassociated and the tensile flow rule is associated. The Mohr–Coulomb model takes into consideration the plasticity of the rock mass, which is the nonlinearity of its stress–strain characteristics. The plastic flow formulation in FLAC2D rests on basic assumptions from plasticity theory that the total strain increment may be decomposed into elastic and plastic parts, with only the elastic part contributing to the stress increment by means of an elastic law. In addition, both plastic and elastic strain increments are taken to be coaxial with the current principal axes of the stresses [36].

The factor-of-safety calculation can be performed for stability analyses in FLAC2D. This calculationis basedupon the “strength reduction method” to determine a factor of safety. The strength reduction method is an increasingly popular numerical method for evaluating the factor of safety in geomechanics. This method is typically applied in factor-of-safety calculations by progressively reducing the shear strength of the material to bring the slope to a state of limiting equilibrium. The method is commonly applied with the Mohr–Coulomb failure criterion. A series of simulations are made using trial values of the factor _{trial} to reduce the cohesion, ^{t}

The strength reduction method can be applied to essentially any material failure model to evaluate a factor of safety based upon the reduction of a specified strength property or property group. The method has been used extensively in the context of Mohr–Coulomb material and, principally, the simultaneous reduction of cohesion and frictional strength. In FLAC2D, in addition to Mohr–Coulomb strength properties, the method is also applied automatically to ubiquitous-joint strength properties and to Hoek–Brown strength properties. The strength reduction method can also be applied for interface strength properties, friction, and cohesion (e.g., faults) and to grout shear strength of the structural elements (e.g., cable, pile, and rockbolt elements) [36].

In FLAC2D, one of the many indicators that can be used to assess the state of the numerical model is the plasticity indicator. It determines the possibility of the failure of individual rock mass points as a result of tensile or shear stresses. Moreover, displacements, velocity vectors, stresses, strains, and so on can be plotted if necessary. These parameters were used to analyze the impact of the Khe Cham II–IV UG operation on the Cao Son OP.

Model includes a 3200 m long (x-axis) and 700–1100 m high (y-axis), presenting a 400 m open pit highwall, the coal seam No. 9 with an average thickness of 3 m located under the OP and two major faults crossed the pit slope. All models were conducted following typical boundary conditions for slope stability analysis: the bottom was fixed, a roller was applied to the sides of the model, and the top surface was set free (Figure 4). The modeled site study consisted of approx. 40,000 elements of four-node quadrilateral. The model was originally developed as an elastic model to achieve the initial stress state. Then, the displacement and velocity vectors were reset. In the next step, the “null” model was assigned to the zones corresponding to the extracted coal seam No. 9, and the Cao Son overburden and the model were recalculated.

The rock mass was modelled as a Mohr–Coulomb material using built-in constitutive model available in FLAC2D. A summary of the mechanical parameters for numerical modeling is presented in Tables 1 and 2.

Determining the geometry of the caved zone was the subject of a number of studies that were later used to develop various empirical formulas (Table 3).

Hypotheses for calculating thickness of the caved zone.

Peng and Chiang, 1984 [38] | (2–10) |

Bai et al., 1995 [39] | 100_{1}_{2}) |

Mazurkiewicz et al., 1997 [40] | _{r} |

Heasley, 2004 [41] | (10–18) |

Biliński, 2005 (simplified) [42] | (_{s}t_{c}^{0.5}+0.02) |

Wang et al., 2017 [43] | (3–4) |

c_{1}, c_{2} – constants dependent on the compressive strength of roof rocks,

_{r}

_{s}_{s}

_{c}

The height of the caved zone tends to be proportional to the thickness of the exploited seam. The variety of hypotheses is large due to the different mining conditions in which the tests were conducted. However, only Bai et al. [39] and Biliński [42] considered the compressive strength of roof rocks in order to estimate the height of the caved zone. Therefore, these hypotheses are considered to be the most reliable for determining the geometry of caved zone induced by longwall mining.

Access to the caved zone is limited. Therefore, it is very difficult to find equivalent mechanical properties that are reliable to express the heterogeneity of these materials.

Various studies were carried out to determine the equivalent mechanical properties of the caved zone. The elastic modulus of this zone can be calculated as a function of the compressive strength of undisturbed roof rocks and the bulking factor [44]. Tajduś [45] used the back analysis method to determine the value of parameters of a disturbed rock mass caused by mining. He found that the elastic modulus in the horizontal and vertical directions is very low and ranges from 50 to 150 MPa. Cheng et al. [46] and Jiang et al. [47] assumed that the elastic modulus and Poisson's ratio in the caved zone are 190 MPa and 0.25, respectively. Ahmed et al. [48] suggested that the elastic modulus of the caved zone is about 2.1% of the elastic modulus of roof rocks.

In order to assess slope stability for the 2025 stage and 2030 stage of the Cao Son pit slope, the factor of safety was implemented with the Shear Strength Reduction technique. Next, numerical calculations were carried out to show the impact of UG Khe Cham II–IV operation in coal seam No. 9 on the pit slope in 2025 and 2030. According to the UG operation plan [33], about 20% of coal seam no. 9 under the Cao Son pit slope will be mined by 2025 (two longwall panels which are each 160 m long) and approximately 50% of the coal seam will be mined by 2030 (five longwall panels, each 160 m long). Three different scenarios were considered for the investigation of the UG impact on the 2025 slope stability: I – UG operation starts from the left-hand side; II – UG operation starts from the center; and III – UG operation starts from the right-hand side (Figure 5). Then, the best option of the extraction location was selected. From the selected option, two other scenarios were considered in order to examine the impact of further UG extraction on the 2030 slope stability: IV – extracting the further longwall panels near the mined panels, V – extracting the further longwall panels in another side of the mined panels.

The results of the numerical calculations are presented in the form of maps of the plasticity indicator and maps of the FoS (Factor of Safety) contour, which enable assessment of the impact of UG on slope stability. Due to the large number of results obtained, only selected results for certain calculation variations are shown in this manuscript.

Plasticity indicators were also calculated for scenarios I, II, and III, as shown in Figure 8. The failure zone will reach the slope surface in scenario II and scenario III by 2025 (Figure 8b and c). These failures could negatively affect OP operations such as transportation machinery and equipments. Additionally, they could form the path-inrush water from the pit floor to further UG exploitation in the coal seam No. 9. The failure zone in scenario I did not reach the slope surface (Figure 8a). Hence, scenario I is considered to be better than scenarios II and III in terms of UG impact on slope stability.

The calculated FoS values for the three scenarios (I, II, III) with the impact of UG operation are shown in Figure 9. In general, the size of the slope failure surface was larger and deeper inside the slope body in all scenarios. The global lowest FoS value descreased in scenario II (from a value of over 2.0 to less than 1.75 on the left wall) (Figures 6b and 9b) and scenario III (from a value of over 2.5 to less than 2.25 on the right wall) (Figures 6b and 9c). In scenario I, the global lowest FoS value was the same (over 2.0) (Figures 6b and 9a). A summary of the FoS value and the size of the slope failure surface by 2025 is shown in Table 4. FoS values also confirm that scenario I is a better option than scenarios II and III in terms of the impact of UG on slope stability.

Change of FoS value and size of the slope failure surface with different scenarios of UG operation by 2025.

Left slope wall | 1.75–2.0 | 1.75–2.0 |
1.5–1.75 |
1.75–2.0 |

Right slope wall | >2.5 | >2.5 |
>2.5 |
2.0–2.25 |

In scenario IV, the impact of UG operation could be eliminated with time by the overburden extraction of the left side of the slope in the next stage (2038), and the right side part could be independently extracted while the pit slope above ends its operation. In this case, scenario IV is considered to be the better extraction option, rather than scenario V.

In general, longwall subsidence increases the size of slope failure surface inside the slope body in all calculation scenarios. For an OP–UG transition case study at the Cao Son mine, the calculated values of FoS were reduced by approximately 10%, depending on the location of the longwall panel. Despite reduction in the FoS values, the 2025 pit slopes are considered stable due to high FoS values of the planned slopes (much greater than 1.5). In case of slopes with low FoS values (e.g., less than 1.3), such a reduction of FoS values (10%) could have a significant impact on slope stability or even lead to loss of stability. Increase in the length of longwall panel enlarges the longwall subsidence. As a consequence, the size of slope failure surface increases and the FoS values tendentiously reduce. In a simultaneous operation, the impact of longwall subsidence on the previous slope phase can be reduced in the next slope phase by modifying the slope geometry. Geomechnical assessments should be carried out on an ongoing basis for each phase advance of both OP and UG operations. Such an action aims to indetify the possible slope failures and minimize the impact of longwall mining on slope stability.

The obtained results confirm the impact of UG operation on slope stability. UG extraction decreases the slope FoS and increases the size of slope failure surface inside slope body. The grade of impact is dependent on the UG extraction arrangements (extraction location and length). The UG-induced deformation occurring on the slope surface can negatively impact on OP operations, for example, transport or mining equipment. It also may initiate the slope failure surface and lead to pit slope failure with further UG advance. In the case of the Cao Son and the Khe Cham II–IV operations, the recommended option for UG operation in coal seam No. 9 is to start the UG operations from the left side (scenario I) and expand it to the other side (scenario IV). This operation arrangement minimizes the UG-induced impact on the pit slopes at the 2025 stage and 2030 stage.

The impact of UG operation on the OP slope is a continuous process related to the advances of both the OP and UG exploitations. During simultaneous operations, such as Cao Son and Khe Cham II–IV, changes in slope geometry can also reduce the impact of UG exploitation. Hence, it requires an appropriate schedule and the coordination of both OP and UG extraction.

From field practices, some suggestions are drawn: (1) applying full/partly hydraulic backfill as a method of ground control can reduce the impact of UG on slope stability and (2) applying slope regrading aims to prevent slope landslide if needed.

The analysis conducted provides the preliminary assessments of impact of UG on slope stability based on actual geo-mining conditions at the Cao Son OP mine. While these results are valuable, monitoring programs are recommended in order to control uncertainty and risk in the long term.

Further investigations should focus on the development of 3D geological and topographic models for 3D analysis, which would allow the direction and location of longwall panels in single coal seams and multi-coal seams to be taken into consideration. The outcomes would deliver better understanding and accuracy of OP–UG interaction.

#### Change of FoS value and size of the slope failure surface with different scenarios of UG operation by 2025.

Left slope wall | 1.75–2.0 | 1.75–2.0 |
1.5–1.75 |
1.75–2.0 |

Right slope wall | >2.5 | >2.5 |
>2.5 |
2.0–2.25 |

#### Mechanical properties of rock mass in the studied region [33].

^{t} |
^{3}) |
|||||
---|---|---|---|---|---|---|

Mudstone | 2.31 | 1.54 | 1.87 | 30 | 1.34 | 2620 |

Claystone | 2.23 | 1.34 | 2.05 | 26 | 1.14 | 2600 |

Anthracite | 2.17 | 1.36 | 2.14 | 27 | 1.22 | 1500 |

Conglomerate | 4.76 | 3.57 | 3.23 | 28 | 2.27 | 2510 |

Sandstone | 3.91 | 2.46 | 3.56 | 28 | 1.96 | 2600 |

#### Strength parameters of the fault in the studied region [33].

Value | 12–20 | 3.4–4.4 | 0 |

#### Hypotheses for calculating thickness of the caved zone.

Peng and Chiang, 1984 [ |
(2–10) |

Bai et al., 1995 [ |
100_{1}_{2}) |

Mazurkiewicz et al., 1997 [ |
_{r} |

Heasley, 2004 [ |
(10–18) |

Biliński, 2005 (simplified) [ |
(_{s}t_{c}^{0.5}+0.02) |

Wang et al., 2017 [ |
(3–4) |

Bearing Capacity Evaluation of Shallow Foundations on Stabilized Layered Soil using ABAQUS Evaluating the Effect of Environment Acidity on Stabilized Expansive Clay Overstrength and ductility factors of XBF structures with pinned and fixed supports Comparative Analysis of Single Pile with Embedded Beam Row and Volume Pile Modeling under Seismic Load Comment On Energy-Efficient Alternative for Different Types of Traditional Soil Binders Vertical and horizontal dynamic response of suction caisson foundations