GIS-based Landslide Susceptibility Evaluation Using Analytical Hierarchy Process (AHP) Approach: The Case of Tarmaber District, Ethiopia

Landslides are one of the natural threats that often result in great loss of life and destruction of property in Ethiopia. One of the areas that is affected by landslides of different types and sizes is the Tarmaber district in the rift margin in the central part of Ethiopia. Keeping in view the cause and effect relationship and mitigation, Analytical Hierarchy Process (AHP) approach is used in the present case to understand the possible causes for landslides. Based on AHP, landslide susceptibility map is produced for Tarmaber using field survey data, remote sensing data, and geographic information system tools. The factors considered in the present case that can influence landslides are lithology, proximity to fault, land use, proximity to drainage, slope gradient, aspect and elevation. The results are validated with the inventory of landslide occurrences. The landslide susceptibility index (LSI) is calculated using the weighted -linear combination (WLC) technique based on the assigned weight and rating given by the AHP method. The accuracy of the results verified using the existing landslide locations is about 88.6%.The landslide susceptibility zonation map has identified our classes/zones: very high (29%), high (44%), moderate (20.0%) and low (7%). Keywords:  Landslide susceptibility; Analytical Hierarchy Process (AHP); Tarmaber; Ethiopia.


INTRODUCTION
Earth being dynamic, the surface of the earth is always under change due to various geoprocesses causing natural hazards and threat to mankind. Landslides are one of the destructive natural hazards that commonly lead to serious problems, especially in hilly and mountainous Tarmaber district, one of the most landslide-affected districts in southern Afar Rift Margin of Ethiopia is chosen study (Fig 1) with the aim of (1) mapping the landslide occurrences, (2) identifying the causative factors, (3) correlating and evaluating the contribution of the factors to the slope failures, and (4) developing landslide susceptibility maps. The study area covers about 220km 2 , with an average rainfall of 1947 mm and an elevation between 1370m © CNCS, Mekelle University 16 ISSN: 2220-184X and 3100m above mean sea level. Less elevated areas are densely populated and intensely cultivated while the elevated have steep slopes with good vegetation cover. The area is prone to high seismic activity and the flat terraces and cliffs are attributed to the rift margin faults. The drainage system of the area is dense and most of the rivers are affected by gully erosion resulting in deeply -cut rivers and associated mass movements. Large-scale and complex landslide occurred on 13 th September 2005 in Tarmaber district has resulted in the loss of >1500 hectares of arable lands, displacement of 4049 people, destruction of >1200 dwelling houses (Tukul settlements) (Fig 2), and >75% crop harvest failure in many localities (Izaba and Shotel Amba) (Gebresilassie, 2007;Woldearegay, 2008;Abay and Barbieri, 2012;Abay, 2012). New or reactivated different intensity landslides are also a common phenomenon in this district in every rainy season or during seismic tremor events or combination of both (Woldearegay, 2008;Abay and Barbieri, 2012). Analysis (AHP) with help of GIS is used for correlating landslide locations and the chosen causative factors in the study area.

Methodology
There are various methods to study and evaluate the landslide phenomena and its causative and triggering factors. The elements that affect slope stability are numerous and varied, and interact in complex ways (Varnes, 1984). In the present case, seven causative factors are chosen as inputs for the landslide hazard evaluation which include: lithology, proximity to fault, land use, slope steepness, aspect, elevation and proximity to drainage. Data were extracted and collected using remote sensing and field based surveys to prepare the preliminary lithological, structural, land use, landslide inventory maps from the interpretation of Aster image of various years, Landsat images of 2001, time series of historical Google Earth(from 2006 to 2017). This was supported by the existing regional geological and topographic maps. The images have been interpreted using ENVI4.5 by applying different enhancing and band composition techniques to visualize various features and were digitized using ArcGIS10.1. The topographic factors like elevation, slope steepness and slope aspect have been prepared from 30m DEM (Digital Elevation Model).
Then, detailed field surveys have been conducted and the preliminary input maps have been crosschecked and updated into more detailed maps (Fig 4).The procedures followed to assign weights to various causative factors and produce landslide susceptibility maps is given below.

The AHP Model
The AHP Model as developed by Saaty (1980), is a decision-aiding tool for dealing with complex and multi-criteria decisions. It is one of the very popular multi-criteria decision making methods with a wide application in many fields such as site selection and suitability analysis (Bantayan and Bishop, 1998;Mahsa et al., 2011;Anagnostopoulos and Vavatsikos, 2012), regional planning (Jankowski, 1989), urban landuse planning (Dai et al., 2001;Feng and Chan, 2004), and environmental impact assessment (Ramanathan, 2001;Gregory et al., 2005), and design and engineering (Hambali et al., 2009;András, 2010). One of its wide applicability in recent years is in the field of landslide study. Several landslide studies have been published using the AHP approach (Yagi, 2003;Ayalew et al., 2004;Woldearegay, 2005;Bachri and Shresta, 2010;Narumon and Songkpt, 2010;Long and De Smedt, 2011;Mezughi et al., 2012;Moradi et al., 2012;Khodadad and Jang, 2015). AHP has an advantage of permitting a hierarchical structure of the criteria which provides users with a better focus on the specific criteria (factor) and sub-criteria (classes) when allocating the weight (Ishizaka and Labib, 2009 decision makers typically use their expert knowledge (judgments) about the elements' relative meaning and importance. Then, the method converts these evaluations into numerical values that can be processed and compared over the entire range of the problem. A numerical weight or priority is derived for each element of the hierarchy, allowing diverse and often numerous elements to be compared to one another in a rational and consistent way. The AHP helps to overcome the problems with arbitrary weights and scores approaches, by its ability to enable decision-makers to derive ratio scale priorities or weights as opposed to arbitrarily assign them (Yalcin, 2007).
To make a decision in an organized way and generate priorities, we need to decompose the decision into steps. Application of AHP to a decision problem considers four steps (Saaty, 1980;Zahedi, 1986;Saaty, 2008): (1) Structuring of the decision problem in a hierarchical approach, (2) Making pair-wise comparisons and obtaining the judgmental matrix, (3) Computing local weights and consistency of comparisons, and (4) Weights aggregations.
Step-1: Structuring of the decision problem in a hierarchical approach This step allows a complex decision to be structured into a hierarchy descending from an overall objective to various 'criteria', 'sub-criteria', and so on until the lowest level. The overall goal of the decision is represented at the top level of the hierarchy while the criteria and sub-criteria contributing to the decision are represented at the intermediate levels. Finally, the decision alternatives are positioned at the last level of the hierarchy. Although, there is no clear set of procedures for generating the levels to be included in the hierarch, Saaty (2000) indicated that a hierarchy can be constructed by creative thinking, recollection and using people's perspectives.
In this study the hierarchy of the landslide causative factors (criteria or sub-criteria) is adopted from the correlations studied by applying frequency ration model (Abay and Barbieri, 2012) besides to the field observations.
Step-2: Constructing pair-wise comparisons and obtaining the judgmental matrix Once the hierarchy has been structured, the next step is to construct the pair-wise comparison matrix as proposed by Saaty (1980). The input data for problem consists of matrices of pair-wise comparisons of elements of one level that contribute to achieve the objectives of the next higher level. That is, the elements of a particular level are compared with respect to a specific element in the immediate upper level. Once the matrix is created, elements are compared pair-wise to determine their relative importance in terms of each criterion (factors) based on the scale introduced by Saaty (1980). According to this scale, the available values for the pair-wise comparisons are members of the set: {1/9, 1/8, 1/7, 1/6, 1/5, 1/4, 1/3,1/2, 1,2,3,4,5,6,7,8,9 } (Table 1). Using this scale, the verbal judgments for each pair wise elements is transformed into numerical quantities. Usually, an element receiving higher rating is considered as superior (or more influential) compared to another one that receives a lower rating.  (Saaty, 2008). Extreme importance The evidence favoring one activity over another is of the highest possible order of affirmation Reciprocals of above nonzero

Intensity of Importance
If activity i has one of the above nonzero numbers assigned to it when compared with activity j, then j has the reciprocal value when compared with i.
A reasonable assumption 1.1-1.9 If the activities are very Close May be difficult to assign the best value but when compared with other contrasting activities the size of the small numbers would not be too noticeable, yet they can still indicate the relative importance of the activities.

Step-3: Computing local weights and consistency of comparisons
The aim of this step is to find a set of priorities or local weights, which is the normalized Eigen vector of the elements of the matrix. Although different methods have been used to derive the priorities, we employed the Eigen-value technique for computing the weights under AHP, which is one of the common methods developed by Saaty (1980). The steps followed to compute the criterion weights (Eigen vector) of a reciprocal matrix involved the following operations: (i) Summing of the values in each column of the reciprocal matrix; (ii) Divide each element in the matrix by its column total (the resulting matrix is referred to as the normalized pair wise comparison matrix and the sum of each column is 1); and (iii) Compute the average of the elements in each row of the normalized matrix, that is, divide the sum of normalized scores for each row by the number of criteria. These averages provide an estimate of the relative weights of the criteria being compared. The priority vector shows relative weights among the things that we compare. The higher the weight is the more important the criteria. Priorities or relative weights make sense only if derived from a consistent or near consistent matrices. Thus, checking the consistency of reciprocal matrix is vital. To do that, we need what is called Principal Eigen value ( ) which is an important validating parameter in AHP. It is used as a reference index to screen information by calculating the Consistency Ratio (CR) (Saaty, 2000) of the estimated vector in order to validate whether the reciprocal matrix provides a completely consistent evaluation. The CR is calculated as per the following steps: i) Calculate the eigenvector and largest Eigen value ( ) for each matrix of order n. Principal Eigen value ( ) is obtained from the summation of products between each element of Eigen vector and the sum of columns of the reciprocal matrix.
ii) Compute the consistency index (CI) for each matrix of order 'n' by the formulae: iii) Calculate the CR using the formulae: Where, RI is random consistency index obtained from a large number of simulation runs and varies depending upon the order of matrix (Table 2).  (Saaty, 1980(Saaty, , 2000 CR range varies according to the size of matrix i.e. 0.05 for a 3 by 3 matrix, 0.08 for a 4 by 4 matrix and 0.1 for all larger matrices, n ≥ 5 (Saaty, 2000;Cheng and Li, 2001). The value of CR ≤ 10%, indicates a good level of consistency in the comparative judgments represented in that matrix and is acceptable. CR value with>10% has resulted in inconsistency of judgments within that matrix. This suggests that the evaluation process needs to be reviewed and improved.
Step-4: Weights Aggregations Final priorities of the alternatives can be obtained by aggregating the local priorities of elements of different levels, which are obtained in the above steps (steps 1-3). The AHP approach adopts an additive aggregation (equation 3) with normalization of the sum of the local priorities to unity (Ishizaka and Labib, 2009) Where, Pi: over all final priority of the alternative i (LSI in our case); Lij = local priority; wj; weight of the criterion j.
The specific application of the AHP method is presented below in modeling the landslide susceptibility of the study area.

Input Data Preparation
The factors that cause landslides are varied, and interact in complex ways (Varnes, 1984). There is no standard method for selecting these factors (Yalcin, 2008). In the present case, seven factors are considered for landslide hazard evaluation based on the landslide occurrences and their characteristics in Tarmaber, and published information (Abay and Barbieri, 2012; Meten et al., 2015). These include lithology, proximity to fault, land use, slope steepness, aspect, elevation and proximity to drainage.

Landslide Inventory Mapping
Landslide inventory maps depict spatial distribution of past landslides and it is an essential information for producing landslide zoning (susceptibility, risk and hazard zoning) (Moradi et al., 2012). The landslide distribution is determined using image interpretation (Aster image of various years, Landsat images, and Google Earth) and field survey, and then digitized directly into inventory map using GIS. Accordingly, 165 landslides of different types and sizes are identified in the study area (Table 3, Fig 3). The inventory map is also used an input map for verification of the susceptibility of landslide occurrences prediction for the study area.

Landslide Causative Factors i). Slope Gradient
Slope gradient is very regularly used in landslide susceptibility studies since land sliding is directly related to the slope angle (Dai et al., 2001;Lee, 2005;Woldearegay, 2005;Yalcin, 2007;Long et al., 2011). For the study area, the slope is derived from 30m DEM using the slope function of the spatial analyst of ArcGIS 10.1. The slope map is in the form of a raster having the same pixel size as the DEM. A map of slope classes is generated by grouping the slope angles into six different classes (Fig 4a): (1)

ii). Aspect
Aspect is considered as a landslide controlling factor by several other studies (Saha et al., 2005;Yalcin, 2011). Aspect in general refers to the orientation to which a mountain slope faces. The aspect of a slope can make very significant influences on its local climatic factors such as amount of rainfall which in turn influences the occurrences of landslides. Aspect related parameters such as exposure to sunlight, drying winds, rainfall (degree of saturation), and discontinuities may control the occurrence of landslides (Dai and Lee, 2002). The Aspect map of the study area was derived from the DEM of the study area. During mapping, the slope aspect was grouped into eight main directions (Fig 4b). removal of materials that provided support at the toe. For this reason the drainage proximity is considered as one causative factor in the landslide susceptibility study. The study area is divided into six different drainage proximity zones (Fig 4c) and about 52.7% of the study area is found within 300m distance from the drainage and about 70% of the landslides occur in these areas.

iv). Elevation
Elevation is another factor that plays an important role in landslide susceptibility assessments considered in several researchers (Dai et al., 2002;Yalcin et al., 2011;Dou et al., 2015). The variation in elevation may be related to different environmental settings such as vegetation types and rainfall (Catani et al., 2013). The elevation map of the area of interest has been prepared from the 30m DEM using the spatial analyst of ArcGIS10.1 and categorized into 5-class ranges (Fig 4d). The minimum elevation is 1,368 m, while the maximum is 3,200 m. Elevation classes in the range of 1,500-2,000m and 2,000-2,500m covers 44.5% and 32.5% of the total area respectively. No major landslides were recorded at elevation class with greater than 2500m because this portion of the study area is covered by dense and deep rooted trees and sound rock mass which eventually reduced the potential land sliding process as it is covered by deep rooted forest and sound rock exposures. About 86% of the landslides occurred at elevation range between 1,500 and 2,500 m.

v). Proximity to Faults
Faults have been considered as a critical factor in activating landslide in tectonically active areas (Tien et al., 2011). The main Ethiopian rift system and its escarpments (where the study area is located) are characterized by normal faults. Fault map for the study area is prepared and compiled from various images (Aster and Landsat images) and previous works. The images have been interpreted using ENVI4.5 by applying different enhancing and band composition techniques, identifying the various features and digitizing them using ArcGIS 10.1. Faults are very important factors for landslide initiations due to the fact that they are not only weak zones, but also mostly characterized with: (a) deeper weathering, resulting in greater thickness of soil masses; (b) higher potential for concentrated groundwater flow, which can act as lubricant and also can produce water pressures causing landslides. Keeping this in view, fault proximity maps are produced to evaluate its contribution to the landslide (Fig 4e).The major fault trends in the study area can be grouped into the NNE-SSW, NNW-SSE and East-west trending fault systems.

vi). Land Use
The land use pattern is often has a great influence in the landslide occurrences because they relate to the anthropogenic interference on hill slopes (Pradhan et al., 2010). Accordingly, Aster and Landsat image, Google Earth, topographic maps of 1:50,000 and field surveys are used in preparing the land use map (Fig 4f). The major land use pattern in the area comprises influence. People in this region are still actively involved in agriculture and are moving into steep slopes and cultivating without constructing proper terraces. This is leading to deforestation and land degradation which contributes to landslide initiation/reactivation in the area.

vii). Lithology
Lithology is one of the most important parameters in landslide studies because different lithological units have different susceptibility (Dai et al., 2001;Yalcin, 2007

RESULTS AND DISCUSSION
Seven landslide causative factors mentioned above, depending on their relative influence, each factor is further classified into a number of classes. In putting priorities, weighing factors, and determining relative influence of the factors-(a) field-based expert judgment, and (b) review of published data (Abay and Barbieri, 2012) are used. Values ranging from 9 (extremely) to 1 (equally) and 1/9 (opposite extremely) are assigned based on table 1 to each pair of parameters resulting a square reciprocal matrix by rating rows relative to columns, shown in tables 4 and 5.
Pair-wise comparisons and obtaining the judgmental matrix, consistency checks are done.
Although all the mentioned factors induces landslide, their individual relative influences on the slope instabilities are different. To weigh the relative importance of the above mentioned causative factors and their subdivisions (classes) quantitatively on the initiation of the landslide, a pair-wise comparison and a judgment matrix was made based on the proposal of Saaty (1980Saaty ( , 2000. When comparing two attributes (layer classes or parameters in a layer), the above stated numerical relational scale is used (Table 1).
Once the comparisons matrices are made, the priorities or relative weights or Eigen vectors, as well as the Principal Eigen value ( ) are calculated following the procedures stated in step-3 above. Then, consistency index (CI) and consistency ratio (CR) values of each developed factor or class matrices have been determined using equations 1 and 2 above while values of Random consistency index (RI) are referred from table 2. The average eigen vectors (relative weight) for each factor, in the columns are initially calculated following the steps stated under step-2 above.
As can be noted from table 4 and figure 6, lithology is the major parameter contributing to the landslide of the Tarmaber district followed by the faulting and drainage, which is also true from the prospective of the field observation. The area is characterized by the fragile types of lithology i.e. thick colluvium-eluvium, highly fractured and deeply weathered rhyolitic ignimbrites, welded tuffs, and basalts. Slope aspect and slope gradient are the next influential parameters, with more or less equal influences on the landslide occurrences.   Table 5. Pair-wise comparison matrixes, principal eigenvectors (relative weights) and consistency ratios of classes within the various parameters (causative factors) and the data layers. Finally, the aggregation and integration of the various weights of the factors and classes to a single landslide susceptibility index (LSI) is calculated on the basis of weighted -linear sum (WLS) (Voogd, 1983)    From the calculation, it was found that the LSI had a minimum value of 0.492 (low susceptibility), and a maximum value of 2.294 (very high susceptibility). Finally all the identified values of LSI were re-classified and categorized it to four classes by adopting manual classifier method as suggested by Galang (2004) to signify five classes of the Landslide susceptibility zone of the area. Susceptibility evaluation (using the AHP method) shows that 29% and 44% of the study area is categorized into very high and high levels of susceptibility respectively (Figs 7 & 8).
Verification of the result was executed by comparing known landslide location data (i.e. the landslide inventory map) with the landslide susceptibility map. The rate curve was created and its areas of the under curve was calculated for all cases.
The rate explains how well the model and factor predict the landslide. Thus, the area under curve can measure the prediction accuracy qualitatively. To create the validation curve, the computed index values of all cells in the study area were arranged in descending order and divided into 100equal classes ranging from very highly susceptible classes to low susceptible classes. Then, 100 classes were overlaid and intersected with known landslides to establish the percentage of landslide incidences in each susceptible class. The rate verification results appear as a line in figure 9. The fitness the rate curve can be judged by the fact that more percentage of landslides must occur in very high susceptibility zone as compared to other zones. For example, Remondo et al. (2003) states that: (1) a hypothetical "validation curve" coinciding with a diagonal from 0 to 1 would be equivalent to totally random prediction; and the further up and away the validation curve from that diagonal the better the predictive value of the model, (2) similarly, the greater the gradient in the first part of the curve the greater its predictive capability.
Landslide susceptibility index (LSI) indicates 29% class of the study area (Fig 8)

CONCLUSION
Landslide susceptibility map is prepared using Analytical Hierarchical Process (AHP) for Tarmaber district located in the rift margin of Ethiopia. Seven landslide causative factors were analyzed and included in the study to prepare the susceptibility map of the area. The relative importance of causative factors/classes on the initiation of landslide and susceptibility of the study area is evaluated based on the relative weights, and consistencies are checked.
© CNCS, Mekelle University 32 ISSN: 2220-184X The areas highly prone to landslide in the Tarmaber district are those areas covered by: The landslide map generated using AHP method has identified four susceptibility zones: very high (29%), high (44%), moderate (20%) and low (7%), and accordingly suggest that 73% of the area is prone to landslides.

ACKNOWLEDGMENTS
Authors duly acknowledge "The Developing Countries Cooperation Board of the Italian Ministry of Foreign Affairs" for funding the research work; "Tekeze Deep Water Wells Drilling P.L.C" for providing logistic support; "The USGS-Northern California Earthquake Data Centre and USGS/NASA" for allowing to use the earthquake data and Aster images of the study area; and "Mekelle University, College of Natural and Computational Sciences" for covering part of the field expenses. Also thank reviewers for their valuable comments on the paper.