Application of Limit Equilibrium Analysis and Numerical Modeling in a Case of Slope Instability

: The application of limit equilibrium analysis and numerical simulation in case of slope instability is described. The purpose of the study was to use both limit equilibrium methods (LEMs) and numerical simulations (ﬁnite element method (FEM)) to understanding the common factor imposing the selected slope into slope instabilities. Field observations, toppling analysis, rotational analysis, and numerical simulations were performed. The results of the study showed that the selected unstable slopes were associated with the sliding types of toppling; it was observed that the slopes were governed by tension cracks and layered soil mass and dominated with approximately two joints sets throughout. The simulated factor of safety (FoS) of the slopes composed of clay soil was denoted to be prone to slope instability while others were categorized as moderately stable. The simulated FoS of the slopes correlated very well with the visual observations; however, it is anticipated that properties of soil mass and other characteristics of the slopes contributed largely to the simulated FoS. The sensitivity of the model was further tested by looking into the e ﬀ ect of the slope angle on the stability of the slope. The results of the simulations showed that the steeper the slope, the more they become prone to instability. Lastly, Phase 2 numerical simulation (FEM) showed that volumetric strain, shear stress, shear strain, total displacement, and σ 1 and σ 3 components of the slope increase with the stages of the road construction. It was concluded improper road construction, steepness of the slope, slope properties (soil types), and multiple geological features cutting across are the common mechanisms behind the slope instability.


Introduction
The history of landslides in a form of slope instability has been documented to be one of the global threats to humans and the environment, as well as infrastructures [1][2][3][4][5][6][7]. Although landslides are considered to be a global challenge, Africa is one of the continents that has been identified to be leading among others [2,8]; however, it appears that there are several landslides that occurred in Africa but few have been reported and studied [5,[9][10][11][12]. It is considered a challenge in most African countries to purely study the mechanism as well as a prediction of the recurrence of the landslides. One would say failure to perform detailed studies or report landslide is not due to ignorance or lack of technological advancement but rather a lack of access to some remote areas and also an insufficient number of specialists in the field of study. Studies such as those of Hong et al., [13], Kirschbaum et al., [14], Nadim et al., [11], and Stanley and Kirschbaum [7] have contested that the global landslides susceptibility maps of Africa are governed by very little or limited data; such limited information affects the calibration and reliability of the maps [15]. Furthermore, the adequateness of the global landslide susceptibility maps of Africa are considered to be dramatically affected due to an increase in the population [13]. Additionally, this implies that there are several developments such as road construction along mountainous environments, disturbing the topography of several slopes; as a result, both micro and large landslides are expected to occur. For argument's sake, studies such as those of Gariano and Guzzetti, [2] and Maes et al. [10] have argued that there is a gap of knowledge on landslides research in Africa, as such the gap knowledge has to be bridged by furthering the studies in landslides' occurrence and prediction.
The ultimate aim of this paper is to apply the limit equilibrium and numerical modeling in assessing factors contributing to the slope instability of the R523 road in Thohoyandou, Limpopo Province, South Africa. In addition, we aim to identify safety factors, stresses, strain, and displacement that have occurred during and after road construction. Post and active slope instability along the regional road (R523) in the province were identified as the case studies.

Locality of the Study Area
The scope of the study was limited to a national road located in the Thulamela Municipality: the R523. The Thulamela Municipality is a Category B municipality situated within the Vhembe district in the far north of the Limpopo province as shown in Figure 1. The Kruger National Park forms its eastern boundary while the southern and south-western side is bordered by Makhado Municipality. As the smallest of the four municipalities in the district, Thulamela is spread over about 2642 km 2 . It is, however, the largest in the province in terms of population size [16].
Sustainability 2020, 12, x FOR PEER REVIEW 2 of 33 Furthermore, the adequateness of the global landslide susceptibility maps of Africa are considered to be dramatically affected due to an increase in the population [13]. Additionally, this implies that there are several developments such as road construction along mountainous environments, disturbing the topography of several slopes; as a result, both micro and large landslides are expected to occur. For argument's sake, studies such as those of Gariano and Guzzetti, [2] and Maes et al. [10] have argued that there is a gap of knowledge on landslides research in Africa, as such the gap knowledge has to be bridged by furthering the studies in landslides' occurrence and prediction.
The ultimate aim of this paper is to apply the limit equilibrium and numerical modeling in assessing factors contributing to the slope instability of the R523 road in Thohoyandou, Limpopo Province, South Africa. In addition, we aim to identify safety factors, stresses, strain, and displacement that have occurred during and after road construction. Post and active slope instability along the regional road (R523) in the province were identified as the case studies.

Locality of the Study Area
The scope of the study was limited to a national road located in the Thulamela Municipality: the R523. The Thulamela Municipality is a Category B municipality situated within the Vhembe district in the far north of the Limpopo province as shown in Figure 1. The Kruger National Park forms its eastern boundary while the southern and south-western side is bordered by Makhado Municipality. As the smallest of the four municipalities in the district, Thulamela is spread over about 2642 km 2 . It is, however, the largest in the province in terms of population size [16].

Geological Setting of the Study Area
The R523 road identified for this research is situated in the rugged topography of the Soutpansberg Group. This geological group is partially buried beneath sedimentary and volcanic rocks of different thicknesses. The group is dominated by large faults and joints, thereby creating a discrete and blocky rock mass especially in Thulamela [17][18][19][20]. Furthermore, conjugate faults and

Geological Setting of the Study Area
The R523 road identified for this research is situated in the rugged topography of the Soutpansberg Group. This geological group is partially buried beneath sedimentary and volcanic rocks of different thicknesses. The group is dominated by large faults and joints, thereby creating a discrete and blocky Sustainability 2020, 12, 8870 3 of 33 rock mass especially in Thulamela [17][18][19][20]. Furthermore, conjugate faults and joint sets have been evidenced in the northwest-southeast and northeast-southwest directions with quartz vein filling the brecciated joints and faults (see Figure 2). It is believed that these geological features have led to high rockfalls and slope stability problems on the roadside. An attempt is made to understand these phenomena with the help of computer-based simulation tools built around the finite element method (FEM) framework and limit equilibrium.
Sustainability 2020, 12, x FOR PEER REVIEW 3 of 33 joint sets have been evidenced in the northwest-southeast and northeast-southwest directions with quartz vein filling the brecciated joints and faults (see Figure 2). It is believed that these geological features have led to high rockfalls and slope stability problems on the roadside. An attempt is made to understand these phenomena with the help of computer-based simulation tools built around the finite element method (FEM) framework and limit equilibrium. These roads are the bloodline between villages and surrounding townships. Safety is, therefore, critical when one considers the fact that they were excavated by the blasting of very sloppy and unstable terrain [21]. These roads also pass through several rivers [22]. Intensive summer rainfall and human activities in the study area also make roads susceptible to frequent slope failure.
Several remedial measures have been suggested and implemented but rockfalls and slope instability are still reported at an alarming rate [21]. For example, Mutanamba [22] analyzed the stability of cut-slope along the R523 road between Thathe Vondo and Khalavha. First, geological structures, properties of the soil, and groundwater regime were surveyed. Then, conventional analytical methods were used to classify the type of soil and produce slope stability charts. Finally, two types of failures occurring close to the road were identified from the findings: rotational slump failure and rockfall in the form of toppling [21]. The qualitative findings, however, could not be used to comprehensively describe the underlying mechanism of the active and post landslide, occurrence, and deposition. That is what the present research is attempting to explore.

Materials and Methods
Detailed methodologies were followed in this paper. These methods include toppling analysis, rotational analysis, transitional analysis, and finite element methods (Phase 2). Thulamela Municipality road (R523) was used as the case study. Historically, limit equilibrium methods (LEM or LE analysis) have been extensively used to analyze landslides [23][24][25][26][27]. The majority of reported cases have been based on the toppling and rotational analyses of landslides.
Kinematic analysis is an alternative set of techniques that can also be used for the analysis of rock slope stability. The term generically refers to any empirical method dealing with landslides and rockfalls from the premise of their specific geological structures, physical properties of the material involved, and groundwater regime, as well as slope geometry [28].
Kinematic analysis is commonly used to identify the translational failures that occur due to the formation of wedges or planes among others [29,30]. Kinematic techniques can be regarded as useful tools for the simplified analysis of rock slope stability. However, their relevance is limited to desktop studies beyond which limit equilibrium methods are better suited. These roads are the bloodline between villages and surrounding townships. Safety is, therefore, critical when one considers the fact that they were excavated by the blasting of very sloppy and unstable terrain [21]. These roads also pass through several rivers [22]. Intensive summer rainfall and human activities in the study area also make roads susceptible to frequent slope failure.
Several remedial measures have been suggested and implemented but rockfalls and slope instability are still reported at an alarming rate [21]. For example, Mutanamba [22] analyzed the stability of cut-slope along the R523 road between Thathe Vondo and Khalavha. First, geological structures, properties of the soil, and groundwater regime were surveyed. Then, conventional analytical methods were used to classify the type of soil and produce slope stability charts. Finally, two types of failures occurring close to the road were identified from the findings: rotational slump failure and rockfall in the form of toppling [21]. The qualitative findings, however, could not be used to comprehensively describe the underlying mechanism of the active and post landslide, occurrence, and deposition. That is what the present research is attempting to explore.

Materials and Methods
Detailed methodologies were followed in this paper. These methods include toppling analysis, rotational analysis, transitional analysis, and finite element methods (Phase 2). Thulamela Municipality road (R523) was used as the case study. Historically, limit equilibrium methods (LEM or LE analysis) have been extensively used to analyze landslides [23][24][25][26][27]. The majority of reported cases have been based on the toppling and rotational analyses of landslides.
Kinematic analysis is an alternative set of techniques that can also be used for the analysis of rock slope stability. The term generically refers to any empirical method dealing with landslides and rockfalls from the premise of their specific geological structures, physical properties of the material involved, and groundwater regime, as well as slope geometry [28].
Kinematic analysis is commonly used to identify the translational failures that occur due to the formation of wedges or planes among others [29,30]. Kinematic techniques can be regarded as useful tools for the simplified analysis of rock slope stability. However, their relevance is limited to desktop studies beyond which limit equilibrium methods are better suited.
Limit equilibrium analysis finds relevance in the study of the complex failure of the heavily fractured and weathered soil and rock mass. Detailed methodologies of the study are fully explained below.

Toppling Analysis
Toppling analysis is one of the crucial analyses when studying the stability of slope in mountainous or highway areas. The appreciation of whether soil or rock is likely to slide or topple along a slope is generally based on a chart known as the toppling and sliding chart. Illustrated in Figure 3, the chart enables one to categorize slope instability as sliding, toppling, or both sliding and toppling. Limit equilibrium analysis finds relevance in the study of the complex failure of the heavily fractured and weathered soil and rock mass. Detailed methodologies of the study are fully explained below.

Toppling Analysis
Toppling analysis is one of the crucial analyses when studying the stability of slope in mountainous or highway areas. The appreciation of whether soil or rock is likely to slide or topple along a slope is generally based on a chart known as the toppling and sliding chart. Illustrated in Figure 3, the chart enables one to categorize slope instability as sliding, toppling, or both sliding and toppling. Hoek and Bray [31] who developed the chart argued that direct toppling occurs when the center of gravity of the discrete soil or rock mass lies outside of the overall base of the block. The conventional approach to soil or rock toppling and sliding is based on the above premises. In this study, all selected landslides were also analyzed for sliding and toppling possibility. The purpose of performing such analysis was to categorize the selected landslides into sliding or toppling and then suggest which types of failure should be expected within the study area. This analysis also provides a broad view of a mechanism that could be associated with slope failure.

Rotational Analysis
The rotational analysis is an LE technique that differs from the transitional analysis covered earlier. This analysis technique specifically applies to circular rock slope failure. Circular landslides usually occur within the weak rock or soil material of strengths of a magnitude similar to that of the induced stresses [24,[32][33][34][35][36]. In such cases, a simplifying assumption is made to exclude geological structures. This then leads to the production of landslides that occur along a circular or rotational failure profile. A graphical summary of this rotational analysis is provided in Figure 4.
The simulation was focused on the safety factor of the defined slope. The two-dimensional (2D) numerical software called SLIDE is based on the rotational analysis in Figure 4. SLIDE finds use in the analysis of slip surfaces and the determination of safety factors along landslides involving materials with soil-like properties. To simplify the results of the study, the first analysis was to simulate the safety factor of silt clay soil, clay soil, and, lastly, clay loam soil. Therefore, the safety factor of each slope was then simulated. Hoek and Bray [31] who developed the chart argued that direct toppling occurs when the center of gravity of the discrete soil or rock mass lies outside of the overall base of the block. The conventional approach to soil or rock toppling and sliding is based on the above premises. In this study, all selected landslides were also analyzed for sliding and toppling possibility. The purpose of performing such analysis was to categorize the selected landslides into sliding or toppling and then suggest which types of failure should be expected within the study area. This analysis also provides a broad view of a mechanism that could be associated with slope failure.

Rotational Analysis
The rotational analysis is an LE technique that differs from the transitional analysis covered earlier. This analysis technique specifically applies to circular rock slope failure. Circular landslides usually occur within the weak rock or soil material of strengths of a magnitude similar to that of the induced stresses [24,[32][33][34][35][36]. In such cases, a simplifying assumption is made to exclude geological structures. This then leads to the production of landslides that occur along a circular or rotational failure profile. A graphical summary of this rotational analysis is provided in Figure 4.
The simulation was focused on the safety factor of the defined slope. The two-dimensional (2D) numerical software called SLIDE is based on the rotational analysis in Figure 4. SLIDE finds use in the analysis of slip surfaces and the determination of safety factors along landslides involving materials with soil-like properties. To simplify the results of the study, the first analysis was to simulate the Sustainability 2020, 12, 8870 5 of 33 safety factor of silt clay soil, clay soil, and, lastly, clay loam soil. Therefore, the safety factor of each slope was then simulated.

Transitional Analysis
Like toppling, the transitional analysis was also introduced by Hoek and Bray [31]. The technique assumes that there is transitional sliding of a riding body along the plane of weakness as illustrated in Figure 5. It is crucial to indicate that all the forces pass through the centroid of the block because sliding does not experience rotation. In addition to this, it is always assumed in the limit equilibrium solution that "all points along the sliding plane(s) are on the verge of failure" [38]. This assumption makes the problem statically determinate and allows for the smooth calculation of the factor of safety (FoS). When analyzing the potential of a certain slope for failure, the criterion is always given to the critical slip surface and the determination of the safety factor. This criterion helps to select or identify the potential of the instability of the slope; however, the slope is mostly divided into slices with the consideration of forced and movement equilibrium acting on each slice. A 2D numerical software called SLIDE was utilized to analyze the slip surface as well as the determination of factor safety along with the selected landslides material using a variety of soil properties found within the selected section. In this study, eight methods were used to quantify the results; those methods are Bishop simplified [39], Janbu simplified and Janbu corrected [40], Spencer [41], Corps of Engineers Number One and Corps of Engineers Number Two [42], Lowe-Karafiath [43], and GLE/Morgenstern-Price [44]. All these methods were filtered within the model to produce reliable results of the slope safety factor. To simplify the results of the study, the first analysis was to simulate the safety factor of silt clay soil, clay soil, and, lastly, clay loam soil.

Transitional Analysis
Like toppling, the transitional analysis was also introduced by Hoek and Bray [31]. The technique assumes that there is transitional sliding of a riding body along the plane of weakness as illustrated in Figure 5. It is crucial to indicate that all the forces pass through the centroid of the block because sliding does not experience rotation. In addition to this, it is always assumed in the limit equilibrium solution that "all points along the sliding plane(s) are on the verge of failure" [38]. This assumption makes the problem statically determinate and allows for the smooth calculation of the factor of safety (FoS). When analyzing the potential of a certain slope for failure, the criterion is always given to the critical slip surface and the determination of the safety factor. This criterion helps to select or identify the potential of the instability of the slope; however, the slope is mostly divided into slices with the consideration of forced and movement equilibrium acting on each slice. A 2D numerical software called SLIDE was utilized to analyze the slip surface as well as the determination of factor safety along with the selected landslides material using a variety of soil properties found within the selected section. In this study, eight methods were used to quantify the results; those methods are Bishop simplified [39], Janbu simplified and Janbu corrected [40], Spencer [41], Corps of Engineers Number One and Corps of Engineers Number Two [42], Lowe-Karafiath [43], and GLE/Morgenstern-Price [44]. All these methods were filtered within the model to produce reliable results of the slope safety factor. To simplify the results of the study, the first analysis was to simulate the safety factor of silt clay soil, clay soil, and, lastly, clay loam soil.

Numerical Simulation Procedures for SLIDES
It is, however, noteworthy to state here that SLIDES was utilized to estimate the FoS values of the slope for several scenarios. The FoS values can be computed following various models (Bishop's, among others). The different outcomes can be compared at once. In terms of model building, the exercise starts with the creation of a new project. Similar to other modeling platforms, the new project required the delimitation of the model limitation in XY coordinates. For that, various X and Y coordinates defining the region were entered. The ultimate goal of this step was to draw the model of the region. Upon generating the boundaries of the model, the actual initial conditions of the simulation were defined next for the project. Inputs such as the statistics associated with groundwater conditions, the computational methods, and the failure directions were captured. selected section. In this study, eight methods were used to quantify the results; those methods are Bishop simplified [39], Janbu simplified and Janbu corrected [40], Spencer [41], Corps of Engineers Number One and Corps of Engineers Number Two [42], Lowe-Karafiath [43], and GLE/Morgenstern-Price [44]. All these methods were filtered within the model to produce reliable results of the slope safety factor. To simplify the results of the study, the first analysis was to simulate the safety factor of silt clay soil, clay soil, and, lastly, clay loam soil.  The correct values of input data and the appropriate selection of procedural approaches were of paramount importance. This is because the exercise determined how realistic the computer modeling of the actual rockmass associated with the study area would be. For a thorough interpretation of the results, several methods (Bishop simplified, Corps of Engineers #1 and #2, GLE, Janbu simplified and corrected, Lowe-Karafiath, Ordinary, and Spencer) of analysis were considered. The motivation was to compare various approaches to the problem. However, corresponding results were interpreted based on the performance of each approach.
The other important factor was the groundwater regime set at 9.81 KN per cubic meter. This value was estimated in the laboratory after performing several tests on soil material. Statistics of the yearly rainfall in the area were also considered in estimating the groundwater regime. It is important to stress that theoretical models underpinning all these computer programs are generally based on simplifying assumptions. These depart from the real behavior of soil and rock masses and diminish the value of output results. To circumvent the limitation, a degree of randomness and unpredictability is supposed to be allocated to the modeling frameworks. In this paper, the exercise was implemented by resorting to the probabilistic analysis built upon the Monte Carlo sampling method. Indeed, the Monte Carlo method refers to the allocation of randomness to a phenomenon by means of random number generators.
The next step was the definition of the boundaries of the excavation upon which the roads were constructed. The coordinates entered were based on field surveying and existing maps of the region. The inputs are rendered automatically and diagrammatically as the excavation model (see Figure 6). Finally, computational routines are selected and defined in terms of grid spacing of the limit equilibrium analysis of the excavation model. Other relevant input parameters are also incorporated before the model is finally solved. The convergence solution to the problem is exemplified in the results sections with the most likely FoS defined in the inset. Various scenarios were simulated and their output results were recorded for later analysis.

Numerical Simulation Procedures for RocPlane
RocPlane followed the same basic workflow as DIPS; however, the difference is the input parameters required. Here, the dips and dip directions of the set of joints, the upper face, and the slope face were entered. In addition to this, the dip, dip direction, cohesion, and trace length of tension cracks were incorporated. From this, the analysis was run as per theoretical models while the options for displaying results can be either probabilistic or deterministic (see Figure 7). Furthermore, using this method, the expected deposition of the falling wedge can also be predicted (see example in Figure 8).
The inputs are rendered automatically and diagrammatically as the excavation model (see Figure 6). Finally, computational routines are selected and defined in terms of grid spacing of the limit equilibrium analysis of the excavation model. Other relevant input parameters are also incorporated before the model is finally solved. The convergence solution to the problem is exemplified in the results sections with the most likely FoS defined in the inset. Various scenarios were simulated and their output results were recorded for later analysis.

Numerical Simulation Procedures for RocPlane
RocPlane followed the same basic workflow as DIPS; however, the difference is the input parameters required. Here, the dips and dip directions of the set of joints, the upper face, and the slope face were entered. In addition to this, the dip, dip direction, cohesion, and trace length of tension cracks were incorporated. From this, the analysis was run as per theoretical models while the options for displaying results can be either probabilistic or deterministic (see Figure 7). Furthermore, using this method, the expected deposition of the falling wedge can also be predicted (see example in Figure 8). The deterministic option presents the results as a force field of the slope with a wedge in it. The benefit is the possibility of seeing how the wedge would rotate in time steps once the stability conditions of the area are changed. The probabilistic option renders the RocPlane results in two ways: a histogram plot of FoS versus relative frequency or a scatter plot of FoS versus wedge weight, wedge length, a tension crack friction, and joint strength among others. Monte Carlo simulation has been used for this study; the setting for Monte Carlo simulation settings are shown in Figure 9.

Numerical Simulation Procedures for RocPlane
RocPlane followed the same basic workflow as DIPS; however, the difference is the input parameters required. Here, the dips and dip directions of the set of joints, the upper face, and the slope face were entered. In addition to this, the dip, dip direction, cohesion, and trace length of tension cracks were incorporated. From this, the analysis was run as per theoretical models while the options for displaying results can be either probabilistic or deterministic (see Figure 7). Furthermore, using this method, the expected deposition of the falling wedge can also be predicted (see example in Figure 8). The deterministic option presents the results as a force field of the slope with a wedge in it. The benefit is the possibility of seeing how the wedge would rotate in time steps once the stability conditions of the area are changed. The probabilistic option renders the RocPlane results in two ways: a histogram plot of FoS versus relative frequency or a scatter plot of FoS versus wedge weight, wedge length, a tension crack friction, and joint strength among others. Monte Carlo simulation has been used for this study; the setting for Monte Carlo simulation settings are shown in Figure 9.  The deterministic option presents the results as a force field of the slope with a wedge in it. The benefit is the possibility of seeing how the wedge would rotate in time steps once the stability conditions of the area are changed. The probabilistic option renders the RocPlane results in two ways: a histogram plot of FoS versus relative frequency or a scatter plot of FoS versus wedge weight, wedge length, a tension crack friction, and joint strength among others. Monte Carlo simulation has been used for this study; the setting for Monte Carlo simulation settings are shown in Figure 9.

Numerical Simulation of Slope Stability Using Phase 2
A Phase 2 numerical model was utilized to simulate stress, strain, and shearing of the material as the road has been developed. Three stages of road development were generated and the abovementioned parameters were simulated and compared as the stages progress. The model involves boundaries construction, mesh generation, boundary conditions, adding traction, field stress, excavating each stage, and, lastly, simulation. Like with all the previous tools, a new project was launched for each simulation scenario.
The next step was to generate the limits of the road excavation and rock face slopes using existing coordinates of the cross-sectional view of the road. The required coordinates were inputted based on information collected on the field and in the various archival maps of the region. Next was to define the project settings, that is, the initial and boundary conditions of the model to be used once the Phase 2 solver was launched. In addition to this, the number of iterations and the convergence criteria of the numerical computation were decided upon. Upon setting up the Phase 2 solver, it was necessary to enter the boundary geometry and conditions. This critical step was also the opportunity to define the boundary functions of the excavation. Boundary functions specifically refer to the loading conditions of the site. These may include the vibration created on the road surface as vehicles pass by and the effects of rainfall and groundwater amongst others.
At this point, model discretization and meshing were to follow. Based on the level of details required of the FEM solution, graded meshing was used so that the vicinity of the excavation had finer meshes. The meshing element was a 3-noded triangle as this type tends to be numerically stable. The next step was to define the field stresses at play around and within the rock mass model. Since the model is a surface excavation, the gravity stress field needed to be defined. Then, the vertical and horizontal distribution of the stress field was set at 1.5 throughout. It is well established that surface excavation experiences more high horizontal stresses than vertical stresses, therefore, the K-ratio (horizontal:vertical stresses) was set as 1.5. The remaining entries were captured based on the measured properties of the corresponding type of material.
The last step before the simulation model could be run was to assign appropriate properties corresponding to the various layers making up the road (see Figure 10 and Figure 11). A similar assignment could be done for the various stages of road construction. Finally, a choice of the type of outputs to be made available for later interpretation was made as indicated. Once the solver had converged, the results of the model could be analyzed, refined, and reported for the various stages of road construction. Some of the input parameters used are indicated in Table 1 below. The soil properties indicated in the table differ greatly because the area of study has different soil types, as such soil properties tend to differ.

Numerical Simulation of Slope Stability Using Phase 2
A Phase 2 numerical model was utilized to simulate stress, strain, and shearing of the material as the road has been developed. Three stages of road development were generated and the above-mentioned parameters were simulated and compared as the stages progress. The model involves boundaries construction, mesh generation, boundary conditions, adding traction, field stress, excavating each stage, and, lastly, simulation. Like with all the previous tools, a new project was launched for each simulation scenario.
The next step was to generate the limits of the road excavation and rock face slopes using existing coordinates of the cross-sectional view of the road. The required coordinates were inputted based on information collected on the field and in the various archival maps of the region. Next was to define the project settings, that is, the initial and boundary conditions of the model to be used once the Phase 2 solver was launched. In addition to this, the number of iterations and the convergence criteria of the numerical computation were decided upon. Upon setting up the Phase 2 solver, it was necessary to enter the boundary geometry and conditions. This critical step was also the opportunity to define the boundary functions of the excavation. Boundary functions specifically refer to the loading conditions of the site. These may include the vibration created on the road surface as vehicles pass by and the effects of rainfall and groundwater amongst others.
At this point, model discretization and meshing were to follow. Based on the level of details required of the FEM solution, graded meshing was used so that the vicinity of the excavation had finer meshes. The meshing element was a 3-noded triangle as this type tends to be numerically stable. The next step was to define the field stresses at play around and within the rock mass model. Since the model is a surface excavation, the gravity stress field needed to be defined. Then, the vertical and horizontal distribution of the stress field was set at 1.5 throughout. It is well established that surface excavation experiences more high horizontal stresses than vertical stresses, therefore, the K-ratio (horizontal:vertical stresses) was set as 1.5. The remaining entries were captured based on the measured properties of the corresponding type of material.
The last step before the simulation model could be run was to assign appropriate properties corresponding to the various layers making up the road (see Figures 10 and 11). A similar assignment could be done for the various stages of road construction. Finally, a choice of the type of outputs to be made available for later interpretation was made as indicated. Once the solver had converged, the results of the model could be analyzed, refined, and reported for the various stages of road construction. Some of the input parameters used are indicated in Table 1 below. The soil properties indicated in the table differ greatly because the area of study has different soil types, as such soil properties tend to differ.   There are two densities of the material (saturated and unsaturated density) shown in Table 1; these densities are based on two conditions (rainy and sunny conditions). Due to the fact the focus is on sunny conditions, saturated density was used. In the meantime, when measuring the densities,   There are two densities of the material (saturated and unsaturated density) shown in Table 1; these densities are based on two conditions (rainy and sunny conditions). Due to the fact the focus is  There are two densities of the material (saturated and unsaturated density) shown in Table 1; these densities are based on two conditions (rainy and sunny conditions). Due to the fact the focus is on sunny conditions, saturated density was used. In the meantime, when measuring the densities, for example, unsaturated densities were based on the volume of rainfall usually received which was, therefore, injected within the soil mass and the bulk density was measured. In saturated conditions, similar soil mass was dried up for 24 h then the bulk density was measured. Finally, the tension strength of the soil mass was considered zero throughout the entire clay soil types. The tension strength of the soil mass has less impact in this study which is why it was assumed constant.

Results and Discussions
Results of toppling analysis, rotational analysis, transitional analysis, and advanced numerical modeling are presented in this section. The results were produced by using a limit equilibrium (LE) analysis together with the FEM model.

Toppling Analysis
The concepts of toppling and sliding of blocks have been applied in six identified areas associated with slope instability. The analysis revealed that all the selected areas can be categorized as sliding rather than toppling as illustrated in . With reference to Table 2, the slope dimensions are illustrated as well as the ratio determined from the analysis; it was calculated that all slopes were experiencing sliding only. This could be due to the composition of the material that is mainly a clay soil and to the tension cracks surrounding the blocks. The results of this section correlate very well with some of the detailed studies on the application of toppling failure in soil mass slope; such studies include those of Goodman and Bray [45], Zanback [46], Bobet [47], Alejano et al. [48], and Amini et al. [49]. The above-mentioned authors have evidence of the occurrence of toppling failure in soil slope with most of the slope experiencing so-called complex toppling failure (toppling with circular failure). This type of failure is governed by soil layers that usually fail in a form of layers, however, tension cracks and the influence of rainfall have been documented to control such failure.
Similarly to the current case study, tension cracks were denoted to be the most important aspects that initiate the sliding of the blocks. However, the initiation was also noted to be controlled by the amount of rainfall (see Figure 18) which in turn reduces the resistance of the material against sliding. Although one analysis could not be used to draw a broad conclusion from the initial analysis (toppling analysis) it appeared that tension cracks and together with the climate change of the location had an influence.             Another aspect of toppling analysis is the common joint patterns around the area. These have been noted to influence the occurrence of the sliding of blocks since joints usually create wedges. Similarly, previous studies such as those of Alejano et al. [48], Amini et al. [49] have indicated joint sets also define the toppling failure of the soil mass. To be specific, Alejano et al. [48] have documented numerous case of various areas including open pit slope and mountainous slope among others; however, the common feature which could be denoted across the cases were joint sets that    Another aspect of toppling analysis is the common joint patterns around the area. These have been noted to influence the occurrence of the sliding of blocks since joints usually create wedges. Similarly, previous studies such as those of Alejano et al. [48], Amini et al. [49] have indicated joint sets also define the toppling failure of the soil mass. To be specific, Alejano et al. [48] have documented numerous case of various areas including open pit slope and mountainous slope among others; however, the common feature which could be denoted across the cases were joint sets that In summary, toppling analysis has shown that all active landslides are associated with sliding and to a certain extent a combination of toppling and sliding. Nevertheless, toppling analysis alone cannot give a broad view of the mechanisms associated with slope instability of the study area. That is why the next section presents the safety factors of six selected landslides estimated by simulation using SLIDES and RocPlane models. While the SLIDES software was used for rotational analysis, the RocPlane model took care of the transitional analysis. The purpose was to compare the impact of the slope angle on the recurrence of the slope instability.

Rotational Analysis
The rotational analysis covered in this section was done on three soil types found around the six landslides within the study area. The soil properties of the layers used for simulation are denoted in Table 1, wherein clay soil is indicated as upper layer A, silt clay soil indicated as upper layer B, and, lastly, loam clay soil properties are denoted in upper layer C. The above-mentioned three types of  Another aspect of toppling analysis is the common joint patterns around the area. These have been noted to influence the occurrence of the sliding of blocks since joints usually create wedges. Similarly, previous studies such as those of Alejano et al. [48], Amini et al. [49] have indicated joint sets also define the toppling failure of the soil mass. To be specific, Alejano et al. [48] have documented numerous case of various areas including open pit slope and mountainous slope among others; however, the common feature which could be denoted across the cases were joint sets that created a failure. The above-mentioned cases also correlate very well with the observations made in this study.
In summary, toppling analysis has shown that all active landslides are associated with sliding and to a certain extent a combination of toppling and sliding. Nevertheless, toppling analysis alone cannot give a broad view of the mechanisms associated with slope instability of the study area. That is why the next section presents the safety factors of six selected landslides estimated by simulation using SLIDES and RocPlane models. While the SLIDES software was used for rotational analysis, the RocPlane model took care of the transitional analysis. The purpose was to compare the impact of the slope angle on the recurrence of the slope instability.

Rotational Analysis
The rotational analysis covered in this section was done on three soil types found around the six landslides within the study area. The soil properties of the layers used for simulation are denoted in Table 1, wherein clay soil is indicated as upper layer A, silt clay soil indicated as upper layer B, and, lastly, loam clay soil properties are denoted in upper layer C. The above-mentioned three types of clay soil used for simulation differ slightly in their composition and properties; the details are denoted in Table 1. The simulation work underpinning the analysis focuses on the safety factor of several slices representing sections of the landslides. The results of the FoS of the clay soil slope and subsequent discussion are presented in the next paragraphs.
With reference to Figures 19-22, FoS values of the six locations with clay soil were found to be between 0.612 and 0.894 for all computation methods considered. The simulated minimum and maximum FoS values rate slopes in the area as unstable or prone to failure. It is clear that the material properties have played a major role in generating such low FoS. This is supported by field observations whereby the extent of slope failure was very wide in clay soils in comparison to other soils. In a sense, simulated results mimicked empirical analysis and observations. This provides a certain degree of reliability in the simulation model produced in this study.          In Figures 23-26, eight methods are compared for consistency. All have shown that clay loam soil slopes are expected to collapse if mitigating solutions are inexistent. This is made worse by hydrometric streams that usually play a major role in triggering landslides [51]. The area of study was noted to have a number of periodical streams cutting across. Streams and geological features have the ability to weaken the material strength by reducing the resistance of the material. It is, therefore, arguable that extreme rainfalls reported within the study area may have also contributed to the observed cases of slope instability in the six locations.
The second type of soil simulated is silty clay. Results of the silty clay soil slope show that all areas are stable to a large extent. Indeed, estimated FoS values range from 1.424 to 1.717 as summarized in Figures 23-26.   In Figures 23-26, eight methods are compared for consistency. All have shown that clay loam soil slopes are expected to collapse if mitigating solutions are inexistent. This is made worse by hydrometric streams that usually play a major role in triggering landslides [51]. The area of study was noted to have a number of periodical streams cutting across. Streams and geological features have the ability to weaken the material strength by reducing the resistance of the material. It is, therefore, arguable that extreme rainfalls reported within the study area may have also contributed to the observed cases of slope instability in the six locations. In Figures 23-26, eight methods are compared for consistency. All have shown that clay loam soil slopes are expected to collapse if mitigating solutions are inexistent. This is made worse by hydrometric streams that usually play a major role in triggering landslides [51]. The area of study was noted to have a number of periodical streams cutting across. Streams and geological features have the ability to weaken the material strength by reducing the resistance of the material. It is, therefore, arguable that extreme rainfalls reported within the study area may have also contributed to the observed cases of slope instability in the six locations.
The second type of soil simulated is silty clay. Results of the silty clay soil slope show that all areas are stable to a large extent. Indeed, estimated FoS values range from 1.424 to 1.717 as summarized in Figures 23-26.  In Figures 23-26, eight methods are compared for consistency. All have shown that clay loam soil slopes are expected to collapse if mitigating solutions are inexistent. This is made worse by hydrometric streams that usually play a major role in triggering landslides [51]. The area of study was noted to have a number of periodical streams cutting across. Streams and geological features have the ability to weaken the material strength by reducing the resistance of the material. It is, therefore, arguable that extreme rainfalls reported within the study area may have also contributed to the observed cases of slope instability in the six locations.
The second type of soil simulated is silty clay. Results of the silty clay soil slope show that all areas are stable to a large extent. Indeed, estimated FoS values range from 1.424 to 1.717 as summarized in Figures 23-26.     From the field observations and the empirical analysis, the slopes appeared to be less stable under normal conditions. However, slope stability is expected to be compromised or reduced when rainfall occurs. Here also, conclusions similar to those applicable to the clay soil slope can be drawn. For example, the material properties affect the FoS. Simulation results also make sense when compared to empirical observations. However, the extent of slope failure is moderate compared to pure clay soils. Finally, it is crucial to indicate that rainfall reduces the stability of the slope. The soil phase quickly changes during heavy rainfall and the soil material becomes liquid. At that stage, the shear strength of the material rapidly drops and later the inter-particle binding weakens and allows deformation to occur at any given time.
The last simulations were undertaken to verify the performance of the clay loam soil material. Results show that slopes were stable for all simulation methods considered (see Figures 27-30).  From the field observations and the empirical analysis, the slopes appeared to be less stable under normal conditions. However, slope stability is expected to be compromised or reduced when rainfall occurs. Here also, conclusions similar to those applicable to the clay soil slope can be drawn. For example, the material properties affect the FoS. Simulation results also make sense when compared to empirical observations. However, the extent of slope failure is moderate compared to pure clay soils. Finally, it is crucial to indicate that rainfall reduces the stability of the slope. The soil phase quickly changes during heavy rainfall and the soil material becomes liquid. At that stage, the shear strength of the material rapidly drops and later the inter-particle binding weakens and allows deformation to occur at any given time.
The last simulations were undertaken to verify the performance of the clay loam soil material. Results show that slopes were stable for all simulation methods considered (see Figures 27-30). From the field observations and the empirical analysis, the slopes appeared to be less stable under normal conditions. However, slope stability is expected to be compromised or reduced when rainfall occurs. Here also, conclusions similar to those applicable to the clay soil slope can be drawn. For example, the material properties affect the FoS. Simulation results also make sense when compared to empirical observations. However, the extent of slope failure is moderate compared to pure clay soils. Finally, it is crucial to indicate that rainfall reduces the stability of the slope. The soil phase quickly changes during heavy rainfall and the soil material becomes liquid. At that stage, the shear strength of the material rapidly drops and later the inter-particle binding weakens and allows deformation to occur at any given time.
The last simulations were undertaken to verify the performance of the clay loam soil material. Results show that slopes were stable for all simulation methods considered (see . Figures 23-26 suggest that the selected slope is not prone to failure provided it is not exposed to extreme rainfalls. It is also important to highlight that the selected slopes are stable but in a very steep environment. So, geological structures and streams cutting across are expected to affect the simulated FoS which is mostly controlled by soil properties. pure clay soils. Finally, it is crucial to indicate that rainfall reduces the stability of the slope. The soil phase quickly changes during heavy rainfall and the soil material becomes liquid. At that stage, the shear strength of the material rapidly drops and later the inter-particle binding weakens and allows deformation to occur at any given time.
The last simulations were undertaken to verify the performance of the clay loam soil material. Results show that slopes were stable for all simulation methods considered (see .         Figures 23-26 suggest that the selected slope is not prone to failure provided it is not exposed to extreme rainfalls. It is also important to highlight that the selected slopes are stable but in a very steep environment. So, geological structures and streams cutting across are expected to affect the simulated FoS which is mostly controlled by soil properties.

Transitional Analysis
The transitional analysis is well established for the study of slope instability in highways and surface mining. The method is used in this section to simulate the safety factor, the slope behavior, and the effect of the slope angle. It was implemented using the RocPlane program. Simulation results showed that an increase in slope angle leads to slope instability (see   Table 3). From Table 3, it has been observed that the steeper the slope the lower the FoS or the slope becomes prone to instability in the study area. The Monte simulations results correlate very well with the observations made across the study that most active and post slope instability are dominated along steep slopes rather than flat slopes.

Transitional Analysis
The transitional analysis is well established for the study of slope instability in highways and surface mining. The method is used in this section to simulate the safety factor, the slope behavior, and the effect of the slope angle. It was implemented using the RocPlane program. Simulation results showed that an increase in slope angle leads to slope instability (see   Table 3). From Table 3, it has been observed that the steeper the slope the lower the FoS or the slope becomes prone to instability in the study area. The Monte simulations results correlate very well with the observations made across the study that most active and post slope instability are dominated along steep slopes rather than flat slopes.     It was also observed that the slope angle and water pressure increase the magnitude of the driving force. This eventually slides the solid material to the slope toe as evidenced in Figures 34-39. These results suggest that the presence of geological features such as joints, heavy rainfall, steepness of the slope, and material properties are the common factors influencing the slope instability along the study area.    45 40 It was also observed that the slope angle and water pressure increase the magnitude of the driving force. This eventually slides the solid material to the slope toe as evidenced in Figures 34-39. These results suggest that the presence of geological features such as joints, heavy rainfall, steepness of the slope, and material properties are the common factors influencing the slope instability along the study area.          In summary, the results obtained from SLIDES and RocPlane were observed to agree well with each other. It was noted that material properties affect the stability of the slope in all situations. The steepness also plays a role in the stability of the slope with the simulated steep slope being prone to instability. The greatest challenge is that the wedge formation process cannot be simulated by the abovementioned technique. However, it can be inferred from the simulations that deformation or displacement within the Thulamela material is time dependent. This is because it is activated by heavy rainfall and intersecting geological features.   In summary, the results obtained from SLIDES and RocPlane were observed to agree well with each other. It was noted that material properties affect the stability of the slope in all situations. The steepness also plays a role in the stability of the slope with the simulated steep slope being prone to instability. The greatest challenge is that the wedge formation process cannot be simulated by the abovementioned technique. However, it can be inferred from the simulations that deformation or displacement within the Thulamela material is time dependent. This is because it is activated by heavy rainfall and intersecting geological features. In summary, the results obtained from SLIDES and RocPlane were observed to agree well with each other. It was noted that material properties affect the stability of the slope in all situations. The steepness also plays a role in the stability of the slope with the simulated steep slope being prone to instability. The greatest challenge is that the wedge formation process cannot be simulated by the abovementioned technique. However, it can be inferred from the simulations that deformation or displacement within the Thulamela material is time dependent. This is because it is activated by heavy rainfall and intersecting geological features.

Advanced Numerical Simulation of Slope Stability
The conventional techniques covered in the previous sections are generally limited to simplistic problems. However, many rock slope stability problems involve complexities such as in situ stresses, seismic loading, and material anisotropy. Because of this, advanced numerical simulation becomes a viable alternative. In this section, stages that are normally followed during road construction were simulated using the Phase 2 software. The aim was to look into the change in stress, strain, and shear as road construction progresses. The FEM-based results generated with Phase 2 are presented below.

Volumetric Strain of Slope in Phase 2
Simulation results rendered in Figure 40 for example show that during the first stage of road construction, the volumetric strain of the left and right slopes was minimal. This translated into little slope instability expected to occur from the upper part of the slope. Slope instability was nonetheless noted to gradually increase as the second stage of road construction was introduced (see Figure 41). However, it appears that much deformation was expected to occur on the sides of the road walls during stage two. This is because the slope angle and height had rapidly increased which implied that the slope stability was to be affected largely. Simulation results rendered in Figure 40 for example show that during the first stage of road construction, the volumetric strain of the left and right slopes was minimal. This translated into little slope instability expected to occur from the upper part of the slope. Slope instability was nonetheless noted to gradually increase as the second stage of road construction was introduced (see Figure 41). However, it appears that much deformation was expected to occur on the sides of the road walls during stage two. This is because the slope angle and height had rapidly increased which implied that the slope stability was to be affected largely.  Unlike other simulation tools used so far, Figures 40,41 indicate that a portion of the slope head is expected to slide into the excavation due to increased deformation. In reality, the model somehow agrees with reconnaissance and field observations as well as commonly accepted behaviors of slopes. The final stage of road construction was also simulated as shown in Figure 42. It can be seen that the Simulation results rendered in Figure 40 for example show that during the first stage of road construction, the volumetric strain of the left and right slopes was minimal. This translated into little slope instability expected to occur from the upper part of the slope. Slope instability was nonetheless noted to gradually increase as the second stage of road construction was introduced (see Figure 41). However, it appears that much deformation was expected to occur on the sides of the road walls during stage two. This is because the slope angle and height had rapidly increased which implied that the slope stability was to be affected largely.  Unlike other simulation tools used so far, Figures 40,41 indicate that a portion of the slope head is expected to slide into the excavation due to increased deformation. In reality, the model somehow agrees with reconnaissance and field observations as well as commonly accepted behaviors of slopes. The final stage of road construction was also simulated as shown in Figure 42. It can be seen that the depth of the road, the slope height, and the slope angle have increased, causing a proportionate Unlike other simulation tools used so far, Figures 40 and 41 indicate that a portion of the slope head is expected to slide into the excavation due to increased deformation. In reality, the model somehow agrees with reconnaissance and field observations as well as commonly accepted behaviors of slopes. The final stage of road construction was also simulated as shown in Figure 42. It can be seen that the depth of the road, the slope height, and the slope angle have increased, causing a proportionate increase in the volumetric strain of the slope. Ultimately, simulation results suggest that road construction has been disturbing the behavior and morphology of the slope. This has redistributed the total volumetric strain of the rockmass as blasting is used for excavation.

Maximum Shear Strain of Slope in Phase 2
The same simulation outputs used to produce Figures 40-42 were used to extract the maximum shear strain. The results of the Phase 2 simulation model show that the strain due to the shearing of the material has been increasing from the first stage to the final stage of road construction. The output results summarized in Figures 43-45 suggest that the resistance of the material had been gradually dropping with each stage. This has affected the shear strength of the material as a result. Such inference could explain not only the tension cracks observed along the slopes but also the manner in which the slope failure has been happening along the road.  Ultimately, simulation results suggest that road construction has been disturbing the behavior and morphology of the slope. This has redistributed the total volumetric strain of the rockmass as blasting is used for excavation.

Maximum Shear Strain of Slope in Phase 2
The same simulation outputs used to produce Figures 40-42 were used to extract the maximum shear strain. The results of the Phase 2 simulation model show that the strain due to the shearing of the material has been increasing from the first stage to the final stage of road construction. The output results summarized in Figures 43-45 suggest that the resistance of the material had been gradually dropping with each stage. This has affected the shear strength of the material as a result. Such inference could explain not only the tension cracks observed along the slopes but also the manner in which the slope failure has been happening along the road. Ultimately, simulation results suggest that road construction has been disturbing the behavior and morphology of the slope. This has redistributed the total volumetric strain of the rockmass as blasting is used for excavation.

Maximum Shear Strain of Slope in Phase 2
The same simulation outputs used to produce  were used to extract the maximum shear strain. The results of the Phase 2 simulation model show that the strain due to the shearing of the material has been increasing from the first stage to the final stage of road construction. The output results summarized in Figures 43-45 suggest that the resistance of the material had been gradually dropping with each stage. This has affected the shear strength of the material as a result. Such inference could explain not only the tension cracks observed along the slopes but also the manner in which the slope failure has been happening along the road.   One peculiar remark is the concentration of maximum shear strain around the middle part of the slope regardless of the stage of the road construction. This implies that material shearing is commonly expected around the middle section of the slope face. It is still not clear why the Phase 2 simulation model produced such an outcome.

Maximum Shear Stress of Slope in Phase 2
Shear stress and strain were also analyzed from the Phase 2 simulation outputs. By comparing  at the similar construction stages, it can be that the distribution of shear stress relates well with that of shear strain. However, the noticeable difference is that shear stress in Figures 46-48 steadily increases with the stages of road construction. This could be attributed to the fact that shear stress is generally proportional to the depth and width of the road [52][53][54].  One peculiar remark is the concentration of maximum shear strain around the middle part of the slope regardless of the stage of the road construction. This implies that material shearing is commonly expected around the middle section of the slope face. It is still not clear why the Phase 2 simulation model produced such an outcome.

Maximum Shear Stress of Slope in Phase 2
Shear stress and strain were also analyzed from the Phase 2 simulation outputs. By comparing  at the similar construction stages, it can be that the distribution of shear stress relates well with that of shear strain. However, the noticeable difference is that shear stress in Figures 46-48 steadily increases with the stages of road construction. This could be attributed to the fact that shear stress is generally proportional to the depth and width of the road [52][53][54]. One peculiar remark is the concentration of maximum shear strain around the middle part of the slope regardless of the stage of the road construction. This implies that material shearing is commonly expected around the middle section of the slope face. It is still not clear why the Phase 2 simulation model produced such an outcome.

Maximum Shear Stress of Slope in Phase 2
Shear stress and strain were also analyzed from the Phase 2 simulation outputs. By comparing  at the similar construction stages, it can be that the distribution of shear stress relates well with that of shear strain. However, the noticeable difference is that shear stress in Figures 46-48 steadily increases with the stages of road construction. This could be attributed to the fact that shear stress is generally proportional to the depth and width of the road [52][53][54].
Finally, shear stress was also found to be concentrated along the middle section of the slope face (see . This simply means that the failure of the slope is expected to be initiated along that stress zone. In fact, most slope failures were observed not to occur from the mountain top but rather from the central area of the slope as shown in Figure 49. Sustainability 2020, 12, x FOR PEER REVIEW 25 of 33 Figure 46. Simulated maximum shear stress in the rockmass at the first stage of road construction. Note that maximum shear stress simulated ranges between (0.00 × 10 −1 to 2.240 × 10 1 ).

Figure 47.
Simulated maximum shear stress in the rockmass at the second stage of road construction. Note that maximum shear stress simulated ranges between (2.00 × 10 0 to 5.00 × 10 1 ).

Figure 48.
Simulated maximum shear stress in the rockmass at the final stage of road construction. Note that maximum shear stress simulated ranges between (3.00 × 10 0 to 7.50 × 10 1 ).
Finally, shear stress was also found to be concentrated along the middle section of the slope face (see . This simply means that the failure of the slope is expected to be initiated along that stress zone. In fact, most slope failures were observed not to occur from the mountain top but rather from the central area of the slope as shown in Figure 49.  Note that maximum shear stress simulated ranges between (2.00 × 10 0 to 5.00 × 10 1 ).

Figure 48.
Simulated maximum shear stress in the rockmass at the final stage of road construction. Note that maximum shear stress simulated ranges between (3.00 × 10 0 to 7.50 × 10 1 ).
Finally, shear stress was also found to be concentrated along the middle section of the slope face (see . This simply means that the failure of the slope is expected to be initiated along that stress zone. In fact, most slope failures were observed not to occur from the mountain top but rather from the central area of the slope as shown in Figure 49.   Note that maximum shear stress simulated ranges between (3.00 × 10 0 to 7.50 × 10 1 ).
Finally, shear stress was also found to be concentrated along the middle section of the slope face (see . This simply means that the failure of the slope is expected to be initiated along that stress zone. In fact, most slope failures were observed not to occur from the mountain top but rather from the central area of the slope as shown in Figure 49.

Total Displacement of Slope in Phase 2
The total displacement of a slope is an important aspect to simulate as it gives an idea of the expected deformation as a result of rock excavation. The simulation results are summarized in Figures 50-52 for the three stages of road construction. Much development in the distribution of displacement can be seen throughout the stages. What is evident from Figures 50-52 is that displacement steadily increases from the first stage to the final stage of road construction.

Total Displacement of Slope in Phase 2
The total displacement of a slope is an important aspect to simulate as it gives an idea of the expected deformation as a result of rock excavation. The simulation results are summarized in Figures 50-52 for the three stages of road construction. Much development in the distribution of displacement can be seen throughout the stages. What is evident from Figures 50-52 is that displacement steadily increases from the first stage to the final stage of road construction.

Total Displacement of Slope in Phase 2
The total displacement of a slope is an important aspect to simulate as it gives an idea of the expected deformation as a result of rock excavation. The simulation results are summarized in Figures 50-52 for the three stages of road construction. Much development in the distribution of displacement can be seen throughout the stages. What is evident from Figures 50-52 is that displacement steadily increases from the first stage to the final stage of road construction.  Note that total displacement simulated ranges between (0.00 × 10 0 to 2.04 × 10 −1 ).

Figure 52.
Simulated total displacement in the rockmass at the final stage of road construction. Note that total displacement simulated ranges between (0.00 × 10 0 to 3.60 × 10 −1 ).
Lastly, the total displacement of the slope changes as the stages of construction progress. However, total displacement remains the highest on the road base and the foot of the slope at all stages. This has serious implications on the stability of the road foundation that requires a solid base and stabilizing gabions to curb the anticipated deformations and cracks.

Principal Stress of Slope in Phase 2
The last simulation outputs pertain to the minor and major principal stresses developed within the rockmass defining the slope. Figures 53-55 summarized the simulation outputs of the minor principal stress (symbolically known as σ1) at different construction stages. It can be seen from the stress distribution of the entire road excavation that the stress cloud is concentrated in the lower section of the slope. The findings are somewhat similar to the shear stress ones (see Figures 46-48 for comparison). Note that total displacement simulated ranges between (0.00 × 10 0 to 2.04 × 10 −1 ).

Figure 52.
Simulated total displacement in the rockmass at the final stage of road construction. Note that total displacement simulated ranges between (0.00 × 10 0 to 3.60 × 10 −1 ).
Lastly, the total displacement of the slope changes as the stages of construction progress. However, total displacement remains the highest on the road base and the foot of the slope at all stages. This has serious implications on the stability of the road foundation that requires a solid base and stabilizing gabions to curb the anticipated deformations and cracks.

Principal Stress of Slope in Phase 2
The last simulation outputs pertain to the minor and major principal stresses developed within the rockmass defining the slope. Figures 53-55 summarized the simulation outputs of the minor principal stress (symbolically known as σ1) at different construction stages. It can be seen from the stress distribution of the entire road excavation that the stress cloud is concentrated in the lower section of the slope. The findings are somewhat similar to the shear stress ones (see Figures 46-48 for comparison). Lastly, the total displacement of the slope changes as the stages of construction progress. However, total displacement remains the highest on the road base and the foot of the slope at all stages. This has serious implications on the stability of the road foundation that requires a solid base and stabilizing gabions to curb the anticipated deformations and cracks.

Principal Stress of Slope in Phase 2
The last simulation outputs pertain to the minor and major principal stresses developed within the rockmass defining the slope. Figures 53-55 summarized the simulation outputs of the minor principal stress (symbolically known as σ 1 ) at different construction stages. It can be seen from the stress distribution of the entire road excavation that the stress cloud is concentrated in the lower section of the slope. The findings are somewhat similar to the shear stress ones (see   Note that the concentration of principal stress σ1 quickly changes between construction stages. Indeed, stress levels around the surface and for shallow excavation are minimal, then, develop and provide some indication of slope weakening with the deepening road excavation.
The last note is to state that the major principal stress σ1 increases with the rate of excavation and blasting as evidenced in Figures 56-58. This implies that the strength of the road walls is compromised by the excavation.   Note that the concentration of principal stress σ1 quickly changes between construction stages. Indeed, stress levels around the surface and for shallow excavation are minimal, then, develop and provide some indication of slope weakening with the deepening road excavation.
The last note is to state that the major principal stress σ1 increases with the rate of excavation and blasting as evidenced in Figures 56-58. This implies that the strength of the road walls is compromised by the excavation.   Note that the concentration of principal stress σ1 quickly changes between construction stages. Indeed, stress levels around the surface and for shallow excavation are minimal, then, develop and provide some indication of slope weakening with the deepening road excavation.
The last note is to state that the major principal stress σ1 increases with the rate of excavation and blasting as evidenced in Figures 56-58. This implies that the strength of the road walls is compromised by the excavation. Note that the concentration of principal stress σ 1 quickly changes between construction stages. Indeed, stress levels around the surface and for shallow excavation are minimal, then, develop and provide some indication of slope weakening with the deepening road excavation.
The last note is to state that the major principal stress σ 1 increases with the rate of excavation and blasting as evidenced in Figures 56-58. This implies that the strength of the road walls is compromised by the excavation. Sustainability 2020, 12, x FOR PEER REVIEW 29 of 33 . Figure 56. Simulated major principal stress σ3 at the first stage of road construction.
. Figure 57. Simulated major principal stress σ3 at the second stage of road construction.

Conclusions
The concepts of toppling and sliding of blocks have been utilized in this study to identify the mechanisms behind the slope failure in the area of study. To this end, rotational and translational analyses were performed on selected slopes. The rotational analysis was performed using classical theoretical equations while transitional analysis relied on the RocPlane numerical software.

Conclusions
The concepts of toppling and sliding of blocks have been utilized in this study to identify the mechanisms behind the slope failure in the area of study. To this end, rotational and translational analyses were performed on selected slopes. The rotational analysis was performed using classical theoretical equations while transitional analysis relied on the RocPlane numerical software.

Conclusions
The concepts of toppling and sliding of blocks have been utilized in this study to identify the mechanisms behind the slope failure in the area of study. To this end, rotational and translational analyses were performed on selected slopes. The rotational analysis was performed using classical theoretical equations while transitional analysis relied on the RocPlane numerical software.

Conclusions
The concepts of toppling and sliding of blocks have been utilized in this study to identify the mechanisms behind the slope failure in the area of study. To this end, rotational and translational analyses were performed on selected slopes. The rotational analysis was performed using classical theoretical equations while transitional analysis relied on the RocPlane numerical software.
It was found that the failure mechanism of all the selected areas conforms to sliding rather than toppling. In other words, the selected areas appear to be dominated by sliding failure due to the fact that the slope material is mainly clay soil. Tension cracks were also found to be responsible for the initiation of the sliding motion of blocks. However, the initiation was noted to be controlled by the amount of rainfall which reduces the resistance of the material against sliding. Slope failure occurred for two reasons: first, wedges failed to anchor the ground due to the lubricating effect of rainwater and, second, tension cracks develop and eventually sliding takes place.
The rotational analysis was also performed in three different types of soil: clay, silt clay, and loam clay. The purpose was to confirm the mechanism behind the slope instability. The first simulation using clay soil slope yielded FoS values ranging from 0.612 to 0.894. The second set of simulations on silt clay produced 1.424 < FoS < 1.717. Similar results were produced with the third set of simulations. The extent of slope failure is expected to be moderate in silty clay areas compared to purely clay soils. From the outputs, it can be said that the slopes appear to be less stable in normal conditions but when extreme rainfall occurs, their stability gets compromised further.
Finally, the results of the transitional analysis using RocPlane indicate a more unstable ground with an increase in slope angle. This basically means that as the slope gets steeper, FoS drops. The other finding from the RocPlane simulations was that the slope angle and the driving force of the water pressure are the common factors that can initiate sliding. In conclusion, the transitional analysis also concurs with other methods. This is because here also joints, heavy rainfall, slope steepness, and material properties were found to be the main contributors to slope instability along the study area.
In the meantime, it was also observed that the results from the Phase 2 numerical simulation have shown that volumetric strain, shear stress, shear strain, total displacement, and sigma 1 and 3 components of the slope increase with the stages of road construction and, as such, they compromise the stability of the slope. Looking into the volumetric analysis of the slope, results showed that some minimal instability is to be expected during the first stage of road construction. The situation gradually worsens as road construction progresses while much deformation is expected to occur on the sides of road walls. These simulated results were noted to agree with other techniques such as rock fracturing on the sidewall of the excavation and other numerical simulations.
The maximum shear strain was also constructed using the Phase 2 modeling tool. A gradual increase of shear strain from the first stage to the final stage of road construction was noted. Similar conclusions were reached with shear stress as well as minor and major principal stresses. This result implies that the resistance and integrity of the material were gradually reduced with each stage of road construction, thereby affecting the shear strength of the material. It is concluded that the improper road construction, steepness of the slope and slope properties (soil types), multiple geological features, and streams cutting across the area are the common mechanisms behind the slope instability around the selected study area. It is anticipated that sophisticated methods such as fuzzy model can be implemented in the future to verify and further analysis, such methods are well documented in several studies [55][56][57].