Slope stability assessment of an open pit using lattice-spring-based synthetic rock mass (LS-SRM) modeling approach

Subash Bastola, Ming Cai,c,*, Branko Damjanac

a MIRARCO, Laurentian University, Sudbury, Canada

b Bharti School of Engineering, Laurentian University, Sudbury, Canada

c Key Laboratory of Ministry of Education for Safe Mining of Deep Metal Mines, Northeastern University, Shenyang, China

d Itasca Consulting Group, Inc., Minneapolis, MN, USA

Keywords:Lattice-spring-based synthetic rock mass(LS-SRM) modeling Non-planar discontinuities Slope stability Slope model Discrete fracture network (DFN) modeling

ABSTRACT Discontinuity waviness is one of the most important properties that influence shear strength of jointed rock masses, and it should be incorporated into numerical models for slope stability assessment.However, in most existing numerical modeling tools, discontinuities are often simplified into planar surfaces. Discrete fracture network modeling tools such as MoFrac allow the simulation of non-planar discontinuities which can be incorporated into lattice-spring-based geomechanical software such as Slope Model for slope stability assessment. In this study, the slope failure of the south wall at Cadia Hill open pit mine is simulated using the lattice-spring-based synthetic rock mass (LS-SRM) modeling approach. First, the slope model is calibrated using field displacement monitoring data, and then the influence of different discontinuity configurations on the stability of the slope is investigated. The modeling results show that the slope with non-planar discontinuities is comparatively more stable than the ones with planar discontinuities. In addition, the slope becomes increasingly unstable with the increases of discontinuity intensity and size. At greater pit depth with higher in situ stress, both the slope models with planar and non-planar discontinuities experience localized failures due to very high stress concentrations,and the slope model with planar discontinuities is more deformable and less stable than that with non-planar discontinuities.

1. Introduction

Slope stability is an important geotechnical issue in surface mining operation, because slope failure can result in production halt, losses of mining machineries, and human casualties. Reliable slope design in civil and mining projects can not only improve safety but also avoid unexpected construction cost overrun due to slope failure. In assessment of slope stability, the most important factors that need to be considered are the geometry of the slope and the mechanical properties of the jointed rock mass(intact rock and discontinuities).Structural discontinuities control the stability of a slope,because these are the weakest component of the rock mass.Logically, the slope becomes increasingly unstable with the increase of unfavorably oriented discontinuities.Hence,it is vital to have a reliable estimate of the mechanical properties of intact rock and discontinuities to increase the reliability of rock slope engineering design. Mechanical properties of rock and discontinuities can be estimated experimentally and empirically, and mechanical properties of jointed rock masses can be estimated using analytical,empirical,and numerical methods(Cai et al.,2004,2007;Cai,2010;Hoek and Brown,2019).These estimated mechanical properties can be used in slope design.

Slope design can be conducted using methods such as limit equilibrium method,empirical design method,numerical modeling method, physical model testing, and probabilistic-based design method (Sjoberg,1996). Out of these design methods, the numerical modeling approach is the most popular in recent years.Among the numerical modeling methods, the discontinuum method is better suited than the continuum method to simulate slope stability problems, because the continuum method cannot simulate crack development,large deformation,and rotation of pre-existing discontinuities. Discrete element methods (DEMs) such as the bonded particle method (BPM) and the bonded block method(BBM) are widely used.

Synthetic rock mass (SRM) modeling based on BPM and BBM have been used for stability assessment,because they can explicitly represent pre-existing discontinuities and crack evolution from these discontinuities. Scholtès and Donzé (2012) identified the critical failure surface on a jointed rock slope using the BPM-based SRM approach in Yade (Kozicki and Donzé, 2009). Camones et al.(2013) investigated the step-path failure mechanism in a jointed rock slope using the BPM-based SRM approach in PFC (Potyondy and Cundall, 2004). Based on the BBM-based SRM modeling analysis in 3DEC,Brideau et al.(2011)and Havaej et al.(2016)concluded that slope stability can be influenced by slope topography,discontinuity friction angle, and discontinuity persistence.Kulatilake and Shu (2015) simulated the slope deformation of an open pit mine using the BBM-based SRM approach and indicated that the lateral to the vertical stress ratio influences the model displacements significantly. Using an BBM-based SRM model in 3DEC, Sainsbury et al. (2016) demonstrated that intermediate discontinuity structures have a large influence on the stability of mine-scale slopes.

In recent years, the lattice-spring-based SRM (LS-SRM)modeling approach has been used to study brittle rock fracturing and stability of rock slopes (Herrero, 2015; Cundall et al., 2016;Havaej and Stead, 2016), using the SRMTools-Slope Model (Itasca,2016). Herrero (2015) found that both the stability and the failure mode of slopes are influenced by the in situ stress in deep open pit mines. Havaej and Stead (2016) quantified the rock damage intensity due to brittle rock fracturing and highlighted the close relationships among kinematics, failure surface, and rock damage.Based on insights from the numerical modeling, Cundall et al.(2016) argued that slope failure is the outcome of the brittle fracturing of the intact rock bridges and that the assumption of equivalent shear strength of the slope is not justified.

In all the previous studies, discontinuities are simplified into planar surfaces for convenience, which is clearly not the case in natural geological settings.There is a large influence of non-planar crack geometry on the cracking mechanism and mechanical properties (strength and deformability) of laboratory-scale precracked rocks (Bastola and Cai, 2019). Hence, it is hypothesized in this study that the large-scale discontinuity waviness can have a large influence on slope stability and deformability.

To validate the hypothesis,the south wall slope failure at Cadia Hill open pit mine in Australia is simulated using the LS-SRM modeling approach in the Slope Model. LS-SRM models are built using Slope Model as it allows explicit representation of preexisting discontinuities (both planar and non-planar). In the LSSRM, intact rock is represented as a random assembly of nodes that are connected to each other with massless springs in threedimensional space. The mechanical formulation and contact logic of the lattice spring model in Slope Model are similar to the BPM in PFC (Potyondy and Cundall, 2004) except that the particles are represented as concentrated point masses at the lattice nodes and contacts are represented by springs in both the normal and shear directions. Large-scale discontinuities are represented as an assembly of smooth joints which have the same smooth joint logic in PFC. The lattice spring modeling approach is efficient in terms of computational time,because it is formulated in small strain which does not require contact detection and updating as in other DEM tools. Unlike other DEM modeling approaches,calibration of some micro-property parameters is not required in the lattice code,because calibration factors are built-in and most of the desired macroscopic properties can be provided as direct model inputs(Damjanac et al.,2016).Furthermore,the lattice code is efficient in handling discontinuities by avoiding the formation of very small blocks due to ill-conditioned meshing.Details on the lattice spring method and SRMTools can be found in some publications(Cundall,2011; Damjanac et al., 2016; Bastola and Cai, 2019).

In the following discussion, an overview of the Cadia Hill open pit mine site and previous studies conducted on the open pit slope stability assessment is summarized in Section 2. Discrete fracture network (DFN) and the slope model for the slope are developed,and the mechanical properties of intact and jointed rocks are calibrated in Section 3. Using the calibrated slope model parameters,the influences of different discontinuity configurations and in situ stresses on the stability of the slope are investigated in Section 4.The conclusions of the modeling results are presented in Section 5.

2. A brief overview of Cadia Hill open pit

Cadia Hill open pit gold mine is operated by Cadia Valley Operations,a subsidiary of Newcrest Mining Limited(Sainsbury et al.,2007).The mine is located within a complex geological setting near Orange City in New South Wales,Australia,approximately 250 km west of Sydney. Mining of the pit began in 1997. The pit width is about 1500 m and the final pit depth as per the design in 2005 is 580-720 m. The designed inter-ramp slope angle is 35°-65°,depending on the strength of rock mass (Fig.1).

Fig.1. Cross-section of Cadia Hill open pit south slope (Franz, 2009).

Fig. 2. (a, b) Failure location and (c) conceptual model of the failure mechanism on the south wall of Cadia Hill open pit mine (Sainsbury et al., 2007).

In 2006, the mine experienced a multi-bench failure in the center of the south wall from 535 to 625 m ramp levels (RLs)(Domain 18) while mining at 505 m RL (Fig. 2a and b). The rock mass at the failure location was composed of monzonite. The failure mechanism involved was a combination of both geological structure and rock mass failures. A large shear structure subparallel to the pit face at an orientation of 56°/004°formed a basal sliding failure plane (Fig. 2b). The failed structure does not form a daylighting wedge; hence, the failure cannot be backanalyzed using traditional limit equilibrium or simple continuum models (Sainsbury et al., 2007). The conceptual model for the slope failure mechanism is displayed in Fig. 2c, showing a combined failure of sliding along the shear structure and tensile and shear failures of the rock mass (Sainsbury et al., 2007;Franz, 2009).

The prevalent lithology of the south wall is monzonite volcanics and Silurian sediments associated with faulting. The rock mass at the south wall is characterized by moderate-to-high rock quality designation (RQD) values and high intact rock strength (120-160 MPa)(Sainsbury et al.,2007).Based on the drillhole and bench face mapping data, geological strength index (GSI) values for the rock mass at Domain 18 and the BE-Fault zone are 60-70 and 40-50,respectively.The structural geology of the south wall is complex because of four deformation episodes that have resulted in curvilinear discontinuity structures.

Sainsbury et al. (2007) simulated the failure of the south wall using 3DEC. Although the model predicted slope displacements agreed,to some extent,with the field measured data,the model did not explicitly represent the discontinuities and the mechanical properties of the rock were downgraded using the GSI approach.Franz(2009)conducted a parametric study with various conceptual slope models of Cadia Hill open pit using both continuum and discontinuum modeling approaches and indicated that the slope stability is influenced more by discontinuity continuity (persistence) than by discontinuity spacing (intensity). Neither of these studies modeled discontinuities in consideration of non-planar features.

3. Model development and calibration

In this study, the Slope Model software is used to simulate the main failure of the south wall at Cadia Hill open pit mine. Discontinuities are represented as a DFN,which is generated using MoFrac that can represent both planar and non-planar discontinuities.

Fig.3. Calibration approach for upscaling the mechanical properties of laboratory-scale specimens to the ones that are equivalent for massive rock masses.E-Young’s modulus of intact rock,σc-compressive strength of intact rock,c-cohesion of intact rock,φ-angle of internal friction of intact rock,knj-joint normal stiffness,ksj-joint shear stiffness,σm- compressive strength of rock mass, and Em - rock mass deformation modulus.

First, the intact rock model is created, and the lattice microparameters are adjusted to match the reduced macro-mechanical properties for the intact rock (Section 3.1). Then the DFN model(Section 3.2)that is representative of the discontinuity structures in Domain 18 is generated in MoFrac and incorporated into the intact rock model to generate the LS-SRM models (Section 3.3). The calibration methodology previously established for the laboratoryscale jointed rock models (Bastola and Cai, 2019) is adjusted for calibration of the mechanical properties of large-scale rock mass models (Fig. 3). The target mechanical properties of rock mass for calibration of the LS-SRM model are derived from Sainsbury et al.(2007). The slope is excavated in the calibrated LS-SRM model using Slope Model and the result is validated using the field displacement monitoring data (Section 3.4).

Fig.4. (a)Axial stress-axial strain relation for intact monzonite(solid line)and BE-Fault material(dotted line)under unconfined compression;and(b)Strength envelopes for intact monzonite (solid line) and BE-Fault material (dotted line) under triaxial compression.

Table 1Input lattice parameters used for intact rock property calibration.

3.1. Intact rock model

The intact rock model is generated with dimensions of 5 m×5 m×10 m(width×length×height).The intact rock model is calibrated to two sets of lattice parameters, one for the intact monzonite, and the other for the slightly weaker rock within the BE-Fault zone. The size of the intact rock model that is to be calibrated is significantly larger than the size of conventional laboratory specimens. Cai and Kaiser (2014) suggested that the value of the peak in situ strength of the massive intact rock should be about 80% of uniaxial compressive strength (σc) that is measured in the laboratory-scale specimens. Thus, the calibration targets of σcfor both the intact monzonite and the BE-Fault material are reduced to 80% of their respective laboratory values of σc. Voronoi-shaped lattice is used, because the lattice resolution of Voronoi-shaped lattice does not exhibit scale effect on σcand E for massive intact and jointed rocks (see Section 3.3). For initial model calibration, a lattice size of 18.75 cm is used,which results in about 26 lattices in the shortest dimension(5 m)of the synthetic specimen.A loading rate of 0.02 m/s is used along with the flat-joint contact model in the simulations.

Calibration of the massive intact rocks is performed by carrying out unconfined and confined compression tests (σ3≤6 MPa) and comparing the simulation results with the scale-adjusted laboratory test results. The confinement limit of 6 MPa is sufficient,because the in situ stress in the pit is roughly in the range of 0-6 MPa σcvalues of the laboratory-scale intact Monzonite and the weaker BE-Fault material are 120-160 MPa and 80-120 MPa,respectively, as reported by Sainsbury et al. (2007) and Franz(2009). The uniaxial compressive strengths of the intact monzonite and the weaker intact BE-Fault material are calibrated to σcof 116.5 MPa and 85.3 MPa,respectively(Fig.4a).Elastic moduli(E)of the intact monzonite and the weaker BE-Fault material are calibrated to 63.1 GPa and 43.6 GPa, respectively. Similarly, based on the triaxial compression test data of the intact monzonite and the weaker BE-Fault material,c and φ values are calibrated to 19.4 MPa and 51.9°and 16.6 MPa and 46.5°,respectively for the two materials(Fig.4b).The calibrated lattice model parameters for both rocks are shown in Table 1.

3.2. DFN model

The dimensions of the DFN model generated in this study are 250 m×200 m×150 m.DFNs are generated using MoFrac,because it can generate realistic non-planar discontinuities (Junkin et al.,2017). To ensure that the spatial locations of all the discontinuities in the DFN with planar and non-planar discontinuities are the same, the DFN with non-planar discontinuities is exported to the Fracman file format (.fab) in MoFrac to generate planar discontinuities. Joint parameters such as joint orientation, intensity, size,and waviness factor are used to generate DFNs. The joint parameters for the DFN model generation are adapted from Sainsbury et al.(2007), Wilson (2003) and Franz (2009). Discontinuities are stochastically represented as discontinuity sets based on orientation distribution.The discontinuity orientation data from Wilson(2003)and Franz (2009) were plotted in DIPS (Rocscience, 2016) and six discontinuity sets were identified(Fig.5).The dominant joint set is south dipping with an average dip of 68°.All discontinuity sets are assigned a Fisher distribution value(Fisher constant)of K=700 to reduce the dispersion of orientation distribution.

Fig. 5. Stereonet for all six fracture groups.

The fracture intensity parameter, cumulative area distribution(CAD),is derived from the discontinuity length data(Wilson,2003;Franz, 2009). In order to convert fracture length into area, a discontinuity having a length ratio of 2 is assumed along the strike and the dip direction. The number of discontinuities per unit volume or fracture density(P30)is extracted for each discontinuity set by weighting the discontinuity sets for the number and size of discontinuities within the experimental volume. Non-planar discontinuities are generated using the slope deviation parameter(SD)in MoFrac which allows the random deviation of the discontinuity along the orientation plane. The DFN parameters used to generate the base case DFN model are shown in Table 2.The DFN model with non-planar discontinuities generated using the field measured discontinuity parameters of the Cadia Hill south wall is referred to as the base case DFN. The waviness of the non-planar discontinuities corresponds to the SD value of 0.5.The range of this parameter is 0-1.The mean value of 0.5 is estimated by back-analysis through trial-and-error method.The close resemblance of the input and the output CAD parameters for all six discontinuity sets in Fig. 6 validates the generated base case DFN model.The base case DFN model has a volumetric intensity (P32) value of 0.4 m-1.

In addition to the base case DFN model (Fig. 7a), DFN models with smaller discontinuity size (Fig. 7b) and higher discontinuity intensity (Fig. 7c) are also generated. The DFN with smaller discontinuities is generated by reducing the value of maximum fracture sample area (Amax) in Table 2 to 100 m2while keeping all the other parameters the same. The DFN with a higher discontinuity intensity (P32= 0.72 m-1) is generated by changing the values of fracture areas along the slope of CAD curve (A1, A2), and Amaxin Table 2 to 1 m2, 10,000 m2, and 20,000 m2respectively while keeping all the other DFN parameters the same. The generated DFNs are also exported to the.fab file format to produce DFNs with planar discontinuities in the same spatial locations.It is noted that the DFNs generated in this study are purely stochastic and each realization is different from the others.Hence,it is recommended to constrain the DFN to the known spatial locations of the field mapped discontinuity trace data, if available.

3.3. Jointed rock mass model

The dimensions of the pit slope model are 250 m ×200 m×150 m,which are the same as the DFN model.Jointed rock mass models are generated by integrating the DFN models(Section 3.2)with the calibrated intact rock model(Section 3.1).Six jointed rock mass models are generated using DFNs with different discontinuity intensities, sizes, and waviness. The stability and deformability of the slope excavated in these rock mass models are assessed using Slope Model in Sections 3.4 and 4.

Fig. 6. Comparison of input and output cumulative area distributions (CAD) for the base case DFN.

The generated base case SRM model is calibrated by performing uniaxial compression tests on the representative elementary volume (REV) samples and comparing the corresponding strength and deformation modulus to the ones back-calculated from the empirical relations (Sainsbury et al., 2007). A study on the scaleeffect of the mechanical properties of rock mass suggests that the rock mass sample size of 10 m×10 m×20 m is the REV size of the rock mass.The base case jointed rock mass model with non-planar discontinuities is calibrated to the peak compressive strength(σm)and the deformation modulus (Em) of 3.2 MPa and 8.1 GPa,respectively using a joint friction angle of φj= 30°, joint normal stiffness of knj=1 GPa/m,and joint shear stiffness of ksj=0.1 GPa/m, as these reported by other researchers (Sainsbury et al., 2007;Franz, 2009). All other joint parameters such as tensile strength,cohesion,and dilation angle are set to zero.

3.4. Jointed slope models

The lattice model parameters of the calibrated SRM model are used in the slope stability assessment in Sections 3.4 and 4. As mentioned above, the calibrated lattice parameters for the intact rock are for Voronoi-shaped lattices with a lattice size of 18.75 cm and model dimensions of 5 m×5 m×10 m.The large-scale slope model with the dimensions of 250 m×200 m ×150 m cannot be simulated using this small lattice size with the currently available PCs to the authors.To resolve this issue,the lattice size needs to be upscaled.To justify the use of the upscaled lattice size,a sensitivity study of Voronoi-shaped lattice structure on strength and deformation modulus of rock mass is conducted. A sensitivity study on the laboratory-scale intact Zhenping marble shows that there is a small influence of lattice resolution on σcand E (Bastola and Cai,2018). A similar sensitivity study is conducted herein to investigate the influence of lattice resolution of Voronoi-shaped lattice onstrength and deformation modulus of both massive intact rock and rock mass that have dimensions of 40 m×40 m×80 m,which is 64 times the size of the rock mass REV.

Table 2DFN parameters for Cadia Hill open pit.

Fig. 7. DFN models with non-planar discontinuities: (a) Base case, (b) Reduced discontinuity size, and (c) Increased discontinuity intensity.

Fig. 8. (a) Influence of lattice size (LS) on normalized σc and normalized E; and (b) Axial stress-axial strain plots for massive intact rock and rock mass with dimensions of 40 m × 40 m × 80 m and different lattice sizes using Voronoi-shaped lattice structure.

Fig. 9. (a) Influence of lattice resolution on peak shear strength; and (b) Shear stress-shear displacement plots of a non-planar discontinuity under a constant normal stress of 0.5 MPa.

Fig.10. Shear stress-shear displacement and shear displacement-normal displacement curves illustrating the developments of microcracks along the wavy section of the nonplanar discontinuity during direct shear under the normal stress of 0.5 MPa and lattice sizes of (a) 200 cm, (b) 150 cm, (c) 100 cm, and (d) 50 cm.

The results of the sensitivity study shown in Fig.8 suggest that the Voronoi-shaped lattice does not have a significant influence on the normalized strength and deformation modulus of both the massive intact rock and the large-scale rock mass (dimensions of 40 m × 40 m × 80 m). However, the post-peak deformation modulus increases with the increase of lattice size,suggesting that the brittleness of the intact rock decreases with the decrease of lattice size. The strength and deformation modulus of the largescale rock mass are normalized using the corresponding intact rock values for the lattice size range of 80-200 cm.

The base case SRM model, which includes a large-scale release structure, is used for calibration and validation of the base case slope model.The release structure forms the basal sliding surface of the slope failure(Sainsbury et al.,2007).Lattice model parameters of the calibrated jointed rock mass model(Section 3.3)are used for calibration of the base case slope model except that the lattice size is 150 cm.

Fig.11. Slope model for Cadia Hill open pit illustrating 24 excavation stages(a)with non-planar discontinuities and(b)without discontinuities.The BE-Fault zone and displacement monitoring points are shown as sub-horizontal planes and cubes, respectively.

Fig.12. Comparison of displacements recorded on the monitoring prisms(dotted line)in the field and in the slope model with non-planar discontinuities (solid line) during excavation. The time of excavation corresponding to the pit depth (RL) is also shown.

Five direct shear tests are conducted under the normal stress of 0.5 MPa using the lattice size of 250 cm, 200 cm,150 cm,100 cm,and 50 cm on rock mass models of dimensions of 150 m × 150 m × 15 m with non-planar discontinuity. This is to determine the practical lattice resolution that is sufficient to correctly represent the mechanics of sliding and cracking in the walls of wavy rough fractures but is computationally possible for the large-scale slope model. The direct shear test results (Fig. 9)show that the lattice resolution corresponding to the lattice size of 150 cm is sufficient to simulate the cracking mechanism in the wall of the rough discontinuity.With further increase of the lattice size beyond 150 cm,there is an increase of the peak shear strength even though the amount of increase is small with the increase in lattice size. Fig. 10 shows the microcrack development along the wavy section of the non-planar discontinuity.

Fig.14. Cross-section along north-south (NS) plane at the center of the slope model illustrating rock failure (cluster of black discs) induced by the propagation and coalescence of cracks (discs) along pre-existing non-planar discontinuities.

Based on the results of sensitivity analyses of compression tests(Fig. 8) and direct shear tests (Figs. 9 and 10), the Voronoi-shaped lattice structure with lattice size of 150 cm is used in slope model. Although direct stress measurements were not conducted at the Cadia Hill Open Pit site, in situ stresses were adopted from the stress measurements at Ridgeway underground mine,which is situated within the same thrust fault system approximately 3 km away from the open pit(Li et al.,2003).The bottom of the open pit is under the influence of confining stresses to the order of σZ=4.05 MPa,σX=6.90 MPa,and σY=5.06 MPa such that the ratio of the maximum horizontal to the vertical in situ stress(k=σX/σZ)is 1.7 and the top of the open pit is free of any confining pressure(Fig.11a).The bottom boundary of the slope model is fixed,the top surface is free, and roller boundary conditions are applied to the surfaces perpendicular to the X- and Y-axis. Failure occurs under the influence of gravity and the in situ stress as a result of excavation.

The slope models are simulated in 25 steps. In the first step,initial equilibrium is established under the influence of in situ stress and gravitational force without any excavation. Displacements incurred during the first step are reset to zero and the induced microcracks are healed. Then, excavations are made to simulate the mining sequence in 24 stages with a 5 m bench increment, as shown in Fig. 11. Both vertical and horizontal displacements are recorded using 18 monitoring points shown as cubes in Fig.11b. The numbers of microcracks formed during each stage are also recorded. The bands of the weaker rocks within the BE-Fault zone are shown as sub-horizontal planes with a dip of 20°and are located at the lower section of the pit(Fig.11b).The largescale release structure with a dip of 56°and a joint friction angle of 20°is also defined which forms the basal sliding surface of the slope failure (Fig.11b). The extent of the release structure is about 90 m along the vertical direction of the slope which is the approximate vertical extent of slope failure.

Fig.13. Displacement contours along (a) vertical Z direction and (b) horizontal X direction in the base case slope model. Displacement unit in m.

Fig.15. Comparison of displacements recorded in the slope models with planar discontinuities (dotted lines) and non-planar (solid lines) discontinuities.

The slope model with the base case DFN,release structure,and non-planar discontinuities is used for model validation.The model is calibrated using the displacements recorded along X, Y, and Z directions at the south wall of the open pit.Displacements recorded in the field on the monitoring prism are compared with the displacements recorded in the slope model during excavation at the same location.Excavation stages 20-24 are located in the region of weaker rock mass pertaining to the BE-Fault zone. The increase in slope movement has also been reported to be associated with rainfall and blasting events (Sainsbury et al., 2007). Therefore, the blasting and rainfall during excavation of the weak rock mass could have had a combined influence on the sudden changes in X-and Zdisplacement at the later excavation stages. It is seen from Fig.12 that the slope model displacements agree reasonably well with the field monitoring data and the displacement contours of the simulated slope model are concentrated at or near the location of the release structure (Fig.13). Compared with Z-displacement, the higher magnitude of X-displacement is attributed to the higher horizontal in situ stress in X direction which corresponds to the respective displacement contours. In addition, the displacement contours of the slope model show a very good correlation with the vertical extent of failure of 90 m as observed in the field.This leads to a calibrated slope model. The SRM model for rock slope simulation can represent slope failure more realistically, because the failure is the result of crack propagation and coalescence associated with the pre-existing discontinuities (Fig. 14) along the release structure as investigated by Sainsbury et al. (2007). Cracking is predominant on the edges of the excavations(benches)and at the location of the weaker band of rock mass within the BE-Fault zone.The lattice parameters of the calibrated slope model are used in all the other models in Sections 4.1-4.4.

4. Influence of different parameters on slope stability

Slope deformation and failure depend on many factors such as lithology, slope and discontinuity geometry, geomechanical properties of intact rock and discontinuities,hydrogeology,in situ stress and dynamic loadings (Kulatilake and Shu, 2015). In this section,using the lattice model parameters of the calibrated slope model,the influence of discontinuity waviness(Section 4.1),discontinuity intensity (Section 4.2), discontinuity size (Section 4.3), and in situ stress state(Section 4.4)on slope displacements and factor of safety are investigated.The factor of safety is computed using the strength reduction method. The release structure is not included in the models in Sections 4.2 and 4.3, because the relative size of the release structure is too small in comparison with the size of the discontinuities in DFN with a higher P32of 0.72 m-1, while the relative size of the release structure is significantly larger than the size of the discontinuities in DFN with a lower P32of 0.1 m-1.

Fig.16. Contours of factor of safety for the slope excavated in base case SRM model with (a) planar discontinuities and (b) non-planar discontinuities.

Fig.17. Cross-sections along NS plane at the center of the slope model illustrating contours of factor of safety for the slope excavated in base case SRM model with (a) planar discontinuities and (b) non-planar discontinuities.

4.1. Influence of discontinuity waviness

It has been demonstrated that non-planar crack surface has a large influence on the mechanical properties and crack evolution of the laboratory-scale cracked marble under compression (Bastola and Cai, 2019). In this section, parameters of the calibrated largescale slope model with the release structure are used in conjunction with the jointed rock mass models generated in Section 3.3 to investigate the influence of discontinuity waviness on slope stability, except that the slope model with planar discontinuities is assigned a higher friction angle of 45°(= 30°+ 15°). This 15°increment in the basic friction angle would account for the additional frictional component contributed by the large-scale discontinuity waviness and small-scale roughness(Cai et al., 2004).

Displacements and contours of factor of safety after each excavation stage are compared between the jointed slope models with planar and non-planar discontinuities.The modeling results shown in Fig.15 suggest that the magnitudes of displacements along the X,Y, and Z directions in the slope model with planar discontinuities(dotted lines) are higher than those in the slope model with nonplanar discontinuities (solid lines). In addition, the contours of factor of safety suggest that the slope model with planar discontinuities consists of larger unstable zones than the slope model with non-planar discontinuities (Fig.16). Cross-sectional analysis along the NS plane at the center of the slope suggests that the slope model with planar discontinuities constitutes unstable zones deeper into the slope than the slope with non-planar discontinuities(Fig.17).This is because the slope excavated in the SRM model with non-planar discontinuities is relatively stiffer than the one with planar discontinuities because of block interlocking along the wavy surfaces. As a result, fewer cracks are observed on the slope with non-planar discontinuities (Fig.14) than the one with planar discontinuities (Fig. 18). The results show that the displacements recorded in the slope model with the non-planar discontinuities are closer to the field measured displacements than the ones recorded in the slope model with planar discontinuities. This suggests that the adjusted discontinuity friction angle to account for the roughness component is not enough to simulate the physics of the slope with non-planar discontinuities using the slope model with planar discontinuities. Therefore, it is recommended to employ the realistic geometry of the discontinuities to solve the geomechanical problems in rock masses with wavy discontinuities.

Fig.18. Cross-section along NS plane at the center of the slope model illustrating rock failure(cluster of black discs)induced by propagation and coalescence of cracks(discs)along the pre-existing planar discontinuities.

Fig.19. Comparison of displacements recorded in the slope models with non-planar discontinuities and different P32 values.

Fig. 20. Contours of factor of safety for slopes excavated in SRM models with P32 = 0.72 m-1: (a) Planar discontinuities, and (b) Non-planar discontinuities.

Fig. 21. Comparison of displacements recorded in the slope models excavated in SRM models with non-planar discontinuities and different discontinuity sizes.

4.2. Influence of discontinuity intensity

Strength and deformation moduli of rock mass decrease with the increase of the number of discontinuities.Hence,there is a very high likelihood that a slope excavated in a rock mass with higher discontinuity intensity (P32) would exhibit larger slope deformation.To investigate the influence of P32on slope displacements and associated factor of safety, the DFN in the base case slope model with a P32value of 0.4 m-1is replaced with a DFN model that has a P32value of 0.72 m-1.

The modeling results presented in Fig. 19 show that the displacements along X, Y, and Z directions in the slope model excavated in the SRM model with P32=0.72 m-1(solid lines)are higher than the ones with P32= 0.4 m-1(dotted lines). The slopes excavated in the jointed rock mass with both planar and non-planar discontinuities and a higher P32are highly unstable. The contours of factor of safety suggest that the slope models excavated in the SRM models with P32=0.72 m-1(Fig.20)consist of larger unstable zones than the models with P32=0.4 m-1(Fig.16),for both planar and non-planar discontinuities. This is because that the slopes excavated in the SRM models with P32= 0.72 m-1are heavily jointed with large discontinuities,which decrease the strength and stiffness of the rock mass.

Fig.22. Contours of factor of safety for the slope models excavated in the base case SRM model with Amax=100 m2:(a)Planar discontinuities,and(b)Non-planar discontinuities.

Fig. 23. Comparison of displacements recorded on the slope models with (a) planar discontinuities and (b) non-planar discontinuities under the influence of different horizontal stresses.

Fig. 24. Crack development in the slope models with planar discontinuities under the influence of (a) k = 0.85, and (b) k = 2.55.

Fig. 25. Crack development in the slope models with non-planar discontinuities under the influence of (a) k = 0.85, and (b) k = 2.55.

Fig. 26. Contours of safety factor for the slope models with planar discontinuities under the influence of (a) k = 0.85, and (b) k = 2.55.

4.3. Influence of discontinuity size

The strength and deformation moduli of rock mass decrease with the increase of the discontinuity size relative to the scale of the slope. To investigate the influence of the maximum discontinuity area(Amax)on the slope displacements and the factor of safety,the DFN in the base case slope model (with Amax= 10,000 m2) is replaced with a DFN that has Amax= 100 m2.

The modeling results presented in Fig. 21 suggest that the displacements along X, Y, and Z directions in the slope excavated in the SRM model with Amax= 100 m2(solid lines) are smaller than the ones for the slope excavated in the base case SRM model with Amax= 10,000 m2(dotted lines). The slope excavated in the SRM model with Amax= 100 m2is very stable. The contours of factor of safety suggest that slope excavated in the base case SRM model is composed of larger unstable zones(Fig.16)but the slope excavated in the SRM model with Amax= 100 m2(Fig. 22) has no unstable zone. This is because that the SRM model with Amax= 100 m2consists of massive rock blocks formed by nonpersistent discontinuities which increase the strength and stiffness of the rock mass. The slopes excavated in the jointed rock mass with both planar and non-planar discontinuities for Amax= 100 m2are stable (Fig. 22).

4.4. Influence of in situ stress

With the increase of the pit depth, the local in situ stress state changes due to the change in regional geology and tectonics. High horizontal stresses tend to promote the extension and coalescence of pre-existing discontinuities, resulting in the formation of macroscopic failure planes. To investigate the influence of in situ stress on slope displacements and factor of safety for models with planar and non-planar discontinuities, the magnitude of in situ stress state in the slope models with both planar and non-planar discontinuities in Section 4.1 are altered such that slope models with maximum horizontal to vertical stress ratios (k) of 0.85 and 2.55 are simulated,in addition to the base case slope model with k of 1.7.

Fig. 27. Contours of safety factor for the slope models with non-planar discontinuities under the influence of (a) k = 0.85, and (b) k = 2.55.

Fig. 28. Cross-sections along NS plane at the center of the slope models with planar discontinuities under the influence of (a) k = 0.85, and (b) k = 2.55.

The modeling results suggest that the magnitudes of displacements along X, Y, and Z directions in the slope model excavated in the SRM models with both planar discontinuities (Fig. 23a) and non-planar discontinuities (Fig. 23b) increase with the increase of k. As witnessed in Section 4.1, the displacements on the slope models with planar discontinuities are higher than those on the slope models with non-planar discontinuities.This is because of the higher rate of extension of the pre-existing discontinuities due to the increase in induced microcracks, leading to the localized instability of the pit benches (Figs.14, 24 and 25). The rate of increase of extensile cracks with the increase in ‘k’ is larger for the slope model with planar discontinuities.

The instability is also illustrated by the contours of safety factor and cross-sections. The slope models with both planar discontinuities(Figs.16a and 26) and non-planar discontinuities(Figs.16b and 27) with higher horizontal stresses are composed of larger unstable zones.The cross-sectional plots along the NS plane at the center of the slope suggest that the slope model with higher horizontal stresses constitutes unstable zones deep into the slope(Figs. 17, 28 and 29). The observations are in line with the ones reported by Herrero (2015) in the slope model excavated in a hypothetical rock mass with planar discontinuities under the influence of different magnitudes of horizontal stresses. As expected,the extent of unstable zone is larger and deeper for the slope model with planar discontinuities.

5. Conclusions

Slope deformation and failure on the south wall of Cadia Hill open pit mine are successfully simulated using the LS-SRM modeling approach in Slope Model. Rock failure and larger displacements are localized around the release structure. Slope stability and displacements are highly dependent on the discontinuity characteristics such as waviness,intensity,size,and the magnitude of in situ stress.

The slope model excavated in the rock mass with planar discontinuities is more deformable and less stable than the one excavated in the rock mass with non-planar discontinuities. In addition, the deformability and instability increase for the slopes excavated in the rock mass with higher discontinuity intensity and larger discontinuity sizes. At greater pit depth with higher in situ stress, both slope models with planar and non-planar discontinuities experience localized failures due to very high stress concentrations, and the slope model with planar discontinuities exhibits larger deformation and a wider unstable zone than that with nonplanar discontinuities.

LS-SRM models used in this study have opened new avenues for slope stability assessment using realistic non-planar discontinuities. Moreover, these models are computationally more efficient than other DEM tools;they can also consider intersecting and nonplanar discontinuities in slope stability analysis. Due to limited computer resources, the dimensions of the analyzed slope models in this study are intermediate. The model dimensions can be extended to full pit scale using more powerful computing services such as a supercomputer and cloud computing.

Declaration of competing interest

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

Acknowledgments

The authors would like to thank Ontario Trillium Scholarship for supporting the doctorate program at Laurentian University.Financial supports from the Natural Sciences and Engineering Research Council of Canada (NSERC CRD 470490-14) of Canada,Nuclear Waste Management Organization (NWMO), and Rio Tinto for this work are gratefully acknowledged. Special thanks are extended to Dr. Varun Varun and Maurilio Torres from Itasca for providing the academic license and technical support for the Slope Model software courtesy of Itasca Education Partnership. The authors would also like to thank the editor and the two anonymous reviewers for their valuable feedbacks.