Numerical simulation study on safety thickness of boundary pillar and mining-induced stability during open-pit to underground transition

2.1. Hydrogeological conditions

The regional hydrogeological survey area is characterized by denudation-accumulation topography within a low mountainous hilly region. The terrain slopes from north to south, located at the northern margin of the Taozao Basin hydrogeological unit, within the recharge area of a minor watershed. The southern low-lying area exhibits progressively decreasing southward gradients, forming a NW-SE trending hilly zone. The northern regional elevation peak (+324.90 m) is situated at Ximashan, while the lowest point (+70.00 m) occurs near Mengzhuang. The southern Quaternary coverage area gradually descends southward, with the Shiliquan spring group located 11 km from the mining area and the lowest discharge point at +55.00 m, representing the regional hydrogeological base level.

The drainage system is moderately developed, with four seasonal tributaries of the Guoji River (a branch of the Dasha River) distributed east-west across the mining area. These ephemeral streams remain dry in spring but convey water during rainy seasons, eventually meandering southward into the Dasha River.

The primary water-filling source for underground mining is directly recharged by bedrock fracture aquifer formations. The orebody’s roof and floor strata consist of biotite granulite, which is characterized by high rock strength, relatively intact structure, poorly developed fractures, low water-bearing capacity, and weak hydraulic conductivity. These properties ensure minimal hydrogeological impacts on deep mine workings and water inflow.

2.2. Engineering geological conditions

The orebody occurs within metamorphic strata of the Taishan Group Shancaoyu Formation. The roof and floor strata are primarily composed of biotite granulite and quartz-hornblende schist, maintaining conformable contact with the orebody.

Biotite granulite, gray-black in color, is distributed throughout the mining area, exhibiting fine-grained xenoblastic texture and banded structure. Surface exposures are intensely weathered, with most weathered materials eroded. The weathering depth ranges 0.00-49.50 m: strongly weathered zones display granular disintegration and low strength (easily crumbled by hand), while moderately weathered rocks occur as blocky fragments with higher strength but discontinuous structure. The stratum thickness decreases east-westward and thickens centrally. Shallow rock masses are fractured with developed joints and low mechanical strength, indicating poor engineering geological conditions. Deep fresh bedrock demonstrates relative integrity, high strength, and limited joint development (joint frequency: 2-3 joints/m), exhibiting weak water-bearing capacity. The rock exhibits favorable physico-mechanical properties, with saturated compressive strength of 82.7-157.3 MPa and shear strength of 7.97-12.5 MPa, classifying as hard rock.

Quartz-hornblende schist, dark green and unevenly distributed, locally occurs in orebody roof/floor strata with 1-2 m thickness. It displays granoblastic texture, schistose structure, and pronounced foliation. This rock type exhibits poor compressive/shear strength but similar hydro-mechanical behavior to biotite granulite.

The steeply dipping stratified/semi-stratified orebody comprises magnetite-hornblende quartzite, magnetite-quartz hornblendite, and biotite-hornblende schist. These ore-bearing rocks exhibit granoblastic or schistose textures with striated, banded, or foliated structures. Their hydro-mechanical properties mirror surrounding strata, showing weak water-bearing capacity. The discontinuous orebody, influenced by faulting, maintains relative integrity with undeveloped fractures and high mechanical strength (saturated uniaxial compressive strength: 263.7 MPa; shear strength: 15.8 MPa), classified as extremely hard rock.

2.3. Rock mechanics parameters

Specimens for rock mechanics testing were prepared from field-collected rock blocks. Standard specimens were manufactured from three rock types: surface waste rock, surrounding rock, and orebody, all maintained in their natural moisture state. Tested parameters included unit weight, elastic modulus, Poisson’s ratio, compressive strength, tensile strength, cohesion, and internal friction angle.

Given scale discrepancies between laboratory specimens and in situ rock masses, rock mechanics parameters were converted to engineering-scale values through empirical methods including RQD classification, CSIR and Q-system classifications, GSI classification, and geotechnical code-based approaches [20]. Comprehensive reduction factors were applied to derive underground rock mass parameters. Backfill mechanical parameters were obtained through backfill testing, while fault mechanical parameters were determined via laboratory experiments, as detailed in Table 1.

3.1. Open-pit status

Two surface pits (No. 1 and No. 2 open-pits) are present. No. 1 open-pit has been backfilled using surrounding waste rock piles, tailings sand generated during mining operations, and tailings from the tailings pond. After backfilling, the pit mouth was leveled with crushed stone soil. A layer of cultivated soil was subsequently applied to the leveled surface. While No. 1 open-pit has been fully rehabilitated, No. 2 open-pit remains unbackfilled (Fig. 1).

Fig. 1Current situation of unfilled No. 2 open-pit (Wang Feifei, taking photos in Shandong, China on July 15, 2025)

3.2. Underground mining status

The mine is divided into three sectors according to the design plan: Eastern Mining Area (Line –56 to –70, strike length ~1,400 m). Western Mining Area (Line –70 to –106, strike length ~3,600 m). Western Periphery Mining Area (Line –106 to –123, strike length ~1,700 m). Construction in the Eastern Mining Area has been suspended due to incomplete relocation of Hengshanqian Village on the surface. Only one vertical shaft has been completed in this sector. The Western Mining Area and Western Periphery Mining Area are currently in normal production following safety facility completion acceptance.

The underground workings are divided into five levels from top to bottom: ±0 m, -55 m, –110 m, –165 m, and –220 m, with each level interval measuring 55 m. The -220 m level is currently the active production horizon. The IV-2 and III-1 orebodies are being extracted using two methods: Shallow-Hole Shrinkage Stoping with Subsequent Backfill. VCR (Vertical Crater Retreat) Method with Subsequent Backfill. Six production stopes (762, 772, 782, 852, 902, and A2) are operational at the -220 m level. Among these, stopes 762, 772, 782, 852, and A2 employ shallow-hole shrinkage stoping with subsequent backfill, while stope 902 utilizes the VCR method with subsequent backfill.

4.1. Fine modeling method for multi-source data coupling

The mining area is characterized by denudation-accumulation topography within a hilly region, with higher elevations in the north and gentle slopes in the south. The ground elevation ranges from +77.50 m to +271.00 m, creating a relative relief of 193.50 m. The open-pit area is primarily situated on plains surrounded by farmland. Based on mine-provided drawings, a surface model was established in MIDAS-GTS/NX software. Three-dimensional solid models of surface open pits, ore bodies, stopes, and pillars were created through extrusion commands. The established surface model was extended into a solid model, while previously generated void entities were displayed in the software. The model origin and boundary dimensions were defined according to computational spatial requirements. Mining zone entities were embedded into the global model, forming internal mining districts as shown in Fig. 2. The model measures 3,130 m in length, 1,450 m in width, and 595 m in height.

Fig. 2Process of establishing the overall model of the research area (unit: m): a) surface topographic map; b) surface model established in software; c) design drawing for underground mining; d) underground mining ore body model; e) three-dimensional numerical model

Stope dimensions (50 m width × orebody thickness) were modeled at 1:1 scale using the mine’s latest design drawings. The model incorporates surface open pits, ore bodies, and various rock masses, totaling 233 entities. Materials including surface waste rock, biotite granulite, ore-rock, and tailings backfill were included. Surface soil layers were not separately modeled due to minimal thickness and negligible impact.

Three-dimensional mechanical analysis required mesh generation. Given the model’s large scale, non-critical regions were assigned coarse meshes (20 m density) to optimize computational efficiency, while critical components (ore bodies, backfill, crown pillars) received refined meshes (2 m density with 1 m node spacing). The final mesh contained 691,107 nodes and 593,762 elements.

The global model was exported from MIDAS-GTS/NX as a neutral file containing node/element information [21]. Through the MIDAS-FLAC converter, this data was transformed into FLAC3D-compatible *.f3grid format. The FLAC3D model (Fig. 3) was generated via the import command in the Grid module.

Fig. 3FLAC3D computational model

Boundary conditions were applied as follows: $X$-direction constraints on all nodes in $X$-plane (sliding hinges). $Y$-direction constraints on all nodes in $Y$-plane (sliding hinges). Full fixation ($X$-$Y$-$Z$) on nodes in the lowest $Z$-plane (fixed hinges). Initial geostatic stress calculation considered gravitational effects exclusively, given the maximum mining depth (~220 m) and simple geological structure [22-24].

Fig. 4FLAC3D calculation model stress balance process

A mesh sensitivity analysis was conducted to ensure the numerical results were independent of element size (Fig. 4). Three global mesh densities with average element sizes of 25 m (coarse), 10 m (medium), and 5 m (fine) were tested, with proportional refinement in the critical regions (ore body, stopes, and boundary pillar). The key output parameters, namely the maximum vertical displacement at the crest of Open-Pit No. 2 and the maximum principal stress within the boundary pillar, were monitored. The results showed that the change in these values was less than 2 % when the mesh was refined from the 5 m (fine) configuration to an even finer 3 m configuration. Therefore, the 5 m global mesh (with 2 m local refinement in critical zones) was selected for all subsequent simulations, as it provided a converged solution with optimal computational efficiency. The final mesh contained 691,107 nodes and 593,762 elements.

4.2. Simulated mining scheme

Based on the content and objectives of this study, the main tasks include: 3D modeling and numerical simulation of the original open-pit slope mining and backfilling; 3D modeling and numerical simulation analysis of the isolation top pillar under the current mining conditions; stability analysis of different thickness schemes for the isolation top pillar; and determination of reasonable parameters for the isolation top pillar, underground mining sequence, and stope structure. To achieve the above research objectives, a 3D numerical simulation method will be employed to simulate the mining methods and sequence, reproducing the geo-stress and surrounding rock deformation of the mine. Based on the numerical simulation results, the designed mining plan will be evaluated.

The numerical simulation plan is as follows: Step 1: Geo-stress equilibrium → Step 2: Mining No. 1 open-pit → Step 3: Mining No. 2 open-pit → Step 4: Backfilling No. 1 open-pit → Step 5: Underground stope mining to the current status. Given the isolation layer thickness of 45 m designed for the mine, the mining process of the ore body in the underground stope will be simulated to investigate the deformation, displacement, and stress of the isolation layer, thereby determining whether the thickness of the isolation layer is reasonable.

The selection of this specific sequence was based on two primary criteria: (1) Historical Accuracy: to faithfully replicate the actual chronological order of excavation and remediation activities at the mine, ensuring the stress evolution in the model mirrors the field conditions. (2) Mechanical Criticality: this sequence creates a conservative scenario where the boundary pillar is first unloaded by open-pit excavation and then loaded by subsequent underground mining, allowing us to evaluate its stability under the most demanding conditions. The mining steps simulate the bulk extraction of entire stopes or large sections, a standard approach in elastoplastic models to analyze the final equilibrium state after each major phase of excavation. The model's equilibrium was confirmed after each step by a significant reduction in the maximum unbalanced force.

Given the isolation layer thickness of 45 m designed for the mine, the mining process of the ore body in the underground stope was simulated to investigate the deformation, displacement, and stress of the isolation layer, thereby determining whether the thickness of the isolation layer is reasonable. The impact of this key parameter (45 m thickness) on overall mining stability is presented and discussed in detail in Section 4 (Analysis of results), focusing on its effectiveness in controlling displacement, stress redistribution, and plastic zone development.

4.3. Stability assessment criteria

According to the allowable values for displacement and deformation required by different protection levels specified in relevant regulations [25-26], the criteria for determining the impact of underground mining on surface stability in mining areas are deformation $\epsilon =$ ± 2 mm/m, inclination $i=$ ± 3 mm/m, and curvature $k=$ ± 0.2×10^-3/m.

The criterion for determining the deformation and displacement of underground rock formations based on the analogy of underground mining experience is that displacements below 20 mm have little effect on the stability of the rock mass. A displacement of 20-50 mm can maintain the stability of the rock mass. Displacement on the order of 50-100 mm poses potential stability issues for the rock mass 100 mm belongs to the problem of large deformation displacement, and there is a phenomenon of rock mass failure, which may also occur on a large scale. In addition to using deformation displacement as a criterion, comprehensive judgment is assisted by stress magnitude and plastic zone distribution [27-28]. This numerical simulation used indicators such as maximum unbalanced force, rock stress, surface displacement, and plastic zone distribution to determine mining stability.

These quantitative criteria provide the benchmark for our stability assessment. In the following sections, the results from the numerical simulation – including the magnitude and distribution of displacements, the values of induced stresses compared to rock mass strengths, and the extent of plastic zones – will be rigorously evaluated against these standards to conclusively determine the stability of the open-pit slopes, underground stopes, and the boundary pillar.

5.1. Stability of existing surface open pits

The mining process was simulated numerically to reproduce open-pit and underground extraction, yielding deformation displacements, rock stress distributions, and plastic zones (Fig. 5).

Fig. 5Stability calculation results of open pit mining site status

a) Deformation displacement

b) Rock stress

c) Distribution of plastic zone

Displacement nephograms indicate minimal deformation at No. 1 open-pit due to backfilling, with maximum displacements (~1 cm) localized in backfilled waste rock areas, attributed to compaction. No. 2 open-pit slopes exhibited displacements of 1.0-1.2 cm, below the 20 mm threshold for significant stability impacts. This confirms No. 2 open-pit slopes remain stable under current conditions, with no large-scale landslides or collapses observed.

Stress nephograms reveal compressive stress gradients increasing with depth. Stress redistribution and concentration occurred near slope toes and adjacent stope boundaries. Maximum compressive stresses in open-pit areas (1.5-3 MPa) and underground stopes (2-4.5 MPa) remain well below the compressive strengths of biotite granulite (9.6 MPa) and orebody (16.8 MPa), indicating low compressive failure risks. Tensile stresses (0.24-0.34 MPa) localized at slope faces and stope roofs are substantially lower than the tensile strengths of biotite granulite (0.64 MPa) and orebody (1.12 MPa), suggesting negligible tensile failure potential. Combined principal stress analysis confirms No. 2 open-pit slopes remain stable without underground-induced collapse risks.

Plastic zone distributions show no development at No. 2 open-pit slopes. Localized plastic zones in No. 1 open-pit backfill result from granular material compression. Underground backfill exhibits isolated plastic zones without inter-stope connectivity, demonstrating minimal upward stability impacts.

5.2. Stability of existing underground stop

Underground mining processes were numerically reconstructed to assess deformation, stress, and plastic zones in overlying strata (Fig. 6).

Maximum roof displacements of 7-9 mm was recorded, below the 20 mm stability threshold. This confirms intact rock mass stability without large-scale collapses. Compressive stress distributions follow depth-dependent gradients, with maximum values (2-4.5 MPa) in underground stopes remaining below rock strength thresholds. Tensile stresses (0.24-0.34 MPa) at stope roofs are subcritical relative to material tensile strengths. Principal stress integration verifies current underground stability. Plastic zones are confined to backfill materials without propagation to adjacent stopes or surface operations, confirming effective stability control.

Fig. 6Stability calculation results of underground mining status

a) Deformation displacement

b) Rock stress

c) Distribution of plastic zone

5.3. Stability of boundary pillar

Three-dimensional simulations evaluated the designed 45m boundary pillar during phased mining: geostatic equilibrium → No. 1 open-pit extraction → No. 2 open-pit extraction → No. 1 open-pit backfilling → underground mining to current status → full underground depletion. The spatial relationship between the isolation layer and the open-pit and underground mining areas is shown in Fig. 7. Four monitoring points in the boundary pillar are shown in Fig. 8.

Fig. 7Spatial relationship between the isolation layer and the open-pit and underground mining areas

Fig. 8Four monitoring points in the boundary pillar

Displacement monitoring (Fig. 9) shows maximum displacements of 1.0 cm at eastern No. 2 open-pit (M1/M2), with other points below critical thresholds. The 45 m boundary pillar effectively decouples underground disturbances from surface operations, maintaining displacements under 20 mm.

The minimal deformations and the absence of failure mechanisms within the 45 m boundary pillar have significant practical engineering implications. Firstly, they validate the current design and provide confidence that underground mining activities will not induce catastrophic collapse or landslides into the upper open-pit, ensuring the safety of surface infrastructure and equipment. Secondly, the model demonstrates that the pillar successfully acts as a stress-cap, transferring mining-induced loads away from the pit bottom and preventing the development of a continuous failure path to the surface. This effectively mitigates the risk of one of the most critical hazards in transition mining: sudden and unpredictable subsidence.

Fig. 9Displacement of monitoring points in the model

a) Monitoring point M1

b) Monitoring point M2

c) Monitoring point M3

d) Monitoring point M4

5.4. Model validation

A crucial step in numerical simulation is to validate the model’s accuracy against real-world field measurements. To ensure the reliability of the constructed model in predicting the mechanical behavior of the mining area, a validation process was carried out by comparing the simulation results with on-site monitoring data.

Field monitoring included GNSS surface displacement stations installed around the perimeter of the No. 2 open-pit and inclinometers for measuring subsurface movement. The simulated displacements after the excavation of the open-pit were extracted at locations precisely corresponding to the monitoring points.

The comparison revealed a strong correlation between the simulated and measured data. For instance, the maximum simulated horizontal displacement at the slope crest was approximately 1.1 cm, while the monitored value averaged 1.15 cm. The spatial distribution of the displacement field from the simulation also matched the trends observed in the monitoring data. The minor discrepancies (within an acceptable margin of error of ~5 %) can be attributed to the inherent heterogeneity of the rock mass and simplifications in the model. The close agreement confirms that the calibrated numerical model accurately represents the geotechnical conditions and can be reliably used for subsequent analysis of the isolation pillar's stability.