Landslide Susceptibility Assessment of a Part of the Western Ghats (India) Employing the AHP and F-AHP Models and Comparison with Existing Susceptibility Maps

Abstract

Landslides are prevalent in the Western Ghats, and the incidences that happened in 2021 in the Koottickal area of the Kottayam district (Western Ghats) resulted in the loss of 10 lives. The objectives of this study are to assess the landslide susceptibility of the high-range local self-governments (LSGs) in the Kottayam district using the analytical hierarchy process (AHP) and fuzzy-AHP (F-AHP) models and to compare the performance of existing landslide susceptible maps. This area never witnessed any massive landslides of this dimension, which warrants the necessity of relooking into the existing landslide-susceptible models. For AHP and F-AHP modeling, ten conditioning factors were selected: slope, soil texture, land use/land cover (LULC), geomorphology, road buffer, lithology, and satellite image-derived indices such as the normalized difference road landslide index (NDRLI), the normalized difference water index (NDWI), the normalized burn ratio (NBR), and the soil-adjusted vegetation index (SAVI). The landslide-susceptible zones were categorized into three: low, moderate, and high. The validation of the maps created using the receiver operating characteristic (ROC) technique ascertained the performances of the AHP, F-AHP, and TISSA maps as excellent, with an area under the ROC curve (AUC) value above 0.80, and the NCESS map as acceptable, with an AUC value above 0.70. Though the difference is negligible, the map prepared using the TISSA model has better performance (AUC = 0.889) than the F-AHP (AUC = 0.872), AHP (AUC = 0.867), and NCESS (AUC = 0.789) models. The validation of maps employing other matrices such as accuracy, mean absolute error (MAE), and root mean square error (RMSE) also confirmed that the TISSA model (0.869, 0.226, and 0.122, respectively) has better performance, followed by the F-AHP (0.856, 0.243, and 0.147, respectively), AHP (0.855, 0.249, and 0.159, respectively), and NCESS (0.770, 0.309, and 0.177, respectively) models. The most landslide-inducing factors in this area that were identified through this study are slope, soil texture, LULC, geomorphology, and NDRLI. Koottickal, Poonjar-Thekkekara, Moonnilavu, Thalanad, and Koruthodu are the LSGs that are highly susceptible to landslides. The identification of landslide-susceptible areas using diversified techniques will aid decision-makers in identifying critical infrastructure at risk and alternate routes for emergency evacuation of people to safer terrain during an exigency.

1. Introduction

Landslides are one of the most destructive disasters occurring worldwide, with an estimated loss of USD 20 billion [1] and causing more than 4000 deaths annually worldwide [2]. The recent global dataset of landslide fatalities [2] confirmed 55,997 deaths (from 4862 distinct occurrences) between January 2004 and December 2016. According to Li et al. [3], a total of 9223 landslides have been recorded worldwide between 2000 and 2020. Froude and Petley [2] reported that 75% of the landslides occurred in Asia and listed India as the country with the highest number of incidences. Additionally, Broeckx et al. [4] reported that every year, landslides mobilize roughly 56 billion m³ of hillslope material, of which Asia accounts for 68% (38 billion m³) and the Central and Southeast Asian mountains for 50% (28 billion m³). The Himalayas, the Western Ghats, and the Eastern Ghats of India were listed as the most landslide-prone regions by the National Disaster Management Authority of India [5].

In Kerala’s Western Ghats region, landslides are the most prevalent hazard, which is temporally restricted to the monsoon season. In the year 2018, a total of 4728 landslides occurred in Kerala [6,7]. The catastrophic landslide disasters that occurred in Kerala in the recent past include Puthumala (which resulted in the deaths of 17 people), Kavalappara (which killed 59 people) in the year 2019 [8], and the Pettimudi disaster (which killed 70 people) in the year 2020 [9]. Two such major landslides that were reported on 16 October 2021, were in the Koottickal village of the Kottayam district (Figure 1). These landslides resulted in the deaths of ten people (four at Plappally and six at Kavali) [10]. Such catastrophic landslides demand the necessity of relooking into the existing landslide susceptibility map and warrant additional models if required. This will help the government agencies that deal with landslide hazards find the elements at risk and also help minimize the effects of imminent landslides. Thus, this study proposed two models, namely the AHP and F-AHP, and compared the results with existing landslide-susceptibility maps of Kerala, namely the GIS Tool for Infinite Slope Stability Analysis (GIS-TISSA) map [11,12] and the map created by the National Centre for Earth Science Studies (NCESS) [13].

Worldwide, researchers employed methods such as the AHP [14,15], F-AHP [15], fuzzy logic [16], frequency ratio [17], analytical network process [18], FROC [19], and Mora–Vahrson–Mora (MVM) [20,21] for the mapping of landslide-susceptible zones. Researchers also applied artificial intelligence (AI)/machine learning (ML) models such as support vector machine (SVM) [22,23,24], Naïve Bayes (NB) [25,26,27], decision tree [28,29,30], K-nearest neighbor [31,32], random forest (RF) [33,34], adaptive neuro-fuzzy inference system (ANFIS) [35], convolutional neural network (CNN) [36,37,38], artificial neural network [39,40], logistic regression [41,42,43], support vector regression [44,45], recurrent neural network [36,37], Adaptive Boosting [46], extreme gradient boosting [35], Random Subspace (RSS) [47], Reduced Error Pruning Tree (REPTree) [48], etc., for landslide susceptibility modeling. The mixed ensemble models of ML techniques such as multi-layer perceptron [49], ANFIS [50], genetic algorithm [51], RSS [52], Credal Decision Tree [53], boosted regression tree [54], SVM [55], RF [56,57], REPTree [58], deep learning neural network [59], Multi-layer Perceptron Neural Network Classifier [60], NB [61], logistic model tree [62], Bagging [63], CNN [64], and Particle Swarm Optimization [65] have also been employed for landslide susceptibility modeling. The effectiveness and efficiency of an ML solution are determined by the quality and type of data as well as the performance of the learning algorithms [66].

Amidst the plethora of models available, we preferred the AHP and F-AHP, as the AHP simplifies complex issues into specific sub-issues using pair comparisons, and the factors are prioritized once the hierarchical structure is constructed [67]. Because there are just two options to evaluate, the decision-maker can only prioritize those that are being compared [67]. The F-AHP is used to solve the AHP’s inability to handle imprecision and subjectivity in assessments [68]. In addition, the AHP and F-AHP will provide better results in data-scarce conditions when compared to ML methods [69,70]. Hence, the AHP and F-AHP models are suitable when a large number of landslides are not available to train and validate the results. Abdi et al. [71] and Akshaya et al. [15] employed both the AHP and F-AHP models for demarcating landslide-susceptible zones and found both models effective. Furthermore, the use of the AHP and F-AHP on the one hand, and comparison with existing susceptible maps on the other, distinguishes this work. The susceptibility maps were created using the AHP and F-AHP models and the existing NCESS landslide susceptibility map, which is also a product of multi-criteria decision-making (MCDM) analysis. However, the existing GIS-TISSA map for Kerala state was not created using MCDM analysis. Furthermore, this study also has the distinction of comparing two MCDM models with an existing MCDM-derived map and a non-MCDM-derived susceptibility map.

To fulfill the aim of this study, an ideal study area is warranted. With numerous recent landslides in the high-range local self-government (LSG; LSG is the smallest administrative unit in India) bodies in the Kottayam district, resting in the Western Ghats, thus became the obvious choice. Thus, this study applied the AHP and F-AHP models to the aftermath of the 2021 Koottickal landslide disaster, compared these models with the GIS-TISSA [11,12] and NCESS [13] models, and assessed the influence of the conditioning factors (slope, soil texture, land use/land cover, geomorphology, road buffer, normalized difference road landslide index, normalized difference water index, normalized burnt ratio, soil-adjusted vegetation index, and lithology) on landslide initiation. The intention of considering LSG rather than a geographic unit (such as a river basin) is because a geographic unit might pass through different administrative units, and hence, implementing management practices becomes an administrative hurdle.

2. Materials and Methods

2.1. Study Area

The area selected for this study lies between 9°20′ and 9°55′north latitudes, and 76°30′ and 77°0′east longitudes. This area comprises 20 LSGs (19 panchayat (panchayat is the LSG in a hamlet) and one municipality (municipality is the LSG in a town)), namely, Bharananganam, Erumely, Kadanad, Karoor, Koottickal, Koruthodu, Melukavu, Moonnilavu, Mundakkayam, Parathodu, Poonjar, Poonjar-Thekkekara, Ramapuram, Teekoy, Thalanad, Thalappalam, Thidanad, Uzhavoor, Veliyannoor, and Erattupetta municipality. The study area spans around 737.77 km² (Figure 2). A total of 73 landslides occurred in this region during the 2018 monsoon season [6].

2.2. Major Steps and Methodology

The modeling process involves the following eight steps (Figure 3):

• The data were generated from different sources as depicted in

  Table 1. Data sources.

  Data	Source	Thematic Layers Derived	Data Type	Spatial Resolution	Scale
  ASTER GDEM	https://earthexplorer.usgs.gov/ (accessed on 22 August 2022)	Slope angle	Continuous	30 m
  Landsat 8 OLI image	https://earthexplorer.usgs.gov/ (accessed on 15 November 2022)	NBR NDWI SAVI LULC NDRLI Geomorphology	Continuous	30 m
  Soil map	National Bureau of Soil Survey and Land Use Planning	Soil texture	Discrete		1:250,000
  Geological map	Geological Survey of India	Lithology	Discrete		1:50,000
  Topographical map	Survey of India	Road	Discrete		1:50,000
  Google Earth Pro	https://www.google.com/intl/en_in/earth/versions/ (accessed on 24 November 2022)	Road (updated)	Discrete	15 cm to 15 m

  . ArcGIS Pro and ERDAS Imagine 10.0 were utilized to generate different layers.
• The landslide incidence locations for the year 2018 acquired from Hao et al. [6] were utilized for training and validation.
• The continuous factor layers (slope, normalized difference road landslide index (NDRLI), normalized difference water index (NDWI), normalized burnt ratio (NBR), and soil-adjusted vegetation index (SAVI)) were categorized into five classes using the natural breaks method [72,73,74,75].
• The multi-collinearity analysis of the conditioning factors was carried out using R 4.2.2 software.
• To validate the susceptibility maps, the incidence locations were split into two: training (60%) and validation datasets (40%).
• The AHP and F-AHP models were utilized to create maps of landslide susceptibility. The weights for these models were computed using Microsoft Excel and FisPro 3.8 [76,77] software.
• Two existing landslide susceptibility maps, the GIS-TISSA model [11,12] and the NCESS model [13], were used for comparing the prediction capability of these four maps.
• The ROC, accuracy, MAE, and RMSE techniques were used to validate the four susceptibility maps. The R 4.2.2 software was used to create the ROC curves and to compute accuracy, MAE, and RMSE values, and IBM SPSS Statistics 25 was utilized to compute the Kappa index.

2.3. Conditioning Factors

The slope was computed from the ASTER GDEM using ArcGIS surface (spatial analyst) tools. The road cuttings were derived from the topographic maps and Google Earth Pro, and a 100 m buffer was computed utilizing the ArcGIS proximity tool. Soil and lithology were extracted from the National Bureau of Soil Survey and Land Use Planning (NBSS & LUP) soil map and Geological Survey of India (GSI) geological maps, respectively, using ArcGIS tools. The ERDAS Imagine software was employed to extract the LULC types from the Landsat 8 images (acquired on 16 January 2022, Path: 144, Row: 053, Landsat scene identifier: LC81440532022016LGN00). In the ERDAS Imagine software, various LULC types found in this area were classified using the maximum likelihood supervised classification approach [78,79]. The geomorphic classes were extracted by visual interpretation of Landsat 8 images aided by field knowledge. The topographical map and high-resolution Google Earth imagery were used to verify this. The NDRLI, NBR, NDWI, and SAVI were derived from the Landsat images utilizing the ArcGIS map algebra tool. Equations (1) [80], (2) [81,82], (3) [83], and (4) [84] were used to compute the NDRLI, NBR, NDWI, and SAVI, respectively:

$$\mathrm{NDRLI}=\frac{\left(\mathrm{SWIR}1-\mathrm{Blue}\right)}{\left(\mathrm{SWIR}1+\mathrm{Blue}\right)}$$

(1)

$$\mathrm{NBR}=\frac{\left(\mathrm{NIR}-\mathrm{SWIR}\right)}{\left(\mathrm{NIR}+\mathrm{SWIR}\right)}$$

(2)

$$\mathrm{NDWI}=\frac{\left(\mathrm{Green}-\mathrm{NIR}\right)}{\left(\mathrm{Green}+\mathrm{NIR}\right)}$$

(3)

$$\mathrm{SAVI}=\left(\frac{\mathrm{NIR}-\mathrm{Red}}{\mathrm{NIR}+\mathrm{Red}+\mathrm{L}}\right)\times \left(1+\mathrm{L}\right)$$

(4)

where SWIR, Blue, NIR, Green, and Red stand for spectral reflectance in short wave infrared, blue, near-infrared, green, and red bands, respectively, and L stands for the soil brightness correction factor (0.5).

2.4. Accuracy Assessment Using Cohen’s Kappa Coefficient

Cohen’s Kappa coefficient [85] was used to assess the LULC classification. SPSS software was utilized to compute the Kappa coefficient using Equation (5) [86]. A total of 110 locations have been collected with the aid of a handheld Global Positioning System for the accuracy assessment. A Kappa coefficient value between 0.8 and 0.9, and between 0.9 and 1.0 represent strong and almost perfect agreements, respectively [87].

$$\mathrm{k}=\frac{{\mathrm{P}}_{\mathrm{o}}-{\mathrm{P}}_{\mathrm{c}}}{1-{\mathrm{P}}_{\mathrm{c}}}$$

(5)

where ${\mathrm{P}}_{\mathrm{o}}$= proportion of observed similarity and ${\mathrm{P}}_{\mathrm{c}}$ = proportion of similarities by chance

2.5. Multi-Collinearity Analysis

Multicollinearity occurs when one or more independent variables have a significant correlation with one another [88]. The tolerance and variance expansion factor (VIF) can be used to assess multicollinearity [89]. A factor is believed to have multicollinearity when the VIF is greater than 10 and the tolerance is less than 0.1 [90]. The VIF was computed using Equation (6) [91]. Tolerance is the reciprocal of the VIF and was determined using Equation (7) [92]:

$$VI{F}_j=\frac{1}{1-R_j^2}$$

(6)

$$Tol=1-R_j{}^2$$

(7)

where ${\mathrm{R}}_{\mathrm{j}}^2$ is the multiple coefficients of determination of the jth factor.

2.6. AHP Modeling

The AHP is an MCDM method established by Thomas L. Saaty [93]. This approach deconstructs complex problems into a hierarchy and identifies the best solution for the goal [94]. The most important steps are making the matrix for pair-wise comparisons and finding the eigenvector, the weighting coefficient (

Table 2. Pairwise comparison matrix.

	SA	Soil	LULC	Geom.	RB	NDRLI	NDWI	NBR	SAVI	Litho.	Vp	Cp
SA	1	2	3	4	5	6	7	8	9	10	4.529	0.290
Soil	1/2	1	2	3	4	5	6	7	8	9	3.392	0.218
LULC	1/3	1/2	1	2	3	4	5	6	7	8	2.414	0.155
Geom.	1/4	1/3	1/2	1	2	3	4	5	6	7	1.707	0.109
RB	1/5	1/4	1/3	1/2	1	2	3	4	5	6	1.196	0.077
NDRLI	1/6	1/5	1/4	1/3	1/2	1	2	3	4	5	0.836	0.054
NDWI	1/7	1/6	1/5	1/4	1/3	1/2	1	2	3	4	0.586	0.038
NBR	1/8	1/7	1/6	1/5	1/4	1/3	1/2	1	2	3	0.414	0.027
SAVI	1/9	1/8	1/7	1/6	1/5	1/4	1/3	1/2	1	2	0.298	0.019
Litho.	1/10	1/9	1/8	1/7	1/6	1/5	1/4	1/3	1/2	1	0.221	0.014
∑	2.93	4.83	7.72	11.59	16.45	22.28	29.08	36.83	45.50	55.00	15.59	1.00

SA = slope angle, Geom. = geomorphology, RB = road buffer, Litho. = lithology.

), and the consistency ratio (CR) (

Table 3. Normalized matrix.

	SA	Soil	LULC	Geom.	RB	NDRLI	NDWI	NBR	SAVI	Litho.	∑ Rank	[C]	[D] = [A]*[C]	[E] = [D]/[C]	λmax	CI	CR
SA	0.34	0.41	0.39	0.35	0.30	0.27	0.24	0.22	0.20	0.18	2.90	0.290	3.17	10.95	10.54	0.06	0.04 (4.1%)
Soil	0.17	0.21	0.26	0.26	0.24	0.22	0.21	0.19	0.18	0.16	2.10	0.210	2.32	11.06
LULC	0.11	0.10	0.13	0.17	0.18	0.18	0.17	0.16	0.15	0.15	1.52	0.152	1.66	11.00
Geom.	0.09	0.07	0.06	0.09	0.12	0.13	0.14	0.14	0.13	0.13	1.09	0.109	1.18	10.83
RB	0.07	0.05	0.04	0.04	0.06	0.09	0.10	0.11	0.11	0.11	0.79	0.079	0.83	10.58
NDRLI	0.06	0.04	0.03	0.03	0.03	0.04	0.07	0.08	0.09	0.09	0.56	0.056	0.58	10.34
NDWI	0.05	0.03	0.03	0.02	0.02	0.02	0.03	0.05	0.07	0.07	0.40	0.040	0.40	10.15
NBR	0.04	0.03	0.02	0.02	0.02	0.01	0.02	0.03	0.04	0.05	0.28	0.028	0.28	10.08
SAVI	0.04	0.03	0.02	0.01	0.01	0.01	0.01	0.01	0.02	0.04	0.20	0.020	0.20	10.14
Litho.	0.01	0.02	0.02	0.01	0.01	0.01	0.01	0.01	0.01	0.02	0.15	0.015	0.15	10.29
∑	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	1.00	10.00	1.00		105.46

) [73].

Equations (8) and (9) were used to compute the eigenvector (Vp) and weighting coefficient (Cp) [73].

$$\mathrm{Vp}=\sqrt[\mathrm{k}]{\mathrm{W}1\times \dots \mathrm{Wk}}$$

(8)

where k = no. of factors and W = ratings.

$$\mathrm{Cp}=\frac{\mathrm{Vp}}{\mathrm{Vp}1+\dots \mathrm{Vpk}}$$

(9)

As described by Danumah et al. [95], the normalized matrix, priority vector (C), overall priority (D), and rational priority (E) were determined.

Equations (10)–(12) were used to compute the eigenvalue (λmax), consistency index (CI), and consistency ratio [73]:

$${\lambda }_{\max }=\frac{\left[\mathrm{E}\right]}{\mathrm{k}}$$

(10)

$$\mathrm{CI}=\left(\lambda \max -\mathrm{k}\right)/\left(\mathrm{k}-1\right)$$

(11)

$$\mathrm{CR}=\frac{\mathrm{CI}}{\mathrm{RI}}$$

(12)

where the (RI) random index = 1.49 [93].

A consistency ratio value of less than 0.1 is admissible [93]. The analysis is repeated if the CR is greater than 0.1. This model yields an acceptable CR (0.041). As a result, the outcomes are trustworthy.

The AHP weights are depicted in Equation (13).

$$\begin{array}{l}\mathrm{LSZ}=\left(0.290\times \mathrm{SA}\right)+\left(0.218\times \mathrm{Soil}\right)+\left(0.155\times \mathrm{LULC}\right)+\left(0.109\times \mathrm{Geom}.\right)+\\ \left(0.077\times \mathrm{RB}\right)+\left(0.054\times \mathrm{NDRLI}\right)+\left(0.038\times \mathrm{NDWI}\right)+\left(0.027\times \mathrm{NBR}\right)+\\ \left(0.019\times \mathrm{SAVI}\right)+\left(0.014\times \mathrm{Litho}.\right)\end{array}$$

(13)

2.7. F-AHP Modeling

The F-AHP, which is an ensemble of the AHP and fuzzy logic models, was used to weigh the factors [96]. This model solves the limitations of the AHP model by allowing decision-makers to evaluate their preferences within an acceptable interval [97]. Buckley [98] presented a method for comparing fuzzy ratios (triangular membership functions), and that method was employed in this study. The major procedures involved are the establishment of a pair-wise comparison matrix (

Table 4. Pair-wise comparisons of factors.

	SA	Soil	LULC	Geom.	RB	NDRLI	NDWI	NBR	SAVI	Litho.
SA	(1,1,1)	(1,2,3)	(2,3,4)	(3,4,5)	(4,5,6)	(5,6,7)	(6,7,8)	(7,8,9)	(8,9,10)	(10,10,10)
Soil	(1/3,1/2,1)	(1,1,1)	(1,2,3)	(2,3,4)	(3,4,5)	(4,5,6)	(5,6,7)	(6,7,8)	(7,8,9)	(8,9,10)
LULC	(1/4,1/3,1/2)	(1/3,1/2,1)	(1,1,1)	(1,2,3)	(2,3,4)	(3,4,5)	(4,5,6)	(5,6,7)	(6,7,8)	(7,8,9)
Geom	(1/5,1/4,1/3)	(1/4,1/3,1/2)	(1/3,1/2,1)	(1,1,1)	(1,2,3)	(2,3,4)	(3,4,5)	(4,5,6)	(5,6,7)	(6,7,8)
RB	(1/6,1/5,1/4)	(1/5,1/4,1/3)	(1/4,1/3,1/2)	(1/3,1/2,1)	(1,1,1)	(1,2,3)	(2,3,4)	(3,4,5)	(4,5,6)	(5,6,7)
NDRLI	(1/7,1/6,1/5)	(1/6,1/5,1/4)	(1/5,1/4,1/3)	(1/4,1/3,1/2)	(1/3,1/2,1)	(1,1,1)	(1,2,3)	(2,3,4)	(3,4,5)	(4,5,6)
NDWI	(1/8,1/7,1/6)	(1/7,1/6,1/5)	(1/6,1/5,1/4)	(1/5,1/4,1/3)	(1/4,1/3,1/2)	(1/3,1/2,1)	(1,1,1)	(1,2,3)	(2,3,4)	(3,4,5)
NBR	(1/9,1/8,1/7)	(1/8,1/7,1/6)	(1/7,1/6,1/5)	(1/6,1/5,1/4)	(1/5,1/4,1/3)	(1/4,1/3,1/2)	(1/3,1/2,1)	(1,1,1)	(1,2,3)	(2,3,4)
SAVI	(1/10,1/9,1/8)	(1/9,1/8,1/7)	(1/8,1/7,1/6)	(1/7,1/6,1/5)	(1/6,1/5,1/4)	(1/5,1/4,1/3)	(1/4,1/3,1/2)	(1/3,1/2,1)	(1,1,1)	(1,2,3)
Litho	(1/10,1/10,1/10)	(1/10,1/9,1/8)	(1/9,1/8,1/7)	(1/8,1/7,1/6)	(1/7,1/6,1/5)	(1/6,1/5,1/4)	(1/5,1/4,1/3)	(1/4,1/3,1/2)	(1/3,1/2,1)	(1,1,1)

), the computation of the geometric mean (

Table 5. Geometric means of fuzzy comparison values.

	Fuzzy Geometric Mean Value (${\tilde{\mathbf{r}}}_{\mathbf{i}}$)
SA	3.64	4.53	5.32
Soil	2.59	3.36	4.23
LULC	1.83	2.41	3.13
Geom.	1.28	1.71	2.25
RB	0.90	1.20	1.59
NDRLI	0.63	0.84	1.12
NDWI	0.44	0.59	0.78
NBR	0.32	0.41	0.55
SAVI	0.24	0.30	0.39
Litho.	0.19	0.22	0.28

), the determination of fuzzy weights (see

Table 6. Relative fuzzy weights of each factor.

	Fuzzy Weight (${\tilde{\mathbf{W}}}_{\mathit{i}}$)
SA	0.19	0.29	0.44
Soil	0.13	0.22	0.35
LULC	0.09	0.16	0.26
Geom.	0.07	0.11	0.19
RB	0.05	0.08	0.13
NDRLI	0.03	0.05	0.09
NDWI	0.02	0.04	0.06
NBR	0.02	0.03	0.05
SAVI	0.01	0.02	0.03
Litho.	0.01	0.01	0.02
$\sum {\tilde{r}}_i$	12.046	15.557	19.626
${\left(\sum {\tilde{r}}_i\right)}^{-1}$	0.051	0.064	0.083

), and the computation of averaged and normalized relative weights (

Table 7. Averaged and normalized relative weights of factor.

	Weight (${\mathbf{M}}_{\mathbf{i}}$)	Normalized Weight (${\mathbf{N}}_{\mathbf{i}}$)
SA	0.306	0.283
Soil	0.233	0.215
LULC	0.169	0.157
Geom.	0.121	0.112
RB	0.085	0.079
NDRLI	0.059	0.055
NDWI	0.042	0.039
NBR	0.029	0.027
SAVI	0.021	0.020
Litho.	0.016	0.01
∑	1.08	1.00

) [73]. The following are the steps in F-AHP modeling:

Step 1:

The factors were compared.

The fuzzy triangular scale (1/4, 1/3, 1/2) will be used when factor 1 (P1) is less essential than factor 2 (P2). For the comparison matrix, the scale will be (1/4, 1/3, 1/2) [99].

Equation (14) [73] depicts the matrix:

$${\tilde{\mathrm{A}}}^{\mathrm{k}}=\left[\begin{array}{ccc}{\tilde{\mathrm{d}}}_{11}^{\mathrm{k}}& {\tilde{\mathrm{d}}}_{12}^{\mathrm{k}}& \begin{array}{cc}\dots & {\tilde{\mathrm{d}}}_{1\mathrm{n}}^{\mathrm{k}}\end{array}\\ \begin{array}{c}{\tilde{\mathrm{d}}}_{21}^{\mathrm{k}}\\ \dots \end{array}& \begin{array}{c}\dots \\ \dots \end{array}& \begin{array}{c}\begin{array}{cc}\dots & {\tilde{\mathrm{d}}}_{2\mathrm{n}}^{\mathrm{k}}\end{array}\\ \begin{array}{cc}\dots & \dots \end{array}\end{array}\\ {\tilde{\mathrm{d}}}_{\mathrm{n}1}^{\mathrm{k}}& {\tilde{\mathrm{d}}}_{\mathrm{n}2}^{\mathrm{k}}& \begin{array}{cc}\dots & {\tilde{\mathrm{d}}}_{\mathrm{mn}}^{\mathrm{k}}\end{array}\end{array}\right]$$

(14)

where ${\tilde{\mathrm{d}}}_{\mathrm{ij}}^{\mathrm{k}}$ reflects the kth decision maker’s priority for the ith factor over the jth factor [99].

Step 2:

${\tilde{\mathrm{d}}}_{\mathrm{ij}}$ was computed using Equation (15) [73].

$${\tilde{\mathrm{d}}}_{\mathrm{ij}}=\frac{{{\displaystyle \sum }}_{\mathrm{k}=1}^{\mathrm{K}}{\tilde{\mathrm{d}}}_{\mathrm{ij}}^{\mathrm{k}}}{\mathrm{K}}$$

(15)

Step 3:

Equation (16) [73] was used to modify the matrix.

$$\tilde{\mathrm{A}}=\left[\begin{array}{ccc}{\tilde{\mathrm{d}}}_{11}& \cdots & {\tilde{\mathrm{d}}}_{1\mathrm{n}}\\ ⋮& \ddots & ⋮\\ {\tilde{\mathrm{d}}}_{\mathrm{n}1}& \cdots & {\tilde{\mathrm{d}}}_{\mathrm{nn}}\end{array}\right]$$

(16)

Step 4:

Equation (17) [98] was used to compute the geometric average.

$${\tilde{\mathrm{r}}}_{\mathrm{i}}={\left({\displaystyle \prod }_{\mathrm{j}=1}^{\mathrm{n}}{\tilde{\mathrm{d}}}_{\mathrm{ij}}\right)}^{\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$\mathrm{n}$}\right.},\mathrm{i}=1,2,\dots ,\mathrm{n}$$

(17)

where ${\tilde{\mathrm{r}}}_{\mathrm{i}}$ = triangular values.

Step 5:

The next three sub-processes (5a, 5b, and 5c) were used to compute the fuzzy weight:

Step 5a:

The vector summation was computed;

Step 5b:

The fuzzy triangular number was replaced to arrange it in ascending order after computing the summation vector’s (−1) power;

Step 5c:

Each ${\tilde{\mathrm{r}}}_{\mathrm{i}}$ was multiplied by the reverse vector to determine the fuzzy weight as shown in Equations (18) and (19).

$${\tilde{\mathrm{w}}}_{\mathrm{i}}={\tilde{\mathrm{r}}}_{\mathrm{i}}\times {\left({\tilde{\mathrm{r}}}_1+{\tilde{\mathrm{r}}}_2+\dots +{\tilde{\mathrm{r}}}_{\mathrm{n}}\right)}^{-1}$$

(18)

$${\tilde{\mathrm{w}}}_{\mathrm{i}}=\left({\mathrm{lw}}_{\mathrm{i}},{\mathrm{mw}}_{\mathrm{i}},{\mathrm{uw}}_{\mathrm{i}}\right)$$

(19)

Step 6:

The de-fuzzification of fuzzy weights using Equation (20) [100].

$${\mathrm{M}}_{\mathrm{i}}=\frac{{\mathrm{lw}}_{\mathrm{i}},{\mathrm{mw}}_{\mathrm{i}},{\mathrm{uw}}_{\mathrm{i}}}{3}$$

(20)

Step 7:

For the standardization of ${\mathrm{M}}_{\mathrm{i}}$, Equation (21) [73] was used.

$${\mathrm{N}}_{\mathrm{i}}=\frac{{\mathrm{M}}_{\mathrm{i}}}{{{\displaystyle \sum }}_{\mathrm{i}=1}^{\mathrm{n}}{\mathrm{M}}_{\mathrm{i}}}$$

(21)

The F-AHP weights are depicted in Equation (22).

$$\begin{array}{l}\mathrm{LSZ}=\left(0.283\times \mathrm{SA}\right)+\left(0.215\times \mathrm{Soil}\right)+\left(0.157\times \mathrm{LULC}\right)+\left(0.112\times \mathrm{Geom}.\right)+\\ \left(0.079\times \mathrm{RB}\right)+\left(0.055\times \mathrm{NDRLI}\right)+\left(0.039\times \mathrm{NDWI}\right)+\left(0.027\times \mathrm{NBR}\right)+\\ \left(0.020\times \mathrm{SAVI}\right)+\left(0.010\times \mathrm{Litho}.\right)\end{array}$$

(22)

2.8. GIS-TISSA Model

The GIS-TISSA [11], a python-based model, offers an alternative way of deriving a factor of safety (FS) (slope stability) in a GIS domain. This model represents slope stability as an FS; the slope is stable when the FS is greater than 1 and unstable when the FS is less than 1 [101]. The model estimates where a landslide is most likely to start, but the spreading process, in which the mass of the failing slope breaks away from the ground and goes downslope, must be assessed individually using different methods.

2.9. NCESS Model

The landslide susceptibility map, which was created by NCESS and made available to the public by the Kerala State Disaster Management Authority [13], is the map that is included in the disaster management plans of the state and districts [102].

2.10. Validation of the Maps

2.10.1. ROC Technique

A two-dimensional metric of classification performance is the ROC technique [103]. The model has outstanding performance when the AUC value is within the range of 0.9–1.0; excellent performance when the value is within the range of 0.8–0.9; and acceptable performance when the value is within the range of 0.7–0.8 [75]. The AUC values ranged between 0.6 and 0.7, and values of 0.5 and 0.6 depicted poor performance and failure, respectively [104]. Utilizing the training and validation datasets, the AUC value can reflect the model’s success and prediction rate [105]. In this study, the incidence data were randomly divided into training and validation datasets of 60% and 40% [106], respectively, and the AUC values were calculated with SPSS software. The 60–40% split is acceptable when the dataset is small [106].

2.10.2. Accuracy

The ratio of true instances retrieved (both positive and negative) among all instances retrieved is the accuracy [107]. Accuracy was computed using Equation (23) [108]:

$$\mathrm{Accuracy}=\frac{\mathrm{TP}+\mathrm{TN}}{\mathrm{TP}+\mathrm{TN}+\mathrm{FP}+\mathrm{FN}}$$

(23)

where TP = true positive, TN = true negative, FP = false positive, and FN = false negative.

2.10.3. MAE

The average of the absolute error values is referred to as the MAE [109]. The MAE was computed using Equation (24) [110]:

$$\mathrm{MAE}=\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}\left|{\mathrm{Y}}_{\mathrm{i}}-{\tilde{\mathrm{Y}}}_{\mathrm{i}}\right|$$

(24)

where ${\mathrm{Y}}_{\mathrm{i}}$ is the predicted value, $\tilde{\mathrm{Y}}{ }_{\mathrm{i}}$ is the observed value, and n is the number of incidences/points.

2.10.4. RMSE

The square root of the mean of the square of all the errors is known as the RMSE [111]. The RMSE was computed using Equation (25) [112]:

$$\mathrm{RMSE}=\sqrt{\frac{1}{\mathrm{n}}{\displaystyle \sum }_{\mathrm{i}=1}^{\mathrm{n}}{\left({\mathrm{Y}}_{\mathrm{i}}-{\tilde{\mathrm{Y}}}_{\mathrm{i}}\right)}^2}$$

(25)

3. Results

3.1. Multicollinearity

The collinearity statistics are shown in

Table 8. Collinearity statistics.

Factors	Collinearity Statistics
	Tolerance	VIF
Slope	0.829	1.206
Soil	0.853	1.173
LULC	0.932	1.073
Geomorphology	0.853	1.173
Road buffer	0.920	1.088
NDRLI	0.208	4.808
NDWI	0.206	4.854
NBR	0.215	4.651
SAVI	0.217	4.608
Lithology	0.868	1.152

. From the analysis, it was found that LULC (0.932) has the highest tolerance value, followed by road buffer (0.920), whereas NDWI (0.206) has the lowest tolerance value, followed by NDRLI (0.208). The VIF values range from 1.073 (LULC) to 4.854 (NDWI). This ascertained that there are no collinearity issues among the factors since the tolerance value is above 0.1 and the VIF value is below 10 for all the conditioning factors.

3.2. Conditioning Factors

3.2.1. Slope

One of the most relevant conditioning factors impacting the probability of landslides is the slope gradient [113]. Topography with a higher slope is more likely to fail [17,114]. The shear stress in the soil rises as the slope angle rises, causing earth materials to move downslope [17]. The slope, which ranges from 0 to 74.47°, is depicted in Figure 4a. The steeper slopes are found along the eastern portion of this area, where a considerable number of landslide incidences have occurred. It was found that 32 (43.83%) and 30 (41.09%) landslides happened in the classes with slopes between 22.48° and 31.83° and 31.83° and 74.47°, respectively.

3.2.2. Soil

In steep terrain, the clay content in the soil forms a possible slip zone, which can induce slope failure and landslides [114]. As clayey soil has low permeability and high porosity, it holds more water than other soil types [17]. The shear strength of the soil is reduced when pore–water pressure rises, resulting in slope failure [115]. The soil types observed in this region include loam, gravelly loam, clay, and gravelly clay (Figure 4b). A total of 63 (86.30%) landslides were recorded in locations with clayey soil, which confirms the above statement.

3.2.3. Land Use/Land Cover (LULC)

LULC reflects natural and human-made influences that affect slope stability in a given location, either directly or indirectly [114]. Landslides are less likely in areas with vegetation because plant roots provide mechanical support, which enhances the soil’s shearing strength [116]. Shu et al. [117] observed that landslide density was substantially lower in areas covered by forests and shrubs. According to Abedin et al. [118], built-up areas are one of the most landslide-prone land use types. This is due to the modifications made to mountainous terrain for various development activities. Landslides are more likely to occur on barren soil [119] due to the weaker soil strength, which leads to instability. The different LULC types observed in this region are water bodies, forests, agricultural land, built-up areas, and wasteland (Figure 5a). It was found that 35 (47.94%) landslides occurred on wasteland, followed by 30 (41.09%) landslides on agricultural land. Only five (6.84%) landslides occurred in the forest, which substantiates the study by Shu et al. [117]. The Kappa coefficient for the LULC classification is 0.818, which shows strong agreement (

Table 9. Theoretical error matrix of LULC classification.

		LULC classes				Total
		Forest	Built-Up Area	Agricultural Land	Wasteland
GPS points	Forest	5	0	0	1	6
	Built-up area	0	3	0	1	4
	Agricultural land	2	1	82	0	85
	Wasteland	2	1	0	12	15
Total		9	5	82	14	110

and

Table 10. Cohen’s Kappa measures.

	Value	Asymptotic Standard Error a	Approximate T b	Approximate Significance
Measure of agreement–Kappa	0.818	0.057	12.225	0.000
No. of valid cases	110	-	-	-

^a Not assuming the null hypothesis. ^b Employing the asymptotic standard error assuming the null hypothesis.

).

3.2.4. Geomorphology

Wang et al. [120] describe geomorphology as a significant factor that influences the initiation of landslides. The geomorphic classes in this area include water bodies, floodplains, plateaus, and denudational hills (Figure 5b). The elevated areas of the study area are occupied by denudational hills. All the 73 landslides reported are confined to denudational hills.

3.2.5. Road Buffer

The excavation of slopes for road construction and widening is common in mountainous terrain, and it has a significant impact on slope stability [121,122]. A total of nine landslides (12.32%) were recorded within the 100-m buffer distance of roads (Figure 5c). This shows that road cuttings did not have much of an effect on the number of landslides that occurred.

3.2.6. Normalized Difference Road Landslide Index (NDRLI)

Road-induced landslides are extracted using the NDRLI [80]. Higher NDRLI values imply sites of unscientifically modified slopes for road cuttings, which might facilitate landslides and make them more vulnerable to future landslides. The NDRLI values of this area range from −0.25 to 0.50 (Figure 5d), which is grouped into five classes (−0.25 to 0.04, 0.04 to 0.09, 0.09 to 0.13, 0.13 to 0.18, and 0.18 to 0.50). According to the analysis, 29 (39.72%) landslides were recorded in the class with an NDRLI value range of 0.18–0.50, followed by 21 (28.76%) landslides occurring in the class with a value range of 0.13–0.18. This confirms that the NDRLI has a significant influence on landslide occurrence.

3.2.7. Normalized Difference Water Index (NDWI)

Water bodies are extracted using the NDWI [123,124]. The NDWI scale ranges from −1 to 1 [125], with values ranging from 0.2 to 1 representing water surface, 0.0 to 0.2 representing flooded areas or areas with higher humidity, −0.3 to 0.0 representing moderate drought, non-aqueous surfaces, and −1.0 to −0.3 representing drought, non-aqueous surfaces [126]. Landslides are more likely in places with NDWI values between 0.0 and 0.2, as the soil moisture in this area is higher. As pore pressure increases and soil shear strength decreases, soil moisture plays a crucial role in causing landslides [127]. This area’s NDWI values range from −0.49 to 0.19 (Figure 6a) and are classified into five classes (−0.49–−0.22, −0.22–−0.17, −0.17–−0.14, −0.14–−0.09, and −0.09–0.19). From the analysis, it is found that only five (6.84%) landslides occurred in areas with positive NDWI values. This confirmed that the NDWI does not have much influence on the initiation of landslides.

3.2.8. Normalized Burnt Ratio (NBR)

The NBR is a factor for identifying burned areas and estimating burn severity [128]. Higher NBR values indicate burnt areas, which are more vulnerable to landslides because fire lowers soil infiltration capacity and increases soil erodibility by obliterating plant and ground cover, making the area more prone to landslides [129]. The NBR values of this area range from −0.56 to 0.56 (Figure 6b), which are grouped into five classes (−0.56 to 0.14, 0.14 to 0.25, 0.25 to 0.33, 0.33 to 0.40, and 0.40 to 0.56). From the spatial analysis, it is proven that five (6.84%) landslides were recorded in the class with an NBR value range of 0.40–0.56, followed by 31 (42.46%) landslides in the class with a value range of 0.33–0.40, and 22 (27.40%) landslides in the class with a value range of 0.25–0.33, respectively. Thus, it is proved that the NBR has only a considerable role in the initiation of landslides.

3.2.9. Soil-Adjusted Vegetation Index (SAVI)

The SAVI values vary from −1.0 to 1.0, with low values indicating a lack of green vegetation [130]. This index is similar to the NDVI and is appropriate for places with little vegetation cover and open soil surfaces [130,131]. The probability of landslides is high in places with lower SAVI values. The SAVI values of this area range from −0.09 to 0.88 (Figure 6c), which is grouped into five classes (−0.09 to 0.37, 0.37 to 0.48, 0.48 to 0.57, 0.57 to 0.65, and 0.65 to 0.88). The spatial analysis found that 41 (56.16%) landslides were recorded in areas with SAVI values ranging from 0.57 to 0.88, and only 14 (19.17%) landslides were reported in areas with SAVI values between −0.09 and 0.48. This confirmed that the SAVI does not have much influence on the initiation of landslides.

3.2.10. Lithology

Rock strength and permeability are determined by lithological units [132]. Landslides often happen in rock strata with low shear strength and permeability [114]. The rock types observed in this region are charnockite, granite, and hornblende gneiss (Figure 6d). All 73 landslides occurred on charnockite rocks. Though all the landslides occurred in areas with the charnockite rock type, this factor cannot be confirmed as a landslide-initiating factor in the study area. This is because a major portion (99.5%) of this region is occupied by this rock type. Moreover, this is a massive, high-grade Precambrian metamorphic rock, which has little influence on landslides. This indicates that lithology does not have much impact on slope failures.

3.3. Landslide Susceptible Zones

The landslide susceptibility is high in the eastern portion of this region. The susceptibility maps are depicted in Figure 7a–d. The area and percentage of different susceptible zones are depicted in

Table 11. Area and percentage of susceptible zones.

	AHP		F-AHP		GIS-TISSA		NCESS
Susceptible Zones	Area in sq. km.	Percentage of Area	Area in sq. km.	Percentage of Area	Area in sq. km.	Percentage of Area	Area in sq. km.	Percentage of Area
Low	320.34	43.42	317.32	43.01	411.67	55.80	50.54	6.85
Moderate	280.35	38.00	282.71	38.32	82.19	11.14	133.68	18.12
High	137.08	18.58	137.74	18.67	111.11	15.06	61.61	8.35
Total	737.77	100	737.77	100	604.97 *	82 *	245.83 **	33.32 **

* 18% of the study area belongs to the no-susceptible zone. ** 66.68% of the study area belongs to the no-susceptible zone.

. A total of 18.58%, 18.67%, 15.06%, and 8.35% of the study area belong to highly susceptible zones for the AHP, F-AHP, GIS-TISSA, and NCESS models, respectively. It was found that the highly susceptible zone has steeper slopes (22.48–74.47°), clayey soil, LULC types such as wasteland and agricultural land, denudational hills, and NDRLI values between 0.13 and 0.50. This study confirmed that the NBR, NDWI, SAVI, and road buffers do not have much influence on landslide initiation. This ascertained that the landslide occurrence is due to natural as well as human-induced reasons. Koottickal, Poonjar-Thekkekara, Moonnilavu, Thalanad, and Koruthodu are the LSGs that are highly susceptible to landslides in the Kottayam district.

3.4. Validation of the Maps

The AUC values for the AHP, F-AHP, TISSA, and NCESS models for the training dataset are 0.838 (83.8%), 0.849 (84.9%), 0.868 (86.8%), and 0.759 (75.9%), whereas the AUC values for the validation dataset are 0.867 (86.7%), 0.872 (87.2%), 0.889 (88.9%), and 0.789 (78.9%)(Figure 8 and Figure 9). From the results, it is clear that the AHP, F-AHP, and TISSA models have an excellent AUC value (AUC value above 0.80), whereas the NCESS model has an acceptable AUC value (AUC value above 0.70). When comparing all these models, the TISSA model is found as the most efficient, followed by the F-AHP, AHP, and NCESS models.

Other validation matrices such as accuracy, MAE, and RMSE also provided better results (higher accuracy and lower MAE and RMSE values) for the TISSA model, followed by the F-AHP, AHP, and NCESS models. The accuracy, MAE, and RMSE of the TISSA model for the training dataset are 0.863, 0.274, and 0.196, respectively, whereas, for the validation dataset, the values are 0.869, 0.226, and 0.122, respectively (

Table 12. Validation matrices and values of the training and validation datasets.

Matrices	Training Dataset				Validation Dataset
	AHP	F-AHP	TISSA	NCESS	AHP	F-AHP	TISSA	NCESS
Accuracy	0.831	0.842	0.863	0.752	0.855	0.856	0.869	0.770
MAE	0.292	0.278	0.274	0.422	0.249	0.243	0.226	0.309
RMSE	0.218	0.214	0.196	0.298	0.159	0.147	0.122	0.177

).

4. Discussion

This modeling study observed that the highly susceptible zone has slopes that are steeper (22.48–74.47°), clayey soil, wasteland and agricultural land, denudational hills, and NDRLI values between 0.13 and 0.50. The study by Khan et al. [133] found 89.76% of the landslide occurrences on slopes above 20°. The research by Oh et al. [89] observed a sharp increase in the number of landslides for slopes greater than 26°, whereas Getachew and Meten [119] found this trend for slopes above 25°. The landslide susceptibility modeling by Thomas et al. [17] observed 91.35% of landslides on clayey soil, whereas Babitha et al. [72] found 65.39% of landslides on the same type of soil. The study by Ajin et al. [58] also found higher incidences on clayey soil. The increase in water content in clay resulting from the percolation of rainfall can reduce its stability (increase stress and decrease strength) [134,135]. Similar to the findings of this study, Khan et al. [133] found higher landslide incidences (66.51%) on barren land, and Thomas et al. [17] observed higher incidences (56.19%) on plantation areas. Abedin et al. [118] in their study observed that massive landslides occurred mostly in plantation agricultural fields as the trees were not able to prevent erosion. Ajin et al. [10] also found that the replacement of vegetation with plantations is one of the reasons for landslide initiation. Denudational hills are the result of erosional and weathering processes [136], and this geomorphic feature is susceptible to landslides [137,138]. Vijith et al. [139] in their study identified the majority of landslides (49.49%) on denudational hills. The factor that represents road-induced landslides, the NDRLI, shows a significant influence in the study area. Though this factor was not much used in landslide modeling studies, road density or road buffers have been employed by many researchers [15,72]. However, the road buffer does not show much influence since it is a single-class factor and the majority of the sliding occurred beyond this limit. The significant influence of the NDRLI confirmed that modification of slopes for road widening or construction is one of the human-induced reasons for landslides. Akshaya et al. [15] in their research on susceptibility modeling found a majority of sliding (69.56%) along road cuttings. The landslide susceptibility zonation conducted by Ajin et al. [132] in the Meenachil and Kanjirappally taluks of the Kottayam district found slope and elevation as the major influencing factors, followed by drainage density, road density, and LULC. Even though several studies have mentioned the strong link between landslides and factors such as NBR, NDWI, and SAVI [127,140,141,142], this study did not find a strong link between landslide occurrences and the influence of factors such as the NBR, NDWI, and SAVI.

This study found TISSA as the best model, followed by the F-AHP, AHP, and NCESS models. The AUC analysis for the validation dataset confirmed that the TISSA model has the highest performance (an AUC value of 88.9%), which shows a 1.7% increase over the F-AHP model, a 2.2% increase over the AHP model, and a 10.0% increase over the NCESS model. For the training dataset, the AUC value of the TISSA model confirmed better performance with a 10.9% increase over the NCESS model, a 3.0% increase over the AHP model, and a 1.9% increase over the F-AHP model. The performance of one or more models can be compared using the ROC curve, which can also be used to evaluate a model’s overall performance [143]. The accuracy indices generated with the ROC analysis are not affected by fluctuations brought on by the application of arbitrary decision criteria or cut-offs, which is the technique’s most desirable characteristic [144]. This study also found that the TISSA model has the highest accuracy value of 0.869 for the validation dataset, which confirmed an increase of 0.099 over the NCESS model’s accuracy value. The accuracy value highlights how many predictions were correct. In addition, the MAE and RMSE values of the TISSA model confirmed a decrease in error values of 0.083 and 0.055 over the NCESS model. The accuracy, MAE, and RMSE values also reiterate that the TISSA model is the best, followed by the F-AHP, AHP, and NCESS models. With the units of the original data, the MAE and RMSE both express the average magnitude of the prediction error [145]. While the MAE is unambiguous [146], the RMSE has the advantage of penalizing large errors more strongly than the MAE [145]. The studies by Escobar-Wolf et al. [11] and Pikul [101] both found that the GIS-TISSA model works well for identifying areas that are likely to slide. Of the two MCDM models (AHP and F-AHP) selected for this study, the F-AHP model is found better than the AHP model. The studies by Akshaya et al. [15], Meshram et al. [147], Tripathi et al. [148], and Vilasan and Kapse [149] also found the F-AHP model to be more efficient than the AHP model.

Furthermore, except for the NCESS model, all other models have excellent AUC values. The reason for the lower value is that the NCESS susceptibility map is outdated. The NCESS map is 12 years older and is still in use. van Westen et al. [150] assert that every landslide susceptibility map should be continuously updated; the update frequency may range from one to ten years, depending on the dynamics of land use change, since land use determines both the elements at risk and the occurrence of new slides. Considering the high population density in Kerala and the intense land use change [151,152,153,154], this map needs an immediate update to ensure that it is capable of addressing the larger number of landslides that occurred as a result of an extreme event such as the 2018 Kerala floods.

Due to the lack of an optimal number of rain gauges or automatic weather stations, rainfall, the sole triggering factor, is not considered in this model. The global precipitation datasets and India Meteorological Department (IMD) gridded rainfall data have a low spatial resolution, and hence, their integration may result in erroneous outputs. In an analysis of global precipitation datasets, Sun et al. [155] observed considerable disparities in the variability and magnitude of precipitation estimates. Similarly, Prakash et al. [156] analyzed four precipitation products and reported that all of them underestimate precipitation significantly. They found that the number and spatial coverage of rain gauge stations, data assimilation models, and satellite algorithms all limit the reliability of precipitation datasets. In addition, the inconsistencies stated to limit the product’s ability to monitor climate change and validate models. Jena et al. [157] investigated the effectiveness of IMD gridded data (0.25° × 0.25°) and satellite estimates for identifying cloudburst occurrences over the northwest Himalayas, finding that only 6 out of 18 events could be recognized using the IMD data when they compared different precipitation datasets (including IMD gridded data).

5. Conclusions

Landslides are a common hazard in Kerala’s Western Ghats region. The Ghat region of the Kottayam district has also witnessed many landslides in the past, of which the Koottickal disasters of 2021 were the most devastating. The present study demarcated the susceptibility of high-range LSGs in the Kottayam district into three zones using the AHP and F-AHP models and compared the results with the existing GIS-TISSA and NCESS susceptibility maps. The spatial analysis revealed that three maps (TISSA, F-AHP, and AHP) have excellent performance with AUC values above 0.80 (80%). Though the difference in prediction accuracy is negligible, of these three maps, the one produced using the TISSA model yields the best result. However, the NCESS map has only acceptable performance (an AUC value of 0.789). This is because the NCESS map is outdated (published in 2010), and hence, the change in land use/land cover types for the last decade was not revealed. According to the TISSA map, a total of 15.06% of this region is highly landslide susceptible. Slope, soil texture, LULC types, geomorphology, and NDRLI were found to be the most influential factors, confirming the influence of both natural as well as anthropogenic factors on landslide occurrence. Though rainfall distribution (the triggering factor) is crucial for landslide susceptibility modeling, it is not utilized in this research due to a lack of an adequate number of rain gauges and automated weather stations. In addition, the IMD gridded rainfall data and other satellite data are of poor spatial resolution, which is not suitable for a smaller area.

The created maps will be extremely useful to agencies dealing with landslide hazards and implementing effective mitigation measures such as landslide early warning systems, disaster risk reduction guidelines, and sustainable construction practices. This study also proposes to develop effective guidelines for agricultural practices on hilly terrains, including crop selection, irrigation methods, and harvesting methods. Hydroseeding of vegetative species (such as vetiver grass or lemon grass), which can stabilize the slope, should be planted, and nature-based solutions should be promoted in high- and very-high susceptible zones. This can minimize the impact of landslides.