Nationwide Susceptibility Mapping of Landslides in Kenya Using the Fuzzy Analytic Hierarchy Process Model

Landslide susceptibility mapping (LSM) is a cost-effective tool for landslide hazard mitigation. To date, no nationwide landslide susceptibility maps have been produced for the entire Kenyan territory. Hence, this work aimed to develop a landslide susceptibility map at the national level in Kenya using the fuzzy analytic hierarchy process method. First, a hierarchical evaluation index system containing 10 landslide contributing factors and their subclasses was established to produce a susceptibility map. Then, the weights of these indexes were determined through pairwise comparisons, in which triangular fuzzy numbers (TFNs) were employed to scale the relative importance based on the opinions of experts. Ultimately, these weights were merged in a hierarchical order to obtain the final landslide susceptibility map. The entire Kenyan territory was divided into five susceptibility levels. Areas with very low susceptibility covered 5.53% of the Kenyan territory, areas with low susceptibility covered 20.58%, areas with the moderate susceptibility covered 29.29%, areas with high susceptibility covered 29.16%, and areas with extremely high susceptibility covered 15.44% of Kenya. The resulting map was validated using an inventory of 425 historical landslides in Kenya. The results indicated that the TFN-AHP model showed a significantly improved performance (AUC = 0.86) compared with the conventional AHP (AUC = 0.72) in LSM for the study area. In total, 31.53% and 29.88% of known landslides occurred within the “extremely high” and “high” susceptibility zones, respectively. Only 8.24% and 1.65% of known landslides fell within the “low” and “very low” susceptibility zones, respectively. The map obtained as a result of this study is beneficial to inform planning and land resource management in Kenya.


Introduction
Every year, landslides cause a large number of deaths and enormous property losses in mountainous areas [1]. Landslides are a prehistoric issue. It is currently receiving considerable attention since damage induced by landslides has risen in recent years. Estimated fatalities from landslides reached 32,322 between 2004 and 2010, though this value is likely underestimated [2]. The situation varies in different countries. Landslides are concentrated in developing countries or regions, such as the Himalayan region and its surrounding areas in China and African countries. Hence, more effort is required to reduce landslide risks within those countries. Within this topic, the preemptive identification of landslide-prone areas through landslide susceptibility mapping (LSM) is a very promising hazard mitigation approach.

Study Area
Kenya is an east African country ( Figure 1). Kenya has a territorial area of 582,646 km 2 and a population of 41.8 million. The elevation of Kenya stretches from sea level in the coastal regions to over 5000 m above sea level (a.s.l.) at Mount Kenya. Geomorphologically, the landforms in Kenya are dominated by highlands in the central and the west regions, plains in the northeast and the coastal regions [20]. The Great Rift Valley (GRV) cuts through the western territory of Kenya from south to north, separating the highlands into two parts, the western highland and the eastern highland. The narrowest part of the rift valley basin is about 30 km near the Naivasha Lake, while the broadest of that is about 300 km in width near the Turkana Lake ( Figure 1). The elevation difference in the GRV zone ranges from 500 to 1000 m.

Study Area
Kenya is an east African country ( Figure 1). Kenya has a territorial area of 582,646 km 2 and a population of 41.8 million. The elevation of Kenya stretches from sea level in the coastal regions to over 5000 m above sea level (a.s.l.) at Mount Kenya. Geomorphologically, the landforms in Kenya are dominated by highlands in the central and the west regions, plains in the northeast and the coastal regions [20]. The Great Rift Valley (GRV) cuts through the western territory of Kenya from south to north, separating the highlands into two parts, the western highland and the eastern highland. The narrowest part of the rift valley basin is about 30 km near the Naivasha Lake, while the broadest of that is about 300 km in width near the Turkana Lake ( Figure 1). The elevation difference in the GRV zone ranges from 500 to 1000 m.  Geologically, Kenya is mainly constituted of five parts [20]: (1) the Archean rocks, of which the major rock types are: shales, mudstones, greywackes, phyllites, and conglomerates; (2) the Proterozoic rocks, of which the major rock types are: rhyolites, basalts, quartzites, and conglomerates; (3) the Paleozoic and Mesozoic sediments, which were dominant by rock types of granites, granodiorites, and leucogranites; (4) the Tertiary/Quaternary volcanic rocks and sediments, of which the major rock types are: sands, marls, clays, conglomerates, and limestones; (5) the Pleistocene to Recent deposits, in which clays, diatomite, shales, and silts are major rock types.
Kenya has a mild climate with annual temperatures ranging from 16 to 26 degrees Celsius. The mean annual precipitation (MAP) in Kenya ranges from <200 mm to 2500 mm. The precipitation of Kenya is characterized by its nonuniform distribution in both time and space. From the historical meteorological data of the Jomo Kenyatta weather station (Figure 2), two distinct rainy seasons can be observed. One is from March to May (the heavy rainy season), and the other is from October to December (the light rainy season). The rest is the dry season. The spatial distribution of MAP is illustrated in Figure 3a. Because of the sudden elevation changes in the GRV zone (from highland to valley then to highland again), there is a sharp transition between wet and dry regions across the GRV zone in southwestern Kenya. Both sides of the GRV in this region had the highest value of MAP, while the northeastern and northern parts of Kenya had the lowest MAP of less than 800 mm. Kenya has experienced a series of geohazards arising from floods, storms, landslides, and debris/mudflows. Most of these geohazards are related to climate extremes [18,20]. Rainfall and human activities (farming, devegetation, construction, etc.) have triggered the majority of the landslides that have occurred in Kenya.
Land 2020, 9, x FOR PEER REVIEW 4 of 22 Geologically, Kenya is mainly constituted of five parts [20]: (1) the Archean rocks, of which the major rock types are: shales, mudstones, greywackes, phyllites, and conglomerates; (2) the Proterozoic rocks, of which the major rock types are: rhyolites, basalts, quartzites, and conglomerates; (3) the Paleozoic and Mesozoic sediments, which were dominant by rock types of granites, granodiorites, and leucogranites; (4) the Tertiary/Quaternary volcanic rocks and sediments, of which the major rock types are: sands, marls, clays, conglomerates, and limestones; (5) the Pleistocene to Recent deposits, in which clays, diatomite, shales, and silts are major rock types.
Kenya has a mild climate with annual temperatures ranging from 16 to 26 degrees Celsius. The mean annual precipitation (MAP) in Kenya ranges from <200 mm to 2500 mm. The precipitation of Kenya is characterized by its nonuniform distribution in both time and space. From the historical meteorological data of the Jomo Kenyatta weather station (Figure 2), two distinct rainy seasons can be observed. One is from March to May (the heavy rainy season), and the other is from October to December (the light rainy season). The rest is the dry season. The spatial distribution of MAP is illustrated in Figure 3a. Because of the sudden elevation changes in the GRV zone (from highland to valley then to highland again), there is a sharp transition between wet and dry regions across the GRV zone in southwestern Kenya. Both sides of the GRV in this region had the highest value of MAP, while the northeastern and northern parts of Kenya had the lowest MAP of less than 800 mm. Kenya has experienced a series of geohazards arising from floods, storms, landslides, and debris/mudflows. Most of these geohazards are related to climate extremes [18,20]. Rainfall and human activities (farming, devegetation, construction, etc.) have triggered the majority of the landslides that have occurred in Kenya.

Historical Landslides in Kenya
A landslide inventory usually portrays the date, location, cause, type, and geometry of landslides. A detailed landslide inventory is a mandatory input for LSM using quantitative methods. Landslide inventories can be produced through field surveys, review of relevant documents, and remote sensing (satellite images, aerial photos, etc.). Since compiling a detailed landslide inventory remains a time-and labor-consuming task there has been no reliable landslide inventory with national coverage available in Kenya until now.
Considering the conditions mentioned above, the landslide inventory utilized in the present study was composed of three subsets from different sources, as illustrated in Figure 1. The first subinventory (LS1) was obtained from the Kenya Open Data of the Ministry of Information, Communications, and Technology of Kenya. LS1 contained 39 historical landslides that occurred in Kenya during the 1999-2013 period. Landslide information regarding the specific longitude/latitude, dates of occurrence, number of people affected, and estimated economic losses were provided in the LS1. The second subinventory (LS2) was derived from the global fatal landslide database (GFLD) (version 2) extracted from Froude and Petley [28]. Similar to LS1, LS2 also provided the location, date of landslide occurrence, induced fatalities, general description, and related reports of landslides. In total, 63 landslides were recorded in LS2. The third subinventory (LS3) was extracted from a landslide inventory of Africa, which was compiled by Broeckx et al. [29]. Within Kenya's territory, 323 landslides were detected in LS3. The majority of landslides in LS3 were mapped through visual interpretation of Google Earth imagery. This was time-consuming work and involved subjective interpretation. In contrast to LS1 and LS2, only the locations of landslides were provided in LS3. For the present study, a total of 425 landslides were stored in the inventory.
As indicated in Figure 2, the landslide occurrences were mostly concentrated in the periods from April to May and from November to December. It should also be noted that the monthly distribution of landslides was consistent with the monthly distribution of precipitation.

Landslide Contributing Factors (LCFs)
From the perspective of geotechnical engineering, landslides are a comprehensive consequence of several contributing factors. Nevertheless, there are no global rules for selecting these factors. In a typical LSM, such LCFs are chosen on the basis of data availability, characteristics and scale of the study area, as well as the expert knowledge or experience. In this study, ten LCFs were utilized in LSM over a nationwide area of Kenya, including four topographic factors (namely, the altitude, aspect, slope, and curvature), two hydrological factors (namely, the topographic wetness index (TWI), stream power index (SPI)), soil texture, precipitation, land use, and landform (as shown in Figure 3). To perform the LSM, all factors were rasterized into 1 × 1 km grids and classified into several classes using ArcGIS software (Version 10.2). In what follows, a brief description of each landslide contributing factor is given.

Mean Annual Precipitation (MAP)
Rainfall is among the dominant inducing factors for landslides not only in Kenya but also in many countries because it increases the soil mass and decreases the soil shear strength. As shown in Figure 2, there are strong correlations between landslide occurrences and precipitation in Kenya. For this study, the MAP data were adopted and categorized into nine levels using a 400 mm interval, as shown in Figure 3a. To initially obtain the MAP factor, monthly total rainfall data from 87 meteorological stations in Kenya were collected and filtered. Only data from stations operated in all types of weather were kept. Then, the filtered monthly rainfall data were summed and averaged annually for each station. Finally, through the inverse distance weighting (IDW) of data for each station, the MAP map of Kenya was derived.

Topographic Factors
Topographic factors are most frequently used in LSM. In this study, several topographic factors were created using a 12.5-m digital elevation model (DEM) obtained from the Alaska Satellite Facility (ASF) [30], which included altitude, slope gradient, slope aspect, and curvature information. Because of variations in temperature, humidity, and vegetation, the degree and type of weathering also varied with altitude. Therefore, altitude has been employed as an LCF in previous studies [8,13,15]. Considering the setting of geomorphology and landforms present in Kenya, the altitude factor was categorized as (1) <50 m, (2) 50-200 m, (3) 200-500 m, (4) 500-1000 m, (5) 1000-2000 m, and (6) >2000 m (Figure 3b). Steepness directly affects slope stability because slopes become more susceptible to landslides as the slope gradient increases. Reviews of LSM studies have suggested that the slope gradient factors is usually categorized using a 5 • interval. Thus, the slope LCF was categorized as follows:  (Figure 3c). Since parameters such as sunlight, precipitation, and vegetation cover vary, slope aspects may also have an effect on landslide occurrences. Consistent with previous studies [15,16,31], the aspect factor was classified into 9 subclasses, as shown in Figure 3d. By controlling the water flow and erosion type in curved terrains, curvature is also a commonly used topographic factor that is associated with landslides. A positive curvature value indicates an upwardly convex terrain, while a negative value indicates an upwardly concave terrain. Terrains with values of zero for the curvature factor were classified as flat ( Figure 3e).

Hydrological Factors
Hydrological factors played a determinant role in affecting slope stability. For this study, the TWI and the SPI were utilized as predictors of landslides. Both the TWI and the SPI are secondary attributes of the DEM and can be calculated using Equations (1) and (2), respectively. The TWI measures the topographic control on groundwater flow and accumulation. Terrains with a higher TWI values are more likely to become saturated during rainfall events. The SPI quantifies the erosion power of flowing water. Gullies are more likely to form at locations with high SPI values. As indicated in Figure 3f,g, the TWI of the study area ranged from 6.80 to 34.72, and the SPI ranged from −2.41 to 25.13. The factors of TWI and SPI were finally divided into five categories through the "natural breaks" function in GIS.
where A is the unit upstream accumulation area, β is the slope gradient, and ln is the natural logarithm.

Environmental Factors
Given the poor data availability in Kenya, three environmental LCFs were considered in this study to produce LS maps. These two factors were soil texture and land use. Soil texture indicates the proportional composition of sand, silt, and clay content in the soil. Because high clay soils usually contain high organic matter content, which is favorable for soil resistance against detachment, soils characterized by high sand or loam content are more susceptible to land sliding than clayey soils [32]. A nationwide soil property GIS database developed by the International Soil Reference and Information Centre (ISRIC) [33] was utilized in this study. Types of soil texture were classified as follows: (1) very clayed, (2) clayed, (3) loamy, and (4) sandy ( Figure 3h). Land use and landform type contribute significantly to slope stability. Specifically, vegetation roots may enhance slope stability by altering the cohesive forces and hydrologic properties. The degradation of forests and vegetation increases the degree of susceptibility of the area to landslides [34]. As for the factor of landform, steep and hill/mountain terrains are prone to landslides compared with flat terrains such as plains, valley floor and foot slope. In addition, from the perspective of geomorphology, landslide itself also plays as a driving role in landform evolution [35,36]. The Kenya National Land Use Dataset (KNLD) [37] was utilized in this study to obtain the LCFs of land use and landform type. As indicated in Figure 3i, the KNLD contains ten land use types ( Figure 3i) and eleven landform types (Figure 3j).

Methodology
A hybrid model of the conventional AHP and fuzzy theory was utilized to conduct the LSM in this study. The AHP has shown good capacity in solving a multicriteria decision-making problem by incorporating expert knowledge into quantitative analysis. To reduce subjectifies involved in conventional AHP analysis, the fuzzy set theory was adopted to handle blurry sets or categories. Hence, the hybrid use of fuzzy sets and conventional AHPs effectively addresses the decision-making issues under multiple criteria. The theoretical background of the conventional AHP and fuzzy theory were briefly introduced in Sections 4.1 and 4.2, respectively. After that, the process of incorporating fuzzy theory into AHP was given in Section 4.3.

The Theoretical Background of the Conventional AHP
The AHP, originally developed by Saaty [38], has shown great potential for handling multicriteria decision-making (MCDM) issues. Implemented in GIS, the AHP has been successfully employed in LSM in many previous studies [13,29,31]. A detailed description of the AHP application steps in LSM was introduced by Van et al. [31] and can be summarized as follows: Step 1: Dividing the decision problem into a hierarchical structure In this step, a complex decision problem was decomposed into a hierarchical structure, including an "objective" level on the top, one or more "criterion" level(s) in the middle, and several decision alternatives at the bottom level. Although there are no universal rules to be followed in constructing such a hierarchy, it was suggested by Saaty [39] that the hierarchy be built based on the decision maker's knowledge and experience with the problem.
Step 2: Constructing the pairwise comparison matrix In this step, a comparison matrix was constructed with each element indicating the pairwise comparison between all the decision elements. By asking the decision maker how important alternative A is compared to alternative B, the pairwise comparison results (relative importance) are usually rated using a linguistic variable, such as "Slightly Important", "Moderately Important", or "Extremely Important" (Table 2). Step 3: Calculate the weights of each decision element and check its consistency For each comparison matrix, the relative weights of each decision element were calculated using the eigenvalue method (or some other methods). Weights could be used only if consistency had been satisfied. Step 4: Hierarchically aggregate weights from all "criterion" levels In this step, the score of each alternative with respect to the final goal was calculated by aggregating the weights of decision elements' weights from all "criterion" levels.
The numerical intensity scale for the relative importance between two decision elements, proposed by Saaty [40], has been broadly used in the AHP. Table 2 shows that in this study, the importance of "Equally Important" to "Extremely Important" was scaled from 1 to 9. The numbers 2, 4, 6, and 8 were used to describe intermediate importance. Inverse importance was scaled using the reciprocals of the numbers from 1 to 1/9. The eigenvalue method was adopted to calculate the weights. In this regard, the consistency index (CI) was calculated as follows: where λ max represents the largest eigenvalue of a matrix. For evaluation of the CI, the term consistency ratio (CR) was introduced. The CR was defined as the ratio of a given CI and that of a randomly generated reciprocal matrix (RI). Consistency is satisfied if

Triangular Fuzzy Number (TFN)
In the practical application of the conventional AHP for LSM, the determination of the exact relative importance of two factors (A and B) is more difficult than to identify one factor as being more important to another. Given this, fuzzy theory was employed to extend the conventional AHP by scaling the experts' decisions as fuzzy numbers. Thus, assigning exact ratio values to pairwise comparison results was avoided. There are many types of fuzzy numbers. For this study, triangular fuzzy numbers (TFNs) were used. Concepts for the TFN-AHP are briefly introduced in the following.
Let M be a TFN on R; then, its member function x ∈ M, µ M (x) : R → [0, 1] can be defined as follows: where a, b, and c represent the left, modal, and right values of M, respectively (see Figure 4).

Triangular Fuzzy Number (TFN)
In the practical application of the conventional AHP for LSM, the determination of the exact relative importance of two factors (A and B) is more difficult than to identify one factor as being more important to another. Given this, fuzzy theory was employed to extend the conventional AHP by scaling the experts' decisions as fuzzy numbers. Thus, assigning exact ratio values to pairwise comparison results was avoided. There are many types of fuzzy numbers. For this study, triangular fuzzy numbers (TFNs) were used. Concepts for the TFN-AHP are briefly introduced in the following.
Let M  be a TFN on R; then, its member function x M Î  , can be defined as follows: where a, b, and c represent the left, modal, and right values of M  , respectively (see Figure 4).  A TFN can be denoted by M = (a, b, c). Let M 1 = (a 1 , b 1 , c 1 ) and M 2 = (a 2 , b 2 , c 2 ) be two TFNs, where a 1 , a 2 > 0, b 1 , b 2 > 0 and c 1 , c 1 > 0. The laws of the operations can be defined as follows: • Summation of two TFNs: • Subtraction of a TFN from another TFN: • Multiplication of two TFNs: • Multiplication of a number and a TFN: • Division of a TFN by another TFN: • Reciprocal of a TFN:

Integration of the AHP and TFN
The integration of fuzzy sets with the AHP has shown great potential not only for use in LSM but also in many other multicriteria decision making processes, such as hospital location selection and tourist risk evaluation. Very reliable results have been obtained in these applications. The following sections will describe the TFN-AHP theory.
In TFN-AHP theory, experts' judgments are scaled using a TFN rather than a definite number. Then, the TFN comparison matrix is defined as follows: where k ij = (a ij , b ij , c ij ) denotes the fuzzy comparison value of criterion i to criterion j and , 1 a ij ) denotes the reciprocal value for i, j = 1, 2, · · · , n and i j. If such judgments are made by consulting more than one expert, the element value is calculated by the average of their decisions as follows: where and h is the number of experts. Consistent with Chang (1996) [41], the required steps to compute the weight vector for a TFN comparison matrix can be expressed using the following procedure: Land 2020, 9, 535 11 of 22 First, calculate the fuzzy synthetic extent with respect to the ith alternative by normalization of the row sums of the fuzzy comparison matrix as follows: Then, calculate weight vectors concerning each decision element under a certain criterion using the degree of possibility of M i ≥ M j , which is defined as follows: This equation can be equivalently expressed as follows: Finally, calculate the normalized vector of weights W = ( w 1 , w 2 , · · · , w n ) of the TFN comparison matrix K as follows: The typical LSM is performed based on raster cells. For the TFN-AHP application in LSM, the criteria refer to a series of LCFs, whereas the alternatives refer to the raster cells within the study area. To perform the LSM using the TFH-AHP, a weighted linear combination (WLC) is conducted to calculate the LSI for each raster pixel as follows: where w i is the weight of ith criterion (LCF) i, s k(x,y) i is the weight of the kth subcriteria (subclass for the LCF) in the ith criterion, and k is determined by the spatial location (x, y) of the raster cell.

Accuracy Validation
A review of the literature has shown that the receiver operating characteristic (ROC) curve is a popular method to evaluate the goodness of fit for classification [1,[4][5][6][7]. The area under this curve (AUC) is adopted to measure the generalization performance of the LSM model. The value of AUC usually ranges between 0.5 to 1.0. A higher AUC value, closer to 1.0, indicates a better performance of the classification model.
The first step in performing the ROC analysis was to construct the validation dataset, which contained both landslide and non-landslide events. For this study, 425 known landslides were used for validation. Additionally, 425 non-landslides were randomly chosen for validation within the study area. Then, by setting different threshold LSI values, the dataset was separated into four groups according to the actual label and precited label. As shown in Table 3, the four groups were true positive (TP), true negative (TN), false positive (FP), and false-negative (FN) events. After that, two indexes were employed; one was the true positive rate (TPR) computed using Equation (19), and the other was the false positive rate (FPR) computed using Equation (20). Eventually, the ROC curve was drawn by plotting the FPR and the TPR on the horizontal and vertical axes, respectively.

Flowchart of Conducting the LSM
In general, the process to conduct the LSM using the proposed TFN-AHP method can be summarized as following 5 steps. Firstly, the 10 LCFs (criteria) that were required to perform the LSM were chosen (as described in Section 3). Then, a two-level hierarchical model was developed with 10 criteria and 41 subcriteria. Next, 11 comparison matrices were established to calculate the criteria weights. After that, using a WLC of weights of all levels, an LSI map was created and reclassified. Finally, the accuracy of the obtained map was validated using ROC curve and the known historical landslides.

Weights of LCFs and Their Subclasses
The weighting vector derivation plays a central role in multicriteria decision making. For the present study, weights were assigned to each LCF. For this purpose, the geotechnical experts were called upon to make a pairwise comparison of the LCFs based on their experiences and knowledge. As illustrated in Table 4, the pairwise comparison matrix for the ten LCFs was constructed by considering expert opinion and similar previous studies [21][22][23]26,27]. From the matrix, the weight vector for the criteria was computed using Equation (14) to Equation (17) and is presented in Table 5. After normalization, the weights for each criterion were derived using Equation (18) and are shown in Table 6. The CR was calculated using Equations (3) and (4). When the CR = 0.086 < 0.1, the judgment was deemed to be consistent.
For the subcriteria (subclasses) under a certain uplevel criterion (an LCF), the weights were derived using the same procedures. The sum of subcriteria weights under each corresponding uplevel criterion should be 1.0. Hence, 10 comparison matrices were created. Additionally, the final weights for the subclasses within each LCF were calculated and are shown in Table 6. Before using these calculated weights, a consistency check was conducted for each comparison matrix. Only if CR < 0.1 was the derived weight accepted. Since all CR values were less than 0.1 (Table 6), a consistency check of the 10 matrices indicated that all the judgments were consistent.
From the TFN-AHP analysis, slope gradient (0.1923), MAP (0.1884), and curvature (0.1651) were considered to be the three most important factors contributing to landslide occurrence, whereas the least important factors were land use (0.0220) and aspect (0.0315). For the factor of slope gradient, the terrain steeper than 30 • was most susceptible to landslides (weight for this subclass is 0.364313), while the category of 5-10 • obtained the lowest weights (0.071825) in determining the landslide occurrences. It also can be seen from results that barren land (0.160743), bush land (0.136464), and grassland (0.121685) were most susceptible to landslides compared with other land use types. In case of curvature, both concave and convex terrain were more prone to landsliding than flat area. Convex terrain (subclass weight is 0.570014) was more favorable for landsliding than concave terrains (0.356956).

Landslide Susceptibility Maps
Using Equation (18), the LSI value for each raster cell within Kenya was calculated. As shown in Figure 5, the resultant LSI map was reclassified into five susceptibility levels using the "natural break" function ArcGIS. In total, 15.44% and 29.16% of the Kenyan territory were mapped as extremely high susceptibility zones. A total of 29.16% of the total area was predicted as a high susceptibility zone. Low and very low susceptibility classes covered 20.58% and 5.53% of the study area, respectively. The remaining 29.29% of the study area was determined to be moderately susceptible to landslides ( Figure 6). The distribution of susceptibility classes differed in each province. As illustrated in Figure 7, the Rift Valley Province and Eastern Province had the highest percentages of EH landslide susceptibility coverage (21% and 19%, respectively), while the Central Province and Nyanza Province had the lowest percentages of EH landslide susceptibility (5% and 6%, respectively).

Landslide Susceptibility Maps
Using Equation (18), the LSI value for each raster cell within Kenya was calculated. As shown in Figure 5, the resultant LSI map was reclassified into five susceptibility levels using the "natural break" function ArcGIS. In total, 15.44% and 29.16% of the Kenyan territory were mapped as extremely high susceptibility zones. A total of 29.16% of the total area was predicted as a high susceptibility zone. Low and very low susceptibility classes covered 20.58% and 5.53% of the study area, respectively. The remaining 29.29% of the study area was determined to be moderately susceptible to landslides ( Figure 6). The distribution of susceptibility classes differed in each province. As illustrated in Figure  7, the Rift Valley Province and Eastern Province had the highest percentages of EH landslide susceptibility coverage (21% and 19%, respectively), while the Central Province and Nyanza Province had the lowest percentages of EH landslide susceptibility (5% and 6%, respectively).

Landslide Susceptibility Maps
Using Equation (18), the LSI value for each raster cell within Kenya was calculated. As shown in Figure 5, the resultant LSI map was reclassified into five susceptibility levels using the "natural break" function ArcGIS. In total, 15.44% and 29.16% of the Kenyan territory were mapped as extremely high susceptibility zones. A total of 29.16% of the total area was predicted as a high susceptibility zone. Low and very low susceptibility classes covered 20.58% and 5.53% of the study area, respectively. The remaining 29.29% of the study area was determined to be moderately susceptible to landslides ( Figure 6). The distribution of susceptibility classes differed in each province. As illustrated in Figure  7, the Rift Valley Province and Eastern Province had the highest percentages of EH landslide susceptibility coverage (21% and 19%, respectively), while the Central Province and Nyanza Province had the lowest percentages of EH landslide susceptibility (5% and 6%, respectively).   Figure 8 shows that over 60% of the landslides had occurred in the extremely high (31.53%, 1 425) and high (29.88%, 127 of 425) landslide susceptibility areas, respectively. Less than 10% of th tal landslides occurred in the area mapped as low (8.24%, 35 of 425) and very low (1.65%, 7 of 42 sceptibility levels. In line with the procedures described in Section 4.5, the ROC curve was draw shown in Figure 9, and the AUC value was 0.86.   Figure 8 shows that over 60% of the landslides had occurred in the extremely high (31.53%, 134 of 425) and high (29.88%, 127 of 425) landslide susceptibility areas, respectively. Less than 10% of the total landslides occurred in the area mapped as low (8.24%, 35 of 425) and very low (1.65%, 7 of 425) susceptibility levels. In line with the procedures described in Section 4.5, the ROC curve was drawn as shown in Figure 9, and the AUC value was 0.86.  Figure 8 shows that over 60% of the landslides had occurred in the extremely high (31.53%, 134 of 425) and high (29.88%, 127 of 425) landslide susceptibility areas, respectively. Less than 10% of the total landslides occurred in the area mapped as low (8.24%, 35 of 425) and very low (1.65%, 7 of 425) susceptibility levels. In line with the procedures described in Section 4.5, the ROC curve was drawn as shown in Figure 9, and the AUC value was 0.86.

Discussion
This research was proposed to apply the TFN-AHP method to map landslide susceptibility in Kenya. Figure 5 directly displays the visible landslide susceptibility information for the entire Kenyan territory, indicating the likelihood of potential landslides. As a developing country, such information would greatly benefit Kenyan efforts to minimize landslide-induced losses and develop optimized land management policies. From Figure 6, it was observed that 44.6% of Kenya is classified as highand extremely high-susceptibility zones, whereas 26.11% of Kenya was mapped as having low and very low susceptibility. High and extremely high landslide susceptibility zones predominantly cover the rift valley region and its surrounding areas. This finding can be attributed to plentiful rainfall, steep terrains, and fractured ground. Low and very low landslide susceptibility areas are primarily distributed in the southwestern and coastal regions. The distribution of susceptibility classifications also varies in different provinces (Figure 7). The Rift Valley Province had a majority of the historical landslides. This province has the largest area coverage of extremely high-and high-susceptibility zones.
Dozens of methods have been used in LSM at different scales. For large areas with poor availability of historical landslide inventories, the spatial multicriteria evaluation (SMCE) method has exhibited overwhelming advantages over statistical and the physically based methods [15,42]. As a representative SMCE method, a review of previous studies (as displayed in Table 1) has suggested that the AHP and its fuzzy extensions are one of the favorable methods in LSM for large areas (e.g., for a whole country). One limitation of the application of the statistical method in this study is the incompleteness of the historical landslide inventory, which reduces the reliability of the results. Despite this, the historical landslide inventory can still be used for validation in a better than no sense.
The validation results demonstrate that the adopted TFN-AHP resulted in promising accuracy with an AUC value of 0.86 (Figure 9). Despite no strict rules for the evaluation of this accuracy, the resultant accuracy seems to be good compared to similar studies in different areas [1,3,8]. For the LSM, an ideal result map should include as many historical landslides as possible in "high" or "extremely high" susceptibility regions. Additionally, few historical landslides should occur in the "low" or "very low" susceptibility region. Figure 8 shows that the concentration of known landslides decreases from the extremely high category to the very low category. For decision making under multiple criteria, it is difficult for humans to quantify criteria weights using extract numbers. However, rational decisions can be made by skilled experts through a certain value with some uncertainties to capture human subjectivity. Given this, the TFN-AHP makes the comparison process

Discussion
This research was proposed to apply the TFN-AHP method to map landslide susceptibility in Kenya. Figure 5 directly displays the visible landslide susceptibility information for the entire Kenyan territory, indicating the likelihood of potential landslides. As a developing country, such information would greatly benefit Kenyan efforts to minimize landslide-induced losses and develop optimized land management policies. From Figure 6, it was observed that 44.6% of Kenya is classified as highand extremely high-susceptibility zones, whereas 26.11% of Kenya was mapped as having low and very low susceptibility. High and extremely high landslide susceptibility zones predominantly cover the rift valley region and its surrounding areas. This finding can be attributed to plentiful rainfall, steep terrains, and fractured ground. Low and very low landslide susceptibility areas are primarily distributed in the southwestern and coastal regions. The distribution of susceptibility classifications also varies in different provinces (Figure 7). The Rift Valley Province had a majority of the historical landslides. This province has the largest area coverage of extremely high-and high-susceptibility zones.
Dozens of methods have been used in LSM at different scales. For large areas with poor availability of historical landslide inventories, the spatial multicriteria evaluation (SMCE) method has exhibited overwhelming advantages over statistical and the physically based methods [15,42]. As a representative SMCE method, a review of previous studies (as displayed in Table 1) has suggested that the AHP and its fuzzy extensions are one of the favorable methods in LSM for large areas (e.g., for a whole country). One limitation of the application of the statistical method in this study is the incompleteness of the historical landslide inventory, which reduces the reliability of the results. Despite this, the historical landslide inventory can still be used for validation in a better than no sense.
The validation results demonstrate that the adopted TFN-AHP resulted in promising accuracy with an AUC value of 0.86 (Figure 9). Despite no strict rules for the evaluation of this accuracy, the resultant accuracy seems to be good compared to similar studies in different areas [1,3,8]. For the LSM, an ideal result map should include as many historical landslides as possible in "high" or "extremely high" susceptibility regions. Additionally, few historical landslides should occur in the "low" or "very low" susceptibility region. Figure 8 shows that the concentration of known landslides decreases from the extremely high category to the very low category. For decision making under multiple criteria, it is difficult for humans to quantify criteria weights using extract numbers. However, rational decisions can be made by skilled experts through a certain value with some uncertainties to capture human subjectivity. Given this, the TFN-AHP makes the comparison process more flexible to minimize the objectivities and uncertainties involved in the conventional AHP process. For the purpose of comparison, the map produced using the conventional AHP is illustrated in Figure 10. The ROC analysis, with an AUC value of 0.72, is also plotted in Figure 9. To perform the conventional AHP, the element value of the TFN-AHP comparison matrix is replaced by a single number according to Table 2. The comparison matrix for the conventional AHP analysis is shown in Appendix A. It should be noted that another source of the subjectivities involved in this study and other similar studies using the SMCE methods or statistical methods may originate from the selection of the LCFs. As shown in Table 1, the number of factors used for LSM ranges from 3 to 10. As discussed in many case studies [1,5,9,16,25,27] and more recently reviewed in [13], the selection of LCFs largely depends on conditions such as data availability, scale, and nature of the study area.
Land 2020, 9, x FOR PEER REVIEW 20 of 21 Land 2020, 9, x; doi: FOR PEER REVIEW www.mdpi.com/journal/land more flexible to minimize the objectivities and uncertainties involved in the conventional AHP process. For the purpose of comparison, the map produced using the conventional AHP is illustrated in Figure 10. The ROC analysis, with an AUC value of 0.72, is also plotted in Figure 9. To perform the conventional AHP, the element value of the TFN-AHP comparison matrix is replaced by a single number according to Table 2. The comparison matrix for the conventional AHP analysis is shown in Appendix A. It should be noted that another source of the subjectivities involved in this study and other similar studies using the SMCE methods or statistical methods may originate from the selection of the LCFs. As shown in Table 1, the number of factors used for LSM ranges from 3 to 10. As discussed in many case studies [1,5,9,16,25,27] and more recently reviewed in [13], the selection of LCFs largely depends on conditions such as data availability, scale, and nature of the study area. Even for skilled departments from China and many other developing countries, no universal rules have been proposed. Hence, for a given study area, comparative studies have always been conducted to select the best maps.

Conclusions
In this research, an integrated method of fuzzy theory and conventional AHP analysis was employed for the LSM of Kenya. A two-level hierarchical index system was established to predict landslide susceptibility with a GIS platform. Ten factors contributing to landslide occurrence were included in the first level of the evaluation system. These contributing factors included slope, altitude, aspect, SPI, TWI, curvature, land use, MAP, landform, and soil texture. For the second level, each of these factors was divided into several subclasses. The weights of these factors and their subclasses were determined using the adopted TFN-AHP theory. A nationwide landslide susceptibility map for the entire Kenyan territory was produced with five different levels ranging from extremely high susceptibility to very low susceptibility. Extremely high and high landslide susceptibility zones primarily covered the rift valley and its nearby regions. Validation results using ROC curves indicated that the TFN-AHP method performed well for developing LS maps of the study area. This method resulted in a higher AUC accuracy than the conventional AHP using the same datasets.
This study was the first attempt to identify landslide susceptibility zones in Kenya on a national scale. The produced map can be used as a general indicator of the relative landslide susceptibility for larger areas rather than an accurate susceptibility measure for each specific site. The results would be Even for skilled departments from China and many other developing countries, no universal rules have been proposed. Hence, for a given study area, comparative studies have always been conducted to select the best maps.

Conclusions
In this research, an integrated method of fuzzy theory and conventional AHP analysis was employed for the LSM of Kenya. A two-level hierarchical index system was established to predict landslide susceptibility with a GIS platform. Ten factors contributing to landslide occurrence were included in the first level of the evaluation system. These contributing factors included slope, altitude, aspect, SPI, TWI, curvature, land use, MAP, landform, and soil texture. For the second level, each of these factors was divided into several subclasses. The weights of these factors and their subclasses were determined using the adopted TFN-AHP theory. A nationwide landslide susceptibility map for the entire Kenyan territory was produced with five different levels ranging from extremely high susceptibility to very low susceptibility. Extremely high and high landslide susceptibility zones primarily covered the rift valley and its nearby regions. Validation results using ROC curves indicated that the TFN-AHP method performed well for developing LS maps of the study area. This method resulted in a higher AUC accuracy than the conventional AHP using the same datasets.
This study was the first attempt to identify landslide susceptibility zones in Kenya on a national scale. The produced map can be used as a general indicator of the relative landslide susceptibility for larger areas rather than an accurate susceptibility measure for each specific site. The results would be helpful in various land resources-related fields to inform decision making, such as regional landslide hazard mitigation, land use management, and infrastructure planning.