A semi-empirical solution for estimating the elastic stresses around inclined mine stopes for the Mathews-Potvin stability analysis

The Mathews-Potvin stability method is widely used in the Canadian mining industry as a starting point to determine the maximum dimensions of mine stopes. However, it cannot be applied to inclined (more frequently encountered) mine stopes without conducting numerical modelling to obtain the stress factor A, defined as a function of the ratio of unconfined compressive strength of intact rock to the induced principal stress on the exposed stope walls. The need to conduct numerical modelling significantly limits the application of the Mathews-Potvin method. In addition, given its empirical nature and main application for preliminary design, it is deemed undesirable to conduct numerical modelling, especially elaborate modelling. Alternatively, theoretical methods can provide a much simpler and quicker way to estimate stresses around stopes and the corresponding stress factors. Over the years, a large number of studies have been conducted to estimate stresses around openings excavated with various crosssections. However, theoretical or graphical solutions remain unavailable for mine stopes that typically consist of horizontal floor and roof, and two parallel inclined walls (hangingwall and footwall). To remedy this situation, a series of numerical simulations is first performed for openings with vertical and inclined walls, including typical stopes commonly encountered in underground mines. A group of empirical solutions is then formulated to estimate the induced principal stresses at the roof centre and mid-height of the stope walls. The validity and predictability of the proposed solution have been verified using additional numerical simulations. The proposed solution can thus be used to calculate stresses and the resultant stress factors A around typical mine stopes with any inclination angle and height to width ratio, under any in-situ stress state, without conducting numerical modelling.


Introduction
Ground stability is a challenging issue frequently faced by rock engineers. The trend towards larger and more powerful equipment to improve productivity requires larger underground openings. However, the dimensions of stable underground excavations are finite, limited by field stresses and rock mass conditions. The correct design of underground openings is thus of paramount importance.
The Mathews-Potvin method is a simple and useful tool for mining engineers. It is commonly used as a starting point to determine the dimensions of stopes or design the required support (e.g., Mathews et al. 1981;Potvin 1988;Hutchinson and Diederichs 1996;Li and Ouellet 2009). The Mathews-Potvin method is also used to estimate the unplanned dilution due to the slough that can take place around the hangingwall and footwall during blasting or muck-out of blasted ore (Scoble and Moss 1994;Clark and Pakalnis 1997;Kaiser et al. 1997;Kaiser 1999, Diederichs, Kaiser, andEberhardt, 2004;Papaioanou and Suorineni 2016). Another application of the Mathews-Potvin method is to estimate the minimum span exposures to ensure the cavability (self-collapse) of ore in caving mining methods (Sunwoo, Jung, and Karanam, 2006).
For mine stope design, a major limitation associated with the Mathews-Potvin method is the need to conduct numerical simulations to obtain a key parameter, called the stress factor (A), which depends on the ratio of the unconfined compressive strength of intact rock to the induced principal stress (σ 1 ) on the walls of the opening. When the geometry of the openings is simple, such as a circular cross-section, analytical solutions exist for estimating the stress around such openings (Kirsch, 1898;Hiramatsu, 1962;Logie and van, Tonder 1967;Hiramatsu and Oka, 1968;Li, 1997). More sophisticated analytical solutions are equally available for estimating the elastic stresses around tunnels with a vertical axis of symmetry, including openings with elliptical, rectangular, and arched walls (Logie and van Tonder 1967;Hoek and Brown 1980;Gerçek 1997;Exadaktylos and Stavropoulou 2002;Brady and Brown 2013).
For vertical mine stopes, graphical solutions used to estimate the induced stresses have been elaborated in two dimensions (2D plane strain ;Potvin 1988;Stewart and Forsythe 1995) and three dimensions (Mawdesley, Trueman, and Whiten, 2001;Vallejos, Delonca, and Perez, 2017). In practice, however, orebodies are always inclined to a greater or lesser degree. Few theoretical or graphical solutions are available to assess the induced stresses around inclined stopes with one horizontal floor, one horizontal roof, and two parallel and inclined walls (hangingwall and footwall). Numerical modelling has to be performed to obtain the induced stresses for each specific mining project (Li and Ouellet, 2009). This requires not only the availability of pertinent software and hardware, but also qualified numerical modellers who have a good understanding of field conditions and the behaviour of the rock mass, and know in particular how to obtain reliable numerical outcomes. It is thus desirable to have theoretical solutions that can be used to estimate the stresses around inclined stopes.
In this paper, the Mathews-Potvin method is first briefly recalled for the sake of completeness. Numerical simulation results are then presented by considering inclined mine stopes surrounded by a homogenous, isotropic and linearly elastic rock mass. A large range of wall inclination angles, height to width ratios, and in-situ stresses are considered. Semi-empirical solutions are proposed to estimate the induced principal stresses on stope walls by applying the principle of superposition of linearly elasticity theory through a curve-fitting technique applied to the numerical results. The prediction capability of the proposed semi-empirical solutions is verified with additional numerical simulations. A typical example is also given to illustrate the application of the proposed solution.

The Mathews-Potvin method
The Mathews-Potvin method is an empirical method based on numerous field observations. This method relates the stability of an exposed wall to two factors -the hydraulic radius (HR) and the stability number (N'). The former is defined as (Potvin 1988): [1] The stability number (N') of the exposed wall is defined by the following equation: where Q' is a modified rock tunnelling quality index, A is the rock stress factor, B is the joint orientation adjustment factor, and C is the gravity adjustment factor.
The parameter Q' resulted from a modification on the Rock Tunnelling Index (Q) of Barton, Lien, and Lunde (1974), is defined as follows: where RQD is the rock quality designation, J n is the joint set number, J r is the joint roughness number, and J a is the joint alteration number.
The rock stress factor A is a function of the ratio between the unconfined compressive strength of the intact rock (σ c ) and the induced major principal stress (σ 1 ) on the exposed walls of a stope. The stress factor A can be expressed as follows (Potvin, 1988), which is further illustrated in Figure 1a: Factor B considers the influence of joints on the stability of the studied exposed wall (Potvin 1988). It represents the effect of the angle between the most critical joints and the wall, as shown in Figure 1b. The gravity factor C depends on the individual influence of the inclination of the exposed wall and the inclination of the critical joints, as illustrated in Figure 1c. Once the stability number N' and hydraulic radius HR are determined, the stability of the exposed wall can be evaluated using the chart of Mathews-Potvin, as shown in Figure 1d.
From Equation [4], one notes that the rock stress factor A proposed by Potvin (1988) has some limitations when the rock is submitted to a tensional stress. In this case, the induced principal stress σ 1 is zero (for a 2D model) or non-zero in the third dimension (for a 3D model). Factor A can thus reach its maximum value of 1.0 independent of the tensile stress and tensile strength of the rock, which is unrealistic. Therefore, Equation [4] is not entirely adequate to describe the stability or failure of rock by tension. To overcome this limitation, Li and Ouellet (2009) proposed two approaches. The first is to neglect the tensile strength of the rock, and to fix A = 0.1 for σ 3 ≤ 0 (where σ 3 is the induced tensile stress around the excavation). The second approach is to compare the tensile stress with the tensile strength of the rock, so that A = 0.1125| σt ⁄ σ3 |-0.125 (same form as Equation [4]; where σ t is the tensile strength of the intact rock). Zhang, Hughes, and Mitri (2011) adopted a similar approach to that of Li and Ouellet (2009) when the rock is submitted to tension. Suorineni (2012) concluded that the stress factor for tension (and other factors) needs to be calibrated. Further discussion on the definition of this factor is beyond the scope of the paper, but it is seen that the determination of factor A depends on knowledge of the induced principal stresses on the exposed stope walls. roof, a horizontal base and two parallel and inclined walls. In the figure, W and H are the width and height of the stope, respectively; b is the inclination angle of the stope walls; σ v and σ h on the stress block represent the vertical and horizontal natural in-situ principal stresses, respectively; the out-of-plane stress is another in-situ principal stress. The mid-points have been denoted as U and V on the surfaces of roof and sidewall respectively. The numerical code Plaxis 2D (Brinkgreve and Vermeer 1999), based on the finite element method and commonly adapted for rock mechanics and geotechnical engineering, is used here to evaluate the stresses around mine stopes. The sign convention used by Plaxis 2D considers compression negative (-) and tension positive (+). However, the results presented in this study follow the sign convention commonly used in rock mechanics analysis, where compression is positive (+) and tension is negative (-).

Numerical modelling of the elastic stresses around inclined stope walls
The linearly elastic model of Plaxis 2D was first validated by comparing the simulated stresses against the analytical solutions for a circular opening (Kirsch, 1898;Hiramatsu, 1962;Hiramatsu and Oka, 1968;Li, 1997). Additional validations were made against the graphical and analytical solutions of Hoek and Brown (1980) in the cases of elliptical and square openings. More details are presented in Pagé (2018). Table I presents the program of numerical simulations. Forty-eight stope geometries were considered by combining the stope width (W), height (H), and wall inclination angle (b). Two regimes of natural in-situ stresses were considered: Case 1 with σ v = 30 MPa and σ h = 0; Case 2 with σ v = 0 and σ h = 30 MPa. It should be noted that the consideration of a zero horizontal insitu stress σ h in Case 1 and a zero vertical in-situ stress σ v in Case 2 is necessary for applying the principle of superposition. This does not mean that the numerical models only correspond to zero vertical or horizontal in-situ stress, although the models remain valid for such extreme cases. The principal stresses tangential to the exposed faces at points U and V (on the wall surfaces) are calculated (see Figure 2). Figure 3 shows a numerical model constructed with Plaxis 2D. An enlarged view of the stope with refined meshes around the stope before excavation is presented. The natural in-situ stresses were first initiated over the entire model. The four outer boundaries were then fixed in all directions. Finally, the excavation of the stope was simulated. For each numerical model with a new stope geometry, domain and meshes sensitivity analyses have been done to ensure that the outer boundaries are far enough from the stope and the meshes around the stope are fine enough. A sufficiently large domain is necessary to avoid the boundary effect, while finer meshes around the stope are required to ensure stable numerical results (see more details presented in Pagé, 2018). Figure 4 presents the minor (σ 3 , Figure 4a) and major (σ 1 , Figure 4b) principal stresses contours around a stope with H/W = 2 and b = 75°, obtained from numerical modelling using an insitu stress state of σ v = 30 MPa and σ h = 0 MPa (with σ v = -30 MPa and σ h = 0 MPa as inputs to Plaxis 2D). Note that the major and minor in-plane principal stresses in Plaxis 2D are represented by σ 1 and σ 3 respectively, while the out-of-plane  Figure 4a shows that the critical tangential stresses on the roof are under tension (positive in Plaxis 2D), while Figure 4b indicates that the critical tangential stresses on the walls undergo compression (negative in Plaxis 2D). The minor principal stress at the roof centre is -26.8 MPa (in tension), while the major principal stress at the mid-height of the hangingwall and footwall is 37.4 MPa (in compression). Figure 5 shows the major (σ 1 , Figure 5a) and minor (σ 3 , Figure 5b) principal stress contours around the stope with H/W = 2 and b = 75°, obtained by numerical modelling with a natural in-situ stress state of σ v = 0 MPa and σ h = 30 MPa (with σ v = 0 MPa and σ h = -30 MPa as inputs to Plaxis 2D). In this case, the critical tangential stress on the roof is under compression (negative in Plaxis 2D) based on the major principal stress (σ 1 , Figure 5a), while the critical tangential stresses on the walls are under tension (positive in Plaxis 2D) based on the minor principal stress (σ 3 , Figure 5b). The major principal (compressive) stress at the roof centre is 61.1 MPa, while the minor principal (tensile) stress at the mid-height of the hangingwall and footwall is -23.6 MPa.

Formulation
To formulate a semi-empirical solution for evaluating the elastic stresses around mine stopes, one uses the principle of superposition valid in elasticity theory for homogenous, isotropic, and linearly elastic material. For a given stope geometry, the stresses around the opening are investigated by applying a horizontal natural in-situ stress. The induced stresses at the point of interest on the stope wall are then normalized by the applied horizontal natural in-situ stress. By changing the stope height to width ratio (H/W) and wall inclination angle (b), a relationship based on curve fitting can then be established between the studied stresses at the point of interest on the stope wall and the horizontal natural in-situ stress, stope width to height ratio, and stope wall inclination angle. The same process is repeated for the vertical natural in-situ stress with different stope width to height ratios and stope wall inclination angles. Applying the curve-fitting technique leads to another equation, which describes the induced stress around the stope opening as a function of the vertical natural in-situ stress, stope width to height ratio, and stope wall inclination angle.
By adding the two equations, one obtains an equation that describes the studied stresses at the point of interest on the wall or roof as a function of the horizontal and vertical natural in-situ stresses, stope geometry, and wall inclination. The procedure can be summarized as follows: [5] [6] where f 1 and f 2 are the geometric functions on the critical tangential stress at the roof centre, associated with the vertical and horizontal natural in-situ stresses, respectively; g 1 and g 2 are the geometric functions on the critical tangential stress at the mid-height of the hangingwall and footwall, associated with the vertical and horizontal natural in-situ stresses, respectively.
To obtain the four geometric functions f 1 , f 2 , g 1 , and g 2 , a second-degree polynomial regression fit (for both f 1 and f 2 ), and a combination of a power regression fit and a second-degree polynomial regression fit (for g 1 and g 2 respectively) were applied to the numerical results of the critical induced stresses at the roof centre and at the mid-height of the wall as a function of the H/W ratio, separately for b = 90°, 75°, 60°, and 45°. Figure  A second calibration of these four geometric functions by considering the wall inclination angle leads to the following equations: Equations [5] to [10] constitute the proposed solution for estimating the elastic stresses at the roof centre and midheight of the hangingwall and footwall around typical mine stopes. These equations are independent of the stope depth and rock mass strength. [7] [8] [9]

Figure 5-Isocontours of σ 1 (a) and σ 3 (b) principal stresses around the stope with H/W = 2 and β = 75°, calculated by applying a natural in-situ stress state of σ v = 0 MPa and σ h = 30 MPa in Plaxis 2D (with σ v = 0 MPa and σ h = -30 MPa as inputs tp Plaxis 2D)
A semi-empirical solution for estimating the elastic stresses around inclined mine stopes ▶ 410 AUGUST 2021 VOLUME 121 The Journal of the Southern African Institute of Mining and Metallurgy [10] When the natural in-situ stress state is σ v > 0 (in compression) and σ h = 0, the solution predicts tension (σ roof < 0) acting on the roof and compression (σ wall > 0) on the midheight of the walls. Conversely, when the natural in-situ stress state is σ v = 0 and σ h > 0 (in compression), the solution leads to compression (σ roof > 0) on the roof and tension (σ wall < 0) on the mid-height of the walls. Figure 6 shows that the critical induced stresses at the roof centre and at the mid-height of the hangingwall and footwall

respectively normalized by the applied horizontal (σ h ) and vertical (σ v ) in-situ stresses as a function of the H/W ratio for stope inclination angles (β) of (a) 90°; (b) 75°; (c) 60°; and (d) 45°
The Journal of the Southern African Institute of Mining and Metallurgy VOLUME 121 AUGUST 2021 calculated by the proposed solution (Equations [5] to [10]), and represented by the full lines, correspond well to those obtained by the numerical modelling. This type of comparison between a proposed solution and numerical (or experimental) results, which is used in the calibration or curve fitting to obtain the required parameters, is usually considered as validation or prediction. This is, however, not true. The validity and predictability of the calibrated model (obtained by calibration or curve filling) should be tested against additional and different numerical (or experimental) results.

Validation and predictability
To test the validity and predictability of the proposed solution, additional numerical simulations were performed by considering more stope geometries and virgin in-situ stress states. Figure 7 shows the variation of the induced tangential stresses, obtained by numerical modelling and predicted by the proposed semi-empirical solution, at the roof centre and at midheight of the walls, for an isotropic natural in-situ stress state of 30 MPa (compression) and stopes having wall inclination angles b = 90°, 75°, 60°, and 45° and H/W ratio varying from 0.1 to 10. It is seen that the agreement between these two different approaches is excellent. Figure 8 presents another validation and test of predictability of the proposed semi-empirical solution using additional numerical simulations conducted with anisotropic in-situ stresses  (Figure 8f). The proposed solution can thus be considered as validated. It can then be used to calculate the stresses and the stress factor A for the case of typical mine stopes with any inclination angle and height to width ratio under any in-situ stress state.

Sample application
In the following, a sample calculation is presented to further illustrate the application of the proposed solution (Equations [5] to [10]).
It is planned to mine out a 6 m (W) wide ore vein inclined at 67° (b), with a 30 m (H) high stope located 500 m below the ground surface. The vertical in-situ stress can be estimated based on the overburden depth, while the horizontal in-situ stress is 1.67 times the vertical in-situ stress. These parameters give the following in-situ stress state and stope geometry:

Discussion
Numerical modelling requires the availability of pertinent software and hardware and qualified modellers who know how to correctly conduct numerical modelling. Currently, the availability of computation resources in terms of hardware and software is no longer an issue, and numerical modelling has become a common practice for various research and design projects. However, knowing how to use a numerical code is often considered equivalent of knowing how to correctly perform numerical modelling. This can partly explain the crisis of confidence in numerical modelling and why many modellers do not believe in even their own numerical results. In fact, knowing how to use a numerical code is totally different from knowing how to conduct numerical modelling. The former needs only short training (a couple of hours) while the latter requires much more advanced training and rich experience in order to obtain stable and reliable numerical outcomes (Chapuis et al., 2001; Barbour and Krahn 2004;Cheng, Lansivaara, and Wei, 2007;Diederichs et al., 2007;Krahn 2007;Chapuis, 2012aChapuis, , 2012bDuncan 2013).
This work is partly motivated by a perception that the Mathews-Potvin method was considered useless for stope analysis, and that the stability and the maximum dimensions of the stopes can be directly analysed using numerical models, instead of determining the stress factor A and applying the Mathews-Potvin method. It should be recalled that the empirical Mathews-Potvin method was based on many case study observations. The numerical models performed to determine the stress factor A are very simple, considering only an isolated opening around a homogenous, isotropic, and linearly elastic rock. The effectiveness of the method has been proven, especially when it is used as a starting point for the determination of stope dimensions in the preliminary stage of mining projects. When numerical modelling is conducted for stope stability analysis, the models are usually much more complex in terms of stope geometry, mining sequence, and material parameters. Calibrations using field data/observations can be necessary to find the required (but unknown) parameters. In the preliminary stage of a mining project, little field data and information are available to allow the construction and calibration of such sophisticated numerical models.
All of these considerations indicate that the Mathews-Potvin method is very useful at the beginning of a mining project, where it can provide a quick and preliminary estimation of the dimensions of stopes. The necessity for more sophisticated numerical modelling at the advanced stage of a mining project does not invalidate the Mathews-Potvin method. Rather, the Mathews-Potvin method can be more appealing if theoretical or graphical solutions are available for estimating the induced stresses around inclined mine stopes. To this end, a semiempirical solution has been proposed in which curve-fitting techniques are applied against numerical modelling and the principle of superposition of linearly elasticity theory. The results show that the proposed semi-empirical solution can be used to evaluate the induced principal stresses at the roof centre and mid-height of the wall around typical mine stopes. However, one should keep in mind that the numerical models presented in this study contain several assumptions.
First, a limitation of the numerical models is associated with the 2D plane strain conditions. The numerical results and the proposed semi-empirical solution are valid only when the stope is very long in one horizontal direction. In an actual mine, this is not always the case. Graphical solutions have been presented by Mawdesley, Treman, and Whiten (2001) and Vallejos, Delonca, and Perex (2017) for 3D vertical stopes. Further work is necessary to consider three-dimensional inclined stopes.
The assumption of a linearly elastic rock mass can be held true at relatively shallow depths. At greater depths (deep mines), the behaviour of rocks and rock masses may change to a nonlinear and non-elastic behaviour. Consequently, the validity of the empirical relationships proposed here may be limited to a certain depth. Additional studies could be conducted to formulate similar empirical relationships in nonlinear rock masses.
Another limitation of the proposed semi-empirical relationships is related to the stope geometry. The stopes considered here have a parallel hangingwall and footwall as well as parallel roof and floor. In practice, stopes with nonparallel walls are commonly encountered. More work is needed to propose solutions for estimating the stresses around stopes with nonparallel walls.
In this study, the vertical and horizontal in-situ stresses were considered to be two principal stresses, implicitly assuming that the out-plane in-situ stress is another principal stress. In practice, the vertical in-situ stress and the two horizontal in-situ stresses could form three normal stresses. Further work is thus necessary to develop a more general solution.
Finally, it is very important to point out that the stress factor A defined in the Mathews-Potvin method corresponds to the maximum induced principal stresses on the exposed faces. However, as shown in Figures 4 and 5, the maximum compressive stresses are close to the four corners rather than at the roof centre. We believe that an accurate estimation of the maximum principal compressive stress at stope corners is difficult and unnecessary due to stress concentration -therefore the critical locations in terms of compression should be at the centre, not the stope corners. For tension, as Figure 5b shows, the largest tensile stresses are located near (but somehow distant from) the stope corners, which correspond to the critical locations (rather than the roof or wall centre). In this study, the maximum tensile stress is not considered as its location varies when the stope geometry or natural in-situ stresses change. This renders the formulation very difficult. More work is needed on this aspect. Nonetheless, given the empirical nature of the Mathews-Potvin method and the still limited considerations of the tensile stresses in applying the method, the proposed solutions can provide useful estimation of stresses for application of the Mathews-Potvin method.

Conclusions
The well-known Mathews-Potvin method is an important design tool for mining engineers. However, the application of this method requires the determination of the induced stresses (and stress factor A) around inclined mine stopes using numerical modelling, as few graphical or theoretical solutions are available for such purposes. To overcome this limitation, a semi-empirical solution has been proposed to estimate the induced principal stresses at the roof centre and mid-height of the walls around mine stopes, by applying the superposition principle of linearly elasticity theory and curve-fitting techniques to numerical results. The validity and the predictive capability of the proposed solution have been verified by additional numerical simulations. The proposed semi-empirical solution can thus be used to evaluate the induced tangential stresses at the roof centre and mid-height of the walls around mine stopes with any inclination angles, height to width ratios, and in-situ stress states, as long as the values of H/W are in the range from 0.1 to 10 and b in the range from 45° to 90°. With these empirical expressions, the stress factor A, a key parameter used in the Mathews-Potvin method, can be determined without conducting numerical modelling.