Full Use of Hydrogeological-Geomechanical Parameters for 3d Stability Analysis in an Open-Pit Iron Mine with Complex Geometry

Vidal Félix Navarro Torres

Mining, Vale Institute of Technology, Belo Horizonte, Brazil.

DOI: 10.4236/gep.2025.132010                

Abstract

This study presents a novel method to evaluate the safety of open-pit slopes by means of three-dimensional numerical modeling with the finite difference method. The method presented here uses a block model as a vehicle to carry relevant information from the rock mass and automatically construct the numerical model. The results suggest that the method is promising because of its capacity to accurately incorporate a large amount of high-complexity rock data by considering spatial location and material behavior. It is expected that the innovations in this method will make the design, construction, and operation of open-pit iron mines safer and more economical.

Keywords

Open-Pit Mine, Hydrogeology, Geomechanics, Numerical Modeling, Finite Difference Method, 3D Slope Stability

1. Introduction

Small changes in the global slope angle of a pit can significantly influence the economy of the entire mining operation (Ning et al., 2011); however, any change must ensure the global stability of the slope, ensuring the safe and effective movement of personnel, equipment, and materials during mining operations (Zevgolis et al., 2018; Steffen et al., 2008; Wyllie & Mah, 2004). Because open-pit mines can reach great depths (Hustrulid et al., 2013), global disruption events are expected to occur. The global slope angle and the height of an open pit mine can be related to the failure mechanisms that may occur. Sjöberg (1999) presented the relationship in Figure 1.

Regarding the probability of failure, maximizing the slope angle is also considered

Figure 1. Relationships between the height and angle of pit slopes and associated failure mechanisms (Sjöberg, 1999).

when studying pit optimization because although the annual probability of failure occurrence is relatively high, on the order of 10⁻¹ to 10⁻², the number of fatalities associated with such failures is relatively low, less than 0.1 per event, as shown in Figure 2. In this regard, efforts to predict geotechnical behavior, given the maximization of the global slope angle, are essential to ensure economically profitable and safe operations.

Wyllie and Mah (2004) indicated that the behavior of mining slopes could be evaluated using 1) limit equilibrium analyses, where the slope safety factor is calculated for different failure mechanisms, the most common being planar failure in wedge, circular and toppling analyses, and 2) numerical analyses, where the stresses and deformations developed in the slope are compared against the resistance of the rock, using sensitivity analysis or probabilistic methods to evaluate a wider range of conditions and parameters that describe the strength of the rock.

Singh and Ghose (2006) considered that a better representation of rock behavior is achieved by 1) physical models, where both the dimensions and the strength and deformation properties are reproduced on a reduced scale to real conditions to identify the failure constraints on the extensive slope of the slope; and 2) analytical solutions, where mathematical models of closed solutions are selected to test their accuracy.

Ning et al. (2011) noted that the method chosen for an ideal project should be able to predict the location, conditions and mechanisms of an eventual slope failure, although the limiting equilibrium analysis, numerical analysis, physical models and analytical solutions may be inadequate for some sites or tend to oversimplify the conditions generating excessively conservative solutions.

Figure 2. Risks associated with engineering projects (Wyllie & Mah, 2004).

He et al. (2013) considered that practical problems of rocky slope stability are intricate in terms of the topological geometry, anisotropy and nonlinear behavior of materials, in situ stresses and several processes occurring in a coupled way, suggesting the adoption of numerical modeling as a tool to generate solutions closer to reality. However, these authors admitted that numerical modeling for rock mechanics analysis is underdeveloped in relation to the engineering demand for solving practical problems related to rocky slope stability. Geological-geotechnical parameters can increase the representativeness of the rock massif behavior.

According to Franco (2006), the great advantage of numerical modeling over other methods lies in the discretization process, which reduces a continuous problem, with an infinite number of variables, into a discrete problem, with a finite number of variables and a finite number of points belonging to the domain, making the problem computationally viable. The main numerical methods used to solve these equations are the finite element method, finite difference method, boundary element method and finite volume method (Chapra & Canale, 2016). The choice of numerical method should consider the limitations, advantages and disadvantages of each method. Thus, for the analyses presented in this study, the method chosen was that of finite differences that act directly on the differential equation of motion.

The main parameters used to carry out the design, construction and operation of mine slopes are obtained from activities such as the execution of geological-geotechnical mapping, the geological-geotechnical description of the drill cores, the application of the rock mass classification system, the elaboration of the geomechanical model, the acquisition of geotechnical parameters of intact rocks, the estimation of the stress state, the analysis of slope stability and the monitoring of groundwater and surface deformation of the slopes.

In this work, we propose a method that automatically incorporates relevant information using a mine block model as a vehicle. The results show that the developed method is promising due to its ability to accurately incorporate a large amount of highly complex spatial and behavioral data.

2. Methodology

The finite difference method (FDM) involves solving differential equations, which are generated by studying cases of engineering situations. The MDF approximates each term in the node mathematical model of a mesh through algebraic expressions.

The software for analyzing the behavior of the rock mass based on finite differences was FLAC3D, as it is suitable for solving the vast majority of problems related to geotechnical engineering. The fast Lagrangian analysis of continua in three dimensions (FLAC3D) simulates the three-dimensional behavior of structures built in soil, rock or other possible materials subjected to plastic flow when its failure limits are reached. The materials are represented by polyhedral elements within a three-dimensional grid that is adjusted to be compatible with the object to be modeled (Itasca, 2012).

A system of algebraic equations that relates the displacements and nodal forces is constructed from the functions that describe the behavior of the elements. However, this system of equations is solved by dynamic relaxation, an explicit procedure as a function of time in which the dynamic equations of motion are integrated at each cycle advance. Static solutions are obtained by including amortization factors that gradually remove kinetic energy from the system (Michalowski & Dawson, 2002).

Three-dimensional modeling in FLAC3D is achieved through its internal programming language, FISH, which creates functions capable of simulating the process of excavation of slopes and monitoring and extracting the desired variables. Its formulation can accommodate analyses with large displacements, deformations and nonlinear behavior of the material even if it exceeds the levels of failure, yield or even collapse. It also contains a programming language that makes it possible to take advantage of the program’s capabilities.

The geomechanical model of the rock mass simulated in FLAC3D considers that the resistant behavior is represented by the generalized Hoek‒Brown (HB) criterion. However, it is known that this failure criterion is isotropic and consequently cannot represent the anisotropic behavior imposed by the presence of structural geological features such as foliations. To consider this type of feature, the ubiquitous joints model was used. Since the base criterion is the generalized Hoek‒Brown criterion, the consideration of a ubiquitous joint model requires the estimation of equivalent Mohr-Coulomb parameters for the isotropic continuum and of the Mohr-Coulomb parameters for the foliations. The parameters of the foliations can be found using the Barton-Bandis criterion, which has parameters that can be more intuitively correlated with geological surveys.

The mine block model was used as a database because it contains lithologies (Figure 3), in addition to allowing the placement of relevant information such as density (Figure 4), location of persistent faults (Figure 5), geomechanical classification (Figure 6), and pore pressure (Figure 7).

Figure 3. Mine block model.

Figure 4. Density.

Figure 5. Fault location.

Figure 6. Geomechanical class.

Figure 7. Pore pressure.

The geometric mesh of the FLAC3D numerical model was constructed by emulating the mine block model. In this process, each zone of FLAC3D was constructed with the same block characteristics, i.e., location, size, geomechanical class (slot 1), density (slot 2), lithology (slot 3), failure indicator (slot 4) and structural sector (slot 5). This was performed automatically through a FISH function that reads the data exported by Vulcan into an *.txt file, which was later imported into FLAC3D and built individually. The mapping of a fault was included in the block model through an indicator variable: i_falha = 1 when the block is cut by the fault and i_falha = 0 when it is not cut by the fault (Figure 5). This problem is highly complex and demanding in terms of computational time.

With all the data entered in the block model and the failure criteria defined, a 3D stability analysis of the iron mine was performed.

3. Routines Generated for Model Optimization

3.1. Meshing Routine

An innovative FLAC3D geometric meshing routine was generated. This routine dramatically reduced the import time from 15 days to a few minutes. This is due to the use of programming in the Python language. The routine consists of building the FLAC3D mesh file synthetically, including all the desired characteristics in the form of zone groupings. The object-oriented routine has the zone class as its main object, which has the following attributes: global index, IJK, centroid, vertex coordinates, lithology, fault index, and GSI.

The reading process (read_blockmodel() function) consists of creating a list of zone objects by filling in and calculating the attributes. For this, the block model must be in text format such as .csv or .txt. The spatial distribution of the GSI was determined externally using the SGeMS program (Stanford Geostatistical Modeling Software). The interpolation of GSI values was performed over the grid of the block model prototype; that is, each IJK presented an estimated value.

Subsequently, the zones are classified and grouped according to the various required groupings (group_materials()), such as lithology and classification. These clusters are finally saved in the FLAC3D model as groups in progressive slot tabs according to the number of classification criteria, such as lithology in Table 1.

Table 1. Typical SF and Pf values as acceptance criteria.

Scale of analysis	Consequences of the rupture	Acceptance criteria
		FS min. Static	FS min. Dynamic	Max. P [FS < 1]
Bank	Low-high	1.1	NA	25% - 50%
Interrupt	Low	1.15 - 1.2	1.00	25%
	Mean	1.2	1.00	20%
	High	1.2 - 1.3	1.10	10%
Overall	Low	1.2 - 1.3	1.00	15% - 20%
	Mean	1.3	1.05	10%
	High	1.3 - 1.5	1.10	5%

Source: Read and Stacey (2009).

Finally, the FLAC3D mesh file is synthetically created using the create_grid() function, which creates the mesh considering the global indices of the zones and vertices, as well as the clusters created.

3.2. Hydrogeological Properties

The rock classes evaluated were grouped considering that the friable rocks correspond to class V rocks and that the compact rocks correspond to class I to IV lithotypes. The values used for compact rocks correspond to the averages of the classes considered. All lithologies present in the model were considered, and the values are shown in Figure 8. The values presented correspond to the values of the scalar component of the conductivity tensor corresponding to each family of discontinuities and to the conductivity of the rock matrix.

The routines for calculating the properties required for the model were programmed in the FLAC3D Python console. Features such as the orientations of the

Figure 8. Hydraulic conductivity values by family of discontinuities.

families are incorporated into the model by calculating the anisotropic conductivity tensors. Figure 9 shows the complete process of creating and solving the steady-state flow problem in the context of pit slopes.

Figure 9. Flowchart of the Python routine for creating FLAC3D meshes.

The conductivity tensor was entered into the model through each of the 6 unique and symmetric components (kxx, kyy, kzz, kxy, kxz and kyz) of the FLAC3D anisotropic conductivity model. Figure 10 shows the logK conductivity fields for the mine in its component aligned with the Cartesian X-axis (kxx). Some lithologies can function as natural barriers due to their low conductivity (blue) compared to the other lithologies. On the other hand, it can also be observed that there are lithologies that present high values of conductivity (red), as is the case for friable materials.

The piezometric surface was created via file import and subsequent Delaunay triangulation. Finally, after fixing the load boundary conditions of the model, the model is run until the steady-state flow condition is reached. Figure 11 and Figure 12 show the pore pressure fields initialized and after the calculation of the stationary condition for the pit, respectively.

Figure 10. LogK conductivity field of the pit—kxx component of the conductivity tensor.

Figure 11. Initial pore pressure field of the final pit.

Figure 12. Final pit pore pressure field.

3.3. Geomechanical Properties

The geomechanical quality of the rock mass under study was evaluated in terms of the geological strength index (GSI). The GSI was evaluated in 821 sections corresponding to 83 boreholes (georeferenced). Figure 13 shows the location of the 821 GSI samples in the mine context. These data were used to perform a GSI interpolation to cover the entire domain using the inverse square distance method. This interpolation was performed on the regular mesh of the block model prototype, which has dimensions of 50 × 50 × 10 m. Each subblock within the original block (even IJK) inherited the GSI of the block. The results of this process are shown in Figure 14.

Figure 13. Location of the GSI map.

Figure 14. Spatial distribution of the GSI.

The rock classes evaluated were grouped by considering that the friable rocks correspond to class V rocks and that the compact rocks correspond to class I to IV lithotypes. The values used for compact rocks correspond to the averages of the classes considered.

Three constitutive models were used to represent the rock mass failure criterion, namely, the Mohr-Coulomb (MC), generalized Hoek-Brown (HB) and ubiquitous joint (UJ) models. However, all lithologies present in the model were considered, and the values are shown in Figure 15. The values presented correspond to the geomechanical parameters necessary to directly save or, if necessary, infer the parameters of the specific model. When the criterion used was the BM, the parameters in Figure 15 were used directly. In the case of the HB and UJ criteria, the model parameters were calculated based on the lithotype and the GSI geomechanical quality corresponding to each model block. The joint parameters in the UJ model were used directly.

Figure 15. Values of the geomechanical properties of intact rocks for the HB and UJ, MC massif and UJ joints.

The routines for calculating the properties required for the model were programmed in the FLAC3D Python console. Figure 16 shows the process of defining the geomechanical properties. The first 14 extra variables of the zones were allocated with the numerical magnitudes of the model as follows: lithology with the order number in the input list starting from zero (Figure 3), density (Figure 4), failure indicator (Figure 5), GSI (Figure 14), constitutive model with MC:0 and HB:1. UJ:2 (Figure 17), compressive strength of intact rock (Figure 18), deformation modulus of intact rock (Figure 19), Hoek‒Brown constant m_i (Figure 20), foliation friction angle (Figure 21), foliation cohesion (Figure 22), foliation plunge (Figure 23), foliation plunge direction (Figure 24), massif friction angle (Figure 25), and massif cohesion (Figure 26).

Figure 16. Flowchart of the Python routine for creating basic geomechanical properties.

Figure 17. Distribution of the constitutive model.

Figure 18. UCS distribution of intact rock.

Figure 19. Distribution of the deformation modulus of intact rock.

Figure 20. Distribution of the constant mi of the intact rock.

Figure 21. Distribution of foliation friction angle in the pit.

Figure 22. Distribution of foliation cohesion in the pit.

Figure 23. Distribution of foliation dip in the pit.

Figure 24. Distribution of the dip direction of foliation in the pit.

Figure 25. Distribution of the friction angle of the pit.

Figure 26. Cohesion distribution of the pit.

The following constitutive model parameters were obtained: the Hoek‒Brown constant m_b, the Hoek‒Brown constant s, the Hoek‒Brown constant a, the deformation modulus of the rock mass, the equivalent friction angle HB for the MC to use in the UJ, and the cohesion equivalent HB to the MC to use in the UJ. The generalized Hoek‒Brown criterion parameters were estimated using the formulation presented by Hoek, Carranza-Torres & Corkum (2002), where the envelope is described by Equation (1).

${\text{<msup> σ <mo>′</mo> </msup>}}_1={\text{<msup> σ <mo>′</mo> </msup>}}_3+{\text{σ}}_{ci}{\left(m_b\frac{{\text{<msup> σ <mo>′</mo> </msup>}}_3}{{\text{σ}}_{ci}}+s\right)}^a$ (1)

where σ_ci is the uniaxial compressive strength of the intact rock and

$m_b=m_i\exp \left(\frac{GSI-100}{28-14D}\right)$ (2)

$s=\exp \left(\frac{GSI-100}{9-3D}\right)$ (3)

$a=\frac{1}{2}+\frac{1}{6}\left({\text{e}}^{-GSI/15}-{\text{e}}^{-20/3}\right)$ (4)

where m_i is a constant that describes the resistance behavior of intact rock, GSI is the geological quality index and D is a constant that describes the blasting quality (D = 0 for high-quality blasting and D = 1 for low-quality blasting). For the cases of the massif cut by the fault, the massif was considered to be in its residual state. Thus, the residual value of the GSI was used according to Cai et al. (2007).

$GS{I}_r=GSI\exp \left(-0.0134GSI\right)$ (5)

In the case of the ubiquitous joints model, the equivalence of the continuous portion of the massif was determined through an equivalence between the generalized Hoek‒Brown criterion and the Mohr-Coulomb criterion, as proposed by Hoek, Carranza-Torres and Corkum (2002).

${\varphi }_{rm}={\sin }^{-1}\left[\frac{6a{m}_b{\left(s+m_b{\text{<msup> σ <mo>′</mo> </msup>}}_{3n}\right)}^{a-1}}{2\left(1+a\right)\left(2+a\right)+6a{m}_b{\left(s+m_b{\text{<msup> σ <mo>′</mo> </msup>}}_{3n}\right)}^{a-1}}\right]$ (6)

$c_{rm}=\frac{{\sigma }_{ci}\left[\left(1+2a\right)s+\left(1-a\right)m_b{\text{<msup> σ <mo>′</mo> </msup>}}_{3n}\right]{\left(s+m_b{\text{<msup> σ <mo>′</mo> </msup>}}_{3n}\right)}^{a-1}}{\left(1+a\right)\left(2+a\right)\sqrt{1+\left(6a{m}_b{\left(s+m_b{\text{<msup> σ <mo>′</mo> </msup>}}_{3n}\right)}^{a-1}\right)/\left(\left(1+a\right)\left(2+a\right)\right)}}$ (7)

where ${\text{<msup> σ <mo>′</mo> </msup>}}_{3n}={\text{<msup> σ <mo>′</mo> </msup>}}_{3\max }/{\text{σ}}_{ci}$ .

where the maximum confining stress ${\text{<msup> σ <mo>′</mo> </msup>}}_{3\max }$ , where equivalence is needed, will depend on the specific problem. Hoek, Carranza-Torres and Corkum (2002) also presented specific general and specific recommendations for cases involving tunnels and slopes. In the case of the numerical models reported here, the general equivalence used was due to the large spectrum of equivalence cases for particular cases. The equivalence for general cases is ${\text{<msup> σ <mo>′</mo> </msup>}}_{3\max }$ = σ_ci/4. The deformation modulus of the rock mass was estimated using the method of Hoek and Diederichs (2006).

$E_{rm}=E_i\left(0.02+\frac{1-D/2}{{\text{e}}^{\left(\left(60+15D-GSI\right)/11\right)}}\right)$ (8)

where E_i is the deformability module of intact rock, which in the absence of experimental data can be estimated, depending on the lithotype, as:

$E_i=MR\cdot {\text{σ}}_{ci}$ (9)

where MR is the modulus reduction factor.

Finally, the parameters of the Mohr-Coulomb strength criterion of the foliations used in the ubiquitous joint model were estimated using the equivalence proposed by Prassetyo, Gutierrez and Barton (2017).

${\varphi }_j=\arctan \left(A-B\right)$ (10)

$c_j={\sigma }_{n}B$ (11)

where

$A={\tan }^2\left[JRC{\log }_{10}\left(\frac{JCS}{{\text{σ}}_n}\right)+{\varphi }_r\right]$ (12)

$B=\left(1+A^2\right)\frac{JRC}{\ln 10}$ (13)

where B is in degrees, i.e., it must be converted to radians before use.

Figure 27 presents a flowchart showing the sequence of estimation of the parameters of the rock mass strength and deformability models to be used. Figure 25, Figure 26, Figures 28-31 show the parameters obtained in the 3D space of the pit.

4. Results and Discussion

4.1. Stress Deformation

Three-dimensional stress-strain analyses were performed after the pore pressure fields were obtained via piezometric surface generation. Next, the results for displacements, stresses and the resistance-to-stress ratio (SSR) are shown to present the mechanical behavior of the pit. Two stages were evaluated: initial pit and final pit.

4.1.1. Initial Pit

In the initial pit, the displacements were less than 200 mm, with isolated occurrence of values above this threshold. Figure 32 presents the displacement field in the initial pit and details the location of geotechnical section 1-1’ to show the behavior of the west wall of the pit. This section presents an overall slope 170 m high in total, with a break slope of 80 m and a 40˚ inclination. The cut through section 1-1’ showed a maximum local displacement of 100 mm in the region near the bottom of the pit (Figure 33). However, the displacements generally remained below 80 mm.

The SSR field (Figure 34) shows that there is a considerable concentration of stresses on the inter-ramp slope. However, this stress concentration does not generate a kinematic mechanism properly, but the SSR values certainly show that safety levels can be found under marginal conditions.

4.1.2. Final Pit

In the case of the final pit, in addition to the 1-1’ section, the 2-2’ section is presented (Figure 35) to show the behavior of the northeast wall of the pit. Section 1-1’ exhibited a maximum displacement of approximately 180 mm in the lower section, corresponding to the first break of the western wall of the pit (Figure 36). The slope displacement pattern shows solid behavior, without particular concentrations that may indicate anomalies.

Section 2-2’ (Figure 37) shows the northeast slope, which has a height of 250 m and a global slope of 40˚. In general, the displacements are less than 100 mm in magnitude and less than 20 mm in the horizontal direction, with the exception of the middle region of the slope height, where horizontal movements greater than 200 mm and 350 mm in magnitude can be observed.

Figure 27. Flowchart of the Python routine for creating basic geomechanical properties.

Figure 28. Distribution of the constant mb of the pit.

Figure 29. Distribution of the constant s of the pit.

Figure 30. Distribution of the constant a of the pit.

Figure 31. Distribution of the deformation modulus of the pit.

Figure 32. Surface displacements in the initial pit—section 1-1’.

Figure 33. Displacement field in the initial pit—section 1-1’.

Figure 34. SSR in the initial pit—section 1-1’.

Figure 35. Surface displacements in the final pit.

Figure 36. Surface displacements in the final pit—section 1-1’.

Figure 37. Surface displacements in the final pit—section 2-2’.

The displacements on the surface showed a concentration of displacements in the northeast wall of the pit (Figure 37) and a lifting of the bottom on the order of 100 mm. It can also be observed that there is a concentration of displacements in the region of the fault outcrop (Figure 38).

Figure 39 shows that the SSR distribution on the west slope was highly concentrated only near the toe of the global slope. This is consistent with the mobilization of tensions for the development of a global mechanism. However, the lack of continuity in the stress concentration region suggests that the global slope should have a high safety factor. Occasional exceptions may occur for surface ruptures caused by exposure of the compact Batatal shale located below the current iron formation.

The SSR ratio of the northeast slope indicates the presence of two regions of considerable concentration at the top of the slope and in the middle region of the slope height (Figure 40). The concentration region near the top is located in friable Moeda quartzite, and the intermediate region is found in friable Batatal shale with an extension of approximately 40 m. This same intermediate region is found in the fault outcrop, as shown in Figure 38.

4.2. Safety Assessment

In the case of safety evaluations using three-dimensional numerical models, there are still no direct comparison values in the literature. It is worth remembering that typical recommendations are based on two-dimensional analyses, which provide values of factors of safety that are intrinsically different from those obtained through three-dimensional analyses (Albataineh, 2006; Chen et al., 2003; Lu et al., 2013). In addition, high failure probabilities can occur in locations where deterministic safety factors considered acceptable for open pits are obtained (Torres et al., 2022). However, the values presented in Table 1 were taken as a reference only for comparison purposes. More accurate three-dimensional safety analyses should be based on the target failure probabilities presented in the same table; however, this type of analysis is not part of the scope of this study. The safety factor was estimated only for the final pit.

The final pit had an overall safety factor of 1.78, occurring mainly on the

Figure 38. Details of the displacement vectors and fault locations in the final pit—section 2-2’.

Figure 39. SSR in the final pit—section 1-1’.

Figure 40. SSR in the final pit—2-2’ section.

northeast slope. This is due to the presence of the fault, which cuts the pit and conditions the presence of low-strength materials. Based on this typical behavior at the study site, it can be stated that the safety factor of the west slope is at least greater than those found in these global analyses. The determination of the factor of safety of the west slope should be performed in isolation from the northeast and southeast slopes. However, it should be noted that these results correspond to deep mechanisms and do not exclude the occurrence of small surface events such as thin planar landslides or buckling or falling blocks.

5. Conclusion

Geomechanical models are typically constructed in a three-dimensional manner using information obtained from soundings and interpreted based on an understanding of the geology of the studied site. These models can be sophisticated because they incorporate not only geological-geotechnical data but also mining data such as ore grades. However, simplified models are used during the construction of predictive models of mechanical and safety behavior.

In this study, a framework was presented to take advantage of the availability of sophisticated information to reconstruct, generate information and analyze the three-dimensional behavior of the geomechanical stability of slopes in mining pits. The routines were applied to the study of the geomechanical behavior of an iron pit.

Both the methods and the results of the analyses presented here suggest that the method and set of computational routines were shown to be tools with great potential and applicability for robustly and efficiently addressing highly complex geotechnical problems.

Acknowledgments

The author thanks the Vale Institute of Technology and Vale SA for the material and resources used in the preparation of this study.

Date Availability

Not available.

Funding

The author received no financial support for the research authorship and/or publication of this article

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Albataineh, N. (2006). Slope Stability Analysis Using 2D and 3D Methods. MS Thesis, University of Akron.
[2]	Cai, M., Kaiser, P. K., Tasaka, Y., & Minami, M. (2007). Determination of Residual Strength Parameters of Jointed Rock Masses Using the GSI System. International Journal of Rock Mechanics and Mining Sciences, 44, 247-265. https://doi.org/10.1016/j.ijrmms.2006.07.005
[3]	Chapra, S., & Canale, R. (2016). Métodos Numéricos para Engenharia. McGraw Hill Brasil.
[4]	Chen, Z., Mi, H., Zhang, F., & Wang, X. (2003). A Simplified Method for 3D Slope Stability Analysis. Canadian Geotechnical Journal, 40, 675-683. https://doi.org/10.1139/t03-002
[5]	Franco, N. B. (2006). Cálculo Numérico. Person Prentice Hall.
[6]	He, L., An, X. M., Ma, G. W., & Zhao, Z. Y. (2013). Development of Three-Dimensional Numerical Manifold Method for Jointed Rock Slope Stability Analysis. International Journal of Rock Mechanics and Mining Sciences, 64, 22-35. https://doi.org/10.1016/j.ijrmms.2013.08.015
[7]	Hoek, E., & Diederichs, M. S. (2006). Empirical Estimation of Rock Mass Modulus. International Journal of Rock Mechanics and Mining Sciences, 43, 203-215. https://doi.org/10.1016/j.ijrmms.2005.06.005
[8]	Hoek, E., Carranza-Torres, C., & Corkum, B. (2002). Hoek-Brown Failure Criterion. In Proceedings of the 5th North American Rock Mechanics Symposium (pp. 267-273).
[9]	Hustrulid, W. A., Kuchta, M., & Martin, R. K. (2013). Open Pit Mine Planning and Design, Two Volume Set & CD-ROM Pack (3rd ed.). CRC Press. https://doi.org/10.1201/b15068
[10]	Itasca Consulting Group Inc FLAC3D (2012). Fast Lagrangian Analysis of Continua in 3 Dimensions). v. 5.01. Itasca Consulting Group Inc.
[11]	Lu, H. H., Xu, L. M., Fredlund, M. D., & Fredlund, D. G. (2013). Comparison of 3D Finite Element Slope Stability with 3D Limit Equilibrium Analysis. In Proceedings of the 18th International Conference on Soil Mechanics and Geotechnical Engineering (pp. 759-761). https://www.issmge.org/uploads/publications/1/2/759-761.pdf
[12]	Michalowski, R. L., & Dawson, E. M. (2002). Three-Dimensional Analysis of Limit Loads on Mohr-Coulomb Soil. Fundations of Civil and Environmental Engineering, 1, 137-147.
[13]	Ning, Y. J., An, X. M., & Ma, G. W. (2011). Footwall Slope Stability Analysis with the Numerical Manifold Method. International Journal of Rock Mechanics and Mining Sciences, 48, 964-975. https://doi.org/10.1016/j.ijrmms.2011.06.011
[14]	Prassetyo, S. H., Gutierrez, M., & Barton, N. (2017). Nonlinear Shear Behavior of Rock Joints Using a Linearized Implementation of the Barton-Bandis Model. Journal of Rock Mechanics and Geotechnical Engineering, 9, 671-682. https://doi.org/10.1016/j.jrmge.2017.01.006
[15]	Read, J., & Stacey, P. (2009). Guidelines for Open Pit Slope Design. CSIRO Publishing/CRC Press.
[16]	Singh, R. N., & Ghose, A. K. (2006). Engineered Rock Structures in Mining and Civil Construction. Taylor & Francis/Balkema.
[17]	Sjöberg, J. (1999). Analysis of Large Scale Rock Slopes. Doctoral Thesis, Luleå University of Technology.
[18]	Steffen, O. K. H., Contreras, L. F., Terbrugge, P. J., & Venter, J. (2008). A Risk Evaluation Approach for Pit Slope Design. In ARMA US Rock Mechanics/Geomechanics Symposium (pp. ARMA-08). ARMA International.
[19]	Torres, V. F. N., Dockendorff, R., Sotomayor, J. M. G., Castro, C., & Silva, A. F. (2022). Probability of Failure and Factor of Safety in the Design of Interramp Angles in a Large Open Iron Ore Mine. Journal of the Southern African Institute of Mining and Metallurgy, 122, 1-14. https://doi.org/10.17159/2411-9717/2025/2022
[20]	Wyllie, D. C., & Mah, C. (2004). Rock Slope Engineering (4th ed.). CRC Press. https://doi.org/10.1201/9781315274980
[21]	Zevgolis, I. E., Deliveris, A. V., & Koukouzas, N. C. (2018). Probabilistic Design Optimization and Simplified Geotechnical Risk Analysis for Large Open Pit Excavations. Computers and Geotechnics, 103, 153-164. https://doi.org/10.1016/j.compgeo.2018.07.024