3D Modelling approach to identify parametric configurations for pillar stability in underground metal mine: a case study

Abstract

In underground metal mining, ensuring the stability of pillars considered as a principal apprehension, given the likelihood of stress increase and relief leading to dynamic pillar failure. Through a detailed analysis examining stress concentration factor (SCF), deformation and yielding surrounding crown and rib pillars, this study explores the impact of various rock mass parameters. A detailed investigation utilizing non-linear 3D numerical models, governed by Mohr–Coulomb elastoplastic failure criteria, is conducted by manipulating seven variables: working depth (D), uniaxial compressive strength (σ_ci), modulus of elasticity (E_m), geological strength index (GSI), crown pillar thickness (T_C), rib pillar thickness (T_R), and stope width (W_S). The findings suggest that increasing D from 230 m level to 290 m level exhibits a 6.62% decrease in major SCF, a substantial 26.37% decrease in minor SCF, and a noteworthy 119.9% increases in yield around crown pillars. Also, at lower depths, the rib pillar experiences a 243.31% increase in yield zone. Present study also exhibits possible failure of pillars in terms of factor of safety (FoS). This study offers crucial insight essential for optimizing pillar design and improving safety of the mine which further benefits industry financially by altering additional operational cost.

Keywords:

• Metal mine
• sublevel stoping
• 3D modelling
• Mohr–Coulomb elastoplastic failure criteria
• pillar stability

Introduction

In underground (u/g) hard rock mines, pillars serve a critical function in maintaining the structural integrity of the overlying strata and ensuring a safe, stable working environment for mining activities (Ghasemi et al. 2017). However, the formation of extensive stopes necessitates careful consideration of geotechnical stability due to stress concentration around natural pillars, particularly crucial for safety and productivity in large-scale u/g metal mining operations. The creation of vast open stopes through excavation leads to stress concentration around natural pillars, leading to significant convergence, lateral disintegration, and increased regime of yielding surrounding the corners of pillars and u/g structures (Zhang and Mitri 2008) potentially leading to the collapse of the stope surface (Urli and Esmaieli 2016). Crown pillar often suffers cracks or damage from the immense load of overlying strata and intense horizontal stresses. Thus, prioritizing geotechnical stability stands as a fundamental aspect for safeguarding and enhancing productivity particularly in sizable u/g mines. Rib pillars serve a critical purpose in open stope mining by improving the stability of an excavation zone. Their function is to mitigate rock mass displacement and confine its presence primarily along stope backs and walls. Moreover, these pillars enable greater adaptability executing ore extraction through several operational faces in partitioned stope blocks. The loss of load-bearing capacity in rib pillars linked with the loss of confinement. The outer rock layers provide essential confinement to the central core of the pillar. Initial signs of pillar deterioration include boundary fractures, releasing some of the original confinement. As load increases, the pillar’s lateral expansion grows, leading to potential slabbing and propagation of fractures into the pillar core, diminishing its strength (Potvin et al. 1989). Severe fracturing is typically necessary before block fallout and disintegration occur. These observations align with previous literature on room and pillar mining (Krauland and Soder 1987), highlighting similarities in response mechanisms between biaxially and uniaxially loaded pillars (Brady and Brown 2006).

Mining operations subjected to extensive depths; extensive u/g stopes concentrate stress, leading to significant deformation and yielding around pillar corners. This process induces stress relaxation in the surrounding areas, posing potential risks such as deterioration of the stope front (Urli and Esmaieli 2016). The stability of highly stressed rock masses relies on several factors, particularly the structural and mechanical attributes of the rocks. Pillar stability is notably influenced by factors such as depth of cover, excavation size, horizontal stresses, rock mass properties, backfilling, reinforcement, and the orientation of the orebody. These factors are extensively studied and documented by various researchers (Kvapil and Blake 1973; Tavakoli 1994; Edelbro 2003; Lu et al. 2008; Deb and Das 2011; Esterhuizen et al. 2011). (Chen and Mitri 2021) conducted a parametric study to analyse the influence of fault characteristics and ratio of horizontal stresses on crown pillar stability whereas (Saeidi et al. 2021) demonstrated the effect of failure under compression, tension and brittleness for the rib pillar using FLAC3D 7.0. While a thick crown pillar contributes to hanging wall support and overall stope stability, its thickness may pose challenges in terms of mineral conservation. Hence, optimizing pillar dimensions holds significant importance for metalliferous mines (Bakhtavar et al. 2011) as it affects the ore dilution around the mining stopes (Guggari et al. 2024). Determining the ideal thickness for a crown pillar is intricate and typically relies on applied knowledge along with array of empirical strategies (Tavakoli 1994; Mark 2016). The reliability of pillars within u/g mine relies on the structural integrity of adjacent rock and induced stresses, influenced by excavation shape, and pre-existed in-situ stresses, which are themselves shaped by the geological background of the rock formation (Martin et al. 2003). (Li et al. 2021a; 2021b) considered Geological Strength Index (GSI) as factor of influence for the stability combining it with uniaxial compressive strength (σ_ci) and pillar width and suggested the influence of all three parameters on the pillar stability. (Heidarzadeh et al. 2020) varied geo-mechanical attributes of rock mass for assessing the effect of these parameters on the pillar stability using numerical techniques associated with deterministic model. Thus, prioritizing geotechnical integrity is crucial, especially in managing pillar integrity and optimizing their dimensions to enhance safety and operational efficiency. The assessment of u/g structure stability stands as a fundamental aspect within the domain of u/g mining. Evaluation approaches often use analytical, empirical, and numerical methodologies. However, it is acknowledged that analytical assessments, constrained by basic geometric considerations, might lack the comprehensive scope required to analyse complex geo-mining scenarios (Brady and Lorig 1988).

A parametric study by Guggari et al. (2023) was conducted to address the design of a crown pillar for u/g deposits, considering factors such as characteristics of the rock mass, excavation span and depth. This study considered the depth of working (D), uniaxial compressive strength (σ_ci), pillar thickness (T), geological strength index (GSI), and modulus of elasticity (E_m) as major parameters for the 2 dimensional (2D) numerical analysis. Its objective was to comprehend the pillar dynamics within the region of stress relief resulting from sub-level open stope excavation. 2D finite element (FE) numerical model was prepared to analyse the effect of stresses, deformation and yielding on the gold bearing u/g metal mine. The major drawback of the study is that nonhomogeneous rock mass analysed using 2D numerical analysis underpredicts deformation as the weak zones present in the rock mass structure has prominent effect on the deformation when subjected to the 3D analysis (Dhawan et al. 2002).

To overcome the limitation of the above study a detailed study is conducted considering additional geo-mechanical parameters to access the influence over the u/g structure, 3D finite element models were created using Ansys Workbench version 19.0. These models were selected for their ease of setup and relatively short computational times, enabling evaluation of sensitivity and refinement (Norén-Cosgriff et al. 2020). Nonlinear materials following Mohr–Coulomb (MC) failure criteria (Xu et al. 2006; 2015) were utilized, chosen due to its compatibility with modelling rock response, accounting for rock brittleness and vulnerability along inherent fault lines. This criterion establishes a failure envelope in the domain of stress via two key indicators -cohesion (c) and internal friction angle (ϕ). Beyond this envelope, failure happens adjacent to a plane oriented aligned with the friction angle. MC criteria aligns well with modelling rock response, particularly accounting associated with the recurring fracturing and collapse observed along pre-existing discontinuities within the rock (Labuz et al. 2012). The outcomes are evaluated in terms of the stress concentration factor (SCF), deformation, yielding zones, and the factor of safety (FoS). Each parameter’s influence on pillar stability is assessed by varying factors such as working depth (D), uniaxial compressive strength (σ_ci), crown pillar thickness (T_C), rib pillar thickness (T_R), geological strength index (GSI), stope width (W_S), and modulus of elasticity (Em).

Overview of case study mine

The general strike of the mine formation extends along N30°E-S30°W, characterized by steep dips of approximately 70° toward the northwest. At the 5000 latitudes of the Mine, there is an eastward shift in the strike, transitioning from N52°E-S52°W to N68°E-S68°W. Mineralization primarily occurs along the gradual contact of the Delhi System, with all ore lenses running parallel or sub-parallel to the host rocks. The surface geological map delineates several transverse faults, notably evident on the footwall side within the feldspathic quartzite horizon. These faults exhibit a clear displacement of thin magnetite bands and dip at around 70° towards the north, occasionally towards the south. Although these transverse faults are mapped, they minimally impact the continuity of ore lenses. Any minor displacements observed do not significantly alter the width of the mineralized zone intersected. Consequently, the ore zones maintain continuity in their strike after traversing these faults.

In the mine area, surface exposures are absent, necessitating an exploration of the wall rocks and host rocks solely through drill core analysis. The examination of borehole cores reveals that both the hanging wall and footwall rocks surrounding the mineralization exhibit solid and competent characteristics, serving as sturdy barriers due to their sub-vertical orientation. However, the mineralized ground displays lower competence due to the presence of joints, shears, and fractures. The physical and chemical attributes of the wall rocks align favourably with those conducive to u/g mining, as observed in neighbouring mine areas. Furthermore, negligible effects of weathering have been observed. The mine spans a latitude range from 5000 N to 6200 N. A decline with a gradient of 1 in 9 serves as the entry point for mineral extraction from the mine. The copper orebody, predominantly chalcopyrite, manifests in narrow lodes of varying widths within the mine. Operations occur at intervals of 60 meters (m) between adjacent main levels. Currently, pillar extraction activities are underway at 240 m level, 180 m level, and 120 m level, with ongoing development work planned for subsequent levels (60 m level and 0 m level). Mineral extraction primarily takes place from transverse and longitudinal stopes using the 'Blast Hole Open Stoping’ method.

Numerical modelling

Geotechnical characterization

Rock samples were gathered from different levels within the mine, ranging in volume from 300 m level to 60 m level, along with samples from the host rock below the mine. Testing included 49 core samples for uniaxial compressive strength (σ_ci) and tri-axial compression, and another 49 for biaxial tensile strength. The core samples varied in length-to-diameter (L/D) ratios (2.5, 2, and 0.5, respectively) with a diameter of 54 mm. The orebody exhibited σ_ci values between 109.96 and 118.87 MPa, E_m from 29.60 GPa to 44.98 GPa, GSI ranging from 45 to 66.8, ϕ from 31.77° to 39.81°, and c between 1.04 and 3.38 MPa. These rock properties were considered for modelling, particularly in relation to variables like T_C, T_R, W_S, D, σ_ci, E_m, and GSI. Variations are considered are in range of these values which are represented in Table 1.

Table 1. Key geotechnical factors dataset for numerical simulation.

Properties	Orebody	Host rock	Variations in values
Working depth (D), m Level	–	–	230, 290, 350
UCS (σci), MPa	✓	✗	100, 120
Modulus of elasticity (Em), GPa	✓	✗	30, 37.5, 45
Geological Strength Index (GSI)	✓	✗	45, 55, 65
Crown pillar thickness (TC), m	✓	✗	8, 10
Rib pillar thickness (TR), m	✓	✗	10, 15
Stope width (WS), m	✓	✗	30, 40
Friction angle (Φ), °	✓	✗	39.81–31.77
Cohesion (c), MPa	✓	✗	1.04–3.38
Dilation angle (Ψ), °	✓	✗	9.87–7.94
UCS (σci’), MPa	✗	✓	122.76
Modulus of elasticity (Em’), GPa	✗	✓	30.65
Geological Strength Index (GSI')	✗	✓	60
Friction angle Φ', °	✗	✓	57.44
Cohesion (c‘), MPa	✗	✓	2.43
Dilation angle (Ψ'), °	✗	✓	14.36

Where ° is symbol for degree.

The geo-mechanical attributes in intact rocks were converted into rock mass properties through the implementation of RSData ver.1.007 within Rocscience software. Utilizing the Hoek–Brown criterion, rock mass properties were determined to facilitate numerical simulations. A total of 144 3D FE excavation models were prepared, having all possible iterations, in addition to 18 in-situ models.

Finite element analysis

All 3D models created for the numerical analysis represent the actual dip and orebody thickness variation across working depths. D is 450 m from the surface and the level are mentioned by 240 m level, 180 m level, 120 m level and 60 m level and 0 m level. The level difference is 60 m with a draw-point dimension of 4 m width and 3.2 m of height. For numerical analysis the dimensions of W_S, T_C, and T_R are varied along with the strength properties of rock mass. Figure 1(a,b) represents front view and isometric view of model geometry, respectively. A fine mesh of element size of 1.5 m is created near the excavation field and the element size was gradually increased with the growth rate of 1.2. Figure 2(a) represents the coarse meshing profile for the host rock whereas Figure 2(b) represents fine meshing for the orebody. For the mesh independency, mesh sensitivity analysis was performed (Table 2) and represented in Figure 3. The stress assessments were consistent with 54,74,165 nodes and 38,46,657 elements, reaching a limit where the element count becomes ineffective with respect to the result. A tetrahedral mesh is a four-node based method which uses two algorithms of patch conforming and patch independent method. This analysis incorporates a patch conforming method commencing from the geometry’s edges and faces, the meshing process gradually extends into its volume. The method generates a high-quality clean mesh which supports a better convergence of stress compared to the patch independent method that is more suitable for poor quality geometry.

Figure 1. (a) Front view of model geometry and (b) Isometric view of model geometry.

Figure 2. (a) Coarse meshing profile for the host rock and (b) Finer meshing scheme for the orebody.

Figure 3. Mesh sensitivity analysis for independency of mesh count on results.

Table 2. Percentage difference of stress vs element count.

Elements count	Major principal stress (MPa)	Percentage difference	Normal stress (MPa)	Percentage difference
24638	15.73	–	6.56	–
602452	17.36	9.37	7.24	9.28
1209468	18.14	4.27	7.57	4.42
2809356	18.27	0.71	7.63	0.74
4079068	18.31	0.23	7.65	0.27
4827036	18.31	0	7.65	0

Horizontal gradient pressure (the maximum σ_H and the minimum σ_h) are applied over the two vertical faces of the 3D FE model with stresses varying with depth from 11.58 MPa to 21.48 MPa for σ_H (Figure 4a) and from 8.29 MPa to 13.15 MPa for σ_h (Figure 4b), respectively. Figure 4(c) shows uniform pressure (σ_V) of 4.71 MPa is exerted onto the model’s uppermost face. Displacement is provided on the opposite faces of the geometry which is constrained in the respective direction whereas fixed support is provided in the bottom of the model geometry (Figure 4d).

Figure 4. (a) Maximum horizontal gradient pressure (σ_H), (b) Minimum horizontal gradient pressure (σ_h), (c) Vertical pressure (σ_V) applied on the model geometry, and (d) Loading and boundary condition on model geometry.

For the estimation of σ_H, σ_h, and σ_V the equation proposed by NGRI (National Geotechnical Research Institute) for the mine is given below: (1) $${\mathrm{\sigma }}_{\mathrm{H}}=\left(6.3+0.033\mathrm{Z}\right)$$(1) (2) $${\mathrm{\sigma }}_{\mathrm{h}}=\left(5.7+0.0162\mathrm{Z}\right)$$(2) (3) $${\mathrm{\sigma }}_{\mathrm{V}}=\text{ρ*g*Z}$$(3) where Z is the depth from the surface in metre (m).

FE models depict the post-stope state after ore recovery. Concurrently, in-situ models are generated to evaluate the impact of rock-mass stresses and deformation states. Calibration, following the methodology proposed by Shnorhokian et al. (2014), is conducted to calibrate these models. This calibration follows a comparative analysis between the output of numerical models and in-situ observations documented at specific locations, outlined in Table 3, before any excavation occurs. The working depth of 230 m level below the surface within the country rock used for the calibration process. For analysing ground behaviour through finite element modelling, several key assumptions have been identified and are listed below.

Table 3. Comparison table: FE model outputs and in-situ observations.

Rock type	In-situ measurements		FE model measurements
	σxx (MPa)	σyy (MPa)	σxx (MPa)	σyy (MPa)
Country rock	7.92	6.85	7.87	6.82

• Joint and bedding plane effects are accounted for within the GSI value for rock mass, rather than analysing the distinct impact of each plane separately.

• The model excludes considerations for dynamic loading.

• A zero value is adopted as the cutoff for tensile strength.

Results and discussion

Distribution of principal stresses around the crown and rib pillar

Prior to estimating the factor of safety (FoS) for the crown and rib pillars, it is imperative to ascertain the stress distribution resulting from stope extraction. An FoS below one signifies a risk of complete pillar failure. Given the transient nature of stopes, intended for subsequent sand stowing, a threshold value of 1.0 could be considered acceptable.

Individual discussions of each model’s outcomes are not within the paper’s focus. However, other input parameters remained consistent as T_C = 8 m, T_R = 10 m, W_S = 30 m, D = 230 m, GSI = 55, σ_ci = 100 MPa, and E_m = 30 GPa while focusing on evaluating the impact of specific input parameters on stress profile. Figure 5(a) illustrates the distribution of the major principal stress (σ₁), while Figure 5(b) showcases the distribution of minor principal stress (σ₃) surrounding the model geometry. Regions characterized by low stress levels (σ₃<0), denoting stress relaxation zones, were identified and confirmed in Figure 5(b) surrounding the stope.

Figure 5. (a) Major principal stress (σ₁) distribution and (b) Minor principal stress (σ₃) distribution.

SCF at a specific location is determined by comparing the stress induced by excavation to the in-situ stress at that location (Zhang et al. 2014). Whereas (Zhang and Mitri 2008) suggested that SCF is a ratio of post-mining major principal stresses to average of pre-miming principal stress value. It signifies alterations or redistribution in stresses concerning initial conditions. SCF near one indicates minimal change or redistribution of in-situ stress caused by excavation. Considering the SCF described above the expression for major and minor SCF (Kumar et al. 2016) is provided as follows (4) $${\kappa }_1=\frac{{\sigma }_1\left(\mathit{\text{induced}}\right)}{{\sigma }_1^{\prime }(\mathit{\text{in}}-\mathit{\text{situ}})}\mathrm{ }$$(4) (5) $${\kappa }_3=\frac{{\sigma }_{3\mathrm{ }}\left(\mathit{\text{induced}}\right)}{{\sigma }_3^{\prime }(\mathit{\text{in}}-\mathit{\text{situ}})}\mathrm{ }$$(5) where κ_1, κ₃ is major SCF and minor SCF, respectively; ${\sigma }_1,$ ${\sigma }_{3 }$ is the mining induced major and minor principal stress, respectively, whereas ${\sigma }_1^{\prime },$ ${\sigma }_3^{\prime }$ is the average in-situ mining major and minor principal stress, respectively.

The SCF serves as a guide in identifying where and how extensive high-stress zones form post-excavation. Consistency in SCF values across the mining opening suggests greater stability for that opening (Obert et al. 1960). When lacking the ability for non-linear modelling, the SCF criterion serves as a valuable tool within linear elastic analysis, highlighting potential failure points.

Effect of D on κ₁ and κ₃

As the host rock is more competent compared to the orebody, the centre of the crown pillar exhibits the most prominent stress concentration and along the skin of the side wall of the rib pillar. Figure 6(a,b) and Table 4 demonstrates the spread of σ₁ and σ₃ around the crown pillar and rib pillar for depths of 230 m level, 290 m level, and 350 m level, respectively.

Figure 6. (a) σ₁ and σ₃ Distribution with different working depth (D) around crown pillar and (a) σ₁ and σ₃ distribution with different working depth (D) around rib pillar.

Table 4. κ₁ and κ₃ Values at different D.

Stress distribution around crown pillar						Along the path A-A'
Depth	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(m level)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
230	26.35	21.26	22.92	1.51	–	10.27	1.33	4.89	0.73	–
290	29.36	21.19	22.53	1.41	6.62 ↓	9.40	1.98	4.29	0.57	26.37 ↓
350	27.48	13.26	21.98	1.29	8.51 ↓	9.53	2.56	5.02	0.53	20.95 ↓
Stress distribution around rib pillar						Along the path B-B'
Depth	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(m level)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
230	34.18	14.96	31.34	2.01	–	6.67	1.81	6.64	0.87	–
290	26.90	15.13	22.75	1.48	26.37 ↓	4.97	1.95	3.91	0.52	40.23 ↓
350	22.52	18.08	19.83	1.17	20.95 ↓	4.01	2.83	3.29	0.38	26.92 ↓

For the crown pillar, κ₁ values were obtained as 1.51, 1.41, and 1.29 at depths of 230 m level, 290 m level and 350 m level below surface level, respectively, exhibiting decreases of 6.62% and 8.51%. Correspondingly, κ₃ values at these depths were 0.73, 0.57, and 0.53, demonstrating decreases of 21.92% and 7.02% specifically at the lower D values.

For the rib pillar, derived κ₁ values at depths of 230 m level, 290 m level and 350 m level below the surface were 2.01, 1.48, and 1.17, respectively, indicating decreases of 26.37% and 20.95%. Simultaneously, the κ₃ values for the rib pillar at these depths were 0.87, 0.52, and 0.38, resulting a decrease of 40.23%, which further reduced to 26.92% at the respective lower D values.

The critical κ₁ threshold was identified at 2.21 for a specific iteration within the crown pillar, characterized by parameters: σ_ci = 120 MPa, E_m = 37.5 GPa, GSI = 65, T_C = 10 m, T_R = 10 m, and W_S = 40 m, observed at a depth of 290 m level. Conversely, the rib pillar exhibited an optimal κ₁ value of 2.31 under conditions: σ_ci = 120 MPa, E_m = 37.5 GPa, GSI = 45, T_C = 8 m, T_R = 10 m, and W_S = 40 m, positioned at a depth of 230 m level. In terms of κ₃, the crown pillar reached its maximum value at 1.12, achieved with parameters: σ_ci = 100 MPa, E_m = 37.5 GPa, GSI = 55, T_C = 10 m, T_R = 15 m, and W_S = 30 m, for a D at 350 m level. Additionally, the value of κ₃ amounted to 1.36 for σ_ci = 100 MPa, E_m = 30 GPa, GSI = 45, T_C = 8 m, T_R = 10 m, and W_S = 4 m, at a depth of 230 m level.

Effect of σ_ci on κ₁ and κ₃

Figure 7(a,b) and Table 5 depicts σ₁ and σ₃ dispersion surrounding the crown and rib pillars concerning variations in σ_ci values of 100 and 120 MPa. In the case of the crown pillar, the κ₁ value remained constant at 1.51 for varying σ_ci values, while the κ₃ value also maintained consistency at 0.73 across the respective σ_ci variations.

Figure 7. (a) σ₁ and σ₃ Distribution with σ_ci variations around crown pillar and (b) σ₁ and σ₃ distribution with σ_ci variations around rib pillar.

Table 5. κ₁ and κ₃ Values at different σ_ci.

Stress distribution around crown pillar						Along the path A-A'
σci	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(MPa)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
100	26.35	21.26	22.92	1.51	–	10.27	1.33	4.89	0.73	–
120	26.05	21.06	23.02	1.51	0	10.21	1.38	4.93	0.73	0
Stress distribution around rib pillar						Along the path B-B'
σci	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(MPa)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
100	34.18	14.96	31.34	2.01	–	6.67	1.81	6.64	0.87	–
120	35.08	15.87	31.68	2.03	0.96 ↑	10.83	2.85	6.40	0.83	4.6 ↓

Conversely, for the rib pillar, the κ₁ values were observed as 2.01 and 2.03 for σ_ci values of 100 and 120 MPa, respectively. The κ₃ values for the corresponding σ_ci were measured at 0.87 and 0.83.

The crown pillar exhibits consistent κ₁ and κ₃ values, showing no variation. Similarly, in the rib pillar, the difference remains within 5%, signifying no substantial alteration with varying σ_ci values of the orebody.

Effect of Em on κ₁ and κ₃

E_m defines the rock mass’s resilience against the strain resulting from external forces, significantly influencing the material’s elasto-plastic behaviour. In the crown pillar (Figure 8a), κ₁ values were measured at 1.51, 1.54, and 1.56 for E_m values of 30, 37.5, and 45 GPa, respectively, while κ₃ values for the crown pillar stood at 0.73, 0.75, and 0.73, respectively.

Figure 8. (a) σ₁ and σ₃ Distribution with E_m variations around crown pillar and (b) σ₁ and σ₃ distribution with E_m variations around rib pillar.

For the rib pillar (Figure 8b), κ₁ values were observed as 2.01, 1.95, and 1.88, while corresponding κ₃ values were recorded at 0.43, 0.44, and 0.46 for varying E_m values.

Overall, an increase in elasticity exhibits decreased σ₁ and σ₃ values (Table 6). The change in κ₁ and κ₃ values within crown pillar (thicknesses of 8 and 10 m) and rib pillar (thicknesses of 10 and 15 m), respectively remains limited within 4% when the E_m of the orebody varies from 30 to 45 GPa. However, the change in κ₃ for the rib pillar shows a 17.24% increase when increasing E_m from 30 to 37.5 GPa. From these findings, several conclusions emerge:

Table 6. κ₁ and κ₃ Values at different E_m.

Stress distribution around crown pillar						Along the path A-A'
Em	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(GPa)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
30	26.35	21.26	22.92	1.51	–	10.27	1.33	4.89	0.73	–
37.5	26.08	20.86	23.22	1.54	1.99 ↑	9.96	1.39	4.97	0.75	2.74 ↑
45	25.91	20.58	23.49	1.56	1.3 ↑	9.72	1.45	5.04	0.73	2.67 ↓
Stress distribution around rib pillar						Along the path B-B'
Em	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(GPa)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
30	34.18	14.96	31.34	2.01	–	6.67	1.81	6.64	0.87	–
37.5	26.90	15.13	22.75	1.48	26.37 ↓	4.97	1.95	3.91	0.52	40.23 ↓
45	22.52	18.08	19.83	1.17	20.95 ↓	4.01	2.83	3.29	0.38	26.92 ↓

• The alterations in κ₁ and κ₃ due to E_m variations are more pronounced in the rib pillar compared to the crown pillar.

• Higher magnitudes of κ₁ and κ₃ are observed around the centre of the crown pillar and along the edges of the rib pillar sidewall.

Effect of GSI on κ₁ and κ₃

GSI holds significance in determining the stability of u/g excavations. Variations in stress magnitude align with the variability in the orebody’s GSI, impacting the yielding zone significantly. In the context of the crown pillar (Figure 9a), κ₁ values exhibit a progression from 1.43 to 1.58 as the GSI increases from 45 to 65. Similarly, κ₃ values rise from 0.69 to 0.74 with the corresponding increase in GSI.

Figure 9. (a) σ₁ and σ₃ Distribution with GSI variations around crown pillar and (b) σ₁ and σ₃ distribution with GSI variations around rib pillar.

Conversely, for the rib pillar (Figure 9b), κ₁ values fluctuate between 1.88 and 1.87 with GSI changes from 45 to 65, while κ₃ values demonstrate a more pronounced shift from 1.08 to 0.74 under the same GSI variations (Table 7). Summarizing the findings:

Table 7. κ₁ and κ₃ Values at different GSI.

Stress distribution around crown pillar						Along the path A-A'
	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
GSI	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
45	26.99	19.14	22.23	1.43	5.3 ↓	11.19	1.10	4.48	0.69	5.48 ↓
55	26.35	21.26	22.92	1.51	–	10.27	1.33	4.89	0.73	–
65	25.24	20.18	23.87	1.58	4.64↑	9.6	1.61	5.11	0.74	1.37 ↑
Stress distribution around rib pillar						Along the path B-B'
	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
GSI	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
45	41.11	10.29	27.65	1.88	6.47 ↓	9.91	1.54	6.26	1.08	24.14 ↑
55	34.18	14.96	31.34	2.01	–	6.67	1.81	6.64	0.87	–
65	40.58	19.18	33.56	1.87	6.97 ↓	7.09	1.39	6.73	0.74	14.94 ↓

• An alteration in GSI from 55 to 45 is associated with a 5.3% decrease in κ₁ and a 5.48% reduction in κ₃, particularly within the crown pillar. However, raising the GSI to 65 results in a 4.64% increase in κ₃ without notable changes in κ₁.

• In the rib pillar, a decrease in GSI from 55 to 45 prompts a 6.47% drop in κ₁ and a substantial 24.14% increase in κ₃. However, shifting the GSI to 65 induces a 6.97% decrease in κ₁ and a 14.94% decline in κ₃.

Effect of T_C on κ₁ and κ₃

The depiction of σ₁ and σ₃ distribution on crown and rib pillars in Figure 10(a,b), respectively reflects the impact of T_C. As the pillar thickness increases, there is a corresponding rise in the extent and scale of the principal stress on the crown pillar. For instance, at 8 and 10 m of pillar thickness, the κ₁ values are 1.51 and 1.55, respectively. Meanwhile, the κ₃ values stand at 0.73 and 0.64 for their respective thicknesses (Table 8).

Figure 10. (a) σ₁ and σ₃ Distribution with T_C variations around crown pillar and (b) σ₁ and σ₃ distribution with T_C variations around rib pillar.

Table 8. κ₁ and κ₃ Values at different T_C.

Stress distribution around crown pillar						Along the path A-A'
TC	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(m)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
8	26.35	21.26	22.92	1.51	–	10.27	1.33	4.89	0.73	–
10	38.23	28.71	23.54	1.55	2.65 ↑	9.48	1.25	4.39	0.64	12.33 ↓
Stress distribution around rib pillar						Along the path B-B'
TC	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(m)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
8	34.18	14.96	31.34	2.01	–	6.67	1.81	6.64	0.87	–
10	33.19	14.04	30.17	2.0	0.5 ↓	6.20	1.59	5.74	0.74	14.94 ↓

The observation indicates that thicker pillars coincide with an increase in the σ₃ at the pillar’s corner, potentially leading to shear yielding. At depths of 230 m level, 290 m level and 350 m level, there is a progressive decrease in σ₁, suggesting that thicker pillars might be suitable for these depths. Summarizing the findings:

• Increasing T_C from 8 m to 10 m results in a 12.33% decrease in the κ₃ value.

• For rib pillars, there is a decrease of 0.5% in κ₁ and a 14.94% decrease in κ₃ when comparing values of 2.01 and 2.0 for κ₁, and 0.87 and 0.74 for κ₃, respectively.

• The impact of increasing T_C from 8 m to 10 m on the κ₁ value is not significant, with a change within 3%.

Effect of T_R on κ₁ and κ₃:

Figure 11(a,b) provides insight into the influence of T_R on the distribution of σ₁ and σ₃ around crown and rib pillars. The data as mentioned in Table 9 illustrates that as pillar thickness expands, there’s a simultaneous increase in the magnitude and extent of κ₁, accompanied by a decrease in minor stress values. Specifically, when T_R is increased from 10 to 15 m, the resulting κ₁ values of 1.51 and 1.52 indicate an insignificant 0.66% increment. In contrast, the respective κ₃ values of 0.73 and 0.65 exhibits a significant decrease of 10.96%.

Figure 11. (a) σ₁ and σ₃ Distribution with T_R variations around crown pillar and (b) σ₁ and σ₃ distribution with T_R variations around rib pillar.

Table 9. κ₁ and κ₃ Values at different T_R.

Stress distribution around crown pillar						Along the path A-A'
TR	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(m)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
10	26.35	21.26	22.92	1.51	–	10.27	1.33	4.89	0.73	–
15	37.26	16.17	23.12	1.52	0.66 ↑	9.28	2.14	4.5	0.65	10.96 ↓
Stress distribution around rib pillar						Along the path B-B'
TR	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(m)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
10	34.18	14.96	31.34	2.01	–	6.67	1.81	6.64	0.87	–
15	33.83	12.08	27.55	1.76	12.44 ↓	6.36	1.11	6.19	0.79	9.2 ↓

This investigation highlighted a distinct elevation in σ₃ at the pillar corner with increasing thickness, suggesting potential implications for shear yielding. Within rib pillars of 10 and 15 m thickness, κ₁ values decreased from 2.01 to 1.76, indicating a 12.44% decrease, while κ₃ values shifted from 0.683 to 0.743, representing a 9.2% decrease. These findings underscore the intricate stress behaviours observed across various values of T_R.

Effect of W_S on κ₁ and κ₃

While assessing the impact of widening W_S from 30 to 40 m on the crown pillar (Figure 12a), κ₁ values were measured at 1.51 and 1.47, indicating a 2.65% decrease (Table 10). Similarly, the κ₃ values recorded were 0.73 and 0.59, reflecting a decrease of 19.18% for the respective W_S values. However, in the case of the rib pillar (Figure 12b), κ₁ values were noted at 2.01 and 2.03 for the varied W_S, whereas κ₃ values stood at 0.87 and 0.80 for the respective W_S values.

Figure 12. (a) σ₁ and σ₃ Distribution with W_S variations around crown pillar and (b) σ₁ and σ₃ distribution with W_S variations around rib pillar.

Table 10. κ₁ and κ₃ Values at different W_S.

Stress distribution around crown pillar						Along the path A-A'
WS	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(m)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
30	26.35	21.26	22.92	1.51	–	10.27	1.33	4.89	0.73	–
40	30.11	18.76	22.33	1.47	2.65 ↓	10.76	2.41	4.12	0.59	19.18 ↓
Stress distribution around rib pillar						Along the path B-B'
WS	σ1 (MPa)			κ1	Percentage	σ3 (MPa)			κ3	Percentage
(m)	Max	Min	Avg.	Avg.	difference of κ1	Max	Min	Avg.	Avg.	difference of κ3
30	34.18	14.96	31.34	2.01	–	6.67	1.81	6.64	0.87	–
40	31.05	14.44	31.79	2.03	1 ↑	5.68	1.68	6.23	0.8	8.05 ↓

This analysis suggests that increasing the stope width from 30 to 40 m minimally impacts both κ₁ and κ₃ values for the rib pillar and κ₁ value of the crown pillar, with changes remaining within 3%. However, the κ₃ value of the crown pillar experiences a notable 19.18% decrease under these conditions.

Effect of various geo-mining properties on vertical deformation profile

Deformations exceeding specific thresholds within the crown pillar and hanging wall can instigate collapses. Examining deformation profiles along line A-A’ for the crown pillar and line B-B’ for rib pillars (Figures 13 and 14, respectively) involved comprehensive analysis across different geo-mining conditions. To understand the impact of excavations on vertical deformation, comparisons between in-situ deformations and excavated models along these lines were conducted.

Figure 13. Vertical deformation of various geo-mining properties for crown pillar.

Figure 14. Vertical deformation of various geo-mining properties for rib pillar.

Numerical analysis demonstrates a 3.24% reduction in deformation (Table 11) when D increases from 230 m level to 290 m level for the crown pillar, decreasing from 81.17 to 78.54 mm. Further deepening to 350 m level shows a marginal 0.36% reduction, reaching 78.21 mm. Conversely, the rib pillar consistently experiences decreasing deformation: from 80.08 mm at 230 m level to a 70.88 mm at 290 m level accounting a 0.36% decline, and a 28.67% reduction at 350 m level, reaching 50.56 mm. This emphasizes a substantial decrease in deformation for lower-level rib pillars, aligning with the 243.31% yield zone increase from 290 m level to 350 m level. From the analysis following outcomes can be summarized:

Table 11. Deformation values of different geo-mining properties on crown and rib pillar.

Properties		Along the path A-A' (crown pillar)	Percentage difference	Along the path B-B' (rib pillar)	Percentage difference
		Average vertical deformation (mm)		Average vertical deformation (mm)
D (m level)	230	−81.17	–	−80.08	–
	290	−78.54	3.24 ↓	−70.88	0.36 ↓
	350	−78.26	0.36 ↓	−50.56	28.67 ↓
σci (MPa)	100	−81.17	–	−80.08	–
	120	−81.52	0.43 ↑	−80.43	0.44 ↑
Em (GPa)	30	−81.17	–	−80.08	–
	37.5	−81.87	0.86 ↑	−80.04	0.05 ↓
	45	−82.2	0.40 ↑	−79.83	0.26 ↓
GSI	45	−76.05	7.11 ↓	−77.75	2.91 ↓
	55	−81.17	–	−80.08	–
	65	−82.8	1.14 ↑	−80.46	0.47 ↑
TC (m)	8	−81.17	–	−80.08	–
	10	−81.44	0.33↑	−79.17	1.14 ↓
TR (m)	10	−81.17	–	−80.08	–
	15	−80.72	1.40 ↓	−79.78	0.37 ↓
WS (m)	30	−81.17	–	−80.08	–
	40	−83.21	2.51 ↑	−79.66	0.52 ↓

• Varying σ_ci from 100 MPa to 120 MPa has minimal impact on deformation profiles, with changes within 0.5% for both crown and rib pillars. Deformation ranges between 81.17 and 81.52 mm for the crown, and from 80.08 mm to 80.43 mm for the rib pillar.

• Varying E_m from 30 GPa to 37.5 and 37.5 GPa to 45 GPa results in 0.86% and 0.40% increased deformation for the crown. Conversely, the rib pillar shows marginal decreases of 0.05% and 0.26%, suggesting no substantial influence on vertical deformation with E_m ranging from 30 GPa to 45 GPa.

• Decreasing GSI from 55 to 45 in the crown pillar leads to a 7.11% deformation reduction (from 81.17 to 76.05 mm), while increasing GSI from 55 to 65 causes a 1.14% increase. However, overall deformation changes remain within 3% considering GSI variations.

• Considering varying T_C from 8 to 10 m and T_R from 10 to 15 m shows no significant impact on deformation, with changes staying within 3%. A general apprehension is that as depth increases, vertical deformation in the cases of T_C and T_R also increases (Kahraman et al. 2003; Kumar et al. 2016). However, as illustrated in the model geometry, the dip angle increases with depth. According to Dintwe et al. (2022), a notable increase in displacement occurs at lower dip angles, leading to a decrease in deformation as the dip angle increases. Additionally, considering the model geometry, the thickness of the orebody gradually increases with depth. From an excavation standpoint, the orebody thickness is at its maximum at 350 m level and minimum at 230 m level. During excavation, the void area created induces stress redistribution, further contributing to deformation. As the void area enlarges, deformation also increases. Nevertheless, due to the greater thickness of the orebody at lower depths, the dimensions of the crown pillar and rib pillar are significantly larger compared to those at 230 m level and 290 m level. Consequently, with the incremental size of the pillar dimensions, deformation decreases as the pillars offer better support compared to the previous levels. Due to these reasons vertical deformation for T_C, and T_R is insignificant.

• While analysing the maximum vertical deformation trend of crown pillar for varying W_S in Figure 13, a gradual increase in the deformation profile occurs. When W_S increases from 30 to 40 m, widening stope cause load redistribution, resulting in higher vertical stress concentrations in the roof and floor, which increases vertical deformation of the crown pillar. Additionally, the removal of a larger volume of material reduces confinement and increases stress relaxation, making the remaining rock mass more prone to deformation under the weight of the overlying strata (Brady and Brown 2006). The larger unsupported span in wider stopes leads to more pronounced bending and sagging of the roof and floor heaving, as the effectiveness of natural arching and pillar support diminishes (Hoek and Marinos 2000) resulting further increase in the vertical deformation in crown pillar.

Extent of yielding zone

The initiation of yielding or failure within a particular zone is indicated by the emergence of plastic strain. The utilization of the Mohr-Coulomb failure criterion aims to ascertain the initiation of yielding within the rock mass, elucidated as follows (6) $$\sqrt{{\mathrm{J}}_2}\left[({\mathrm{N}}_{\mathrm{\varphi }}+1)\mathrm{\text{cos}}\hspace{0.17em}\mathrm{\varphi }+\left({\mathrm{N}}_{\mathrm{\varphi }}-1\right)\frac{\mathrm{\text{sin}}\hspace{0.17em}\mathrm{\varphi }}{\sqrt{3}}\right]-{\mathrm{\sigma }}_{\mathrm{m}}\left({\mathrm{N}}_{\mathrm{\varphi }}-1\right)-2\mathrm{c}\sqrt{{\mathrm{N}}_{\mathrm{\varphi }}}=0\mathrm{ }$$(6) where, ϕ = friction angle, c = cohesion (7) $${\mathrm{N}}_{\mathit{\varnothing }}=\frac{1+\mathrm{\text{sin}}\hspace{0.17em}\mathrm{\varphi }}{1-\mathrm{\text{sin}}\hspace{0.17em}\mathrm{\varphi }}=\text{triaxial factor}$$(7) (8) $${\mathrm{\sigma }}_{\mathrm{m}}=\frac{{\mathrm{\sigma }}_1+{\mathrm{\sigma }}_2+{\mathrm{\sigma }}_3}{3}=\text{mean stress}=1^{\mathrm{\text{st}}}\text{invariant of deviatoric stress}$$(8) (9) $${\mathrm{J}}_2=\frac{1}{2}\left[\left({\mathrm{\sigma }}_1-{\mathrm{\sigma }}_{\mathrm{m}}{)}^2\right.\right.+({\mathrm{\sigma }}_2-{\mathrm{\sigma }}_{\mathrm{m}}{)}^2+(\left.{\mathrm{\sigma }}_3-{\mathrm{\sigma }}_{\mathrm{m}}{)}^2\right]=2_{\mathrm{\text{nd}}}\text{invariant of deviatoric stress}$$(9) (10) $${\mathrm{J}}_3=\frac{1}{3}\left[\left({\mathrm{\sigma }}_1-{\mathrm{\sigma }}_{\mathrm{m}}{)}^3\right.\right.+({\mathrm{\sigma }}_2-{\mathrm{\sigma }}_{\mathrm{m}}{)}^3+(\left.{\mathrm{\sigma }}_3-{\mathrm{\sigma }}_{\mathrm{m}}{)}^3\right]=3_{\mathrm{\text{rd}}}\text{invariant of deviatoric stress}$$(10) (11) $$\mathrm{\text{sin}}\hspace{0.17em}3\mathrm{\theta }=-\frac{3\sqrt{3}{\mathrm{J}}_3}{2{{\mathrm{J}}_2}^{3/2}}\mathrm{ }$$(11)

In this context, a value of 0 indicates no yield or failure, while any value other than 0 signifies potential failure occurring in the rib pillars, around the stope, or within the orebody. From Figure 15(a,b) and Table 12, both models show potential yielding in the roofs and sides of the stopes/pillars. This aligns with the σ₁ and σ₃ plots presented earlier. In scenarios with the minimum GSI, there’s a maximum plastic strain of 0.012 observed around the rib, crown pillars, and stopes.

Figure 15. (a) Maximum plastic strain for minimum GSI value (=45) and (b) Maximum plastic strain for average GSI value (=55).

Table 12. Yield zone extent for the various geo-mining conditions.

Properties	Values		Percentage difference in yield Zone	Extent of yield zone
		Extent of yield zone Crown (m)		Rib (L)	Percentage difference	Rib (R)	Percentage difference
D (m level)	230	3.82	–	1.54	–	2.06	–
	290	8.4	119.9 ↑	1.6	3.9 ↑	1.87	9.22 ↓
	350	5.26	35.85 ↓	1.4	12.5 ↓	6.42	243.31 ↑
σci (MPa)	100	3.82	–	1.54	–	2.06	–
	120	3.47	9.16 ↓	1.42	7.79 ↓	1.82	11.65 ↓
Em (GPa)	30	3.82	–	1.54	–	2.06	–
	37.5	4.63	21.20 ↑	1.76	14.29 ↑	3.02	46.60 ↑
	45	4.07	12.1 ↓	1.47	16.48 ↓	2.62	13.25 ↓
GSI	45	6.21	62.57 ↑	3.08	100 ↑	7.7	273.79 ↑
	55	3.82	–	1.54	–	2.06	–
	65	3.37	11.78 ↓	2.63	70.78 ↑	3.95	91.75 ↑
TC (m)	8	3.82	–	1.54	–	2.06	–
	10	4.53	18.59 ↑	1.38	10.39 ↓	2.76	33.98 ↑
TR (m)	10	3.82	–	1.54	–	2.06	–
	15	4.8	25.65 ↑	1.62	5.2 ↑	3.02	46.6 ↑
WS (m)	30	3.82	–	1.54	–	2.06	–
	40	4.46	16.75 ↑	1.78	15.58 ↑	3.34	62.14 ↑

In comparison, the average GSI model exhibits a strain of 0.0065. These results suggest a possible failure in the skins or sides of the pillars while maintaining the integrity of the core under current geo-mining conditions.

Yield zone based on D

Figure 16(a,b) are used to ascertain the maximum extents of the yield zones at depths of 230 m level, 290 m level, and 350 m level for the crown pillar and rib pillar, respectively. In general, the extent of maximum yield zone increases with increases in D value. The results highlighted a distinct pattern: as the D increased from 230 m level to 290 m level, there was a significant rise in the maximum yield zone, expanding from 3.82 to 8.4 m, signifying a substantial 119.9% increase. Whereas shifting from 290 m level to 350 m level exhibits 35.85bmgv % difference in the yield zone. Notably, even under conditions of extremely low GSI, the yield zone continued to grow, representing a 62.51% increase.

Figure 16. (a) Yield zone along the crown pillar for different D (m level) and (b) Yield zone along the rib pillar for different D (m level).

Additionally, D at 290 m level, the rib pillar demonstrated a 9.22% reduction in the plastic zone compared to its presence at 230 m level, decreasing from 2.06 to 1.87 m. This trend persisted with an additional 243.31% increase in the plastic zone as the D increased to 350 m level, showing a shift from 1.87 to 6.42 m.

Yield zone based on σ_ci and W_S

Analysis of Figure 17(a) reveals changes in plastic deformation for the crown pillar, shifting from 3.82 to 3.47 m when σ_ci values vary between 100 and 120 MPa. This alteration reflects a percentage decrease of 9.16. Conversely, for the rib pillar (Figure 17b), the impact of increasing σ_ci from 100 to 120 MPa results in a shift from 2.06 to 1.82 m, with a percentage decrease of 11.65%. The findings demonstrate that larger σ_ci values coincide with a decrease in the dimensions and extent of the yield. Additionally, when increasing W_S from 30 m to 40 m, both crown and rib pillars (Figure 17(c,d), respectively) show an expansion in the plastic zone, with a percentage increase of 16.75 in case of the crown pillar whereas 62.14% increase in rib pillar.

Figure 17. (a) Yield zone along the crown pillar for different σ_ci, (b) Yield zone along the rib pillar for different σ_ci, (c) Yield zone along the crown pillar for different W_S, and (d) Yield zone along the rib pillar for different W_S.

With a pillar thickness of 8 m and a strength of 100 MPa at 350 m level depth, wider yield zones (243.31% rise) increase the likelihood of roof or rockfall incidents. Similar trends occur across different pillar thicknesses, showing variations in D, GSI, and E_m values.

Yield zone based on E_m

The findings suggest that incremental E_m notably impact the yield zone of the crown pillar by 21.2% increase whereas 46.60% increase in case of the rib pillar when E_m changes from 30 GPa to 37.5 GPa. (Figure 18a,b). However, increasing the E_m value from 37.5 GPa to 45 GPa leads to a 12.1% decrease in the crown pillar and a 16.48% decrease in the rib pillar’s plastic zone.

Figure 18. (a) Yield zone along the crown pillar for different E_m and (b) Yield zone along the rib pillar for different E_m.

Yield zone based on GSI

When examining the crown pillar (Figure 19a), a decrease in GSI from 55 to 45 yields a 62.57% increase in the crown pillar and a 273.79% increase in the rib pillar’s (Figure 19b) plastic zone. In contrast, varying GSI from 55 to 65 shows different trends for yield zone. In case of the crown pillar, a decrease of 11.78% exhibits when GSI ranges from 55 to 65 and the rib pillar showcases a 91.75% increase in the plastic zone. The result suggests that minimum GSI leads to a 7.7 m yield zone which can further lead to failure of the sidewall.

Figure 19. (a) Yield zone along the crown pillar for different GSI and (b) Yield zone along the rib pillar for different GSI.

Yield zone based on T_C and T_R

Varying T_C between 8 and 10 m yields significant effects on the yield zone, with a percentage increase of 18.59 and 33.98 for the crown and the rib pillars, respectively (Figure 20(a,b), respectively). Whereas, increasing the T_R value from 10 m to 15 m induces 25.65% and 46.6% increase in the crown and the rib pillar plastic zone, respectively. The incremental T_C and T_R leads to the yield zone value increasing from 3.82 to 4.8 m for the crown pillar and from 2.06 to 3.02 m for the rib pillar, reflecting a significant rise in yield zone.

Figure 20. (a) Yield zone along the crown pillar for different T_C, (b) Yield zone along the rib pillar for different T_C, (c) Yield zone along the crown pillar for different T_R, and (d) Yield zone along the rib pillar for different T_R.

Pillar stability in terms of factor of safety (FoS)

FoS serves as a critical indicator of stability and safety margins within an u/g metal mine employing the long drill blast hole open stoping method. Understanding its effects across varying D, σ_ci, E_m, GSI, T_C, T_R, and W_S is important to understand how these factors influence the overall stability and structural integrity of the mine. Table 13 illustrates how selected attributes affect the stability of both the crown and rib pillar, quantified in FoS.

Table 13. FoS for various geo-mining conditions.

Properties	Values	Crown pillar	Percentage difference	Rib pillar	Percentage difference
		FoS		FoS
D (m level)	230	1.04	–	0.90	–
	290	1.00	4.81 ↓	0.93	3.33 ↑
	350	1.10	11.11 ↑	1.00	6.45 ↑
σci (MPa)	100	1.04	–	0.90	–
	120	1.14	9.62 ↑	0.95	5.55 ↑
Em (GPa)	30	1.04	–	0.90	–
	37.5	1.04	0	0.91	1.11 ↑
	45	1.03	0.96 ↓	0.91	0
GSI	45	0.84	19.23 ↓	0.83	7.78 ↓
	55	1.04	–	0.90	–
	65	1.28	23.08 ↑	1.04	15.56 ↑
TC (m)	8	1.04	–	0.90	–
	10	0.96	7.69 ↓	0.86	4.44 ↓
TR (m)	10	1.04	–	0.90	–
	15	0.99	4.81 ↓	0.99	10 ↑
WS (m)	30	1.04	–	0.90	–
	40	0.98	5.77 ↓	0.86	4.44 ↓

FoS analysis presents a comprehensive overview of the overarching stability levels across crown and rib pillars under diverse scenarios. The findings underscore the importance of maintaining optimal crown pillar thickness, stope width, and material properties to ensure adequate safety margins. Higher FoS values corresponds to configurations with thicker crown pillars and wider stopes, emphasizing their pivotal roles in increasing structural stability and mitigating potential risks associated with pillar failure. Collectively, these factors intricately shape the fa FoS, directly impacting the stability and safety margins of u/g metal mine employing the long drill blast hole open stoping method. Optimizing these parameters is paramount in enhancing the structural stability and safety of the mine, mitigating potential risks associated with structural failure or instability.

Comparative analysis of simulation result with stability graph approach

The Potvin Stability Graph Method (Potvin 1988), developed by Yves Potvin in the late 1980s, has emerged as a prominent technique in geotechnical engineering for evaluating the stability of u/g excavations within rock masses. This method offers a pragmatic framework for assessing excavation stability, particularly in contexts like mining and civil engineering projects where ensuring rock mass stability is critical. Central to the Potvin Stability Graph Method is its reliance on rock mass classification systems, such as the Rock Mass Rating (RMR) or the Q-system. These systems serve to delineate the geological and geo-mechanical attributes of rock masses, using factors such as rock strength, joint orientation, and overall rock quality. By incorporating data derived from these classification systems, engineers can effectively interpretate the overall stability of the rock mass. The graphical representation at the core of the Potvin Stability Graph Method correlates the excavation span, typically measured as the distance between rock supports, with the assigned rock mass rating. This graphical depiction elucidates the intricate relationship between excavation size and the stability of the encompassing rock mass. Notably, the graph typically comprises stability curves or lines, each indicative of varying levels of stability corresponding to distinct rock mass ratings.

Implementing the Potvin Stability Graph Method entails a systematic approach. Initially, engineers conduct a thorough assessment of geological conditions, employing a rock mass classification system to identify pertinent parameters. Subsequently, they locate the appropriate stability curve on the graph based on the assigned rock mass rating. By plotting the excavation span on the graph, engineers can ascertain whether the proposed excavation size falls within the stable region or necessitates additional support measures.

Assessing the reliability of the FE model required comparing it against Potvin’s Stability Graph Method (SGM). Potvin’s method utilizes key parameters such as Hydraulic Radius (HR) and Modified Stability Number (N'). HR is calculated through the ratio of surface area to perimeter, as defined in the Eq. (12)(12) $$\mathit{\text{HR}}=\frac{\mathit{\text{Area of the surface}}}{\mathit{\text{Perimeter of the surface}}}\mathrm{ }$$(12) . (12) $$\mathit{\text{HR}}=\frac{\mathit{\text{Area of the surface}}}{\mathit{\text{Perimeter of the surface}}}\mathrm{ }$$(12)

Meanwhile, N′ combines Rock Quality Index (Q'), rock stress (A), defect orientation (B), and design surface (C) factors in the Eq. (13)(13) $$N^{\prime }=Q^{\prime }\mathrm{*A*B*C }$$(13) . (13) $$N^{\prime }=Q^{\prime }\mathrm{*A*B*C }$$(13)

Past geotechnical studies within mining areas provided the Q, A, B, and C values associated with crown pillars and tabulated parameters for crown pillars (GSI = 45, σ_ci=100 MPa, E_m=45 GPa, T_C=8 m, T_R=10 m, and S_W=30 m) were utilized to plot a stability graph (Figure 21), indicating consistent stability across various depths.

Figure 21. Potvin stability graph for crown pillar stability.

Further stability analysis was conducted by determining the factor of safety (FoS) for crown pillars using FE models, as described in Eq. (14)(14) $$\mathit{\text{FoS}}=\frac{{\sigma }_3+{\sigma }_{\mathit{\text{ci}}}{\left(m_b\frac{{\sigma }_3}{{\sigma }_{\mathit{\text{ci}}}}+S\right)}^a}{{\sigma }_1}\mathrm{ }$$(14) . The resultant FoS values, detailed in Table 14, consistently surpass 1 at all depths, signifying stability (Hoek et al. 2002). A comparative assessment between FE model outcomes and the SGM revealed similar results, validating the efficacy of the FE simulation findings. (14) $$\mathit{\text{FoS}}=\frac{{\sigma }_3+{\sigma }_{\mathit{\text{ci}}}{\left(m_b\frac{{\sigma }_3}{{\sigma }_{\mathit{\text{ci}}}}+S\right)}^a}{{\sigma }_1}\mathrm{ }$$(14) where σ_ci = uniaxial compressive strength of intact rock. σ₁ and σ₃ = effective principal stresses (major and minor, respectively) at failure; m_b, s, and a are Hoek–Brown rock mass parameters.

Table 14. FE simulation analysis versus stability graph approach.

Depth of working, D (m level)	A	B	C	Q'	N'	HR (m)	Stability state	FoS
230	0.16	0.80	2.84	12.33	4.48	3.38	Stable	1.04
290	0.20	0.80	2.19	12.33	4.32	3.92	Stable	1.00
350	0.12	0.80	2.66	12.33	3.15	2.86	Stable	1.10

For validation of the study, a comparative analysis plot between empirical Potvin’s SGM and numerical simulation results for the crown pillar is presented in Figure 21. The stability state of stopes analyzed using both SGM and numerical methods demonstrate similar results, thereby further validating the analysis.

Conclusions

In conclusion, the study provides a comprehensive analysis of the distribution of stresses around crown pillars in a mine utilizing the long drill blast hole open stoping method. The findings emphasize the significance of considering various parameters, including stress concentration factors, deformation patterns, yield behaviour, and factor of safety. The study reveals intricate insights into the stability dynamics of crown and rib pillars, shedding light on their behaviour under different conditions.

The research highlights that stress concentration is predominantly centred at the core of the crown pillar due to the robust nature of the adjacent rock mass and heightened horizontal stresses. As depth increases, the κ₁ and κ₃ values of the crown pillar decrease, attributed to the widened pillar and stress relief within the stope vicinity. Simultaneously, the tensile zone expands, especially under certain geological conditions. Rib pillars, on the other hand, encounter tensile yielding at deeper levels, influenced by factors such as rock mass properties and variations in thickness.

Deformation analysis indicates that while there is not a notable shift in the average deformation profile of the crown pillar, extreme values of certain parameters significantly influence the frequency and extent of deformation. The study identifies a consistent trend where an increase in GSI values from 45 to 55 leads to a significant rise in crown pillar deformation.

The examination of the crown pillar behaviour reveals notable trends, such as the substantial increase in the yield zone with varying depths and the pronounced influence of factors like σ_ci, GSI and T_R. The rib pillar exhibits a similar trend, with a notable rise in the yield zone under specific conditions.

The research concludes by emphasizing the viability of operating at lower depths despite an initial reduction in FoS. It suggests that improvements in the overall thickness of the orebody contribute to enhanced stope stability, making it feasible to operate even at greater depths. The scientific research provides valuable insights into the complex interplay of factors affecting the stability of crown and rib pillars in u/g mining. Furthermore, the study underscores the importance of specifications such as specific values governing the pillar dimensions for better utilization of optimization strategies to ensure safe mining operations.

Acknowledgments

The authors express gratitude to the mine management for supplying the essential data needed for this study. This study forms a part of the ongoing doctoral research of the first author. The support extended by the personnel from various departments at IIT(ISM), Dhanbad, is sincerely appreciated. It is important to note that the opinions expressed in this paper solely belong to the authors and do not necessarily reflect the views of their respective organizations.

Disclosure statement

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this article.

Data availability statement

The datasets analysed during this study are available from the corresponding author upon reasonable request.

Additional information

Funding

This research has not received any funding from any source.