Numerical Modelling of Slope Stability in Open Pit Phosphate Mines, Algeria: A Comparative Study

The aim of our work is to investigate the comparative performance of three different software packages. Slope stability of an open pit phosphate mine was studied using numerical methods, finite element (FEM) and finite difference methods (FDM). More specifically, the widely used geotechnical software FLAC 2D of Itasca Consulting Group Inc., PLAXIS 2D of Plaxis bv and Phase2 of Rocscience Inc. are used by applying the so-called shear strength reduction (SSR) technique under Mohr-Coulomb constitutive models. The scope of the work is to compare the solutions obtained for the exact same problem from the three software packages and to discuss the comparative similarities and differences among them.


Introduction
The geotechnical engineer frequently uses limit equilibrium methods of analysis when studying slope stability problems. These methods calculate a safety factor (SF), which is defined, with respect to the shear strength of the rock (or soil), as the ratio of the available shear strength to the shear stress at equilibrium [1]. Their basic characteristic is the simplicity and the reduced number of parameters they require (slope geometry, geological parameters, geotechnical parameters and hydrogeological conditions). However, stability problems in open pit mines often involve complexities that are not easily treated by limit equilibrium methods.
In this case, using numerical methods, such as the FEM and FDM, provides the engineer with the opportunity to conduct more comprehensive and sound slope stability analyses. The main advantage of numerical methods is that no assumption needs to be made about the shape, direction and location of slip surfaces [2]. Additionally, numerical models are able to generate stress-strain distributions (deformational response), which may be of crucial importance for a robust interpretation of slope behaviour, and also to address and analyse precisely complicated geometries, simulate pit excavation stages and their effect on the stress state, address possible groundwater seepage, etc. Particularly in the context of large-scale mining slopes, obtaining realistic input information and interpreting the results are the most difficult aspects of numerical modelling [3].
The objective of our work is to investigate the comparative performance of three different software packages. Slope stability of an open pit phosphate mine was studied using numerical methods, FEM and FDM. More specifically, the widely used geotechnical software FLAC 2D of Itasca Consulting Group Inc. [4], PLAXIS 2D of Plaxis bv [5] and PHASE 2 of Rocscience Inc [6]. are used by applying the so-  [7]. Duncan (1996) [8] defined the factor of safety of soil slopes as the ratio of actual shear strength to the minimum shear strength required to avoid failure, or the factor by which shear strength must be reduced to bring a slope to failure. This method was applied for rock masses [9]. Yingren and Shangyi (2004) [10] and Hammah et al. (2007) [11] demonstrated the efficiency of SSR method for slopes of soil and rock masses.
The scope of the work is to compare the solutions obtained for the exact same problem from the three software packages and to discuss the comparative similarities and differences among them.

Methods of slope stability analysis
The numerical methods are a powerful tool for solving many engineering problems. The FEM and FDM are among of these methods. In our work, these analyses are performed through different software PLAXIS 2D of Plaxis bv, FLAC 2D of Itasca Consulting Group Inc and PHASE 2 of Rocscience Inc. In these analyses are, the factors of safety were obtained by Shear Strength Reduction (SSR) Method. Over the last 10-15 years, this technique (SSR) has become an integral part of most geotechnical numerical software packages.
This method was also adopted in several other studies [12,13,14]. The advantages of the SSR method over limit equilibrium method are: (1) No assumption is required of the interslice shear force distribution; (2) the critical failure surface in the slope can be found from the shear strain; and (3) this method is suitable for complex slope conditions in order to interpret details of displacements, stresses and water pressures.
The basic essence of strength reduction is that the rock mass strength index c and φ values are divided by a reduction factor Fs to be a group of new value c' and φ', then the new value c' and φ' are used in the trial calculation as new material parameters of finite element calculation. When the slope conforms to the given critical failure state, the corresponding stress reduction factor (SRF) is the safety factor for the slope, and the slip surface in the slope is the potential sliding surface of the slope. The parameters of c * and φ * are derived from formula (1), (2): * * (1) * * (2) Where: C * and φ * are rock reduced strength characteristic, (adhesion and friction angle) in proportion to the real mode (C, φ).

Site geology and mine description
The site (Kef-Essnoun mine) is located on the southern flank of Djebel El-Onk cretaceous anticline. It is about 73.5 km in the south of Tebessa province in northeast of Algeria and about 21 km to the Algerian-Tunisian border (figure 1) [15]. The geology of the site is mainly composed of a series of Ypresian Dolomite Limestone with silex, locally overcome by Lutetium Limestone, Miocene sands, and quaternary alluvium [16].

Methodology
In the context of this work, different characteristics are sought in the analyses but relevant to each other, such as: 1-The potential sliding surfaces of the slope of the open pit, illustrated by the location of the shear stresses, 2-The kinetics of the collapse mechanism, illustrated by the distribution of displacement vectors and 3-the "failure" extent, illustrated through plastic regions propagation.    Figure 4a, on which maximum values are observed in the foot of the slope and at the vertical tensile crack behind the crest of the slope. Concerning the slip surface developed, this seems to be a non-circular surface commonly observed in this type of problems. In Plaxis, the total deviatoric shear strain contours may be used for the representation of the critical slip surface. In this case, shear strain localization and the global slip surface of the mining slope is shown in Figure 4b. Last, the slope failure surface generated by Phase 2 , is shown in Figure 4c by means of maximum shear strain contours. Again the maximum accumulation of shear deformations seems to be developed in the foot of the slope. Overall, it can be said that in terms of shear strain concentrations and the indicated sliding surfaces, a good qualitative agreement occurs between FLAC, Plaxis and Phase 2 . Finally, it is worth mentioning that the three software, to able to generate a vertical tensile crack behind the slope crest.

Safety factors
In terms of SF, all three software packages ended up with similar values: FLAC, Plaxis and Phase 2 resulted in identical values (SF = 0.87, SF = 0.84 and SF=0.87 respectively). For the cross section under study, this difference (between FLAC / Plaxis and Phase 2 ) is negligible, since the slope appears to be in the case of extreme failures. Therefore, the difference between 0.87 and 0.84 will probably not make any difference in the slope design. In such cases, very small differences in SF (resulting from 0.01 to 0.03) do not need to be discussed.

kinetics of the collapse mechanism (Displacement vectors distribution)
Generally, the kinetics and the collapse mechanism of either potential or already-existing failures, provide useful information to the engineer. In the present work, kinetics and collapse mechanisms are discussed by means of vectors that describe the slope motion in the plane at the time of the collapse. These are the velocity vectors, the incremental displacement vectors and the total displacement vectors, as used respectively in FLAC, Plaxis and Phase 2 ( Figure. 5). It shall be noted that although the direction and the size of the displacement/velocity vectors has no quantitative meaning on an SSR type of analysis, however they provide a good qualitative sense of the developing failure mechanism since they reflect the physical pattern of the motion at the edge of collapse. As shown in Figure 5, in all three cases, the direction of displacement/velocity vectors reflect the developing non-circular failure mode, which is typical on rock slopes. As mentioned earlier, in such conditions (where ground can be realistically simulated as a continuous medium), failure depends rather on the shear strength of the material than on the shear strength of the discontinuities. So, it can be said that the representation of the slope motion in failure is practically identical (qualitatively speaking) for all three different software packages.

Failure extent (plastic region propagation)
Plastic regions propagation demonstrates the areas on which the yield criterion (Mohr-Coulomb in the present study) is satisfied. However, plastic regions do not necessarily coincide with the developing shear bands and do not always reflect the sliding surface of the slope. In fact, it is the combination of plastic failure extent with the shear strain concentrations that gives the critical developing sliding surface. Sometimes, the propagation of the plastic area in regions far away from the shear band either implies the possible presence of a secondary (and less significant) failure mechanism, or it might be simply an alarm of inadequate boundary conditions and model limits. In the light of the above, differences among the results of the three models regarding plastic regions within the slope mass and the corresponding failure extents were quite evident ( Figure. 6). More specifically: Figure 6a depicts the plasticity indicators on FLAC (circles denoting failures at tension, asterisks denoting failures at shear). The extent of the plastic points in this model is quite close to the area covered by the shear band (shear strain contours) described earlier. It is noted that the tension cut-off points observed in the foot of the slope and behind the crest of the slope indicate the tensile cracks formed slope failure phenomenon. As shown in Figure 6a, the plastic zone propagates within the marl material of the slope base.
On Plaxis (Figure. 6b, black squares denoting failure at tension, red squares denoting failures at shear), the area covered by plastic points is fairly more extended than the area covered by the shear band (see Figure. 4b), and more extended than the equivalent area covered by plastic points on FLAC ( Figure. 6a). To be more specific, besides the major failure mechanism (represented by the plastic points mostly in the marl layer, but also slightly within the phosphate at the slope's base).

Figure 5. Displacement vectors distribution
In the Figure 6c, white circles denoting failures at tension and asterisks denoting failures at shear. As we note, plastic points (asterisks) seem to extend well along the marl layer (compared to slight intrusion all these points on FLAC and Plaxis). In addition, tensile failure is concentrated up to the vertical boundary behind the crest of the slope, and slight intrusion of shear failure in the foot of the slope indicates the tensile cracks formed slope failure phenomenon. On the other hand, in the Plaxis and FLAC, the tensile failure is more concentrated on the cover mass (phosphate and limestone Ypresian Lutetian).  Figure 6. Plastic regions propagation Based on all the above, it is fair to say that overall Phase 2 seems to be the most consistent in terms of agreement between the localization of maximum shear strains and the plasticity indicators development, since a good fit between them existed in the first place. In addition, Phase 2 did not require an extension of the bottom boundaries. On the other hand, plastic points in Plaxis and FLAC models are quite extensive and exceed significantly the contours of the maximum shear strains.

Conclusions
Slope stability of an open pit phosphate mine (Kef-Essnoun) was studied using three different numerical software: FLAC (finite difference), Plaxis (finite elements) and Phase 2 (finite elements of Rocscience Inc). Analyses were performed in drained, plane strain, conditions. Elastic -perfectly plastic Mohr- Coulomb constitutive models were used for simulating the mechanical behavior of the involved strata. The shear strength reduction technique was used as the analysis scheme. Based on the results of the study, all three programs demonstrated a good agreement in the determination of safety factors and the kinetics of collapse mechanisms (identification through displacement or velocity vectors). Critical sliding surfaces, identified by shear strains localizations, were found to be similar. Last but not least, noteworthy differences appeared among the three programs in the development of plastic regions (failure in shear) and of regions failing in tension. With respect to these, Phase 2 seemed to obtain the more realistic results and to be the most consistent in terms of agreement between the localization of maximum shear strains and the plasticity indicators development. The other two programs, either required more extended limits in order to achieve consistency between maximum shear strains and plasticity indicators, or insisted on producing large (non-realistic) areas of plastic failures. Overall, it can be said that all three programs may address the examined problem satisfactorily.