Slope stability analysis of deep-seated landslides using limit equilibrium and finite element methods in Debre Sina area, Ethiopia

Slope failure is a recurring natural hazard in the western margin of the Main Ethiopian Rift and especially around the Debre Sina area. To minimize the damage caused by failure events, a detailed investigation of landslide-prone areas identified using numerical modelling plays a crucial role. The main aim of this study is to assess the stability of slopes and to evaluate and compare safety factors calculated by the different available numerical methods. Stability analyses of slopes prone to different types of failures were performed with different techniques. The stability was assessed for slopes of complex geometry composed of aphanitic basalt, porphyritic basalt, tuff, and colluvium (poorly sorted clayey sand to silty sand) using the limit equilibrium method and the shear strength reduction method based on finite elements. Furthermore, numerical analysis was done under static and pseudo-static loading using the horizontal seismic coefficient to model their stability during a seismic event. Satellite images were used to select failure-prone slopes based on slope properties and identified past landslides, as well as to derive structural and geological information for the numerical models. The slope stability analysis indicates that the studied slopes are unstable, and any small-scale disturbance will further reduce the factor of safety and cause failure. The slope stability of landslide prone hills in the study area strongly depends on the saturation conditions and the seismic load.


Introduction
Slope failure is a common phenomenon around the world, often severely affecting human lives. Slope failures cause casualties and damage infrastructure as well as private and public property (Keefer 1984;Varnes 1984;Nadim et al. 2006;Clague and Roberts 2012). The International Association of Engineering Geology (IAEG) database estimates that slope failures are responsible for roughly 14% of injuries and deaths caused by natural hazards worldwide (Aleotti and Chowdhury 1999). Due to the hilly and mountainous terrain, the highlands of Ethiopia experienced numerous slope instabilities. The Debre Sina area is one of the largest mountain chains with undulating topography and a bi-modal monsoonal season. The area is one of the largest deep-seated landslides in Ethiopia superimposed by several shallow landslides. The fragile and highly deformed nature of rocks makes this area vulnerable towards different types of slope failures. These failures cause a considerable loss of life and property, lead to displacement of people, and seriously impact agricultural land, dwellers, and infrastructure. The highlands of Ethiopia are the home for almost 60% of the country's population while experiencing high amounts of precipitation. The high relief is accompanied by deep valleys and gorges with active river incision (Ayalew 1999). The rugged morphology is accompanied by a complex geological setting due to the nearby rift system.
The use of remote sensing data and techniques for the identification and analysis of landslides is widespread and manifold. Applications comprise, but are not limited to, monitoring of landslide movements, identification of past landslides, and mapping of failure-prone hillslopes. Often remote sensing techniques are incorporated into multidisciplinary landslide studies also including geophysical and geotechnical field investigations as well as numerical modelling (Abuzied et al. 2016; Abuzied and Alrefaee 2019; Moradi et al. 2021). From the complementary use of remote sensing with other analysis techniques, discrepancies and limitations arise because not all informations are compatible across methodologies. Mechanical analysis of slope stability provides knowledge of parameters controlling landslides, aiming at entirely removing the guesswork (Cruikshank and Johnson 2002). In the 1990s, most rock slopes were evaluated solely with stereographic projections. Kinematic analysis was used to identify geological structures facilitating slope sliding or toppling (Bar et al. 2019). In the 2000s, slope stability modelling techniques evolved, and the more complex two-dimensional limit equilibrium analysis was developed in combination with numerical modelling techniques for isotropic rock masses (Bar et al. 2019). Numerical methods can be classified as continuum, discontinuum, and hybrid methods. Continuum methods include finite element method (FEM), finite difference method (FDM), finite volume method (FVM), boundary element method (BEM), and meshless methods, while discontinuum methods include discrete element method (DEM) and discrete fracture network method (DFM) (Jing and Hudson 2002;Jing 2003). According to Jing (2003), hybrid methods can be applied for coupled hydraulic flow and mechanical stress calculations. Numerical tools such as the limit equilibrium method (LEM), FDM, BEM, and FEM have been used by researchers for the analysis of slope stability problems (Griffiths and Lane 1999;Jing 2003;Cheng et al. 2007;Sarkar and Singh 2008;Liu et al. 2015;Morales-Esteban et al. 2015;Stianson et al. 2015). In this paper, continuum modelling was used as it is probably best suited for weak and jointed rock masses (Pradhan and Siddique 2020). Until today, the development of numerical models describing the dynamics of intact and fractured rock masses during slope deformation and failure is ongoing and continuously improving (Kundu et al. 2016). Due to their simplicity and their straightforward parameterization, LE models remain popular. A LE model mainly requires information regarding slope geometry and topography, as well as the geological setting. Furthermore, the acting static and dynamic loads and the hydrogeologic conditions need to be specified. Geotechnical material parameters are also required (Ayob et al. 2019). Some of those parameters, such as slope geometry, topography, and geological structures can be derived from remote sensing. However, LE models do not consider the ground behavior, and the safety factors are considered constant along the failure surface. Furthermore, nowadays, remote sensing can provide much more detailed information on geomorphology and geological structures, which cannot be incorporated into LE models due to their technical limitations.
In recent years, the ability of the finite element method to represent complex boundary and interface conditions has made FEM a common tool to assess slope stability. FE models are usually considered robust and very capable of solving complex problems without limiting a prior assumptions regarding the inclinations and locations of interslice forces (Hammah et al. 2004). The technique can address complex geometries and material heterogeneity, compute solutions for non-linear deformations, and consider coupled processes, such as pore pressure and seismic loading (Ayob et al. 2019). The application of FEM in slope stability analysis has become more common as it is able to describe progressive failure. There are several approaches to slope stability analysis using the FEM, from which the shear strength reduction (SSR) is one of the most common ones (Griffiths and Lane 1999). In the SSR, the shear strength envelope is systematically reduced, and the simulation is halted once deformations exceed a mesh-dependent limit or the solution does not converge. The factor reducing the shear strength is called the strength reduction factor (SRF). Once the solution of the FEM does not converge, a critical SRF value has been reached. At this point, the slope became unstable (Hammah et al. 2005(Hammah et al. , 2007. The FEM is based on continuum mechanics and divides the modelling domain into finite elements, typically based on triangulation to combine neighboring numerical grid nodes. Computing large deformations exciding the size of an element or even a detachment of elements from each other is not possible within the classical FEM (Kundu et al. 2016). For blocky rock masses, the computed displacement is small until rock failure occurs. Through the implementation of designated fracture elements, which consider spacing, aperture, infilling, and continuity, joints can be incorporated into the FEM (Jing and Hudson 2002). Depending on past geological and tectonic processes, the rocks can have a wide range of strength and deformability parameters, also causing a spatial heterogeneity (Sari 2019). However, the interaction between discontinuities and the intact rock mass mostly defines the stability when the natural equilibrium is disturbed. For a realistic mechanical model, the stress-strain relationship can be analyzed to identify the most probable failure mechanism (Elmo and Stead 2010). Integrating joint patterns into the rock mass is crucial for an adequate representation of realistic rock configurations. The mechanical properties of joints are controlled by the physical properties in discontinuities that affect the mechanical behavior by friction, compressive strength, weathering, and filling. Also, in comparison to the slice method used in LEM, the SSR does not require a predefined failure surface or the search for a minimum failure surface as the failure plane is an output of the SSR method (Sari 2019).
A number of studies have previously compared the results of slope stability analyses using the LE and FE methods (Griffiths and Lane 1999;Hammah et al. 2004;Khabbaz et al. 2012;Vinod et al. 2017;Zein and Karim 2017), achieving a good agreement between LEM and FEM for homogenous material and simple geometries but revealing an overestimation of the slope stability using LEM for complex geometry and heterogeneous material. However, Hammah et al. (2004) recommended adopting the FE method using shear strength reduction factor as an additional robust and powerful tool for design and analysis. As (Sari 2019) mentioned, it is required to compare the results of different methods to assess slope stability, since a single method might result in an incomplete representation of the stability conditions of rock slopes. Review of previous studies in the study area (Schneider, et al. 2008;Woldearegay 2008), revealed that the failure mechanisms of medium to large rock slides were assessed based on the distribution of rock and soil masses and field observations of features that indicated mass movements. To improve the understanding of the triggering of such large rock slides, numerical models of slope stability with LEM and an elasto-plastic FEM using the SSR technique with the Mohr-Coulomb failure criteria were carried out in the Debre Sina area, particularly in slope sections of Shotel Amba, Yizaba, Nib Amba, and Wanza Beret. These areas were severely affected by landslide incidences in recent years. In this research, the SSR technique was used to determine an SRF or factor of safety (FS) value that is associated with a slope at the verge of failure. Remote sensing was used to identify respective slopes and derive their geometry, parameters, and boundary conditions. Based on the model outcome, future slope performances are evaluated, and various combinations of loading conditions in the natural environment that may affect the area in the future are studied. The factor of safety is a very useful index for determining how close or far away a slope is from failure. Stability was revised in static and seismic conditions. In addition, comparisons were made between four different methods of calculating the FS using the LEM.

Description of the study area
Ethiopia is located close to the tectonically active East African Rift system which results in numerous landslides in many parts of the country. The Ethiopian highlands are susceptible to various types and sizes of landslides due to their variable topography and geology. The study area is located in the central-western highlands of Ethiopia north of Addis Ababa, which is part of the northwestern Ethiopian plateaus (Fig. 1). Debre Sina lies within a tectonically active region along the rift escarpment border. The area is characterized by undulating topography and the presence of alternating hills and valleys. The drainage pattern in the area represents parallel to sub-parallel dendritic patterns developed along faults and master joints in the hard rocks (Mebrahtu et al. 2021a). Geologically, the study area comprised mainly aphanitic basalt-porphyritic-agglomerate, ignimbrite-tuffvolcanic ash, porphyritic basalt-scoriaceous agglomerate, Tarmaber basalt and upper ignimbrite of Tertiary age, and unconsolidated deposits (colluvial and alluvial deposits) (Mebrahtu et al. 2020a) (Fig. 2). In particular, the geology of the 2D-model of the slope sections mainly consists of aphanitic basalt, porphyritic basalt, tuff, and colluvium (poorly sorted clayey sand to silty sand). The aphanitic basalt unit characterized by fine to medium-grained, has a vitreous appearance, clearly developed columnar joints, and intense fracturing, whereas the porphyritic basalt unit is medium to coarse-grained, massive, with plagioclase, and olivine phenocrysts (Mebrahtu et al. 2020a). The ignimbrite found in the study area is associated with vitric tuff. It is characterized by both vesicular and massive variety. The rock fragments include pumice, older ignimbrite, vesicular basalt, fine-grained glass material, and volcanic ash. These rocks of the area are of volcanic type and contain geological structures such as faults, lineaments, and fractures/joints. Intense fracturing, columnar jointing, and spheroidal weathering are very common features.
The morphological setting of the area is closely related to the extensive tectonic activity and deep fluvial dissection. It is typically characterized by an irregular hummocky topography, with deeply incised valleys and sharp steep slopes, gentle slopes, cliffs, and sharp escarpments. Moreover, the deep-seated landslides modify the slope morphology of the study area. In particular, the hummocky topography of the Yizaba-Shotel Amba area covered by colluvium materials, locally including pyroclastic sediments (Mebrahtu et al. 2020b). These deposits consist of unsorted to poorly sorted loose soil sediments (clayey sand to silty sand) and matrix-supported rock fragments, with large blocks of basalt toppled from upslope cliff faces (Mebrahtu et al. 2020b). The alluvial deposits are composed of unconsolidated sediments ranging in size from fine clayey sand, sub-rounded to rounded pebbles, cobbles, and boulders. The unit is characterized by highly variable width, thickness, and composition, both along and across the sequence. The most significant deposit occurs along the Dem Aytemashy, Robi, and Shenkorge rivers, northern, northeastern, and southeastern part of the study area where there is a relatively flat topography (Fig. 2).
The main geological structures identified in the present study area are faults, lineaments, and fractures/joints. The study area has been affected by several normal faults and lineaments, which mainly strike N-S, E-W, NE-SW, NNE-SSW, and NW-SE major trends of faulting and lineaments (Fig. 2). The investigated rock slopes are heterogeneous and dissected by 3-4 sets of discontinuities. It was evidenced during the field surveys that failures in such slopes are controlled by unfavorably oriented discontinuities. The slopes are affected by different types of slope failure. According to (Varnes 1978) classification, the most common types of landslides in the study area are rotational slides, translational slides, rockfalls and toppling, debris slides, and debris and earth flows (Figs. 2 and 3). The most prominent landslide phenomena were observed at Shotel Amba, Yizaba, Nib Amba, and Wanza Beret areas (Fig. 2).

Detailed field survey
A detailed field survey was conducted to examine the geology (lithology and structure), geometric, geotechnical, and groundwater conditions of the slopes with respect to slope stability analysis. In a first step, an inventory of past landslides was established using aerial photography from Google Earth imagery. Past landslides were analyzed with respect to their lithology and geological structure as well as their topographic slopes. From Landsat data, the geological structure of the study area was analyzed after the application of image enhancing techniques as well as band filtering using Erdas Imagine 2014. Faults and lineaments were identified by a structural analysis and used to identify the structural predisposition of the last landslides and possible future failure locations. Geomorphological analysis of the study area was based on a digital elevation model (DEM) derived from high-resolution satellite images. Relevant parameters included slope gradient and aspect ratio of the identified landslides. Detailed fieldwork investigations were carried out along selected traverse lines to cover the study area, focusing on the topographical and geological settings. After the identification of potential slope failure locations (Fig. 2), slope profiles were generated based on the geological data and the topographic maps. The geometric attributes of the slope models were extracted using various spatial manipulation tools. In this study, high-resolution satellite imagery such as Landsat images and ALOS PALSAR with 12.5-m resolution were used to map the lineaments and extract geometric profile of slopes along selected lines. The potential landslide locations selected in Figs. 5, 6, 7, and 8 were taken as the domain for the stability analysis. During the field survey, conditions of discontinuities (i.e., joint set infill, aperture, persistence, roughness, and spacing), fracture plane orientations, slope geometry, and intact rock strength measurements were also conducted (Mebrahtu et al. 2020a). Besides, the geomechanical properties of the different rock masses were determined during field survey using a portable instrument called Schmidt hammer on selected slope sections as per the guidelines suggested by ISRM (1981). From each slope, ten Schmidt hammer rebound values for each joint set were taken and the median ( R ) value was used to determine joint wall strength. The empirical relation developed according to Barton and Choubey (1977) was used to perform the uniaxial compressive strength of the rock using (Eq. (1)).
where c is the uniaxial compressive strength (UCS) in MPa, γ is the dry rock densities in kN/m 3 and R is the average Schmidt hammer rebound value.

Sample collection and determination of mechanical parameters
Slope stability assessment was conducted by coupling various parameters and integrating a variety of information from the geotechnical field and laboratory investigations.
(1) log 10 σ c = 0.00088 × × R + 1.01 c Debris/earth slide. d Debris flow. e Earth flow. f A quasi-rotational slide widening retrogressively with ponded spring water at the toe of the Wanza Beret landslide Geomechanical properties of required for numerical modellling were determined during field and laboratory investigations (Table 1). The rock samples from different zones of the selected slopes were collected for determination of the geomechanical parameters of the failure criterion (Fig. 4). The representative rock block samples were collected for each rock type and their respective geotechnical properties were measured experimentally. The rock sampling was carried out systematically in order to represent different rock formations along the slope sections. Laboratory experiments were carried out for determination of geomechanical properties of the intact rock. The rock samples were cored in the laboratory using diamond drilling and subsequently tested to determine the geomechanical parameters according to (ISRM 1981) specifications. The data obtained from the field were integrated with the laboratory results which were used to perform numerical modelling for stability analysis.
The laboratory testing included index property tests, uniaxial compressive strength (UCS), triaxial compressive strength (TCS), and Brazilian tensile strength (BTS) in the rock mechanics laboratory of Ruhr University of Bochum, Germany (Fig. 4). The geological and geotechnical field data were collected for characterizing identified slopes via geological strength index (GSI) system. GSI system is a major input in numerical modelling and has been used to calculate the deformation and strength parameters of the rock mass (Pradhan and Siddique 2020). RocData 5.0 software from Rocscience was used to calculate equivalent shear strength parameters (cohesion and friction) from Generalized Hoek-Brown (GHB) parameters (Hoek and Marinos   . The residual values of cohesion and friction angle were taken as 25% and 20% of the peak values for the slope stability analysis, respectively. The mean value of the tested samples was taken for the numerical solution of the slope stability analysis (Table 1). Minimum three specimens have been tested from each unit.

Failure criterion of rock mass and joints
The failure surface in highly fracture rock masses mainly runs along discontinuities and partly through the intact rock. The major and widely used failure criteria for numerical modelling of rock mass are Mohr-Coulomb criterion (Mohr 1900), Hoek-Brown criterion (Hoek and Brown 1980), GHB criterion (Hoek et al. 1995), and modified Mohr-Coulomb criterion . The Hoek-Brown failure criterion for rock mass is based on several empirical relationships that characterize the stress conditions associated to the failure of intact and rock mass. Unlike the Mohr-Coulomb linear criterion, Hoek-Brown is an empirically derived criterion which relies on the nonlinear increase in peak shear stress with confining pressure (Eberhardt 2012 (2) in an attempt to develop an empirical realtionship that could be scaled in relation to geological data (Hoek and Marinos 2007).
where ′ 1 and ′ 3 are the major and minor effective principal stresses at failure, respectively; ci is the uniaxial compressive (UCS) of the intact rock; and m and s are the material constants. The parameter m is equivalent to the friction of the rock while s is related to the degree of fracturing that relates to the cohesion of the rock mass (Eberhardt 2012). The concept of Generalized Hoek-Brown criterion was introduced (Eq. (3)). It involves the replacement of rock mass rating (RMR) by GSI system. Along with this, several terms m b , s, and a were introduced for good and poor quality rock masses having GSI values > 25 and < 25, respectively. The Generalized Hoek-Brown criterion is given as follows: where m b = m i exp (GSI − 100)∕28 , s = exp (GSI − 100)∕9 and a = 0.5 for GSI > 25 (good quality); s = 0 and a = 0.5 − GSI∕200 for GSI < 25 (poor quality). m i is the material constant that depends upon the type of rock, m b is the reduced value of m i which accounts for the strength by reducing effects of jointed rock mass that relies upon GSI values and disturbance factor ( D ), and s and a are the curve fiting parameters which are being detrmined by using GSI and D values. Several numerical modelling tools do not incorporate Hoek-Brown criterion but apply the more simple Mohr-Coulomb criterion.
To overcome this problem, a window-based program called "RocLab" was also developed to calculate equivalent shear strength parameters (Hoek et al. 2002). The stability of jointed rock mass is strongly influenced by discontinuities, especially at shallow depths or at low stress levels. Shear strength parameters such as cohesion and friction angle are the main factors that define ressiting force acting normal to the failure plane (Barton and Bandis 1990). The parameters in turn depend on the conditions of discontinuities, such as joint set roughness, shape and roughness of asperities, degree of alteration, continuity, aperture, and type and thikness of infilling material (Hoek and Diederichs 2006). The well-known model for estimating the shear strength of discontinuities is the Mohr-Coulomb criterion (Eq. (4)).
where f is the shear strength at failure (kPa), C j is the cohesion of joint (kPa), σ n is the effective normal stress (kPa), and j is the friction angle of joint (°). The effective normal stress acting on a particular discontinuity surface depends upon its orientation, depth, and weight imposed by the overburden and the local hydrological conditions (Hoek and Diederichs 2006). The shear strength of discontinuity surface within a jointed rock mass is the coupled effect of surface irregularities or asperities, strength, normal stress, and shear displacemnt along the potential sliding surface (Wyllie and Mah 2004). According to Barton and Bandis (1990), shear strength along the failure plane can be evaluated using emperical law of basic friction angle first proposed by (Barton 1973) nonlinear equation (Eq. (5)). The Barton's shear strength criterion can be expressed as follows: where τ f is the peak shear strength (kPa), σ n is the effective normal stress (kPa), ∅ b is the basic friction angle (°), JRC is the joint roughness coefficient, and JCS is the joint wall compressive strength (kPa). The joint roughness coefficient was evaluated during the field survey by comparing with the standard profile (Barton and Choubey 1977), whereas the unit weight of the intact rock was determined in the laboratory. Finally, all the above input parameters were imported into the RocData software, and then cohesion and friction angle along the failure plane were calculated as inputs for numerical modelling.
In this study, the rock structure and patterns of the joints are represented by using the Mohr-Coulomb constitutive model. The discontinuities present in the rock mass play a (4) τ f = C j + σ n tanϕ j (5) τ f = σ n × tan ∅ b + JRC × log 10 JCS σ n 1 3 significant role in controlling the strength and deformational characteristics. The behavior of the discontinuities is usually defined in the form of normal and/or shear stiffness. The stiffness of rock joints defines the deformation under both normal and tangential loads. The stress-strain relationship of joints is described by the joint stiffness parameters, which are crucial for a correct representation of the deformation in jointed rocks. Among others, (Barton 1972) suggested the following Eq. (6) for the estimation of the peak normal stiffness (MPa/m). The faults and interfaces also follow the Mohr-Coulomb failure criterion to evaluate the possibility of slipping failure along the faults ( Table 2). The normal stiffness may be defined as the normal stress per unit closure of the joint while shear stiffness of a joint is the ratio of peak shear stress displacement (Barton 1972). The normal stiffness of joints ( K n ) can be estimated from the moduli of rock mass, intact rock, and joint spacing (Eq. (6)), whereas the shear stiffness of joints ( K s ) was taken as Eq. (7). Barton (1972) stated that the measurement of in situ joint stiffness is expensive and time-consuming and therefore recommends a method for determination of normal and shear stiffness of joints based upon the deformation properties of the rock mass and intact rock (Eqs. (6) and (7)). According to Pain et al. (2014), the shear stiffness in the present study was estimated as one-tenth of the normal stiffness. The rock mass modulus for normal stiffness estimated using the Hoek-Brown criterion and the GSI system for different rock types were determined. The normal stiffness of joint is represented as (Barton 1972) follows: where K n is the normal stiffness of joints, E m is the modulus of rock mass, E i is the modulus of intact rock, and L is the mean spacing of joint. The shear stiffness of joint is represented as (Barton 1972) follows: where K s is the shear stiffness of joints, G m is the shear modulus of rock mass, and G i is the shear modulus of intact rock, and L is the mean spacing of joint.

Model generation
A 2D actual geological model of the study area was first built using the geological modelling method. The geometry of the model for each slope was generated by explicitly introducing in the model domain the faults which mostly influence structurally controlled failure mechanisms (i.e., fault 1 and fault 2). The model development for the landslide slope sections includes the selection and definition of the problem geometry as well as boundary conditions. The input parameters used to determine the safety factor of rock slope failures were slope geometry, discontinuity orientation, shear strength parameters, and unit weight of the rocks. Further, a suitable material model and associated properties need to be selected. The LEM is based on the equilibrium of force and moment. The FEM solves the stress-strain relationship to determine response of the rock. For this study, slope analyses were performed using Rocscience SLIDE2 and RS2 for LE and FE methods, respectively. The geometry was implemented for the different slope sections, and each rock strata in every slope section was assigned with the corresponding properties for the rock mass and joints accordingly in the RS2 model as given in Table 1. Then the boundary conditions were assigned to the slope model (Figs. 5a, 6, 7, and 8a). The boundary between the different subsurface materials (lithology and structure) as well as the varying thickness of soil layers was obtained from the contrast of 2D seismic refraction and detail field survey (Mebrahtu et al. 2020b). During the field survey, information on lithology types, geological structures, rock strength, saturation, and groundwater conditions of the slopes was obtained (Mebrahtu et al. 2020a(Mebrahtu et al. , b, 2021a. The plane strain analysis is being performed with metric units and Gaussian elimination solver. Stress analysis is conducted by 500 maximum iterations with a tolerance of 0.001. The boundary at the base of the FE model was fixed for displacement in the x-and y-directions, while the vertical side boundaries were fixed for displacement in the x-direction (Figs. 5b, 6, 7, and 8b). The Mohr-Coulomb criterion is used to define the intact rock and joint strength characteristics for an elasticperfectly plastic behavior. The modelling is being performed under gravity loading by discretization with a six-nodded graded triangular finite elements with increased density near the faults surface as shown in Figs. 5b, 6, 7, and 8b. Gravitational stress field with horizontal to vertical in situ stress ratio of unity was adopted. According to Eberhardt et al. (2003), stress was initialized assuming horizontal to vertical stress ratio of 0.5. The mesh was made up from approximately 434,318 nodes and 216,659  (Figs. 5b, 6, 7, and 8b), though the results have been shown to be independent of the mesh density. From the springs at the base of the slopes, the groundwater level can be derived and used to specify the hydraulic conditions. The hydraulic condition was determined based on the spring locations and assumed to be in a steady state. Groundwater is a crucial factor in landslide initiation. Increasing groundwater levels and groundwater flow are critical for slope failure because it induces high pore-water pressure reducing the frictional strength. Especially over pressurized groundwater flow is associated with flow-type slope failure. The slope stability analysis of the slopes was carried out using LE and FE methods in terms of factor of safety.

Limit equilibrium method
The slope stability analysis was compiled and synthesized from extensive field investigations, rigorous calculations, and laboratory experiments. The slope stability analysis is performed with the limit equilibrium methods based on assumptions regarding the shape of the sliding surface. LE methods are used extensively for slope stability analysis and use the Mohr-Coulomb failure criterion to determine the shear strength along a slip surface (Deng et al. 2015). A state of limit equilibrium exists when the mobilized shear strength is expressed as a fraction of the shear stress. The slope failure is considered to occur along a pre-defined slip surface. At failure, the shear strength is fully mobilized along the critical slip surface. The critical shear stress at slope failure is the shear strength of the soil. In saturated soils, shear strength for an effective stress analysis is often given by the Mohr-Coulomb failure criteria Eq. (8).
where, τ = shear stress (kPa), c′ = effective cohesion (kPa), σ' = effective normal stress on the surface of rupture (kPa), ϕ′ = the effective angle of internal friction (°), and F s = factor of safety. The FS for a slope failure is calculated as the ratio of the available shear strength to the mobilized shear strength (Eq. (9)). The FS is a very common method for evaluating the stability of slopes. In theory, a FS of 1 means the driving and resisting forces are at equilibrium. Limit equilibrium methods are relatively simple in their application, compared to numerical methods (Eberhardt et al. 2003;Tang et al. 2017). SLIDE2 (Rocscience Inc. 2018) is used to calculate the stability of different slip surfaces using vertical slices. The forces acting on the rock mass are computed for discrete slices along the profile constrained by the surface and the failure plane. The ratio of resisting forces to driving forces at equilibrium defines the FS. LE methods have been broadly applied since they produce satisfactory FS results that can be corroborated from its basic ideologies. However, the LE methods are rather basic in their form as they do not fully consider the stress-strain relationship of the soil, which is also essential for slope stability evaluation. Limit equilibrium analysis is used to give an estimate of the FS and does not manifest information regarding the deformations associated with failure.
There are various LE methods to solve the force and moment equilibrium equations for a body sliding on circular and non-circular surfaces. The circular and non-circular LEM only consider force and momentum of the total mass. Internal conditions are not considered. In this study, four well-known techniques, namely Bishop's simplified method (Bishop 1955), Janbu's simplified method (Janbu 1968), Spencer's method (Spencer 1967), and Morgenstern-Price method (Morgenstern and Price 1965) were employed to locate the critical slip surfaces in the heterogeneous rock mass using the grid search and auto refine search tools provided by the program. Each point in the slip center grid represents the center for rotation of a series of slip circles. These methods are commonly used due to relatively adequate accuracy while calculating the FS and for establishing a common platform for conducting the comparative study between LE and FE methods. Hence, FS was calculated in saturated condition under both static and dynamic conditions (Figs. 5, 6, 7, and 8).

Finite element method
Finite element analysis was performed using RS2 v.9.0 software (Rocscience Inc. 2020). The application of FEM can overcome limitations in LEM, because instead of just the FS, the maximum shear strain, total displacement, and yield elements of the slope can be evaluated. The SSR method using FE has been applied regularly for the Mohr-Coulomb failure criteria reducing cohesion c and the angle of internal friction ϕ by the SRF until the solution does not converge anymore (Zhou et al. 1994;Griffiths and Lane 1999;Hammah et al. 2007;Gover and Hammah 2013). Once unresolved forces occur and the displacement at a node of the FE mesh becomes too large, the solution does not converge anymore (Kainthola et al. 2012). The reduction factor that causes the FE model not to converge is called critical reduction factor and associated with slope failure is considered equivalent to the slope's factor of safety. The SRF parameters are given as follows (Eqs. (10) and (11)). The FS is obtained by dividing the base strength by the lowest strength at which the slope is stable.
(10) C r = C SRF where, c is the cohesion, ϕ is the angle of internal friction, C r and ϕ r are the reduced shear strength parameters, and SRF is the shear strength reduction factor.
Numerical techniques are also useful to see the effect of the variation of the input parameters on the overall response of the rock structures. Geology, discontinuities, material properties (e.g., normal stiffness, shear stiffness, shear strength, and deformability), constitutive equations, failure criterion, groundwater pressure, external loads, in situ stresses are taken as input parameters based on the requirements of the methodology deployed. To obtain realistic results from the FE analysis, some researchers strongly suggested to include the effect of discontinuous media in the analysis (Styles et al. 2011;Agliardi et al. 2013;Satici and Ünver 2015). The inclusion of joints in the FE models could produce an altered failure surface (Hammah et al. 2008;Fu and Liao 2010). In low stress environments such as slopes, the joint properties have a greater influence on the rock mass behavior than the properties of the intact rock. The occurrence of discontinuities changes the stress distribution around the rock mass because the discontinuities alter the mechanical behavior of the rock mass as they often occur along areas of stress concentration and high deformation rates (Barton and Choubey 1977). Slope stability analysis of Shotel Amba, Yizaba, Nib Amba, and Wanza Beret slope sections was further evaluated using the FEM in terms of the critical SRF to calculate the stability conditions of the slopes. Hence, the critical SRF was calculated using RS2 software under different anticipated conditions. The landslide stability of the selected slope sections was analyzed for shear strain, slope stability, and total displacement using the FE method. The RS2 software was utilized for constructing the slopes in Shotel Amba, Yizaba, Nib Amba, and Wanza Beret (Rocscience Inc. 2020). The models were divided into two cases: (i) case 1 is the initial state in a static condition before the earthquake and (ii) case 2 is the dynamic slope model incorporating a seismic event to simulate the influence of an earthquake event on the selected slope sections. For the study area, three sources of seismic load can be identified: the Afar depression, the escarpment, and the Ethiopian Rift system (Mammo 2005). All three sources are within the vicinity of the study area, which is situated along the western margin of the tectonically active Main Ethiopia Rift (MER). The earthquake coefficient of peak ground acceleration (PGA) values for the Afar area ranging from 0.16 g (EBCS 1995) to 0.75 g (Group 1999) for the 0.01 annual probability. Seismic shaking is represented by a horizontal force calculated based on the multiplication of the total sliding mass with the so-called horizontal earthquake coefficient αh. This approach is usually called (11) tan( r ) = tan SRF pseudo-static analysis (Pyke 2002). Here, the coefficient αh = 0.3 is chosen for Shotel Amba and Wanza Beret and αh = 0.2 for Yizaba and Nib Amba slope sections.

Limit equilibrium analysis
The model outputs of the SLIDE2 2018 program for the Mohr-Coulomb material types using grid search methods are shown in Figs. 9, 10, and 11. Factor of safety was calculated under different anticipated conditions, i.e., static condition (without earthquake loading) and dynamic condition (with earthquake loading) as well as water pressure. The calculated FS using the LE modelling are presented in Table 3. Between the four LE methods, MPM resulted in the highest and JSM resulted in lowest FS but the differences are marginal. The calculated FS without seismic load ranged between a minimum of 2.16 and a maximum of 11.31 (Table 3) and between a minimum of 0.92 to a maximum of 4.76 with seismic load for the various cases evaluated in this study (Table 4). This clearly confirms that there is possibility of a sliding circular (rotational) slip failure in the studied slopes. From the Shotel Amba to Wanza Beret cross sections (Figs. 5,6,7,and 8), the Shotel Amba cross section attained the lowest FS of 0.92 with seismic load and 2.16 without seismic load in the 2D limit equilibrium analysis as shown in Fig. 8. It can be seen from Table 4 that the conventional method for calculating the FS with seismic load is the JSM, which has a value of 0.92 in Shotel Amba slope section. The highest FS is produced by MPM with a FS value of 1.02. While BSM and Spencer methods produce almost similar FS values and failure surfaces than the MPM, the JSM results in lower FS values and different failure surfaces.
The JSM generates lower FS because of its simplicity compared to the BSM and Spencer methods, which are both more rigorous methods. This leads to some differences in the FS values, though the rigorous methods are often considered more reliable (Duncan and Wright 1980). From Table 3, it is observed that Nib Amba was the only slope section with very high FS using the four LE methods. Among the slope sections, the Nib Amba slope section differs very high from the other sections. The likely reason for that is when computing the FS, the LEM did not consider joints act as initiation points of the failure plane. In the slope section of Nib Amba, there are weakness zones such as faults and joints distributed along the lithologic boundary between the porphyritic basalt and tuff. The FS values computed without seismic load via BSM, JSM, SM, and MPM are all independently stable, while the FS values computed with seismic load indicate slopes close to failure. This shows that the stability of the area becomes unstable when there is a combination of saturation and seismic load. The upper part of the slopes is dominated by colluvial deposits or intensively fractured rocks. In this case, circular failure surfaces are expected.

Finite element analysis
The vulnerable slopes in the study area were evaluated, and stability grade was quantified by means of critical SRF, as shown in Figs. 12, 13, 14, and 15. The FEM solution for Yizaba and Nib Amba does not converge for αh = 0.3. Once the FEM does not converge within a specified number of iterations, the slope is considered unstable using the SSR technique. The results show that the slopes are unstable as the seismic load increases. Using pseudo-static loading and gravitational forces, the slope sections have with the FEMbased RS2 software to model their stability under a seismic event (Figs. 5,6,7,8,12,13,14,and 15). The RS2 program accurately calculated the failure planes which align partly along fault surfaces and partly pass through intact rock masses (Figs. 12,13,14,and 15).
The slope failures in the rock masses in the FEM are more complex than in the LEM as joints play a crucial role in the development of the failure surfaces (Wyllie and Mah 2004). This is in agreement with previous studies (Styles et al. 2011;Agliardi et al. 2013;Satici and Ünver 2015). Results for the dynamic case are presented in Figs. 12, 13, 14, and 15. From the total displacement and the maximum shear strain contours, a potential shear failure surface can be derived. Shear strain contours by Mohr-Coulomb criteria of each slope were represented in Figs. 12, 13, 14, and 15. The critical SRF or FS for the cases evaluated in this study ranged between a minimum of 0.58 and a maximum of 1.03 (Table 4). This indicates that the slope stability is low due to the weak rock mass. On the other hand, for the results from the static case, the critical SRF ranged between a minimum of 1.51 and a maximum of 2.01 (Table 3). The slope at Shotel Amba showed that the SRF of the slope without seismic load is 1.51, whereas the SRF of slope with seismic load is 0.58 (Fig. 12).
The results from this simulation are shown in Figs. 12,13,14,and 15 showing the total displacement and the maximum shear strain. In comparing the resulting slope failure surfaces using the LE and FE analysis, a large difference in the shape   (Figs. 11 and 15). The critical failure surfaces resulting from the LE analyses were perfectly circular due to the selected search criterion (Figs. 9, 10, and 11), the FE method produces a near-circular zone of failure surfaces near the toe of the slope (Figs. 13 and  15). The FE analyses produce a better-defined failure path than the LE analyses. The FEM computes the failure region without a-priori assumptions with respect to the slip surface in opposition to the LEM. Furthermore, using the FE method, it is possible to compute the total displacement of rock and soil from the input data. Figures 12 and 13 show the critical failure surface and shear strain at the time of failure. Areas of maximum shear strain indicate the likely failure pathway that would develop through the modeled rock mass. The maximum shear strain along the critical failure surface is found to be 18.6 in Shotel Amba and 2.35 in Yizaba slope sections. The maximum total displacement in Shotel Amba and Yizaba is found to be 225 m (Fig. 12) and 36.1 m (Fig. 13), respectively. This is in agreement with the findings of Kropáček et al. (2015), who found that topographic profiles show that the maximum estimated thickness of the active main landslide body is about 150-200 m. The maximum total displacement in the Shotel Amba and Yizaba slope sections coincides with the results of the topographic profiles (Figs. 12 and 13). The maximum stress concentration along with the critical FS indicates a middle zone collapse and subsequent failure of the slope. In the slope section of Yizaba, the shear strain is high along the failure plane. The plane originates from the tension zone and continues towards the toe of the slope (Fig. 13). The failure is a block sliding with the displacement caused by the low strength of the colluvial deposits (poorly sorted clayey sand to silty sand) and tuff. In the slopes Shotel Amba and Yizaba, the middle zone of the slope moves to the right, and a crack opens in the middle of the slope, showing that deepfailure mechanism (Figs. 12 and 13). It is also observed that the failure surface and the location estimated by the numerical method match with the field observation (Fig. 3f).
In the field, many tension cracks can be observed indicating possible future slope failures (Fig. 3). A factor of safety of 1.03 has been obtained through the finite element analysis of the Nib Amba slope with seismic load (Fig. 14). The computed FS with seismic load for the cases evaluated in this study is less than 1 except for Nib Amba (Fig. 14) but which is very close to one. So, there is a risk of failure of all slopes during an earthquake. The calculated FS of Nib Amba slope section using FEM becomes less than 1 when the seismic load increases.
The analysis result revealed that at the Shotel Amba, Yizaba, Nib Amba, and Wanza Beret slope sections, the FS is less than 1 during dynamic condition (with earthquake loading) and saturated conditions (Figs. 12,13,14,and 15). The slope was found to be unstable, indicating that seismicity and pore-water pressure significantly contribute to the destabilization of the slope sections as the study area falls in a high seismic zone (Asfaw 1986). Furthermore, the geometry of the slopes and orientation of faults also played a greater role in destabilizing the study area. Furthermore, the model results at the Shotel Amba and Yizaba slope sections indicate that the scarp is continuously unstable due to its high inclination (Figs. 12 and 13). This is consistent with field observations of rupture of large sizes that occur regularly.

Discussion
The limit equilibrium and finite element methods were applied for slope stability assessment of the landslide prone hills in Debre Sina area with a complex geometry and heterogeneity in material properties. The results of LEM analyses revealed that for the Shotel Amba, Yizaba, Nib Amba, and Wanza Beret slope sections, the FS is more than one, and the slopes were found stable during saturated conditions and seismic load. Whereas the results FEM analyses revealed that all slope sections at dynamic and saturated conditions, the critical SRF is less than 1 and the slopes are unstable. This indicates that the pore-water pressure and horizontal seismic acceleration significantly contributed to the instability of slopes in the study area. The findings indicated that differences between the LE and FE methods are evident when a heterogeneous and complex slope geometry is analyzed. Only the FE method can take benefit of the detailed structural information derived from remote sensing data. The calculation of the FS shows a significant difference between the LE and FE approaches. The FS from LEM with seismic load ranged between 0.92 and 4.76, while critical SRF from finite element analysis was between 0.58 and 1.03 for 2D FEM simulations. Among the slope sections, the Nib Amba slope section differs very strongly from the other sections. The slope was evaluated as unstable (critical SRF = 1.03). The results, i.e., critical SRF and shear strain contour distribution within the slopes are being supported by existing field conditions. This result indicated that there is a strong agreement between the finite element analysis result and the types of failure manifested in the field. This area is the nearest to the most seismically active regions in the country. The results therefore indicate that the seismic shaking due to the tectonic unrest is an important predisposing factor, or even trigger, facilitating landslides in the area (Mebrahtu et al. 2020a). The critical SRF from FEM is significantly lower than the FS from LEM in all studied scenarios. The low value of SRF emphatically suggests that the rock slope can be vulnerable to failure under the influence of any triggering force. The results of the FEM showed that from the static and dynamic analysis, it has been concluded that the slope sections have the lowest FS compared to the calculated FS using the LEM. At the critical SRF, the slope material undergoes significant strength degradation as the failure surface is formed by coalescing fractures, which converge with the kinematically feasible release plane. Any further increment in SRF value accelerates the slide mass movement which causes an increased total displacement of slope mass material.
The lithologic units of the study area are anisotropic in their behavior, and their stress-strain behavior is quite variable, due to the presence of volcanic ash. The rock strength test results indicate that the rock slope sections with medium to high strength, rock slope failure were controlled mainly by structural discontinuities. Due to the adverse orientation of discontinuities, structurally controlled failures are prominent in the study area. Moreover, geological (lithology and structure), extensive rainfall, and geotechnical factors are the major factors accounting for frequent along the slopes. The elasto-plasticity behavior of the ash due to the presence of water poses a serious threat during the periods of prolonged and heavy rainy seasons. Rainfall is one of the most significant triggering factors for landslide occurrence in this study area. The effect of rainfall infiltration on slope could result in changing soil suction and positive pore-pressures as well as raising soil unit weight, which leads to a reduction in soil and rock shear strength. Furthermore, springs originating at different parts of the slope were identified, and it revealed that shallow groundwater table acts as slope-destabilizing factors in addition to rainfall, which triggers stability of the slope mainly during the rainy season (Mebrahtu et al. 2021a). The SRF due to a rise in pore-water pressure might lead to slope failure. The results of the LE and FE calculations indicate that the area is highly susceptible to sliding when it gets moist. The build-up of high hydrostatic pressure results in a reduction of the effective normal stresses bringing the slope closer to failure. Moreover, the concave shape of the slope enhances groundwater flow into the landslide body resulting in comparably high water tables (Mebrahtu et al. 2021a). The results of the simulations show that both the heavy rainfall and seismic events in the tectonically highly active region as well as the preconditioning complex litho-structural layers play an important role in the destabilization of the studied slopes. Thus, the studied slopes are unstable under saturated conditions and seismic events, and any small-scale disturbance will further reduce the factor of safety and probably causing failure (Mebrahtu et al. 2021b). The investigated area is dissected by major normal faults, which are organized in a series of parallel east-facing steep faults. The slopes are located in close proximity to the fault systems that produced a complex displacement across and along the escarpment, manifesting oblique continental rifting (Kiros et al. 2018). Due to such reasons, slopes near Shotel Amba, Yizaba, Nib Amba, and Wanza Beret are very prone to failure.
The numerical model simulations reveal that the LE method resulted in a relatively defined critical failure surface, whereas the FE method shows a much wider zone of critical failure. The critical slip surface developed in FE is near the slope face, while in LE, the critical slip surface is far away from the slope face. It was found that the FS results from LE and FE methods are significantly different. Since the LE and FE methods employ completely different numerical techniques, the simulations did not result in similar FS values or failure surfaces. According to the results of the LE analyses, the slopes are classified as stable, but as a result of the FE analysis, the slopes are classified as unstable. From the FE analysis, it can be inferred that the study area is critically unstable, and any small-scale disturbance, such as heavy precipitation or earthquakes, will further reduce the FS and cause failure. Joints and faults in rock masses further reduce the stability of the slope. The result of the stability analysis shows that the slope stability of landslides in the study area strongly depends on the saturation conditions and seismic load. A significant landslide in September 2005 was most likely triggered by a coincidence of high groundwater levels and seismic activity. From the stability analysis results, both the LE and FE methods produce deep-seated failure surfaces through the foundation of colluvial deposits. The failure surfaces obtained from the numerical simulations agree with the field observations. The strength reduction method is well-suited for complex geometry and further representative normal stress distributions, subsequent generating consequential FS values. The planar distribution of shear strain contours from simulation work indicates that structurally controlled failures are prominent in most of the slopes. The FEM generates a FS values less than 1 for almost all the slope sections, although those are over 1 for the LEM. Stability analysis using FE methods revealed that all slope sections were found unstable under seismic load and saturation condition, as porewater pressure develops, which significantly contributes to destabilization of the slopes. This portrays that the slopes were unstable under all stated conditions in general and during rainy seasons and seismic events in particular. In general, the FS computed from the stress-strain behavior of the rock and soil using FE method is more realistic and provides even more accurate FS results when a heterogeneous and complex slope geometry is analyzed. The FEM can incorporate the detailed structural information derived from remote sensing substantially improving the stability analysis. The overall stability of the slopes performed under gravity loading having critical SRF less than 1 except the Nib Amba slope section.

Conclusions
In this paper, a study on both LE and FE methods for slope stability is performed to understand the general conditions of the stability of the slopes in the study area. Detailed fieldwork and laboratory analysis was conducted on representative samples to determine the index and geomechanical properties of the rocks and soils. The vulnerable slopes in the study area have been evaluated, and stability grade has been quantified by means of FS and critical SRF. The results of the stability analysis show that the slope stability of landslides in the study area strongly depends on the saturation conditions and seismic load. In general, field investigations indicate that several different failure mechanisms are superimposed on the deep-seated Debre Sina landslide. The laboratory tests reveal that the lowest value of peak strength is from less compacted tuff and prone to sliding. The tuff layers with low peak strength are initiation points for the sliding surfaces. The slope stability analysis using LEM and FEM allows us to determine the influence of pore-water pressure and seismicity events on the slope stability. The FEM is found more applicable for stability assessment because of the complex geometry, heterogeneous material, and the failuredominating faults in the study area. The studied slopes are initially close to failure, and increased pore pressure or seismic load are very likely triggers. While the saturation conditions were assumed constant in the shown simulations, it is evident that a pore pressure increase would further reduce the FS. The importance of faults in hydraulic gradient of the groundwater and their influence on the slope stability in the study area has previously been discussed extensively (Mebrahtu et al. 2020a(Mebrahtu et al. , b, 2021a. Using satellite images faults has been identified and incorporated into the model. The SRF in combination with poro-elasto-plasticity enables the evolution of physics-based stress concentrations leading to deformation and slope failure. The FS obtained using LE and FE methods without seismic load ranged between a minimum of 2.16 and a maximum of 11.31 and between a minimum of 1.51 to a maximum of 2.01, respectively. Whereas, the calculated FS using LE and FE methods with seismic load ranged between a minimum of 0.92 and a maximum of 4.76 and between a minimum of 0.58 to a maximum of 1.03, respectively. In particular, the slope sections after the earthquake event were considered as unstable with a SRF between 0.58 and 1.03.
From the finding, it is concluded that the studied slope stability evaluation methods should be obtained collectively as part of a larger slope stability analysis to determine the resultant FS. Furthermore, this study states that the differences between the FE and LE methods become significant for the analysis of heterogeneous slopes. Especially when faults are involved and act as initiation points of the failure plane, the LEM tends to overestimate the slope stability significantly. Both methods show a significant decrease in the obtained FS values for the study area when the seismic load is considered but LEM always results in higher FS values. For the static assessment without a seismic effect, the interpretation of LEM and FEM might not vary so much, as both methods indicate the slopes as stable with FS values around 2. Therefore, using different tools for slope stability analysis, especially for complex and heterogeneous rock masses, is recommended. This research address researchers/professionals using LEM need to have confidence without even going to FE methods. The presented work showcases similarities and differences between LEM and FE models for slope stability analysis using real but complex scenarios from the study area including seismic load and increasing pore pressure as possible triggers. The limitations but also the benefits in applicable cases of the LEM method have been identified. Therefore, this study can serve as a justification and example for future application of LEM or FE models in respective scenarios.