Landslide Susceptibility Mapping in Idukki, Western Ghats, India by Integrating Weighted Multi-Criteria Decision Analysis and Machine Learning

Landslide susceptibility mapping (LSM) plays a critical role in managing natural disaster risks. As climate change leads to more frequent and severe landslides, addressing this issue is the need of the hour. Geospatial techniques need to be applied in conjunction with geoenvironmental factors to develop effective disaster management strategies. Assessing the susceptibility of landscapes to landslides has changed significantly with the advent of machine learning (ML), Geographic Information Systems, and high-resolution geospatial datasets, moving away from traditional statistical methodologies [1]. Modern methodologies effectively identify and model the complex relationships and interdependencies among parameters related to landslides, leading to more robust and accurate predictive models. The advanced approaches improve traditional models by combining remote sensing, ML techniques, and neural networks [2]. This change enhances the landslide risk assessments in regions with complex geo-climatic conditions.

Landslides are among the most catastrophic of natural hazards that cause the loss of thousands of lives and billions of dollars’ worth of infrastructure across the world in mountainous regions every year [3]. In India, the Western Ghats are a geomorphic hotspot, particularly the Idukki district of Kerala, which has a significant history of catastrophic slope failures. The Western Ghats consist of steep gradients, intensely weathered lateritic soils, and high-magnitude monsoonal rainfall. The monsoons are responsible for the Pettimudi landslides, which resulted in the loss of more than 60 lives in Idukki in the year 2020. Since Idukki district is socio-economically vulnerable to such extreme monsoons [4], there is a strong need for an accurate yet operationally interpretable LSM to aid proactive risk management [5]. The preliminary insights, such as bivariate statistics and heuristic methods, have provided early and valuable understanding [6]; however, they often fall short in capturing the complex, non-linear relationships inherent in landslide genesis.

With the rapid development of geoinformatics along with cutting-edge ML techniques, the possibilities of synthesis and evaluation of large, heterogeneous, and multi-source geospatial information sets are unprecedented [7, 8]. Nonetheless, the application of holistic frameworks that utilize sophisticated ML approaches, compare them with expert approaches, and spatially validate and explain them, remains scarce [9, 10].

To address the identified gaps in the related works, this study is guided by the following questions:

RQ1: How do the rankings and contributions of landslide conditioning factors differ between knowledge-driven (Weighted Multi-Criteria Decision Analysis (WMCDA)) and ML methods?

RQ2: How do advanced ML classifiers, such as multilayer perceptron (MLP), compare with traditional WMCDA in terms of predictive accuracy and uncertainty reduction?

RQ3: How can density-based spatial validation enhance the reliability of hazard zonation for regional planning?

As such, the following hypotheses were proposed:

H1: The predictive performance of the MLP model will outperform that of the expert-weighted index and other baseline classifiers.

H2: While natural factors like slope and rainfall are critical, anthropogenic factors (such as proximity to roads) will emerge as primary drivers in the ML models due to the region's human-induced landscape modifications.

H3: The final hazard zonation map will confirm its spatial validity, as the high and very high susceptibility zones will contain a large concentration of historical landslide events.

This study presents three main contributions:

(1) This study presents a systematic comparison between expert-driven WMCDA and six ML classifiers using the six primary conditioning factors and a validation framework, with feature-importance analysis to support interpretability.

(2) Unlike conventional studies that majorly use accuracy metrics, this framework integrates probability calibration, bootstrap-based uncertainty estimation, and density-based spatial validation of hazard zones.

(3) An interpretable framework integrating driver ranking, rainfall-slope triggering thresholds, and hazard zonation is presented.

A high-resolution LSM for the Idukki district was generated using an integrated geospatial framework. Figures 1 (a)-(d) present Kerala within India, Idukki district within Kerala, the spatial distribution of recorded landslide events, and terrain slope classification, respectively.

The workflow comprises five components as illustrated in Figure 2, which presents the general methodology: (a) data ingestion, (b) feature engineering, (c) susceptibility modeling, (d) validation, and (e) decision support.

Figure 1. Geographical setting and input data for the study are: (a) Location of Kerala within India; (b) Location of the Idukki district in Kerala; (c) Landslide inventory comprising 40 documented events, shown relative to the drainage network; (d) Slopes of the Idukki district divided into five slope classes

Figure 2. The five steps involved in the integrated landslide susceptibility mapping (LSM) framework: (a) Data ingestion; (b) Data pre-processing; (c) Core modeling; (d) Validation and analytics; (e) Decision support

3.1 Geospatial database and conditioning factors

An integrated geospatial database comprises the processed seven factors, such as 30 m resolution Shuttle Radar Topography Mission Digital Elevation Model (SRTM DEM), mean monsoon rainfall (2015–2020), lithology, terrain curvature, Land Use/Land Cover (LULC), distance to roads, and distance to waterways and landslide points. The rainfall layer represents spatial patterns of mean monsoon rainfall (2015–2020) and is used to characterize long-term hydro-climatic predisposition to landsliding rather than event-scale triggering. The model was trained and validated against a historical landslide inventory of 40 distinct locations compiled from State Disaster Management Authority reports.

3.2 Feature engineering and data harmonization

To maintain spatial coherence, all the acquired data were projected onto a consistent UTM Zone 43N (EPSG:32643) reference system and resampled to a uniform 30m resolution. We performed a feature engineering pipeline using Rasterio, GeoPandas, SciPy, and NumPy in Python to convert raw input data into a format ready for analysis. We then turned vector data, such as lithology, roads, and waterways, into aligned raster grids. We then extracted key terrain attributes, slope, aspect, Topographic Wetness Index (TWI), and curvature from the DEM. All the primary and derived conditioning factors used for analysis are slated in Table 1. Finally, we performed an Euclidean-based distance transform analysis to produce continuous proximity surfaces for roads and waterways. This process resulted in a consolidated multi-band feature dataset for analysis and modeling.

Table 1. Summary of conditioning factors in a multi-source geospatial database

For integrating continuous factors, the susceptibility model required normalizing the layers onto a common [0, 1] scale. Before scaling, data were subjected to Winsorization (clipping at the 5th and 95th percentiles) in an attempt to reduce the impact of extreme outliers. The normalization equation is shown in Eq. (1).

$X_{\text {norm}}=\frac{X-X_{\min}}{X_{\max }-X_{\min}}$     (1)

where, X denotes the original Winsorized raster value, X_min and X_max represent the minimum and maximum values after Winsorization, and X_norm is the normalized value used in subsequent susceptibility modeling.

A hybrid heuristic-statistical approach was adopted to assign factor importance. For the lithology factor, the weights have been quantitatively derived by applying the FR method, whose values have been computed based on the spatial density of historical landslides. On the other hand, the remaining continuous factors have been assigned heuristic weights according to their geomorphological relevance towards slope stability.

The resulting spatial distribution of these harmonized conditioning factors is visualized in the multi-panel plot (Figure 3). The figure details the region's anthropogenic footprint (Figure 3(a)) and geological setting (Figure 3(b)), while juxtaposing the primary hydrological trigger, monsoon intensity (Figure 3(c)), against the geomorphological susceptibility defined by slope gradients (Figure 3(d)).

Figure 3. Spatial distribution of historical landslides (black dots) and key conditioning factors for landslide susceptibility in Idukki district: (a) Terrain relief and infrastructure; (b) Lithological units; (c) Mean monsoon rainfall; (d) Slope classification

Out of the ten conditioning factors, elevation, rainfall, Lithology, Distance to roads, Distance to waterways, and LULC are the primary conditioning factors. Terrain derivatives, such as slope, aspect, curvature, and TWI, were derived from DEM and treated as independent predictors. Aspect was excluded from the WMCDA framework because its directional and circular nature does not satisfy the monotonicity assumptions required for weighted overlay analysis. However, it was included in ML models, which can capture non-linear and interaction-driven effects without imposing directional bias.

3.3 Modeling approaches

Two contrasting modeling approaches were analysed to select the most suitable for this context: a knowledge-driven Weighted Overlay Model and a data-driven ML Benchmark.

3.3.1 Knowledge-driven model: Weighted overlay

The first approach employed was WMCDA, which combines spatial data considering weights. In order to calculate the Critical Success Index (CSI), we used the normalized conditioning factors, C_i, and expert-assigned weights, W_i, as described in Eq (2).

$C S I=\sum_{i=1}^n W_i C_i$     (2)

3.3.2 Data-driven model: Machine learning benchmark

In addition to the heuristic approach, a rigorous benchmarking analysis was conducted with six ML classifiers: RF, LR, XGBoost, SVM, k-nearest neighbors (KNN), and MLP. A training dataset was constructed using six primary conditioning factors, with additional terrain attributes (slope, aspect, and curvature) derived from DEM and used where explicitly stated for exploratory and statistical analyses. The target variable (1 = landslide, 0 = non-landslide) was derived from the 40 historical inventory points and a randomly sampled set of non-landslide points (1:1 ratio). The preprocessing pipeline ensured data robustness by: (1) outlier mitigation using Isolation Forests, (2) feature scaling using RobustScaler, and (3) class balancing using SMOTE. Class imbalance was addressed only within the training set, after the train and test dataset split, to stabilize model learning under limited samples. No synthetic samples were introduced into the test set.

The dataset was partitioned into 70% for training and 30% for a hold-out test set. All models were calibrated using 5-fold cross-validation. The unseen test set was evaluated using metrics such as, Receiver Operating Characteristic - Area Under the Curve (ROC-AUC), Precision-Recall - Area Under the Curve (PR-AUC), Accuracy, and the Brier score (for calibration reliability). To ascertain statistical rigor, 95% confidence intervals (CI) were calculated using a 1,000-iteration bootstrap resampling technique.

3.4 Statistical analysis and feature importance

We applied permutation-based feature importance analysis to the RF model to rank the drivers of slope instability. This method evaluates the degradation in model predictive accuracy when a single feature’s values are randomly shuffled, thereby isolating its predictive contribution relative to other variables.

Also, using appropriate statistical significance, the key drivers were evaluated. Continuous variables (distance to infrastructure and terrain curvature) at the site of the landslide were compared with a random background distribution using the non-parametric Mann-Whitney U test. In the case of categorical variables (lithology), landslide density was determined, and statistically significant differences in susceptibility were assessed using 95% CI of the Wilson score interval for lithologic types.

3.5 Hazard zonation and validation

For practical planning, the final output raster was converted using quantile classification to divide the continuous susceptibility probability index into five Hazard Zones, namely, Very Low, Low, Moderate, High, and Very High. The predictive capability of the zonation map was validated quantitatively by overlaying the historical landslide inventory and calculating the frequency density of events within the high-susceptibility zones (High and Very High).

We performed advanced statistical analyses to validate the model’s performance and quantify the importance of conditioning factors associated with landslide occurrence (Figure 6).

Figure 6. Statistical and machine learning (ML)-based assessment of landslide hazard: (a) Random forests (RF) feature importance; (b) Lithological susceptibility ranking with confidence intervals; (c) Landslide triggering thresholds; (d) Infrastructure distance analysis; (e) Terrain curvature analysis; (f) Projected climate change impact on landslide risk

5.1 Quantitative hierarchy of landslide drivers

We applied an RF model to rank the significance of the conditioning variables specific to Idukki (Figure 6(a)). As per the analysis, proximity to roads is ranked as the most important factor (0.178). It surpasses the influence of rainfall (0.092) and distance to waterways (0.086). Notably, human-made infrastructure appears to play a more decisive role in slope stability than inherent natural features like slope gradient or curvature.

The prominence of infrastructure-related factors is further supported by non-parametric statistical testing. Box plots and Mann–Whitney U tests (Figure 6(d)) demonstrate that landslide locations are significantly closer to both roads and waterways than randomly sampled background locations (p < 0.001). The results highlight a strong spatial association between infrastructure development and landslide occurrence.

Lithological susceptibility analysis with 95% bootstrap confidence intervals (Figure 6(b)) identifies the Migmatite Gneissic Complex (0.80 landslide density per 10,000 pixels) as the most susceptible lithological unit. This is followed by the Charnockite Gneissic Complex, which shows a moderate landslide density (0.52 landslides per 10,000 pixels) accompanied by wider confidence intervals. On the flip side, the Acid Intrusive rocks are basically solid and barely show any landslides at all.

Beyond lithological controls, the contribution of secondary terrain attributes was evaluated. The probability density analysis of terrain curvature (Figure 6(e)) does not reveal any significant difference (p = 0.1126) between curvature values at landslide sites (mean = 1.366 ± 0.648) and background locations (mean = 0.167 ± 0.0957). This suggests that at this regional scale, broad slope characteristics are more dominant than fine-scale topographic convexities or concavities in controlling landslide occurrence.

5.2 Triggering thresholds and future risk projections

To facilitate the conceptual design of rainfall-triggered EWSs, we examined hydro-geomorphological conditions associated with slope failure, rainfall-slope analysis (Figure 6(c)). A bivariate analysis utilizing an ROC-optimized diagnostic approach was implemented to define the intersection of precipitation and topography. According to the results, a mean monsoon rainfall of 392.5 mm (90th percentile of failure conditions), when combined with a slope gradient of 24.2° (80th percentile), serves as a diagnostic benchmark for identifying climatically predisposed failure conditions (ROC AUC = 0.590).

The potential impact of climate change on future landslide risk was explored through multi-threshold simulations (Figure 6(f)). Scenario-based projections corresponding to increased rainfall intensities, Representative Concentration Pathways (RCP) 4.5 [+15%], RCP 6.0 [+30%], RCP 8.5 [+50%], and extreme-event [+100%] scenarios) indicate a systematic expansion of areas classified as High and Very High susceptibility. Even under moderate rainfall increases, high-risk zones expand appreciably, while more extreme scenarios lead to substantial growth of very high-risk areas.

5.3 Model performance and rigorous validation

The predictive power of the susceptibility framework was rigorously evaluated using both threshold-independent metrics and density-based validation to ensure the reliability of the final hazard map (Figure 7).

Figure 7. Model performance and hazard zone validation

5.3.1 Predictive discriminative power

The model's ability to distinguish between landslide-prone and stable terrain was quantified using the ROC and Precision-Recall (PR) curves (Figure 7(a)).

On the held-out test set, the model achieved an ROC–AUC of 0.845, indicating strong discriminative capability.

Given the inherent class imbalance in landslide inventories, the PR curve provides a more stringent test of performance. The high Average Precision (AP = 0.829) confirms that the model minimizes false positives while successfully retrieving the majority of landslide events, a critical requirement for operational EWSs.

5.3.2 Hazard zone reliability

To validate the practical utility of the zonation map, a landslide density analysis was conducted (Figure 7(b)). This analysis correlates the classified hazard zones (x-axis) with the observed density of historical landslides (y-axis).

The results reveal a striking positive correlation (R² = 0.98), where landslide density increases monotonically from the Very Low to Very High zones.

The Very High hazard zone exhibits the highest landslide density (approx. 0.35 landslides per 10,000 pixels), which is significantly distinct from the lower risk zones as indicated by 95% confidence intervals. This confirms the effectiveness of the model in stratifying the landscape into operationally meaningful risk classes, enabling targeted mitigation and resource prioritization.

To identify the most reliable predictive architecture for LSM, we conducted a rigorous comparative benchmarking analysis across six ML classifiers. The evaluation was done using isotonic calibration and 1,000-iteration bootstrap resampling on a held-out test dataset.

7.1 Discriminative power and model hierarchy

The comparative benchmarking analysis established a distinct hierarchy of model suitability, revealing that non-linear algorithms significantly outperform linear counterparts in capturing the complex geospatial drivers of landslides.

The MLP emerged as the dominant classifier, achieving the highest discriminative power with a mean ROC-AUC of 0.904 (95% CI: 0.804 - 0.971) and a PR-AUC of 0.917 (95% CI: 0.806 - 0.977). The RF model also demonstrated robust performance, yielding an ROC-AUC of 0.859 (95% CI: 0.752-0.945). XGBoost demonstrated competitive performance but exhibited slightly higher variance in the bootstrap confidence intervals compared to the MLP and RF models. Simpler models, such as LR (ROC-AUC = 0.666) and SVM, showed markedly lower performance. This confirms that kernel complexity alone does not guarantee improved generalization under limited and imbalanced landslide inventories.

7.2 Probabilistic calibration and reliability

For operational risk management, the reliability of the predicted probability is as critical as classification accuracy. Accordingly, all classifiers were subjected to post-hoc isotonic calibration using five-fold internal cross-validation, and their probabilistic behavior was evaluated using reliability (calibration) curves and Brier scores (Figure 9).

The MLP and RF models demonstrated excellent calibration, with curves that closely track the perfect calibration diagonal. This indicates high probabilistic coherence; for instance, when the MLP model predicts an 80% probability of a landslide, the observed event frequency is approximately 80%. This visual assessment is quantitatively supported by the MLP achieving the lowest Brier score (0.137).

In contrast, models such as SVM and KNN exhibited significant deviation from the diagonal. This indicates a tendency toward overconfidence or underconfidence, rendering their raw probability outputs unreliable for precise hazard zonation.

7.3 Strategic model selection

Based on the intersection of high discriminative accuracy (ROC-AUC and PR-AUC), robustness under bootstrap resampling, and reliable probabilistic calibration, the MLP model was selected as the primary predictive engine for subsequent LSM and operational analyses. This data-driven selection ensures that the final LSM and subsequent operational analyses are grounded in the most robust predictive framework identified.