GIS-multicriteria evaluation using AHP for landslide susceptibility mapping in Oum Er Rbia high basin (Morocco)

High basin of Oum Er Rbia River, which is located in Middle Atlas Mountain region, is prone to landslide problems due to the geological features combined with the climate change and human activities. The present work including inventory mapping was conducted to establish landslide susceptibility map using GIS-based spatial multicriteria approach. Eight landslide-related factors, including land cover, lithology, distance to road, distance to fault, distance to drainage network, elevation, aspect and slope gradient, were selected for the present assessment. Weight for each factor is assigned using Analytic Hierarchy Process (AHP) depending on its influence on the landslide occurrence. The landslide susceptibility map was derived using weighted overlay method and categorized into five susceptible classes namely, very low (VL), low (L), moderate (M), high (H). The results revealed that 30.16% of the study area is at very low risk, 12.66% at low risk, 25.75% of moderate risk, 22.59% of high risk and 9.11% of very high risk area coverage. The very high landslide vulnerability zones are more common within the river valleys on steep side slopes. Most landslides also involve rocks belonging to the Triassic weathered marl and clay-rich formation. Moreover, human activities namely the construction and the expansion of agricultural lands into forests intervene in inducing landslides through altering the slope stability along the river banks. Lastly, effectiveness of these results was checked by computing the area under ROC curve (AUC) that showed a satisfactory result of 76.7%. The landslide susceptibility map of the Oum Er Rbia high basin provides valuable information about present and future landslides, which makes it viable. Such map may be helpful for planners and decision makers for land-use planning and slope management.


Background
Earth has a life cycle of its own, its rocks are born, grow old, die and then re-emerge in the melting heart of our planet. Landslides are part of this natural process. Generally, a landslide is a slow movement but an exceptional natural event (successive torrential rains for example) or anthropic effects accelerate it. As one of the most geological risks in the world, the landslides caused thousands of victims and deaths, hundreds of billion dollars of damages and environmental losses every year (Aleotti and Chowdhury 1999;Gutiérrez et al. 2015). They occur during the gravity displacement of destabilized soils or rock by natural climatic, geomorphological or geological phenomena or by human activities. Landslides occur in fine scree, moraines, or highly fractured and altered rocks, which are particularly sensitive to landslides, such as clays, marls, gypsum or superficial formations of alterites. Numerous natural and/or anthropogenic parameters determine the appearance and development of landslide movements (topography, geology, hydrology, hydrogeology, rapid erosion of the foot of certain slopes, urbanization, etc.). Meteorological phenomena, however, seem to cause the greatest number of events.
Several recent studies have developed several methodologies for assessing susceptibility to landslides showing that damage from natural processes is increased over the last decades (Althuwaynee and Pradhan 2016;Hong et al. 2015, Shahabi et al. 2015. These studies of landslide susceptibility mapping have used deterministic models and probabilistic approaches while taking account of whether future environmental conditions will meet the requirements for a landslide, identified in previous landslides. The qualitative methods are based on individual or group expert opinions (Neaupane and Piantanakulchai 2006). Inventory and historical information have been helped Experts to evaluate landslides, determine the main factors inducing them, and identify sites that have similar geological and geomorphological features. Some qualitative methods become semi-quantitative by incorporating ranking and weighting (Ayalew et al. 2005), as is the case of the Analytic Hierarchy Process (AHP) (Saaty 1980;Barredo et al. 2000;Yalcin, 2008;Kamp et al. 2008) and the weighted linear combination (WLC) (Ayalew et al. 2005;Akgun et al. 2010). The main disadvantage of these methods is the involvement of many subjective judgments and fails to quantify the weight of each factor. Therefore, the results of these approaches are other subjective and rely on knowledge of the experts. Based on inventorying and heuristic analysis, qualitative or semi-quantitative methods define the risk zones in descriptive terms and are often used for small-scale regional studies (Soeters and van Westen 1996;Carrara et al. 1999;Zumpano et al. 2014). The AHP method, suggested by Saaty (1980), has become a popular tool for multi-criteria decision-making. It supports decision-makers to make the best decision, by reducing complex decisions to a series of comparative pairs and synthesizing the results. The AHP disaggregates a complex decision problem into different hierarchical levels. This method allows quantifying opinions and transforming them into a coherent decision model (Saaty 1980). It was widely used by many authors worldwide (Hong et al. 2015;Shahabiet al. 2015;Sangchini et al. 2016;Althuwaynee and Pradhan 2016).
In Morocco, the areas subjected to the ground movements are mostly the Rif (El-Fengour 2016; Prokos et al. 2016;Benzougagh et al. 2017) and to a lesser extent the Middle Atlas, due to the existence of reliefs with much contrasted geological formations (friable clays, flysch units, marl, etc.) and difficult climatic conditions. This study focused on producing landslide susceptibility map of Oum Er Rbia high basin (Morocco) by combining GIS techniques and AHP method. Besides, the result of this approach is compared to the landslide filed inventory to evaluate the accuracy of the final map with suitable candidate landslide sites.

Study area description
The study was conducted in the high-basin of OumEr-Rbia River, in the Khenifra Province, in the southwest of the Middle Atlas (Fig. 1). The studied basin is delimited to the west by the Hercynian Central Massif, to the north by the Causse of Ajdir and south-east by the plain of the High Moulouya. It occupies an area of 3612.21 km 2 , which lies wholly within the mountainous terrain with a diversity of landforms, structural features, closed depressions, ravines, and accumulation forms represented by alluvial terraces. The elevation variation of the area is between 662 and 2400 m.
The Oum Er Rbia high-basin is part of Middle Atlas Mountain and Central Massif. According to the 1:500,000 scale geological maps of Rabat (1976), the geological formations of the region is composed dominantly of Cretaceous subtabular limestone formations, Liasic dolomitic limestones, Triassic doleritic basalts and red clays, as well as Paleozoic flyshs and quartzite. From large surface pedogenic alteration of these different formations in the past, many soil types were formed such as alluvial and colluvial soils (more frequent), vertisols, calcimagnesic soils and very deep isohumic soils concentrated in valleys and flat areas, and fersiallitic soils on ancient terraces or forests.
The presence of deep valleys associated with steep slopes, exposure to high rainfall and the clayey and marly nature of the outcropping rocks constitute important factors monitoring the subsurface soil dynamic in the area such as landslides and erosion. Also, the study area is characterized by a Mediterranean climate known for warm/hot and dry summers and mild/cool and wet winters. Within the wet periods experienced from October to April, minimum temperature values are usually 5°C, while in the dry season experienced from May to September, an average maximum temperature rises to 50°C. The average annual rainfall that is around 666 mm, takes place between November and March (the winter). Thus, bare soils without vegetation cover are exposed to high erosion rates leading to the vulnerability of taluses in the study area to landslides and soil erosion, in particular for the time of year showing torrential rainfall with high intensity (El Bouqdaoui 2007.
The Oum Er Rbia River and its tributaries of Srou and Chbouka, crossing Atlas mountain chain, drain the studied watershed. These rivers are usually with steep courses and gravel beds, and present irregular regime due to the Mediterranean climate prevailing in the region characterized by prominent seasonality. Exceptional rainfall events combined with the presence of steep slopes and unconsolidated rocks (clay and marls) represent the main factors affecting landslide, soil erosion and leaching phenomena that are negative environmental impacts on soil ) and slope stability and river water quality (Barakat et al. 2016(Barakat et al. , 2018.

Materials and methods
As the susceptibility to landslides starts by the step providing the landslide inventory (Pradhan et al. 2010), this study began with the preparation of an inventory map of Oum Er Rbia high-basin from the field investigation along the road network and the interpretation of Google Earth images. Eight landslide affecting factors namely slope, drainage, lithology, land use, slope aspect, roads, faults, and elevation were used for landslide analysis in the present study. These factors selected either intervene in the stability of slopes and rocky massifs or are exposed to a landslide hazard risk. The various thematic layers relative to these factors were generated and then were combined using weights of factors and sub-factors determined by AHP method to generate the landslide susceptibility map. The combination of all thematic layers in agreement with the AHP results was carried out in a GIS environment using the Weighted Linear Combination (WLC) method. After the landslide susceptibility identification in the study area, the accuracy of landslide susceptibility was evaluated by comparing between the landslide inventory map and the landslide susceptibility map via AUC plots.

Preparing landslide factor layers
The main data required for landslide susceptibility and risk assessment in this study were collected from various sources (Table 1), and constructed of a spatial database.
For each of these thematic maps, the incidence of landslides in their classes was evaluated.
After the landslide inventory, according to our field observations, the thematic layers of the selected factors governing landslides, including slope value, slope aspect, elevation, lithology, land use/cover, distance to roads, distance to drainage network, distance to faults were developed. All produced layers were then combined using weights of factors and sub-factors determined by AHP method to generate the landslide susceptibility map. All thematic layers were integrated into GIS environment by the combination of WLC method and AHP results. The thematic layers of distance to drainage network, to faults and to roads were calculated by Euclidean distance tool in spatial analyst tools of ArcGIS 10.2.2. Landsat 8 Oli image with 30 × 30 m spatial resolution acquired on 15/01/2016 was also used to generate the land use/cover map after radiometric and atmospheric corrections. The land use/cover map was produced by supervised classification likelihood of satellite data using ENVI 5.0 software.
Topographic related factors such as slope degree, slope aspect and elevation were derived from a 30-m Digital Elevation Model of the study area. The digital layer of drainage network was also produced from 30-m DEM by hydrology tools in ArcGIS 10.2.2 software. All thematic layers were converted to raster format having a 30 m × 30 m cell resolution, and each raster was classified into several classes to calculate the numbers of landslide and non-landslide pixels. The preparation procedure for each thematic layer is summarized below.

Topographic factors
In the present study, topographic factor data were extracted from the 30-m DEM. Slope aspect referring the direction of the slope face (Poudyal et al. 2010;Pourghasemi et al. 2012) is frequently used as a landslide-conditioning factor. This factor was reclassified into nine directional classes (Fig. 2a). Slope angle considered a main causal factor, is frequently used to map the susceptibility in landslides (Wei Chen et al. 2017;Nourani et al. 2014). In the current study, the slope angle map was reclassified into seven classes with an interval of 5° (Fig.  2b). Elevation thematic map was extracted and generated from 30-m DEM using ArcGIS software. Ranged between 556 and 2400 m, it was reclassified into five classes (Fig. 2c).

Land use
Land use change is influenced by factors relating to population needs, such as converting the agricultural and forest land to the urban areas, conversion of forest to farmland, and reduction of the involuntary or unethical slope for infrastructure developments. Land use is a major factor in mapping landslide susceptibility. In this study land use was generated using Landsat 8 Oli image, by applying a supervised classification Likelihood with ArcGIS 10.2.2 and Envi 5.0 software, and reclassified to six classes, namely forest, agriculture, uncultivated land, bare land, urban, and water body (Fig. 2d).

Distance to drainage
Rivers play a major role in landslide development (Park et al. 2013). They can lead the failure of banks because of the sub-quotation of slopes, and the modification of the ground caused by gully erosion may also influence landslide initiation (Dai and Lee 2002;Bui et al. 2011). In the present study, the drainage network was produced from 30-m DEM by hydrology tools in ArcGIS 10.2.2. Four different buffers were generated using Euclidean distance method to determine the degree to which the streams could affect the bank slopes (Fig. 3a).

Distance to roads
Distance to roads is one of the major anthropogenic factors influencing landslide occurrences (Nourani et al. 2014;Yilmaz 2010). In fact, during the field works, some landslides owing to road construction work are detected. In the current study, four different buffer zones were generated with an interval of 250 m (Fig. 3b).

Distance to faults
Geological fault areas are in general, highly susceptible to landslides because the surrounding rock strength decreases due to tectonic breaks (Chen et al. 2017). In this study, the fault buffers were reclassified into five categories to produce the distance to faults map at a 1000 m interval using Euclidean distance (Fig. 3c), based on the geological map of Rabat.

Lithology
Lithology is a frequently used factor in landslide susceptibility analyses (Althuwaynee and Pradhan 2016). The lithology map extracted from 1:500000-scale geology map of Rabat showed that the study area is covered with various types of lithological units. These units were classified into seven classes as illustrated in Fig. 3d. Most landslides occur on clayey and marly lands along Srou River and Oum Er Rbia River in the north and the middle of the study area.

Susceptibility mapping
The AHP method was used to find the relative weight and priority of individually factor and sub-factor causing landslides in the high basin of Oum Er Rbia. AHP is an approach to decision-making multi-objective multicriterion, which allows the user to arrive at a scale rather pulled off a set of alternative solutions (Saaty, 1980). It helps decision makers to discover the best suits of their objective and their understanding of the problem. This method is widely used in landslide susceptibility analysis. The process is implemented in some consecutive steps such as: develop the hierarchical structure of the project, perform binary (binary) comparisons of criteria against the objective, establish comparative judgment matrix, calculate priority vectors, give the Random Index value (AI), calculate the average of the value λmax), calculate the coherence index (IC), calculate the Coherence Ratio (RC), establish a table of complete comparisons of criteria, establish a sub-criteria in relation to the number of criteria studied, establish the comparison table sub-criteria, determine the performance of the relative value of each criterion by contribution to the value of the criteria, calculate the aggregation of projects, establish the comparison by peers of alternatives of studied, establish the table of complete comparisons alternatives, determine the performance of the relative value of each alternative by report to project aggregation, calculate final project aggregation, and express final decision (Table 2).
In this hierarchical classification approach, it is also possible to check the coherence of our approach by calculating the consistency or consistency ratio (CR) expressed by Eq.
(1). The latter constitutes an acceptance test of the weights of the various criteria (Saaty 1977). This step aims to detect any inconsistencies in the comparison of the importance of each pair of criteria. The CR consistency ratio is approximately a mathematical indicator of the judgment concerning a decision made randomly; it is calculated using Eq. (1).

CR ¼ CI RI ð1Þ
Where RI is the random consistency index, and CI is the consistency index that is expressed as Eq. (2).
Where λ max is the largest or principal Eigen value of the matrix and is calculated from the matrix and n is the order of the matrix. According to Saaty, the coherence ratio must be ≤10% or an imprecision of less than 10%. The principle consists in comparing the judgment with the random weighting of the elements. Finally, the acquisitive weights were integrated the various causative classes in a single landslide susceptibility index using Eq. (3).
Where R i is the rating classes, each layer and W i is the weights for the each of the landslide conditioning factors.
The resulting LSI-map was classified into five classes (very low, low, moderate, high, and very high) based on natural breaks to define the class intervals in the landslide susceptibility map. Fig. 3 Landslide contributing-factor layers produced for the study area: (a) distance to drainage network, (b) distance to roads, (c) distance to faults, (d) Lithology

Susceptibility mapping
After developing the landslide susceptibility map of the Oum Er Rbia high-basin, it was necessary to verify the accuracy of the landslide susceptibility model used in this study. A proper validation was conducted by comparing between the map obtained from the AHP model and the landslide inventory map. It was realized by means of the Receiver Operating Characteristics (ROC) method. The ROC method is widely used to estimate the validity of a model, which predicts the location of the case (occurrence) of a class by comparing an image of adequacy illustrating the probability that this class occurs, and a Boolean image shows where this class exists.
ROC supplies a diagnosis that can be used to distinguish two classes of events and to display the performance of the classifier ROC curves (Sweets 1988;Gorsevski et al. 2006), typically feature true positive rate on the Y-axis, and false positive rate on the X-axis (Soeters and van Westen 1996;Williams et al. 1999). It means that the top-left corner of the plot is the "ideal" point -a positive forgery of zero and has a positive real rate of one. It is not very realistic, where the zone under the curve (AUC) characterizes the quality of a system of forecast by describing the capacity of the system to be correctly anticipated, the occurrence or non-occurrence of pre-defined "events" (Nie et al. 2001). The ideal model displays a curve that has the largest AUC; the AUC varying from 0.5 (random prediction, represented by the diagonal straight line) to 1.0 (Fawcett 2006;Nandi and Shakoor 2010).

Results and discussion
In this study, a GIS-based AHP as a multicriteria evaluation approach was used to identify the potential landslide occurrences in the Oum Er Rbia high basin. Eight landslide conditioning factors i.e. slope degree, aspect, elevation, lithology, land use, distance to drainage network, distance to roads and distance to faults, were employed for susceptibility analysis.
The AHP model is conventionally based on a rating system provided by expert opinion. In fact, expert opinion is very useful in resolving complex problems like landslides. However, to some extent, opinions may change for every individual expert, and thus may be subjected to cognitive limitations with uncertainty and subjectivity. Therefore, it is important to analyze the spatial relationship between the landslide conditioning factors and landslide locations. In this study, spatial analysis of each parameter and field observations were considered for expert judgment. As shown in Tables 3 and 4, pairwise comparison matrix for the factors and sub-factors and their relative weights are processed based on Saaty (2001) methodology. CR was calculated for all the factors and it was less than 0.10, which means that weights attributed were suitable and reliable. As a result, landslide susceptibility map was produced in GIS. The LSI 3 Moderately more important One decision element is moderately more influential than the other.

5
Much more important One decision element has more influence than the other.

7
Very much more important One decision element has significantly more influence over the other. 9 Extremely more important The difference between influences of the two decision elements is extremely significant.
2, 4, 6, 8 Intermediate judgment values Judgment values between equally, moderately, much, very much and extremely. value of the study area ranging between1.88 and 17 was reclassified into five landslide susceptibility classes: very low, low, moderate, high, and very high-using the natural break classifier (Fig. 4). Based on the results of analyses as shown in Table 5, very low, low and moderate susceptible occurrences represent 30.16%, 12.66%, and 25.75% of the total study area, respectively. The high and very high susceptibility areas represent, respectively, 22.59% and 9.11% of the total study area (Table 5).
According to the landslide susceptibility map (Fig. 4), areas with very high landslide vulnerability are in the northern and eastern parts of the watershed along valleys close the drainage network characterized by steep slopes promoting erosion and slides. The western area has very low vulnerability areas having flatter terrain, dense forest cover and sparse forest cover. As well noted in the landslide susceptibility map, landslides mainly affect the clayey formations locating on steep riverbed slopes. The outcropping formations characterized by the higher percentage of clay or marls are deforested and plowed by local farmers. The work of these lands makes them, in addition to their leaching and erosion, less consolidated and more permeable, and therefore vulnerable to landslides. The carbonate and sandstone rocks, which have high mechanical resistance, presented valueless landslide densities of landslides. However, they are moderately involved in landslides when they form an inclined substratum of the clayey and marly formations. Ground truth verification, the location points of landslides were collected using the Global Position System (GPS) device during field visits (Fig. 4).
For instance, distance to faults, to drainage network, slope and lithology are the most important causative factors followed by slope aspect, elevation, while causative factors like distance to roads and land use are less important. However, sometimes some of these less important preconditioning factors could have a triggering effect but under specific conditions (He and Beighley, 2008). For example, new road excavations or further construction on the land susceptible to landslides may activate landslide occurrence, as well as the presence of temporary water  rainfall. Moreover, this is the case for most of the landslides in our study area that are recorded mostly near roads. Besides, the village of El Kbab located on the south side of the Srou River having a steep slope is subject to landslides due to lack of land use planning, construction of houses and deforestation. To validate the precision rate of landslide susceptibility model (map) in the current study, using the AUC method, the total landslides observed in the study area were divided into two groups: 70 and 30% of 50 landslide locations were used for training and validation of models, respectively. Success rate and prediction rate curves were created on the basis of training data and validating data, respectively. Therefore, fifty landslide location points were collected using the Global Position System GPS device during field visits and compared with five levels of susceptibility map. The ROC curve was produced by plotting cumulative percent of LSI in descending order against cumulative percent of landslides on X and Y axis, respectively (Fig. 5). Area under curve (AUC) value of accuracy curve was calculated by simple trapezium method and its value was 0.767. The analysis revealed that the global success rate of the landslide susceptibility zonation map is 76.7%. Certainly, numerous approaches for the mapping of the susceptibility in landslides developed for the last two decades, i.e. that is a mapping based on the inventory, the statistical analysis, the heuristic, semi-quantitative, quantitative, probabilistic and multicriteria decision.
However, none of these approaches is universally accepted for effective analyses of susceptibility in the landslide, because their accuracies are still discussed. In this study, a GIS-multicriteria decision-making process has  been applied of its utmost importance to determine the probable landslide occurrences. This approach involves consideration of several landslide explanatory variables whose identification of the contribution of each of them constitutes a challenge. Therefore, AHP was solicited to derive priority scales of different landslide causative factors and sub-factors, through pairwise comparisons based on the expert judgments. From the spatial effectiveness of the generated landslide susceptibility map checked by AUC (76.7% of accuracy), it is seen that the used model yielded a good result for landslide susceptibility mapping in the study area. Despite these encouraging results and the flexibility of the model, the main issue is that of causality related to landslides in our study area. This issue has been discussed by some scientists in such worldwide cases (Chacon et al. 2006;Van Westen et al. 2006). The lack of high-resolution images and DEMs, of large-scale geological maps, of detailed forest maps, of soil maps, and of landslide data in the study area, constituted the major constraints in selecting the parameter data related to landslides and in time. For these reasons, especially the lack of publicly available landslide information in the Oum Er Rbia high basin, the AHP method was selected because it is still useful, especially for large-scale assessments or for areas with no available landslide inventory (Zhou et al. 2016).

Conclusions
The landslide is a vital natural hazard, and therefore, the recognition of areas at risk of landslides and the mapping of the susceptibility to landslides are the interest of responsible organizations and researchers. Landslide susceptibility analysis can be done under the circumstance of having few existing data about the factors causing landslides using AHP method, which allows fast and practical analysis of landslides based on the collection of data and manipulation and the analysis of the necessary environmental data for landslide susceptibility.
The Oum Er Rbia high basin is prone to landslides due to their geological and geomorphological setting. In this study, the spatial relationship between field landslide occurrences and causative factors, including elevation, slope degree, slope aspect, lithology, land use/cover, distance to drainage network, roads, and faults were assessed using AHP and GIS techniques. The performance of the model was validated using 30% of the data of landslides mapped by field investigations using AUC plots. The landslide susceptibility map was classified according to the natural break method into five classes with an area of 30.16%, 12.66%, 25.75%, 22.59% and 9.11% of the total study area, for very low, low, moderate, high and very high classes, respectively. The very high and high susceptibility zones are shown along the steep banks of the main rivers (Srou and Oum Er Rbia) and their tributaries and along the main faults across the study area. They are also controlled by clayey and marly rocks, which have the highest landslide density. Moreover, human activities namely the road and house construction and the expansion of agricultural lands into forested lands intervene in inducing landslides through altering the slope stability along the river banks. The AUC values of 76.7% supported the good accuracy of the LSI-AHP model in landslide susceptibility assessment in the study area, provided thus that field conditions and characteristics were correctly determined by proper expertise. Moreover, the landslide susceptibility map of the Oum Er Rbia high basin provides more information about present and future landslides, which makes it viable. Such map may be helpful for planners and decision makers for land-use planning and slope management in the study area to provide prevention of landslide risks and to take preventive and suitable security measures.