A SUSTAINABLE AND PRACTICAL MACHINE LEARNING APPROACH USING SCIKIT-LEARN FOR PREDICTING STOPE INSTABILITY: IDENTIFICATION OF CRITICAL GEOTECHNICAL FACTORS

2.1. Geology and Mining Overview of the Study Area

The mine exploits auriferous sulfide mineralization within shear zones hosted by the Archaean Iron Mask Formation, located in the Harare-Bindura-Shamva greenstone belt of the Zimbabwe Craton (>2.5 Ga) (Blenkinsop et al., 1999). The Iron Mask Formation consists of metamorphosed felsic volcanics (dacite, meta-andesite, meta-rhyolite), with gold mineralization primarily associated with arsenopyrite, pyrrhotite, and pyrite. The ore zones, averaging 1m in width and 3.7 g/t Au, dip variably (10°-70°, mean 40°). The mine utilizes sublevel stoping with random pillar support. Stope dimensions range from 1-3m wide and 10-30m high. Detailed discontinuity mapping and rock mass characterization are essential, considering the lithological variability of the Iron Mask Formation and its impact on rock mass strength and deformation. Notably, a period of stope flooding, followed by resumed mining without a documented geotechnical assessment, necessitates a thorough evaluation of potential water-induced instability. This includes assessing water pressure effects on joint strength and potential pore pressure development.

The regional geological setting of the Zimbabwe Craton, highlighting its major lithological units, is illustrated in Figure 1.

Figure 1. Simplified geological map of the Zimbabwe Craton, showing major lithological units. Adapted from Blenkinsop et al. (1999).

The Harare-Bindura-Shamva Greenstone Belt, a significant mineralized zone within the Zimbabwe Craton, is delineated by a red rectangle on the regional geological map shown in Figure 2.

Figure 2. Regional geological map of the Zimbabwe Craton highlighting the Harare-Bindura-Shamva Greenstone Belt. The red rectangle indicates the approximate location of the greenstone belt. Adapted from Blenkinsop et al. (1999).

Detailed mapping within the Harare-Bindura-Shamva Greenstone Belt has revealed several significant geological structures that control mineralization. These key structures are presented in Figure 3. For instance, the major shear zones are evident.

Figure 3. Detailed map illustrating significant structures within the Harare-Bindura-Shamva Greenstone Belt. Adapted from Blenkinsop et al. (1999).

2.2. Rock Strength Determination

Circular cylindrical core samples were prepared for uni-axial and tri-axial tests from areas under study. The samples were cut to 1150mm size using a diamond saw. For the Brazilian tensile strength test, samples were cut to 16mm length, and all specimens were given identity numbers. The samples were supplied for laboratory testing.

2.3. Rock Mass Classification

The ground water conditions, discontinuity direction, thickness, infill type, joint set number, spacing, roughness, and other geotechnical characteristics impacting rock mass performance at certain places were measured and documented in the study. Field estimates were used to estimate the strength of the rock, and data from uniaxial compressive strength (UCS) tests was added. The joint frequency approach suggested by Palstrom (1982) was utilized to determine the RQD rating because the drill core was not available. The volumetric joint count, J_v, was first calculated using the joint spacing per interval as stated in Equation 1.

[TeX:] \mathbf{Jv =}\frac{\mathbf{1}}{\mathbf{j}\mathbf{1}}\mathbf{+}\frac{\mathbf{1}}{\mathbf{j}\mathbf{2}}\mathbf{+}\frac{\mathbf{1}}{\mathbf{j}\mathbf{3}}$𝐉𝐯=\frac{1}{𝐣1}+\frac{1}{𝐣2}+\frac{1}{𝐣3}$ (1)

where joint sets 1, 2, and 3 are represented by J1, J2, and J3 respectively.

Additionally, information was gathered in compliance with the RMR system of rock mass classification specifications (Bieniawski, 1989).

2.4. Evaluation of Stope Stability

Stability graphs, originally introduced by Mathews et al. (1981), are among the most widely adopted empirical methods for evaluating stope stability in underground mining. This approach was developed using established rock mass classification systems, notably the Q-system by Barton et al. (1974) and the Rock Mass Rating (RMR) system by Bieniawski (1973). The Mathews stability graph focuses on identifying the main factors that influence rock mass stability through specially designed charts that correlate various rock mass properties. These charts incorporate parameters such as the rock stress factor (A), joint orientation adjustment factor (B), modified rock tunneling index (Q') and surface orientation factor (C). These factors are used collectively to compute the stability number (N’), a significant metric in the design of stope dimensions and support systems. The stability number serves as a quantitative measure of the rock mass conditions and stope stability and is determined using the relationship in Equation 2.

[TeX:] N^{'} = Q^{'}\text{ABC}$N^{\prime }=Q^{\prime }\mathrm{\text{ABC}}$ (2)

The Mathews' stability graph method proposed by Mathews et al., (1981) was used to evaluate the stope stability of the stopes. The average of the Q-system data generated for every stope was used to calculate the Q' value. Potvin factor analysis was performed at each stope's specific places where the RMR Q-system had been concluded. The modified stability number N' of each and every stope was computed using the modal Potvin factor values, as well as their respective hydraulic radii. The stope dimensions, thus width and height were measured using a distometer so as to determine the hydraulic radii of the stopes. These values were used to determine the stope stability on the stability graph and the maximum tolerable unsupported length. For each stope under investigation, the modified stability number N’ and the hydraulic radius was determined using Eqaution 2.

A vital element in the accurate assessment of stope stability is the shape factor, commonly known as the hydraulic radius (HR), which links the geometric dimensions of the opening. The hydraulic radius is a major parameter that characterizes the shape of the stope and plays a significant role in determining its structural stability. It is typically defined as the ratio of the area of the exposed hanging wall to its perimeter. In the case of inclined stopes - where the excavation is not perfectly vertical - the hanging wall exposure becomes the most important consideration for calculating the HR. The HR calculation incorporates the stope’s span along both the dip (h) and strike (w) directions (Tishkov, 2018) as shown in Equation 3.

[TeX:] Hydraulic\ Radius = \frac{\text{wh}}{2\left( w + h \right)}$HydraulicRadius=\frac{\mathrm{\text{wh}}}{2(w+h)}$ (3)

2.5. Assessment of Blast Damage

The measuring tools Peak Particle Velocity (PPV) and Half Cast Factor were employed to evaluate blast damage. Vibration sensors were installed at strategic locations to measure the PPV generated by the blast. PPV values were recorded at each sensor location during the blast event. The PPV was determined using the formula proposed by Singh (1994) in Equation 4.

[TeX:] \mathbf{Vmax = K}\left( \frac{\mathbf{R}}{\sqrt{\mathbf{Q}}} \right)^{\mathbf{-}}\mathbf{ᵝ}$𝐕𝐦𝐚𝐱=𝐊{\left(\frac{𝐑}{\sqrt{𝐐}}\right)}^{-}ᵝ$ (4)

where Vmax is the PPV, from an explosive charge, Q, at a known distance R.

The half cast factors (Equation 5) were calculated by measuring the length of half cast barrels in the designated areas that remained after the blast:

Half Cast Factor =[TeX:] \frac{\sum_{}^{}\text{Lengths}\text{ofvisibleblastholes}\left( \text{afterexplosion} \right)}{\sum_{}^{}\text{Lengths}\text{ofperimeterblastholes}\left( \text{beforeexplosion} \right)}$\frac{{\sum }_{}^{}\mathrm{\text{Lengths}}\mathrm{\text{ofvisibleblastholes}}\left(\mathrm{\text{afterexplosion}}\right)}{{\sum }_{}^{}\mathrm{\text{Lengths}}\mathrm{\text{ofperimeterblastholes}}\left(\mathrm{\text{beforeexplosion}}\right)}$ (5)

2.6. Data Analysis

The collected data from field investigations, laboratory experiments, and an extensive literature review were subjected to systematic analysis using a combination of established geotechnical and empirical techniques. Rock strength parameters were evaluated through both uniaxial compressive and triaxial compressive strength tests to determine the mechanical behaviour of the intact rock under varying stress conditions. Rock mass quality was assessed using the RQD and RMR classification systems, providing insight into the structural integrity and geomechanical characteristics of the rock mass. Stope stability was analyzed utilizing the Mathews' stability graph method, an empirical tool that relates rock mass conditions, geometry, and support requirements to the potential for stable excavation. Blast-induced damage was assessed by analyzing PPV data and the Half Cast Factor, both of which are indicators of blast performance and fragmentation control. This multifaceted analytical approach ensured a comprehensive understanding of the subsurface conditions and their implications for excavation and support design.

2.7. Factor Analysis

Machine learning techniques were applied to analyze the dataset, encompassing processes such as feature selection, feature importance evaluation, and decision boundary visualization. The analysis was conducted using Scikit-learn, a widely adopted machine learning library in Python, on a dataset integrating geological, geotechnical, and mining-related parameters. Feature selection was carried out using mutual information and recursive feature elimination (RFE) to identify the most relevant variables contributing to model performance and predictive accuracy. Feature importance was assessed using the Random Forest Classifier, which provided insight into the relative significance of each input variable in the classification tasks. To visualize model behaviour and class separation, decision boundary plots were generated using Support Vector Machines (SVMs). Additionally, feature distribution plots were employed to explore the distribution and potential interactions among variables. To further understand the interrelationships among geotechnical factors, correlation coefficients were computed and analyzed. Figure 4 summarizes the whole factor analysis process.

Figure 4. Workflow Diagram of the Machine Learning Process for Stope Instability Prediction

2.7.1. Feature Importance Analysis

As previously described, the Random Forest Classifier was used for feature importance analysis. The core of the Random Forest lies in the aggregation of multiple decision trees. Each tree Ͳ_j in the forest \T₁,T₂,…, T_m} is built on a bootstrap sample of the training data and considers a random subset of features at each split. The prediction of the forest for a given input x is often based on the majority vote of its individual decision trees. This ensemble prediction, ŷ(x), can be expressed as shown in Equation 6:

[TeX:] \widehat{y}\left( x \right) = \text{mode}\text{\{}T_{j}\left( x \right)\text{\}}_{j = 1}^{m}$\hat{y}\left(x\right)=\mathrm{\text{mode}}\mathrm{\text{{}}{T}_j\left(x\right){\mathrm{\text{}}}}_{j=1}^m$ (6)

Where ŷ(x) is the final predicted output of the Random Forest model, mode denotes the statistical mode, which is the value that appears most frequently, T_j(x) is the prediction of the j^th individual decision tree for a given input x and j is the total number of trees in the Random Forest.

The importance of a feature X_k is quantified by observing how much the prediction accuracy (or impurity) decreases when that feature is randomly permuted. A feature is considered important if its permutation leads to a significant drop in model performance. The importance score I(X_k) for a feature X_k can be expressed as shown in Equation 7:

[TeX:] I\left( X_{k} \right) = \frac{1}{m}\sum_{j = 1}^{m}\left( e_{j} - e_{j,k}^{\text{perm}} \right)$I\left(X_k\right)=\frac{1}{m}{\sum }_{j=1}^m(e_j-e_{j,k}^{\mathrm{\text{perm}}})$ (7)

where e_j is the out-of-bag (OOB) error for tree T_j and e_j,k ^permis the OOB error of Ͳ_j after the values of feature X_k have been randomly permuted. The OOB error is calculated on the data points not used to train that specific tree. Categorical variables (water effect and blast damage) were encoded using ordinal encoding, and numerical features (RQD, Rock Strength, Rock Mass Rating) were standardized using the z-score as shown in Equation 8:

[TeX:] z_{i} = \frac{x_{i} - \mu_{X}}{\sigma_{X}}$z_i=\frac{x_i-{\mu }_X}{{\sigma }_X}$ (8)

where z_i is the standardized value of the i^th data point, x_i is the original value of the i^th data point, µ_x is the mean of the feature X and σ_x is the standard deviation of the feature X.

The z-score: ensures that the variance across features is comparable, which can be important for distance-based algorithms (though Random Forest is less sensitive to feature scaling).

2.7.2. Recursive Feature Elimination Analysis

Recursive Feature Elimination (RFE) is a feature selection technique that recursively removes the least important features until a specified number of features is reached. The process can be described as follows:

The importance of each feature in the Recursive Feature Elimination (RFE) process was calculated using Equation 9:

[TeX:] I\left( X_{k} \right) = \left| w_{k} \right|$I\left(X_k\right)=\left|w_k\right|$ (9)

where I(X_k) is the importance of feature X_k and is the w_x weight assigned to feature X_k in the model.

The RFE process can be mathematically represented by Equation 10:

[TeX:] \text{RFE}\left( D,K \right) = X_{1},X_{2},\ldots,X_{K}$\mathrm{\text{RFE}}(D,K)=X_1,X_2,\dots ,X_K$ (10)

where RFE(D, K) is the set of K features selected from the original dataset D.

2.7.3. Decision Boundary Analysis

Decision boundary plots were generated using the Random Forest Classifier to visualize the classifier's predictions across the feature space of the two most significant factors: water effect and blast damage. The decision boundary Ɓ for separating the feature space into regions corresponding to different class predictions is defined as shown in Equation 11:

[TeX:] \mathcal{B} = \text{\{}x \in R^{d} \mid f\left( x \right) = c\text{\}}$ℬ=\mathrm{\text{{}}x\in R^d\mid f\left(x\right)=c\mathrm{\text{}}}$ (11)

where x is a vector of the 'Water Effect' and 'Blast Damage' features (after encoding), d is the dimensionality of this subspace (here, d=2),R^d I the Euclidean space and c is the class boundary (e.g. the point where the probability of 'Stable' equals the probability of 'Unstable'). For probabilistic classifiers like Random Forest (which can output class probabilities), the decision boundary can be defined at a specific probability threshold (e.g. 0.5).

2.7.4. Data Distribution Analysis

Feature distribution plots were used to visualize the distributions of individual geotechnical factors and the relationships between pairs of factors. The relationships between pairs of variables were explored using joint probability distributions and the Pearson correlation coefficient[TeX:] .$.$ The Pearson correlation coefficient, r, between two variables X and Y is defined as shown in Equation 12:

[TeX:] r_{\text{XY}} = \frac{\text{E}\left\lbrack \left( X - \mu_{X} \right)\left( Y - \mu_{Y} \right) \right\rbrack}{\sigma_{X}\sigma_{Y}}$r_{\mathrm{\text{XY}}}=\frac{\mathrm{\text{E}}\left[(X-{\mu }_X)(Y-{\mu }_Y)\right]}{{\sigma }_X{\sigma }_Y}$ (12)

where E is the expectation operator, µ_x is the mean of variable X, µ_x is the mean of variable Y, σ_x is the standard deviation of variable X and σ_y the standard deviation of variable Y.

2.7.5. Correlation Coefficient Analysis

Correlation coefficient results were used to determine linear relationships among the geotechnical factors, specifically focusing on the correlation between water effect, blast damage, and rock quality (represented by RQD and Q). The t-statistic for assessing the significance of a correlation coefficient is given by Equation 13:

[TeX:] t = r\sqrt{\frac{n - 2}{1 - r^{2}}}$t=r\sqrt{\frac{n-2}{1-r^2}}$ (13)

where n is the number of data points, t is the t-statistic value, r is the Pearson correlation coefficient and n is the number of data points.

The resulting p-value indicates the probability of observing such a correlation if there were no true linear relationship between the variables.

2.7.6. Model Evaluation

The performance of the Random Forest Classifier was evaluated using standard classification metrics, including accuracy, precision, recall, and F1-score. For the 'Unstable' class:

The overall accuracy of the model, precision, recall and F1-score was calculated as shown in Equation 14, 15, 16 and 17.

[TeX:] Accuracy = \frac{\text{Number\ of\ Correct\ Predictions}}{\text{Total\ Number\ of\ Predictions}}$Accuracy=\frac{\mathrm{\text{Number of Correct Predictions}}}{\mathrm{\text{Total Number of Predictions}}}$ (14)

[TeX:] Precision = \frac{\text{True\ Positives}}{\text{True\ Positives} + \text{False\ Positives}}$Precision=\frac{\mathrm{\text{True Positives}}}{\mathrm{\text{True Positives}}+\mathrm{\text{False Positives}}}$ (15)

[TeX:] Recall = \frac{\text{True\ Positives}}{\text{True\ Positives} + \text{False\ Negatives}}$Recall=\frac{\mathrm{\text{True Positives}}}{\mathrm{\text{True Positives}}+\mathrm{\text{False Negatives}}}$ (16)

[TeX:] F1 - score = 2 \times \frac{\text{Precision} \times \text{Recall}}{\text{Precision} + \text{Recall}}$F1-score=2\times \frac{\mathrm{\text{Precision}}\times \mathrm{\text{Recall}}}{\mathrm{\text{Precision}}+\mathrm{\text{Recall}}}$ (17)

To ensure the model's robustness, 5-fold cross-validation was performed. The dataset was divided into 5 folds, and the model was trained and evaluated 5 times, each time using a different fold as the validation set. The average performance across the folds provided a more reliable estimate of the model's generalization ability. The average metric across K folds in cross-validation was computed using Equation 18:

[TeX:] \overline{\text{Metric}} = \frac{1}{K}\sum_{i = 1}^{K}\text{Metric}_{i}$\overline{\mathrm{\text{Metric}}}=\frac{1}{K}{\sum }_{i=1}^K{\mathrm{\text{Metric}}}_i$ (18)

where K is the number of folds used in cross-validation, Metric_i is the value of the evaluation metric on the i^th validation fold.

ROC Curve for Logistic Regression

Furthermore, the ROC curve used was to evaluate the performance of a logistic regression model. The Receiver Operating Characteristic (ROC) curve is a graphical representation of the performance of a binary classifier. It plots the True Positive Rate (TPR) against the False Positive Rate (FPR) at different thresholds.

The TPR and FPR can be calculated using Equation 19 and 20:

[TeX:] TPR = \frac{\text{TP}}{TP + FN}$TPR=\frac{\mathrm{\text{TP}}}{TP+FN}$ (19)

[TeX:] FPR = \frac{\text{FP}}{FP + TN}$FPR=\frac{\mathrm{\text{FP}}}{FP+TN}$ (20)

where TP is the number of true positives, FN is the number of false negatives, FP is the number of false positives, and TN is the number of true negatives.

The area under the ROC curve (AUC) was calculated using Equation 21 to quantify the overall performance of a logistic regression model.

[TeX:] AUC = \int_{0}^{1}{\text{TPR}\left( \text{FPR} \right)\text{dFPR}}$AUC={\int }_0^1\mathrm{\text{TPR}}\left(\mathrm{\text{FPR}}\right)\mathrm{\text{dFPR}}$ (21)

The AUC value ranges from 0 to 1, where 1 represents perfect classification and 0.5 represents random guessing.

The logistic regression model can be represented mathematically as shown in Equation 22:

[TeX:] P\left( Y = 1 \middle| X \right) = 1\text{/}\left( 1 + e^{\left( - z \right)} \right)$P(Y=1|X)=1\mathrm{\text{/}}(1+e^{(-z)})$ (22)

where P(Y=1|X) is the probability of the positive class given the input features X, and z is a linear combination of the input features (see Equation 23):

[TeX:] z = w_{0} + w_{1}X_{1} + w_{2}X_{2} + \ldots + w_{n}X_{n}$z=w_0+w_1{X}_1+w_2{X}_2+\dots +w_n{X}_n$ (23)

Where w_i are the weights assigned to each feature X_i.

3.1. Rock Mechanical Properties

The mechanical properties of intact rock are fundamental in assessing the overall stability and deformation behaviour of underground excavations. In this study, laboratory tests were conducted on representative rock samples to determine key strength parameters, including Uniaxial Compressive Strength (UCS), Triaxial Compressive Strength (TCS), and tensile strength, providing insight into the inherent competence of the rock materials.

The UCS results demonstrated a relatively wide but consistently high strength range, with values spanning from 167 MPa in the ore zone to 215 MPa in the granodiorite unit (see Table 1). These values clearly classify both lithologies as strong to very strong rocks based on standard rock strength classification systems (ISRM). Such high compressive strengths imply that, under unconfined loading conditions, the intact rock is unlikely to fail or undergo significant deformation, even when subjected to the stresses typically encountered in underground mining environments.

Table 1. Summary of Rock Mechanical Properties

Note: Triaxial Compressive Strength Test was only done on the ore zone and the Tensile Strength Test only on granodiorite and the ore zone.

Further testing on ore zone samples revealed a TCS of 93.2 MPa, which, while lower than the UCS, reflects the increased confinement typical of in-situ conditions and demonstrates the material's capacity to sustain stress under triaxial loading. This value is within the expected range for strong, brittle rock and supports the conclusion that the ore body itself, when intact, is mechanically robust.

The tensile strength values were also notably high, with the ore zone recording 11.6 MPa and granodiorite slightly higher at 12.2 MPa. Tensile strength is critical in determining the rock's resistance to crack initiation and propagation, especially in the presence of blast vibrations or stress redistributions. These values suggest that the intact rock would not readily fracture under tensile loading, further reinforcing its competent mechanical nature.

However, despite these strong intact rock properties, field observations and empirical assessments indicate instability in certain stopes. This discrepancy points to a key understanding in rock mechanics: while laboratory strength data provides valuable baseline information, rock mass behaviour in situ is dominantly controlled by discontinuities (such as faults, joints, fractures, and bedding planes) rather than the properties of the intact rock alone (Li et al. 2023; Vinay et al. 2023). The presence of geological discontinuities can dramatically reduce the effective strength of the rock mass, acting as planes of weakness along which shear displacement, dilation, or separation may occur. This is particularly relevant in jointed or faulted ground, where intact rock strength becomes a secondary consideration (Qi et al. 2018). As a result, the rock mass strength is governed not only by the intact strength but also by the orientation, persistence, spacing, and surface conditions of these discontinuities.

In the context of the studied stopes, it is reasonable to conclude that failure mechanisms are likely driven by structural weaknesses rather than by the degradation of the intact rock. This assertion is supported by other findings in the study, such as moderate to high levels of fracturing (as evidenced by RQD values) and elevated blast-induced damage (high PPV values). These factors collectively reduce the overall rock mass quality and increase the likelihood of instability, despite the high laboratory-derived strength values.

3.2. Blast Damage Assessment

Blasting is an essential component of underground mining, enabling efficient rock breakage and extraction. However, improper blast design or execution can result in excessive damage to the surrounding rock mass, compromising excavation stability. One of the most reliable indicators of blast-induced damage is PPV - a dynamic measurement of the vibration intensity caused by explosive charges.

In this study, PPV values recorded during stope development ranged from 229.5 mm/s to 3335 mm/s (see Figure 4). This considerable variation in PPV reveals inconsistent control over blast energy dissipation, with values at the higher end far exceeding typical thresholds associated with safe excavation practices. According to various empirical guidelines and case studies, PPV values exceeding 1000 mm/s are often associated with severe stress redistribution, overbreak, and long-term deterioration of rock mass properties, particularly in jointed or fractured ground. The observed upper range of 3335 mm/s is therefore indicative of substantial blast-induced stress and potential rock mass degradation in affected areas.

Further compounding this issue is the mean half-cast factor, a metric used to evaluate the effectiveness of perimeter control blasting. In this case, the half-cast factor was calculated at only 4.92% (see Figure 6) - a value considerably lower than the acceptable range for well-controlled blasting. The half-cast factor essentially reflects the proportion of blast holes that produce smooth, half-round impressions on the final excavation perimeter. Low percentages indicate poor energy confinement, excessive rock spalling, and a lack of controlled fracturing at the boundary, leading to unintended overbreak and increased wall damage.

Figure 6. Half-Cast Factor Measurements for Different Stope Ends

The combination of high PPV values and a low half-cast factor strongly suggests that blast damage is a significant contributor to the deterioration of stope walls, beyond what might be explained solely by geological factors such as naturally blocky ground or fault structures. While the presence of pre-existing discontinuities undoubtedly plays a role in how energy is transmitted and fractures propagate, the primary mechanism of instability in this case appears to be anthropogenic - specifically, suboptimal blasting practices.

It is also important to note that high PPV values have a cumulative effect on the rock mass, especially in zones with repeated blasting cycles. Damage induced by vibration may not be immediately visible but can manifest over time as gradual loosening, slabbing, or rockfall. In addition, blast-induced microcracks reduce the rock's elastic modulus and shear strength, thereby lowering the overall stability of the excavation even in the absence of visible failures.

3.3. Rock Mass Classification

The quality of the rock mass was evaluated using the RQD and the Q-system, both of which are widely accepted empirical classification methods in geotechnical engineering. These assessments provided quantitative measures of the rock mass integrity and structural competence. The results of this analysis are presented in Table 2.

Table 2. Summary of the rock classification results

3.3.1. Evaluation of Rock Mass Quality Using RQD

The Rock Quality Designation (RQD) serves as a key indicator of the degree of fracturing within a rock mass and is widely used in empirical rock mass classification systems, including the Q-system. In this study, RQD values were assessed across four stopes to evaluate the integrity, continuity, and competence of the rock mass, all of which are critical factors for ensuring the stability of underground excavations.

The results indicate that Stopes 2A and 3C exhibited RQD values of 73% and 70% respectively, which classify them as moderately fractured rock masses, as shown in Table 2. These values suggest the presence of frequent jointing and discontinuities, which can compromise the structural behaviour of the stope walls and roofs, making them more susceptible to instability and overbreak. The moderate RQD values in these stopes are symptomatic of less competent rock that may require additional support measures to ensure safe stope development and long-term excavation performance.

The lower RQD values observed in Stopes 2A and 3C can be attributed to a combination of blast-induced damage and hydrogeological effects. Excessive blasting, especially in poorly controlled rounds, can create new fractures and extend existing ones, thereby deteriorating the surrounding rock mass. Moreover, elevated water pressure within joints and fractures can further degrade rock quality by reducing effective stress and contributing to joint dilation and weakening. These factors act synergistically to exacerbate pre-existing geological weaknesses, leading to the observed reduction in RQD.

In contrast, Stope 2B recorded a significantly higher RQD, indicating improved rock mass quality and reduced fracture density. This suggests a more intact and cohesive rock structure, which enhances load-bearing capacity and decreases the likelihood of shear or tensile failure along joint planes. As such, Stope 2B is expected to exhibit better performance in terms of ground control, reduced support requirements, and lower risk of overbreak. Notably, Stope 2B exhibited the highest RQD at 90%, denoting a very competent and unfractured rock mass. Such high values are indicative of minimal discontinuities, with long, continuous core pieces recovered during drilling, reflecting superior geological conditions. This degree of rock mass integrity is typically associated with enhanced stope stability, increased safety, and more cost-effective excavation due to lower support demands.

This interpretation is further supported by field observations and geotechnical logs, which noted signs of water ingress and irregular fracture patterns in these stopes. The presence of water not only promotes mechanical deterioration but can also lead to chemical alteration of joint infill, compounding the loss in structural cohesion. As a result, the stope stability in these areas is compromised, and targeted mitigation strategies (such as improved blast control, pre-drainage techniques, and localized support installation) may be necessary.

3.3.2. Evaluation of the Rock Mass Quality Using the Q-system

In this study, the Q-system analysis yielded distinct variations across the analyzed stopes, revealing significant differences in rock mass conditions and their implications for stope stability. Stope 2B recorded the highest Q-value of 16.9, categorizing it as a good quality rock mass, as shown in Table 2. This high value is primarily attributed to favourable geological conditions, specifically the absence of groundwater and the presence of clean, clay-free joint infill. The lack of water is particularly significant, as water acts as a destabilizing agent by reducing effective stress, lubricating joint surfaces, and weakening rock-bridging elements. Moreover, clay-free joints contribute to better interlocking and higher shear resistance, both of which are essential for maintaining the structural integrity of excavated openings.

The elevated Q-value for Stope 2B directly correlates with its observed stability and larger tolerable unsupported span, suggesting that it may require minimal ground support, thereby offering potential cost and operational efficiency advantages during mining operations. These findings underscore the importance of dry conditions and joint cleanliness as critical contributors to excavation stability.

In contrast, Stopes 2A and 3C exhibited significantly lower Q-values of 3.8 and 4.9, respectively, placing them in the fair to poor rock quality category. These values reflect moderate stability conditions that warrant more conservative excavation designs and likely necessitate reinforced support systems. A key factor influencing the reduced Q-values in these stopes is the presence of groundwater, which adversely affects several Q-system parameters. Water reduces the Jw factor, indicating wet or saturated conditions, and contributes to the alteration of joint infill materials, thereby decreasing the Ja and potentially increasing joint activity and deformation risk. In such environments, joints are more likely to become critically stressed, leading to increased dilation, rock block detachment, and potential stope wall instability.

Interestingly, Stope 3C displayed a slightly higher Q-value than Stope 2A, which can be attributed to its higher joint friction, represented by a greater Jr value. Joint roughness enhances mechanical interlocking between rock blocks, providing increased resistance to shear failure even in the presence of water or other weakening agents. This factor offers marginally better performance in Stope 3C, although not sufficient to significantly alter its classification within the Q-system. The influence of joint friction also highlights the multifactorial nature of rock mass behaviour, where improvements in one parameter may partially offset deficiencies in others.

Overall, the Q-system analysis provided quantitative insight into the geomechanical performance of the various stopes. The results emphasize that Stope 2B offers a more stable mining environment, while Stopes 2A and 3C require additional engineering interventions. These findings reinforce the value of the Q-system in pre-mining assessments and its ability to guide risk-informed stope design, ground support selection, and excavation sequencing in underground mining operations

3.3.3. Use of the Q results to evaluate stope stability

The maximum tolerable unsupported span, derived using the Q-system empirical design method, serves as a critical parameter in assessing the relative stability of underground stopes. This span refers to the largest excavation width that can be safely maintained without the use of additional ground support, based on the rock mass quality and structural conditions. It is directly influenced by the Q-value, which integrates key factors such as RQD, joint set number, joint roughness, groundwater conditions, joint alteration, and stress reduction due to excavation geometry.

In the context of this study, Stope 2B exhibited the highest tolerable unsupported span, reflecting its superior rock mass conditions and overall stability, as shown in Table 3. This result is consistent with its higher Q-value and RQD, which indicate a competent rock mass with fewer fractures, clean joint surfaces, and minimal water ingress. The ability of Stope 2B to support a larger span without reinforcement suggests that the rock mass is well-interlocked and capable of withstanding the induced stresses from mining without immediate risk of collapse or significant deformation.

In contrast, Stopes 2A and 3C demonstrated considerably lower maximum unsupported spans, signaling reduced stability and a higher probability of requiring artificial support systems. These stopes also recorded lower RQD values (73% and 70%, respectively), indicating more intense fracturing and possibly blast-induced damage or degradation due to groundwater pressure. The reduced span capacity in these areas implies that the rock mass cannot reliably sustain large openings, as discontinuities such as joints, fissures, and weakened zones may serve as planes of failure, especially under the influence of gravity and mining-induced stress changes.

The strong correlation between RQD and the calculated unsupported spans emphasizes the significance of rock mass fragmentation and structural integrity in stope design. RQD, as a measure of the degree of fracturing in the core samples, directly impacts the overall Q-value and, subsequently, the recommended excavation span. This highlights the necessity of accurate geotechnical logging and sampling in the early stages of mine design to ensure safe and efficient stope development.

Furthermore, the variation in span tolerances among the stopes suggests that site-specific ground control strategies are essential, rather than relying on a generalized design approach. While Stope 2B may safely accommodate wider openings, Stopes 2A and 3C would benefit from narrower spans or the implementation of support systems such as cable bolts, mesh, or shotcrete to mitigate the risk of overbreak and ensure worker safety.

Table 3. Maximum tolerable unsupported length of the stopes

3.4. Stope Stability Assessment

To assess the stability of the evaluated stopes, the modified stability number (Q′) and hydraulic radius were calculated, as summarized in Table 4. These parameters are central to the Mathews Stability Graph method, where the stability number (Q′) incorporates key rock mass quality factors such as joint condition, groundwater influence, and joint orientation, while the hydraulic radius reflects the geometry of the excavation.

Table 4. The calculated Modified Stability Number and Hydraulic radius

Among the analyzed stopes, Stope 2B demonstrated the highest Q′ value, indicating the most favourable rock mass conditions overall. Interestingly, this high stability number was maintained despite a relatively low joint orientation factor (B). This suggests that the inherent rock quality (reflected in parameters like intact rock strength, joint spacing, and absence of water) played a dominant role in promoting stability. In other words, while unfavourable joint orientation typically reduces stability, the overall robustness of the rock mass in Stope 2B was sufficient to offset the negative influence of joint orientation. This finding highlights the complex interplay between geological and structural factors in influencing excavation performance, where strong, dry, and tightly interlocked rock masses can compensate for less-than-ideal joint alignments.

However, the analysis also reveals a critical point regarding the mechanism of potential failure. In stopes with lower Q′ values and less favourable geometrical or structural conditions, such as 2A and 3C, failure is likely to initiate along joints that intersect the free face at small angles. These low-angle joints act as potential sliding surfaces, especially when oriented sub-parallel to the excavation walls or roof. Such orientations reduce the shear resistance along the joint planes, particularly when combined with stress relief from excavation or the presence of water, which can lower effective stress and lubricate joint surfaces. This failure mechanism is consistent with observed overbreak patterns in similar mining settings and underscores the importance of joint orientation analysis in ground control planning.

Furthermore, the relationship between Q′ and HR on the Mathews Stability Graph places these stopes within the potentially unstable region, suggesting that while outright failure may not be imminent, there is an elevated risk that warrants attention. In practice, this means that additional support measures or modifications to stope dimensions may be required to maintain stability, particularly in areas where joint geometry and excavation layout intersect unfavourably.

3.4.1. Mathews Stability Graph

The results derived from the Mathews Stability Graph (see Figure 7) reveal that all analyzed stopes plot within the “potentially unstable” zone of the graph. This zone, positioned between the empirically defined stable and failed regions, represents a range of hydraulic radius and stability number (N) combinations where stope performance is highly sensitive to local geological, structural, and operational conditions. It indicates that the analyzed stopes, while not necessarily prone to immediate failure, exhibit a heightened risk of instability and would likely require additional support measures or design modifications to maintain structural integrity.

Figure 7. Stop Stability Results on Mathews’ stability graph

Among the evaluated stopes, Stope 2B demonstrates the highest relative stability, as indicated by its location within the potentially unstable zone at a higher Stability Number (N), positioning it closer to the stable region. This observation is supported by its comparatively higher Q-value, which reflects more favourable rock mass conditions, and a greater RQD, suggesting a higher degree of intact rock within the core samples. The Q-value, which incorporates factors such as joint set number, joint roughness, joint alteration, and groundwater conditions, plays a direct role in the computation of the stability number (N). A higher Q-value in Stope 2B therefore results in a higher N, shifting its plot point further to the right on the stability graph and indicating improved stability.

This trend validates the underlying empirical relationships embedded in the Mathews Stability Graph method. The graphical output aligns well with the field data and quantitative input parameters, especially in the case of Stope 2B, where the correlation between good rock mass quality and predicted stability is clearly evident. However, the fact that all stopes fall within the potentially unstable region (even those with moderate to high rock quality) emphasizes the critical role of stope geometry, particularly the hydraulic radius, in influencing stability outcomes.

3.5. Machine Learning Analysis Using Scikit-Learn

Machine learning analysis was conducted using the Scikit-learn library, a robust and widely adopted framework in Python for building and evaluating machine learning models. The workflow began with data preprocessing, which involved handling missing values, encoding categorical variables, and scaling features to ensure uniformity across the dataset. Following preprocessing, the dataset was split into training and testing subsets to facilitate fair model evaluation and eliminate overfitting. Model evaluation metrics including accuracy, precision, recall, F1-score, were calculated to compare algorithm performance and select the most effective model for the task. Visual tools such as confusion matrices and ROC curves were used to interpret classification results.

3.5.1. Feature Importance Analysis

To understand the relative influence of various geotechnical factors on stope instability, the Random Forest Classifier from the Scikit-Learn library was employed to perform feature importance analysis. The choice of Random Forest was motivated by its ability to handle both categorical and numerical data, its robustness against overfitting, and its capacity to provide feature importance scores based on the mean decrease in impurity.

Prior to model training, categorical variables (water effect and blast damage) were encoded using ordinal encoding, assigning numerical values representing the severity levels (e.g. low=0, medium=1, high=2). Numerical features (RQD, Rock Strength, RMR) were standardized using Scikit-Learn's StandardScaler to ensure consistent scaling.

The feature importance scores, as visualized in Figure 8, revealed that water effect was the most influential factor in predicting stope instability, followed by blast damage and RMR. RQD and Rock Strength exhibited comparatively lower importance. This suggests that the presence of water and the extent of blast damage significantly contribute to instability, potentially by exacerbating existing discontinuities and reducing the effective strength of the rock mass.

Figure 8. Feature importance scores obtained from the Random Forest Classifier

It is important to note that the feature importance scores are relative and should be interpreted within the context of the dataset. While Rock Strength showed low importance, it does not imply it's entirely negligible; rather, its variability within the dataset might be less influential compared to other factors.

3.5.2. Recursive Feature Elimination Analysis

The Recursive Feature Elimination (RFE) process was used to evaluate the importance of each feature in predicting stope stability (see Figure 9). The results showed that water effect, blast damage, and RMR were consistently ranked as the top features, indicating their high importance in predicting stope stability. When these features were removed from the model, the performance dropped significantly, further confirming their importance. The RFE analysis provided further insight into the relationships between the geotechnical factors and stope instability, and supported the findings of the feature importance analysis.

Figure 9. The Recursive Feature Elimination (RFE) process

3.5.3. Decision Boundary Analysis

To further explore the interaction between the two most significant factors, water effect and blast damage, we performed a decision boundary analysis. This technique visualizes the classifier's predictions across the feature space, illustrating how different combinations of water effect and blast damage influence stope stability.

Figure 10. Water Effect against blast Damage Decision Boundary Plot

Figure 10 shows the decision boundary generated by the Random Forest Classifier. The yellow area represents the region where the model predicted "Stable" stope conditions, while the purple area represents "Unstable" conditions.

The analysis revealed that:

Water effect is a dominant predictor: when water effect is low (0) or medium (1), the classifier predicts "Stable" regardless of blast damage severity. This suggests that in the absence of significant water influence, even high levels of blast-induced damage are unlikely to compromise stope stability.

High water effect increases instability risk: when water effect is high (2), the classifier predicts "Unstable" for medium (1) and high (2) blast damage levels. This indicates that elevated water influence substantially increases the risk of instability, particularly when combined with blast damage.

Interaction effect: the decision boundary analysis underscores a strong interaction between water effect and blast damage. Specifically, a high-water effect amplifies the destabilizing impact of blast damage, pointing to a compounded influence on stope stability.

3.5.4. Data Distribution Analysis

Feature distribution plots (see Figure 11) reveal a positive correlation between rock strength and RMR, indicating that higher-strength rocks are generally associated with better rock mass quality. Additionally, a positive correlation was observed between water effect and blast damage, suggesting that the presence of water may exacerbate the extent of blast-induced damage. This distribution analysis is critical for understanding the interrelationships among geotechnical factors and plays a vital role in informing the design of appropriate ground support systems.

Figure 11. Feature distribution plots for the five factors under investigation

3.5.5. Correlation Coefficient Analysis

The correlation coefficient analysis shows a strong negative linear relationship between water effect and rock quality, and a moderate negative linear relationship between blast damage and rock quality (see Figure 12). This suggests that the presence of water is associated with a decrease in rock quality, and that blast damage is also associated with a decrease in rock quality. The correlation coefficient analysis is essential in understanding the relationships between the factors and designing appropriate support systems.

Figure 12. Correlation Heatmap

3.5.6. Model Evaluation and Limitations

The performance of the Random Forest Classifier was evaluated by using multiple metrics, achieving an accuracy of 0.83, precision of 0.88, recall of 0.83, and an F1-score of 0.83, as illustrated in Figure 13. To assess the model’s robustness and its ability to generalize to unseen data, a 5-fold cross-validation procedure was conducted. In this process, the dataset was randomly partitioned into five equal subsets; in each iteration, four subsets were used for training and the remaining one for testing, ensuring that each subset served as the test set exactly once. The performance metrics averaged across the five folds were consistent with those obtained from the initial model evaluation. This consistency indicates that the model does not suffer from overfitting and demonstrates reliable generalization capability. Such validation strengthens confidence in the model's predictive reliability when applied to new or unseen geotechnical datasets.

Figure 13. Classification report extract

It is important to acknowledge several limitations associated with this analysis. First, the dataset used for training the machine learning model was relatively limited in size, which may constrain the model’s ability to fully capture the underlying complexity and non-linear interactions among geotechnical features. This limitation could affect both the generalizability and robustness of the predictive outcomes. Second, while categorical encoding was employed to represent variables such as water effect and blast damage, this approach simplifies the inherently complex and continuous nature of these factors. Although effective for initial modelling purposes, such simplification may lead to a loss of nuanced information. Future research could benefit from the use of more advanced encoding methods, such as ordinal encoding, one-hot encoding with domain-specific thresholds, or even continuous variable modelling, to better reflect the true variability and interactions of these parameters.

Furthermore, the data used in this study were collected exclusively from the 2North Section of the mine. As a result, the model's applicability may be limited to this specific geological and operational context, and its predictions may not be directly transferable to other sections of the mine with differing geotechnical or hydrological conditions. To enhance the model’s predictive power and generalizability, future studies should consider expanding the dataset, incorporating additional features such as joint orientation, stress regime, or time-dependent effects, and applying the model across multiple zones of the mine. These improvements would allow for a more comprehensive assessment of stope stability and broader applicability of the model outcomes.

Figure 14. Receiver Operating Characteristic (ROC) Curve for Logistic Regression

Figure 14 displays the Receiver Operating Characteristic (ROC) curve for the trained Logistic Regression model. The ROC curve plots the true positive rate (sensitivity) against the false positive rate (1 - specificity) at various classification thresholds. The Area Under the Curve (AUC) provides a single scalar value summarizing the overall performance of the classifier for different classes. The AUC values for Class Potentially unstable, Class stable, and Class unstable are 0.545, 0.403, and 0.729, respectively.

The ROC curve and AUC values provide insight into the model's ability to balance true positives and false positives, which is crucial for practical applications in mining engineering. The results indicate that the model performs relatively well in predicting "unstable" stopes (AUC = 0.729), suggesting that it can effectively identify instances that are likely to be unstable. However, the model's performance for "stable" stopes is poor (AUC = 0.403), indicating that it struggles to distinguish between stable and other classes. The model's performance for "potentially unstable" stopes is moderate (AUC = 0.545), suggesting that it can identify some instances that are potentially unstable, but with limited accuracy.

The varying AUC values across classes suggest that the model may be biased towards certain classes or that the features used to train the model are not equally informative for all classes. To improve the model's performance, it may be necessary to collect more data, particularly for the classes with lower AUC values, or to explore alternative feature engineering strategies. Additionally, hyperparameter tuning or the use of ensemble methods may help to improve the model's overall performance.

4.1. Stope Failure Analysis

Rock falls with medium and high hazard intensity are dominating stopes 2A and 3C, indicating the failure of support systems. The support system consists of natural in-situ pillars left randomly in stopes, relying on the excellent rock mass quality of the stopes. The adopted support design from 2South Section, had a maximum tolerable unsupported length of 7.5m which deviated from the calculated values. The maximum tolerable unsupported length for stope 2A and 3C deviated negatively with a magnitude of 2.5m and 1.5m respectively. This entails the deterioration of the rock mass quality in these stopes. As the rock mass quality was altered, the stope shape factor became insignificant for the design thus the stope stability decreased. As a result, new support requirements which may include artificial support are now necessary.

4.2. Support Demand and Pattern

The analysis of support demand and support pattern indicates that the proposed support configuration is adequate for maintaining the stability of the stopes, as detailed in Tables 5 and 6. The support demand was quantitatively estimated based on established rock mass classification systems (RMR and Q-system) and evaluated using the Mathews’ stability graph method, which integrates rock quality, stope geometry, and excavation span to determine required support levels.

The proposed support pattern was designed to ensure an even distribution of load and to accommodate localized stress concentrations, particularly around larger load-bearing pillars. This design approach considers both the mechanical properties of the rock mass and the empirical stability zones defined in the stability graph. The pattern aims to transfer and distribute loads effectively, minimize deformation, and prevent progressive failure, thereby ensuring the long-term structural integrity of the stopes under anticipated mining conditions.

Table 5. Fall out heights and support demand of the stopes under investigation

The support demand was calculated using Equation 24.

[TeX:] \text{Support\ Demand\ per\ }m^{2} = \rho\mathsf{\text{gh.}}$\mathrm{\text{Support Demand per}}{m}^2=\rho \mathrm{\text{gh.}}$ (24)

where ρ is the verage density of the overlying rock mass, g is the acceleration due to gravity and h is the Fall-out height.

Table 6. The tributary area for rock fall conditions and the corresponding support spacing