Landslide Susceptibility Zonation and Aggregation-Scale Sensitivity Across Sub-Districts of the Kumaon Himalaya: A Comparative Analysis Using AHP and Fuzzy-AHP Framework

1. Introduction

Landslides are a pervasive geohazard in mountainous regions worldwide, where steep topography, tectonic activity, and intense precipitation converge to destabilize slopes [1,2,3]. Globally, these events claim thousands of lives annually and impose substantial economic burdens on vulnerable communities, with fatalities disproportionately concentrated in tectonically active, densely populated mountain belts across Asia [2,3]. In the Himalayas, these events frequently disrupt transportation networks, damage settlements, and exacerbate vulnerability in densely populated valleys. The Kumaon sector of the Uttarakhand Himalaya exemplifies this risk, with rainfall-triggered landslides dominating owing to the region's monsoon-driven hydrology and fragile geology. Historical records document recurrent failures along highways and in sub-district-level administrative units [4,5,6], underscoring the need for spatially explicit susceptibility assessments to inform disaster risk reduction and land-use planning. Such mapping delineates zones of potential instability, enabling targeted interventions, such as slope stabilization and restricted development [7], thereby enhancing resilience in infrastructure corridors and rural habitations [1,4,8].

Contemporary approaches to landslide susceptibility mapping (LSM) encompass statistical, machine learning, knowledge-driven, and hybrid methodologies, each tailored to data availability and environmental complexity [9]. Data-driven models, including logistic regression and random forests, excel in regions with comprehensive inventories but fail in data-sparse terrains, where training samples are limited or biased toward recent events [9,10]. Conversely, knowledge-driven multi-criteria decision-making (MCDM) techniques, such as the analytical hierarchy process (AHP), integrate expert judgments with geospatial layers to derive susceptibility indices without relying on exhaustive historical data [10,11,12,13] These methods remain relevant in mountainous settings, such as the Himalayas, where sparse monitoring and heterogeneous factors, ranging from lithology to proximity to linear features, necessitate interpretive weighting [5]. Recent applications have affirmed the efficacy of MCDM for regional-scale mapping, particularly when conditioning factors are reclassified into ordinal susceptibility classes [14,15,16,17].

Central to MCDM is the AHP, which structures pairwise comparisons into hierarchical weights that reflect the relative importance of factors. The AHP has been extensively applied in LSM, facilitating transparent integration of slope gradient, rainfall intensity, and land cover into susceptibility zonations [11,13,15]. However, its reliance on crisp numerical judgments introduces limitations: subjective expert elicitations can amplify bias, whereas crisp scales inadequately capture the vagueness inherent in natural thresholds [18,19,20], such as the transition from stable to unstable slopes. Uncertainty propagation through matrix computations further compromises precision, particularly in expert-driven scenarios with inconsistent pairwise ratings [20]. To address these constraints, the fuzzy analytical hierarchy process (FAHP) extends the AHP by incorporating fuzzy set theory, typically through triangular fuzzy numbers, to model linguistic ambiguity in comparisons. This formulation yields interval weights that systematically propagate uncertainty, enhancing robustness without requiring probabilistic inventories. In geohazard contexts, FAHP has demonstrated superior handling of imprecise inputs, producing susceptibility maps with refined discrimination in boundary zones [16,19,21,22]. Comparative studies in analogous terrains have reported marginal gains in predictive skill, attributed to the alignment of fuzzy aggregation with real-world gradations in factor influence [23,24].

While the adoption of administrative units as the spatial scale of susceptibility reporting has gained traction in policy-relevant hazard mapping [5,6,25,26], a critical and largely unexamined methodological question persists: how sensitive are sub-district-level susceptibility classifications to the specific pixel-to-polygon aggregation statistic employed, and to what extent does the choice of aggregation method mask localized high-susceptibility zones within otherwise moderate-susceptibility administrative units? To the authors' knowledge, no study within the Himalayan susceptibility literature has systematically and quantitatively evaluated multiple zonal aggregation approaches at the sub-district scale using inventory-validated performance metrics, inter-method agreement statistics, and an explicit class-divergence framework; prior sub-district and district-level studies in analogous settings have uniformly adopted mean-based aggregation without assessing whether this choice introduces directional risk underestimation [6,26,27]. This gap has direct practical consequences for disaster governance: an aggregation statistic that suppresses the upper tail of the within-unit susceptibility distribution may lead administrators to systematically underestimate localized risk across a majority of sub-districts, with immediate implications for infrastructure planning, emergency preparedness, and the operationalization of local-level risk disaggregation mandated under contemporary disaster risk reduction frameworks.

Against this background, the present study pursues three explicit objectives: (i) to comparatively evaluate AHP and FAHP for sub-district-level landslide susceptibility zonation across all 53 sub-districts of the Kumaon Himalaya at 30 m resolution, quantifying the effect of fuzzy uncertainty incorporation on model predictive performance against a validated landslide inventory; (ii) to systematically assess the sensitivity of sub-district susceptibility classifications to three alternative pixel-to-polygon aggregation statistics (Zonal Mean, Zonal P90, and Zonal Maximum) using inventory-validated ROC–AUC and Spearman rank correlation metrics; and (iii) to quantify intra-sub-district risk masking through a class-divergence framework and identify priority administrative units for targeted reconnaissance and infrastructure regulation. The study's novel contributions include the first systematic, multi-metric aggregation-scale sensitivity analysis at the sub-district administrative scale within the Himalayan susceptibility literature, a governance-ready class-divergence diagnostic layer, and a spatially explicit, inventory-validated susceptibility framework tailored to the hydro-geomorphic complexity of the Kumaon Himalaya.

2. Materials and Methods

2.1. Study Area

The Kumaon Himalaya constitutes the northeastern division of Uttarakhand, India, situated between latitudes 28°44′–30°49′ N and longitudes 78°45′–81°01′ E (Figure 1). Spanning a geographical area of approximately 21,035 km², the region forms a critical physiographic transition zone bounded by Nepal to the east and the Garhwal Himalaya to the west. The topography is characterized by extreme altitudinal heterogeneity, rising abruptly from the Terai foothills (~190 m) to the Greater Himalayan peaks exceeding 7,000 m. This rugged terrain, shaped by active tectonics and deeply incised valleys, creates a landscape inherently prone to slope instability. Climatically, the region is governed by the Indian Summer Monsoon, where sharp elevation gradients induce intense orographic precipitation–a primary trigger for landslide events. Administratively, the study covers six districts: Almora, Bageshwar, Champawat, Nainital, Pithoragarh, and Udham Singh Nagar. To facilitate localised hazard governance, the analysis was conducted across the region's 53 sub-districts, utilizing them as the fundamental unit for susceptibility assessment.

2.2. Data Sources and Preprocessing

Nine landslide-conditioning factors were prepared as gridded layers at 30 m spatial resolution to enable pixel-wise integration. The datasets included a DEM-derived terrain representation (slope and TWI), hydrographic and infrastructure proximity layers (distance to streams and roads), thematic geology and structure (lithology and distance to faults), land cover, monsoon-season rainfall, and soil erosion susceptibility. In addition, a historical landslide inventory comprising 1,715 georeferenced events, compiled by the Geological Survey of India (GSI), was incorporated for model validation and aggregation-scale sensitivity analysis. All datasets and their primary sources are summarized in Table 1.

All vector datasets were converted to raster format, projected to a common reference system (UTM Zone 44N), and clipped to the Kumaon boundary. Continuous rasters were resampled to 30 m using bilinear interpolation, whereas categorical layers (e.g., LULC and lithology) were resampled using nearest-neighbor interpolation to preserve class integrity. Processing was implemented in QGIS 3.34, and model evaluation was performed using Python. The overall methodological framework adopted in this study is summarized in Figure 2, which illustrates the sequential workflow from data acquisition and factor standardization through weight derivation, LSI computation, zonal aggregation, and final validation.

2.3. Landslide Conditioning Factor Selection and Justification

Nine conditioning factors were selected based on their established relevance to rainfall-triggered landslide initiation in the Himalayan setting, spatial data availability at 30 m resolution, and collective representativeness across topographic, geological, hydro-meteorological, and anthropogenic causal pathways [4,6,16].

2.3.1. Topographic Factors

Slope gradient and the topographic wetness index (TWI) form the topographic foundation of the susceptibility framework. Slope directly governs gravitational shear stress on slope materials and has been consistently identified as the dominant predisposing factor in Himalayan landslide studies, where foliated lithologies frequently fail at angles exceeding 25°–35° [14,28]. TWI, computed as $\ln {\displaystyle \frac{\alpha }{\tan \beta }}$, quantifies terrain moisture accumulation by integrating the upslope contributing area and local gradient; high TWI values identify convergent hollows and valley floors where monsoon-driven pore pressure buildup most critically reduces slope stability [4,22].

Figure 2. Methodological flowchart of the landslide susceptibility assessment framework.

Figure 2. Methodological flowchart of the landslide susceptibility assessment framework.

2.3.2. Geological Factors

Lithology and distance to faults capture the structural and material controls on slope instability [4,14]. The Lesser Himalayan sequences of Kumaon are dominated by mechanically weak, highly fissile phyllites and schists, whose foliation planes, frequently sub-parallel to valley slopes, dramatically reduce shear strength under moisture ingress [4,29]. Proximity to the Main Boundary Thrust (MBT) and Main Central Thrust (MCT) introduces additional vulnerability through fault-proximal fracturing, rock mass weakening, and hydrothermal alteration; slopes within 500 m of mapped fault traces exhibit the highest discontinuity densities and failure susceptibility [5,6,22].

2.3.3. Hydro-Meteorological Factors

Monsoon rainfall and distance to streams capture the primary dynamic trigger and its fluvial amplification. The Indian Summer Monsoon delivers cumulative seasonal totals exceeding 1,500 mm across the higher elevations of Kumaon, with orographic intensification generating sustained high-intensity events responsible for the majority of recorded slope failures in the region [4,6]. Distance to streams reflects the destabilizing influence of active channel incision and lateral toe-undercutting; slopes within 200 m of perennial channels experience elevated pore pressures, progressive debuttressing, and accelerated failure in lithological vulnerable materials [14,30].

2.3.4. Anthropogenic and Environmental Factors

Distance to roads, land use/land cover (LULC), and soil erosion susceptibility collectively represent the anthropogenic and surface environmental dimensions of slope instability. Road construction in mountain terrain, involving slope undercutting, spoil deposition, and disruption of natural drainage, creates concentrated zones of elevated failure probability, particularly within 100 m of cuttings [4,15]. LULC governs root cohesion and infiltration capacity; natural forests provide the greatest slope reinforcement, whereas built-up and barren surfaces intensify runoff and reduce material strength [9,16]. Soil erosion susceptibility serves as a regolith-scale proxy for surface material vulnerability, complementing the lithological factor by capturing shallow translational failure potential in the cultivated and deforested mid-elevation zones characteristic of the Kumaon Himalaya [6].

2.3.5. Factor Standardization and Susceptibility Ranking

To enable multi-criteria overlay, each conditioning factor was reclassified into five susceptibility levels (very low to very high) and assigned ordinal ranks (1–5) based on established geomorphic reasoning and the local Himalayan setting (Table 2). For distance-based factors, proximity to the triggering feature (roads, streams, and faults) was assigned higher ranks; for terrain wetness and slope, higher values were ranked as more susceptible; and for land use/land cover (LULC) and lithology, classes associated with lower material strength or higher anthropogenic disturbance were ranked higher. The reclassified layers were visually checked for class consistency, and the final factor maps were organized as a raster stack for weighting and integration. The spatial distribution of the reclassified conditioning factors across the Kumaon Himalaya is illustrated in Figure 3.

2.4. Weight Determination Approaches

2.4.1. Analytical Hierarchy Process (AHP)

The Analytical Hierarchy Process (AHP), a structured technique for organizing complex decision-making problems, was employed to derive crisp weights for the nine landslide conditioning factors [11,31]. The procedure begins with the construction of an n×n pairwise comparison matrix, A, in which each factor is systematically compared with every other factor using Saaty's 1–9 fundamental scale. The relative importance of the i-th factor with respect to the j-th factor is denoted by $a_{ij}$​, subject to the reciprocal condition $a_{ji}={\displaystyle \frac{1}{a_{ij}}}$​, which ensures matrix consistency. The pairwise comparison matrix is expressed as:

$$A=\left[a_{ij}\right]=\left[\begin{array}{cccc}1& a_{12}& \dots & a_{1n}\\ 1/a_{12}& 1& \dots & a_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ 1/a_{1n}& 1/a_{2n}& \dots & 1\end{array}\right]$$

(1)

In the present study, the factor weights were not obtained from a dedicated expert panel. Instead, the pairwise judgments were established through a comprehensive review of recent landslide susceptibility studies conducted in the Kumaon–Garhwal Himalaya and adjacent Lesser Himalayan regions [4,5,6,22,26]. This approach follows a widely accepted practice in data-constrained mountainous environments, where synthesizing published empirical evidence provides a robust and reproducible alternative to formal expert elicitation [15,16,32]. The resulting 9×9 pairwise comparison matrix, developed based on this literature synthesis, is presented in Table 3.

To compute the normalized factor weights ($W_i$), the pairwise comparison matrix was first normalized by dividing each element by the sum of its corresponding column. The priority weight associated with each factor was then obtained by averaging the elements across the rows of the normalized matrix:

$$W_i={\displaystyle \frac{1}{n}}{\displaystyle {\sum }_{j=1}^n}\left({\displaystyle \frac{a_{ij}}{{\displaystyle \sum _{k=1}^n}{a}_{kj}}}\right)$$

(2)

To assess the logical consistency of the pairwise judgments, the Consistency Ratio (CR) was evaluated. This requires estimation of the maximum eigenvalue (${\lambda }_{max}$​), calculated as the average ratio of the weighted sum vector to the priority weight vector:

$${\lambda }_{max}={\displaystyle \frac{1}{n}}{\displaystyle \sum _{i=1}^n}\left({\displaystyle \frac{{\left(AW\right)}_i}{W_i}}\right)$$

(3)

Finally, the Consistency Ratio was determined as:

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

(4)

where RI is the Random Index corresponding to the matrix order (n=9), and the Consistency Index (CI) is defined as:

$$CI={\displaystyle \frac{{\lambda }_{max}-n}{n-1}}$$

(5)

For the present study, ${\lambda }_{max}=9.294$, CI = 0.037, and CR = 0.025, which is well below the acceptable threshold of 0.10, indicating satisfactory consistency of the pairwise comparisons.

2.4.2. Fuzzy Analytical Hierarchy Process (FAHP)

To account for the uncertainty and subjectivity inherent in expert judgment, the conventional AHP framework was extended using the Fuzzy Analytical Hierarchy Process (FAHP) based on Buckley [33] geometric mean approach. In this method, crisp pairwise comparisons are replaced by Triangular Fuzzy Numbers (TFNs), denoted by

$${\tilde{a}}_{ij}=\left(l_{ij},m_{ij},u_{ij}\right),$$

(6)

where l, m, and u represent the lower, modal, and upper estimates of relative importance, respectively. Accordingly, the fuzzy pairwise comparison matrix is expressed as:

$$\tilde{\mathit{A}}=\left[{\tilde{a}}_{ij}\right]=\left[\begin{array}{cccc}\left(\mathrm{1,1},1\right)& \left(l_{12},m_{12},u_{12}\right)& \dots & \left(l_{1n},m_{1n},u_{1n}\right)\\ \left({\displaystyle \frac{1}{u_{12}}},{\displaystyle \frac{1}{m_{12}}},{\displaystyle \frac{1}{l_{12}}}\right)& \left(\mathrm{1,1},1\right)& \dots & \left(l_{2n},m_{2n},u_{2n}\right)\\ ⋮& ⋮& \ddots & ⋮\\ \left({\displaystyle \frac{1}{u_{n2}}},{\displaystyle \frac{1}{m_{n2}}},{\displaystyle \frac{1}{l_{n2}}}\right)& \left({\displaystyle \frac{1}{u_{2n}}},{\displaystyle \frac{1}{m_{2n}}},{\displaystyle \frac{1}{l_{2n}}}\right)& \dots & \left(\mathrm{1,1},1\right)\end{array}\right]$$

(7)

For each criterion, the fuzzy geometric mean ${\tilde{r}}_i$ is calculated as:

$${\tilde{r}}_i={\left({\displaystyle \prod _{i=1}^n}{l}_{ij}\right)}^{{\displaystyle \frac{1}{n}}},{\left({\displaystyle \prod _{i=1}^n}{m}_{ij}\right)}^{{\displaystyle \frac{1}{n}}},{\left({\displaystyle \prod _{i=1}^n}{u}_{ij}\right)}^{{\displaystyle \frac{1}{n}}}$$

(8)

or, equivalently,

$${\tilde{r}}_i=\left[{\left({\displaystyle \prod _{i=1}^n}{l}_{ij}\right)}^{{\displaystyle \frac{1}{n}}},{\left({\displaystyle \prod _{i=1}^n}{m}_{ij}\right)}^{{\displaystyle \frac{1}{n}}},{\left({\displaystyle \prod _{i=1}^n}{u}_{ij}\right)}^{{\displaystyle \frac{1}{n}}}\right].$$

(9)

The fuzzy weight associated with the $i$-th criterion is obtained by multiplying its fuzzy geometric mean by the inverse of the sum of all geometric means:

$${\tilde{w}}_i={\tilde{r}}_i⨂{\left({\tilde{r}}_1\oplus {\tilde{r}}_2\oplus ···\oplus {\tilde{r}}_n\right)}^{-1}$$

(10)

where $\otimes$ and $\oplus$denote fuzzy multiplication and addition, respectively.

To obtain scalar weights suitable for spatial analysis, the Centre of Area (CoA) defuzzification method was employed. The defuzzified weight corresponding to each criterion was calculated as:

$$w_i={\displaystyle \frac{l_i+m_i+u_i}{3}}$$

(11)

Finally, the defuzzified weights were normalized to ensure that their sum equals unity:

$$W_{FAHP,i}={\displaystyle \frac{W_i}{{\displaystyle {\sum }_{i=1}^n}{W}_i}}$$

(12)

The normalized FAHP weights ($W_{FAHP}$) were subsequently used in the weighted overlay analysis.

2.5. Weight Determination Approaches

The Landslide Susceptibility Index (LSI) for each pixel was calculated using the Weighted Linear Combination (WLC) approach. This method integrates the standardized class rating ($R_i$) and the corresponding factor weight ($W_i$) for all nine conditioning factors ($n=9$):

$$LSI={\displaystyle \sum _{i=1}^n}\left(W_i{R}_i\right)$$

(13)

Separate susceptibility maps were generated using the weights derived from the AHP and FAHP frameworks. Unlike conventional approaches that directly classify pixel-level susceptibility values, the present study retained the continuous LSI values for validation and subsequent aggregation. To support sub-district-scale hazard assessment, pixel-level LSI values were aggregated using zonal statistics. Three aggregation measures—Zonal Mean, Zonal P90, and Zonal Maximum—were evaluated, and the optimal metric was selected based on inventory-based validation (Section 2.6).

The aggregated sub-district susceptibility values were subsequently classified into five categories, namely Very Low, Low, Moderate, High, and Very High, using the Natural Breaks (Jenks) classification method [34]. The Jenks algorithm minimizes within-class variance while maximizing between-class variance by identifying natural discontinuities in the data distribution. Owing to the highly skewed distribution of susceptibility values across the 53 sub-districts, this method produced more meaningful and geomorphologically realistic class boundaries than equal-interval classification schemes, particularly for distinguishing High and Very High susceptibility zones in the upper Himalayan terrain [14,22].

2.6. Aggregation-Scale Sensitivity and Intra-Sub-District Variability Analysis

Because a single zonal statistic may conceal localized hotspots within spatially heterogeneous administrative units, the sensitivity of landslide susceptibility estimates to the choice of aggregation metric was systematically evaluated. Three zonal statistics were derived from the FAHP-based pixel-level LSI raster for each of the 53 sub-districts:

i.

Zonal Mean, representing the arithmetic average of all LSI values within a sub-district and characterizing its overall susceptibility level;

ii.

Zonal 90th Percentile (P90), representing the upper tail of the within-unit LSI distribution and providing greater sensitivity to localized high-risk areas while remaining less influenced by extreme outliers than the maximum statistic;

iii.

Zonal Maximum, representing the highest pixel-level LSI value within a sub-district and corresponding to a worst-case susceptibility scenario. Although useful for identifying extreme conditions, this statistic is highly sensitive to isolated anomalous pixels associated with ridge crests, road cuts, or raster artifacts and therefore may not adequately represent the susceptibility characteristics of the entire administrative unit.

Each aggregated dataset was independently classified into five susceptibility classes (Very Low, Low, Moderate, High, and Very High) using the Natural Breaks (Jenks) method. The predictive performance of each aggregation approach was evaluated at the sub-district scale using receiver operating characteristic area under the curve (ROC-AUC) and Spearman's rank correlation coefficient $\left(\rho \right)$ to assess the extent to which the aggregated statistics preserved the predictive information contained in the original pixel-level LSI surface.

To investigate within-unit spatial heterogeneity, a class-divergence analysis was performed by calculating, for each sub-district, the difference between the Jenks-classified susceptibility categories derived from the P90 and mean aggregation approaches. Positive divergence values indicate sub-districts containing localized high-susceptibility zones that are systematically underestimated by mean-based aggregation. The degree of agreement between the two classification schemes was further quantified using Cohen's kappa coefficient ($\kappa$), which provides a prevalence-independent measure of inter-method consistency. The resulting divergence map was used to identify sub-districts where detailed local investigations and targeted hazard mitigation measures are most urgently required. All analyses were implemented in Python using the rasterstats, scipy, and scikit-learn libraries.

2.6. Model Validation

The predictive performance of the AHP and FAHP models was assessed using receiver operating characteristic (ROC) curve analysis. The validation dataset consisted of 1,715 inventoried landslide locations and an equal number of randomly generated non-landslide points, resulting in a balanced dataset for binary classification (Pourghasemi et al., 2020). During raster extraction, one landslide point coincided with a NoData pixel and was therefore excluded, leaving 1,714 landslide samples for model validation.

Model performance was quantified using the area under the ROC curve (AUC), which measures the ability of the susceptibility model to discriminate between landslide and non-landslide locations. To evaluate whether the observed difference in predictive accuracy between the AHP and FAHP models exceeded sampling uncertainty, a stratified bootstrap procedure was implemented following Mas et al. (2013). A total of 10,000 bootstrap samples were generated with replacement, and the AUC values were recalculated for both models in each iteration. The resulting empirical distributions were used to estimate the 95% percentile confidence intervals (CIs) for the AUC values and their difference.

The statistical significance of the observed improvement in predictive performance was further examined using a one-sided permutation test based on 10,000 permutations under the null hypothesis:

$$H_0:\Delta AUC\le 0$$

(14)

where,

$$\Delta AUC=AU{C}_{FAHP}-AU{C}_{AHP}.$$

(15)

The alternative hypothesis,

$$H_1:\Delta AUC>0,$$

(16)

was tested at a significance level of $\alpha =0.05$. All computations were performed in Python with a fixed random seed to ensure complete reproducibility of the results.

3. Results and Discussion

3.1. Analysis of Factor Importance (AHP vs. FAHP Weights)

The relative importance of the nine conditioning factors derived from both the AHP and FAHP approaches is presented in Table 4. In both models, slope was the most dominant governing factor (AHP weight: 0.2650; FAHP weight: 0.2555), accounting for over a quarter of the total model influence. This prioritization aligns with the established literature on Himalayan landslide dynamics, where slope steepness commonly emerges as the highest-weighted factor because of its direct control on shear stress and gravitational instability [6,14,28]. Monsoon rainfall was assigned the second-highest importance (AHP: 0.1972; FAHP: 0.1964), reflecting the hydrometeorological trigger regime of the Kumaon Himalaya, where orographic intensification during the Indian summer monsoon delivers cumulative seasonal totals exceeding 1,500 mm across higher elevations. This weighting corroborates Dikshit et al. [4] and Chauhan et al. [6], who identified high-intensity monsoon events as the proximate cause of most recorded slope failures in Uttarakhand. Geology(lithology) ranked third (AHP: 0.1501; FAHP: 0.1536), consistent with the geomechanically vulnerability of the Lesser Himalayan sequences, where highly fissile phyllites and schists dominate, and whose foliation-parallel failure planes are well documented in the local literature [29].

The weight deviations between AHP and FAHP, although modest in magnitude, follow a directionally consistent pattern that reflects the structural difference between the two methods. The two most influential factors, slope and monsoon rainfall, were marginally down weighted under FAHP (by −0.95% and −0.08%, respectively), whereas secondary structural and anthropogenic factors, namely, lithology, distance to faults, and distance to roads, were correspondingly upweighted (by +0.35%, +0.29%, and +0.28%, respectively). This redistribution arises from the use of triangular fuzzy numbers in the pairwise comparison stage: by replacing crisp numerical judgments with interval-valued estimates that capture the vagueness inherent in comparative importance assessments, Buckley's geometric mean method distributes weight more evenly across the criteria hierarchy, that prevents the most dominant factors from being overweighted at the expense of secondary synergistic drivers. The practical consequence is a susceptibility surface in which the contribution of mid-ranked factors is more faithfully represented, an effect that prior studies in analogous settings have linked to improved discrimination at susceptibility class boundaries, precisely where misclassification is most consequential for infrastructure planning [16,20]. The combined contribution of LULC and TWI, the two lowest-ranked factors, amounts to less than 6% of total model influence under both methods (AHP: 5.84%; FAHP: 5.76%), a weighting that is geomorphically appropriate for the study area. In the low-gradient depositional terrain of the Terai and Bhabar plains, where slope-driven shear stress and terrain moisture convergence are near-absent, the subdued collective influence of these factors is consistent with the near-zero susceptibility that characterizes the southern sub-districts of Udham Singh Nagar, a pattern corroborated by Chauhan et al. [6], who, using multi-criteria, bivariate, and machine-learning models across Uttarakhand, map Udham Singh Nagar almost entirely within very-low to low susceptibility classes with high-to-very-high zones confined exclusively to the Lesser and Higher Himalayan districts.

3.2. Model Predictive Performance and Validation

Both susceptibility models demonstrated excellent discriminative capacity against the 1,715-point validation inventory, confirming strong spatial concordance between the derived zonation and the historical record of slope failures across the Kumaon Himalaya (Figure 4). Stratified bootstrap resampling over 10,000 iterations [35] yielded AUC values of 0.878 (95% CI: 0.867–0.889) for AHP and 0.887 (95% CI: 0.876–0.897) for FAHP. Although the absolute 95% bootstrap confidence limits of the two models share a marginal intersecting range, the one-sided permutation test returned a p-value of < 0.0001 for an observed ΔAUC of +0.0095 (95% CI: +0.008 to +0.011). This unequivocally confirms that the superior discriminative capacity of the FAHP model is statistically significant at the 0.05 level and cannot be attributed to sampling variation in the landslide inventory.

The magnitude of this improvement (ΔAUC = +0.0095, equivalent to approximately one percentage point in discriminative capacity) falls within the range consistently reported for fuzzy-logic enhancements over crisp MCDM methods in analogous tectonically active settings [10,16,21,36] and is most consequential not in the high-confidence zones of clearly stable or clearly unstable terrain, but at susceptibility class transition boundaries, where infrastructure planning decisions are particularly sensitive to misclassification. The necessity of enhancing standalone AHP is underscored by recent comparative assessments; for example, Nontapot and Chao [37] showed that hybrid and uncertainty-aware approaches consistently outperform unmodified AHP-type models in complex terrain, yielding higher predictive accuracy for rainfall-induced landslide hazard mapping. By replacing crisp pairwise judgments with triangular fuzzy numbers, FAHP captures the vagueness inherent in comparative factor importance assessments, yielding a smoother weight distribution that better reflects the continuous and synergistic nature of slope failure processes. This effect is consistent with the findings of Feizizadeh and Blaschke [20], Sur et al. [22] and Bhagya et al. [16], who report that fuzzy-logic enhancements improve discrimination at susceptibility transition zones, including along road corridors in mountainous and Lesser Himalayan settings. A recent comparative evaluation by Majumder and Fatma [38] in the Shimla district similarly found FAHP to outperform crisp AHP for delineating boundary-zone susceptibility classes in the Western Himalaya, further corroborating the present result.

Both models individually qualify as strong under established AUC benchmarks (AUC > 0.85; Fawcett, 2006), a performance level that compares favourably with recently published knowledge-driven susceptibility models in the broader Himalayan region. Recent studies have reported AUC values of 0.849 for a Fuzzy-AHP framework in the Almora district of Uttarakhand [39] and 0.915 for a similar FAHP approach in the Shimla district of the Western Himalaya [38]. Achieving an AUC of 0.887 across 53 administratively heterogeneous sub-districts highlights the distinct advantage of the FAHP framework for policy-facing hazard assessments. It demonstrates that an inventory-validated approach can deliver high predictive accuracy while maintaining the transparent, auditable weight derivation that is essential for operational decision-making. Consequently, the FAHP model was adopted as the basis for all subsequent subdistrict-level susceptibility analyses.

Having established the statistical superiority of the FAHP model, the analysis subsequently turns to a second methodological challenge: determining the most appropriate approach for aggregating the FAHP-derived pixel-level LSI surface to the sub-district administrative scale. In particular, the study examines whether the choice of aggregation statistic systematically obscures localized high-susceptibility zones within administrative units that would otherwise be classified as having moderate overall susceptibility.

3.3. Aggregation-Scale Sensitivity and Intra-Sub-District Variability

The pixel-level FAHP-derived LSI surface, having been validated as the superior susceptibility model, was subsequently subjected to a systematic aggregation scale sensitivity analysis to determine the most appropriate zonal statistic for sub-district-level hazard classification and to quantify the degree to which mean-based aggregation masks localised high-susceptibility zones within individual administrative units.

3.3.1. Comparative Performance of Zonal Aggregation Methods

The zonal mean achieved the highest sub-district AUC of 0.786 and a Spearman rank correlation of ρ = 0.719 (p < 0.0001) with observed landslide density (Figure 5a), confirming that the central tendency of within-sub-district pixel values provides the most inventory-consistent representation of administrative-unit susceptibility. The zonal P90 returned a moderately lower AUC of 0.728 and ρ = 0.659 (p < 0.0001; Figure 5b), whereas the zonal maximum performed substantially worse at AUC = 0.590 and ρ = 0.455 (p = 0.0006) — performance approaching the no-skill threshold and well below the range of practical utility for governance-level classification (Figure 5d). The inter-method agreement between the Jenks-classified Zonal Mean and Zonal P90 surfaces was moderate (Cohen's κ = 0.431), confirming that the two statistics produce meaningfully divergent susceptibility classifications across a non-trivial share of sub-districts despite their broadly similar rank-order performance.

This progressive deterioration is physically interpretable. The maximum statistic is overwhelmingly dominated by isolated extreme LSI pixels arising from localized topographic features, such as steep ridge crests, road-cut slopes, or micro-topographic artefacts in the SRTM-derived slope layer, which register very high conditioning factor values but do not correspond to the sub-district locations where historical landslides have concentrated in the aggregate. Although the P90 is more robust than the maximum, it still captures the upper tail of the within-unit distribution and is correspondingly better at flagging localized high-susceptibility zones than predicting where events have historically clustered across the full unit. The mean, by weighting every pixel equally, best tracks the dominant landscape character of the administrative unit and therefore most faithfully reproduces the spatial pattern of the inventory.

Figure 5. Comparison of sub-district-level landslide susceptibility zonation derived from three zonal aggregation methods applied to the FAHP LSI raster, classified using Natural Breaks (Jenks): (a) Zonal mean (AUC = 0.786); (b) Zonal P90 (AUC = 0.728); (c) class-divergence map showing sub-districts where P90 and Mean methods yield different susceptibility classes; and (D) ROC curves with AUC bar chart inset comparing the discriminative performance of all three methods against the landslide inventory (n = 3428 validation points).

Figure 5. Comparison of sub-district-level landslide susceptibility zonation derived from three zonal aggregation methods applied to the FAHP LSI raster, classified using Natural Breaks (Jenks): (a) Zonal mean (AUC = 0.786); (b) Zonal P90 (AUC = 0.728); (c) class-divergence map showing sub-districts where P90 and Mean methods yield different susceptibility classes; and (D) ROC curves with AUC bar chart inset comparing the discriminative performance of all three methods against the landslide inventory (n = 3428 validation points).

The conceptual basis for treating the pixel-to-polygon aggregation step as methodologically non-trivial lies in the terrain mapping unit framework of Calvello et al. [40], who formally distinguished terrain computational units used for model derivation from terrain zoning units used for administrative reporting, establishing that the transformation between them, including the aggregation algorithm, is consequential for the integrity of the resulting zoning map. Extending this principle empirically, Jacobs et al. [41] demonstrated that mapping unit choice substantially affects susceptibility classification outcomes independent of the underlying algorithm, a finding that applies with equal force to the choice of aggregation statistic within a fixed administrative boundary, as confirmed by the systematic reclassification effects documented here. These results provide robust, multi-metric justification for retaining the zonal mean as the primary aggregation statistic for sub-district-level susceptibility classification, while simultaneously motivating the use of P90 as a complementary diagnostic indicator for identifying localized risk concentrations.

3.3.2. Intra-Sub-District Class Divergence and Risk of Susceptibility Underestimation

The class-divergence analysis reveals a spatially systematic and operationally significant gap between mean-based and P90-based classifications across the study area (Figure 5c). Of the 53 sub-districts, 24 (45.3%) were assigned to a higher Jenks susceptibility class under P90 than under the mean, while the remaining 29 (54.7%) retained identical classifications under both statistics. Critically, not a single sub-district was downgraded by P90 relative to the mean, a directionally consistent result that follows structurally from the right-skewed LSI distributions characteristic of heterogeneous Himalayan terrain, where the 90th percentile is necessarily greater than or equal to the mean.

This unidirectional pattern has direct practical implications: adopting the zonal mean as the sole aggregation statistic introduces a systematic underestimation of localized risk in nearly half of all sub-districts studied. In these 24 units, villages, slope corridors, or infrastructure segments situated in the upper tail of the within-sub-district LSI distribution are being administratively classified at a lower susceptibility category than the pixel-level evidence warrants, an effect directly analogous to the well-known statistical limitation of using the mean to summarize a strongly right-skewed distribution. The divergence is most pronounced in the mid-Himalayan sub-districts of the Bageshwar and Almora districts, where deeply incised valleys alternate with flatter terraced agricultural terrain, generating wide intra-unit LSI ranges that the mean cannot adequately represent.

Of the 24 divergent sub-districts, 22 exhibited a one-class upward shift under P90, while Dharchula and Haldwani exhibited the most extreme two-class divergence. Dharchula (2,694 km²) shifted from moderate to very high under P90, a consequence of the enormous internal heterogeneity of a sub-district spanning glaciated plateaus, deeply incised river gorges, and steep MCT-proximal slopes within a single administrative boundary. Haldwani shifted from low to high, driven by the abrupt topographic transition from the Bhabar plain to the Shivalik foothills. These two sub-districts represent the highest-priority targets for immediate intra-sub-district susceptibility investigation and constitute instructive cases of how mean-based administrative classification can present a misleadingly benign picture of localized hazard conditions. This finding resonates with the observation by Schlögl et al. [42] that uncertainty-transparent susceptibility products, particularly those expressed through confidence intervals or complementary upper-tail statistics, are essential for honest hazard communication to decision-makers in settings where the administrative reporting unit is geomorphically heterogeneous.

The systematic comparison of aggregation methods presented here addresses a methodological blind spot that has persisted in the administrative unit susceptibility mapping literature. To the authors' knowledge, no prior study in the Himalayan landslide susceptibility domain has empirically evaluated alternative zonal aggregation approaches at the sub-district scale using ROC-validated metrics or quantified the resulting intra-unit risk masking through a class-divergence framework. Prior sub-district and district-level studies in analogous settings, including Chauhan et al. [6], Yadav et al. [24], and Bhandari and Wankhade [27] have uniformly adopted mean-based aggregation without assessing whether this choice systematically under classifies risk in spatially heterogeneous units. The present analysis fills this gap across all 53 sub-districts with inventory-validated evidence, demonstrating that the consequences of this choice are spatially systematic and directionally consistent. The zonal mean and class-divergence map are therefore most appropriately used in tandem: the former provides the inventory-validated regional classification for resource allocation, and the latter flags the 24 sub-districts where intra-unit investigation is most urgently warranted.

The finding that the Zonal Mean maximizes predictive accuracy (AUC = 0.786) while obscuring localized high-susceptibility zones in 24 of 53 sub-districts highlights a fundamental trade-off in administrative hazard mapping: statistical representativeness and worst-case risk sensitivity cannot be optimized by a single metric. While the Mean accurately captures aggregate susceptibility for regional inventory validation, the P90 is indispensable for delineating severe, localized hazards critical to worst-case infrastructure and settlement planning.

Operationally, the class-divergence map aligns with the Sendai Framework for Disaster Risk Reduction 2015–2030, which explicitly calls for the disaggregation of risk data to the local level to enable targeted early-warning systems and community-level preparedness planning (UNDRR, 2015). By explicitly flagging divergent units where aggregate statistics present a deceptively benign picture, it provides district administrations in Bageshwar, Almora, and Pithoragarh with an evidence-based prioritization layer. This framework directly enables targeted village-level slope reconnaissance, site-specific geotechnical investigations, and the mandatory integration of pixel-level overlays into Environmental Impact Assessments for proposed infrastructure.

3.4. Spatial Distribution of Landslide Susceptibility

The FAHP-based jenks-classified susceptibility map (Figure 6) reveals a markedly heterogeneous hazard landscape across the Kumaon Himalaya. As summarized in Table 5, nine (16.98 %) sub-districts fall within the very high class and ten (18.87%) within the high class; together, these 19 units account for 54.5% of all 1,715 recorded landslide events, 934 events concentrated in only 35.85% of administrative units, providing strong empirical validation of the model's discriminative power. The very high susceptibility class is dominated by five sub-districts within the Nainital district: Betalghat, Nainital, Kosya Kutauli, Dhari, and Khanshyu.

This spatial concentration reflects the combined influence of two intersecting drivers. At the structural scale, the Main Boundary Thrust traverses the district, introducing a zone of intensely fractured and mechanically weakened rock masses susceptible to moisture-triggered failure; MBT-proximal slopes within Nainital district carry the highest discontinuity densities and fault-amplified pore-pressure responses in the study area, consistent with the structural controls documented by Sur et al. [22] and Pandey et al. [5] for analogous MBT-crossing corridors in the Garhwal Himalaya. At the anthropogenic scale, decades of unregulated tourism infrastructure, slope loading from construction, road widening, and progressive vegetation clearance have dramatically elevated the failure probability of already precarious slopes; the FAHP weight for distance to roads (0.0911) is directly substantiated by this pattern. The Nainital sub-district alone recorded 196 landslide events and a density of 0.265 events/km², the highest absolute count in the study area, consistent with the bivariate susceptibility assessment by Bhandari and Wankhade [27], who identified road corridors and slopes exceeding 65° around Nainital as among the most failure-prone terrain in the Kumaon Himalayas.

A second concentration of high susceptibility occurs in the Champawat district, where Champawat (85 events; 0.161 events/km²) and Barakot (54 events; 0.364 events/km²) are classified at the highest susceptibility level. Barakot recorded the highest landslide density in the entire study area, reflecting the dominance of fragile, highly fissile phyllites and schists under intense monsoon precipitation, a combination identified as the primary predisposing condition for rainfall-triggered slope failures in the Lesser Himalayan sequences [6,29]. The added role of intense monsoon precipitation as a dynamic trigger amplifying this structural predisposition is consistent with rainfall–failure coupling documented across the Lesser Himalayan sequences of Uttarakhand (Dikshit et al., 2020; Chauhan et al., 2025). Two additional sub-districts with very high susceptibility are Bangapani in the Pithoragarh district (47 events; 0.245 events/km²) and Bhikiyasain in the Almora district (18 events; 0.082 events/km²), both characterized by the same lithological and precipitation conditions, confirming that the phyllite–monsoon coupling is a spatially persistent predisposing mechanism across the Lesser Himalayan belt rather than a localized anomaly.

The high susceptibility class encompasses 10 sub-districts arranged in two geographically contiguous belts. Within the Almora district, Almora, Ranikhet, Syaldey, Bhanoli, and Chaulcheena form a continuous arc of infrastructure-intensive, tourist-frequented slopes across the Lesser Himalayan phyllite zone, where road density is among the highest in the region. Bhanoli (90 events; 0.265 events/km²) and Syaldey (59 events; 0.197 events/km²) recorded the most concentrated landslide activity within this belt. The Pithoragarh district contributed Tejam (57 events; 0.264 events/km²), Pati (53 events; 0.173 events/km²), Thal, and Lohaghat, where high rainfall intensity and structural fragility in the outer Himalayan sequences drive susceptibility. This MCDM-derived pattern is consistent with corridor-scale studies from the Lesser Himalaya, where FAHP- and AHP-based susceptibility mapping along the Kalsi–Chakrata road corridor [22] and National Highway 5 [15] likewise identified lithology and road-related variables as among the most influential conditioning factors, with the highest susceptibility zones concentrated on structurally weak, foliated rock units immediately adjacent to major highways. At the state scale, Chauhan et al. [6], using bivariate, MCDM, and machine learning frameworks, consistently identified Nainital and Almora as the highest-risk districts in Uttarakhand, while Badola et al. [43] similarly flagged phyllitic lithology and road density as leading predictors of very high susceptibility in the Lower Himalayan belt, though the latter finding reflects machine learning variable importance rather than MCDM weight assignment and the two are not directly commensurable.

The large northern sub-districts of Munsyari (2,448 km²), Dharchula (2,694 km²), and Kapkot (1,289 km²) are classified as Moderate by the Zonal Mean, an outcome that is physically interpretable but operationally incomplete. Their vast expanses of glaciated plateaus, uninhabited valley floors, and snowfields at high elevation structurally suppress the sub-district-mean LSI well below what the locally steep, MCT-proximal slopes would register at the pixel scale. As established in Section 4.3.2, Dharchula exhibits the most extreme two-class divergence in the dataset and harbours locally very high hazard pockets, with direct relevance to the proposed Tanakpur–Pithoragarh road corridor passing through this sub-district. Munsyari similarly exhibits a one-class upward divergence from Moderate to High under P90, warranting targeted reconnaissance for its alpine road and hydropower corridors.

The Very Low and Low susceptibility classes are almost entirely confined to the southern Terai and Bhabar plains of Udham Singh Nagar district: Bazpur, Jaspur, Kashipur, Kichha, Rudrapur, Sitarganj, Gadarpur, Nanakmatta, Khatima, and Lalkuan, where low-gradient depositional terrain, stable alluvial geology, and minimal slope steepness produce near-zero susceptibility. No landslide events were recorded in any Very Low-classified sub-district, and only 108 events (6.29% of the inventory) fall within the Low category. This near-total absence of slope failures in Udham Singh Nagar is physically expected and provides an implicit consistency check on the model: the same weight framework that identifies Nainital as having very high susceptibility correctly assigns near-zero susceptibility to the Gangetic foreland, confirming that the FAHP weighting scheme is geomorphically coherent across the full topographic range of the study area.

Figure 6. Sub-district-level landslide susceptibility zonation of the Kumaon Himalaya derived from the FAHP-based Landslide Susceptibility Index (LSI) aggregated to the sub-district administrative scale using the Zonal Mean statistic (sub-district AUC = 0.786; Spearman ρ = 0.719, p < 0.0001).

Figure 6. Sub-district-level landslide susceptibility zonation of the Kumaon Himalaya derived from the FAHP-based Landslide Susceptibility Index (LSI) aggregated to the sub-district administrative scale using the Zonal Mean statistic (sub-district AUC = 0.786; Spearman ρ = 0.719, p < 0.0001).

4. Conclusions

This study presented a comprehensive sub-district-level landslide susceptibility assessment of the Kumaon Himalaya, integrating remote sensing data with multi-criteria decision-making frameworks. The comparative analysis between AHP and Fuzzy-AHP approaches provided critical insights into the spatial distribution of slope instability and the role of expert uncertainty in hazard modelling. The following conclusions are drawn:

a)

Slope morphometry and monsoon precipitation emerged as the most significant predisposing factors, collectively governing over 45% of the model influence. This reaffirms that gravitational instability triggered by high-intensity rainfall remains the primary hazard mechanism in the region.

b)

While both models demonstrated excellent predictive capability (AUC > 0.87), the FAHP approach (AUC = 0.887) proved marginally superior to the conventional AHP (AUC = 0.878). The integration of fuzzy logic effectively mitigated the subjectivity and epistemic uncertainty inherent in expert weighting, resulting in a more robust delineation of transition zones between susceptibility classes.

c)

The study introduces the first systematic multi-metric comparison of zonal aggregation approaches (Zonal Mean, Zonal P90, and Zonal Maximum) for sub-district-scale landslide susceptibility mapping in the Himalayas. Validation against a balanced landslide inventory using ROC–AUC, Spearman rank correlation with landslide density, and Cohen's Kappa agreement demonstrated that the Zonal Mean approach (AUC = 0.786; ρ = 0.719, p < 0.0001) provides the strongest predictive performance and is the most suitable metric for sub-district classification. Class-divergence analysis further revealed that 24 of 53 sub-districts (45.3%) were assigned higher susceptibility classes under P90-based aggregation, while none were downgraded, with a moderate inter-method agreement (κ = 0.431). These findings indicate that mean-based aggregation can mask localized high-susceptibility zones within spatially heterogeneous administrative units. The resulting divergence map provides a practical decision-support tool for prioritizing hazard assessment and risk-management efforts in Bageshwar, Almora, and Pithoragarh.

d)

The sub-district-level susceptibility map, derived from the FAHP-LSI surface aggregated using the inventory-validated Zonal Mean statistic, reveals a spatially differentiated hazard landscape across the Kumaon Himalaya. Nine sub-districts (16.98%) are classified as Very High susceptibility and 10 sub-districts (18.87%) as High, together accounting for 54.5% of all recorded landslide events despite comprising only 35.85% of administrative units — a concentration ratio that empirically validates the model's discriminative power. The highest absolute hazard is concentrated in the Nainital district (Betalghat, Nainital, Kosya Kutauli, Dhari, Khanshyu) and Champawat district (Champawat, Barakot), with Barakot recording the highest landslide density in the study area (0.364 events/km²). Critically, the large northern sub-districts of Munsyari and Dharchula, though classified as Moderate in aggregate, harbour locally Very High susceptibility zones identified through the class-divergence analysis, warranting targeted intra-sub-district reconnaissance for their major infrastructure corridors.

The spatially explicit susceptibility framework at the administrative (sub-district) scale is a vital decision-support tool for governance. It is imperative that future infrastructure planning, particularly road widening and hydropower expansion in the identified high-risk sub-districts, incorporates these susceptibility layers into Environmental Impact Assessments (EIA) to minimize disaster risks and economic losses.

5. Limitations and Future Directions

A key limitation of the present study is its exclusive reliance on static conditioning factors, which capture spatial predisposition to failure but do not incorporate dynamic triggering variables, particularly event-level rainfall intensity and antecedent soil moisture, which are the proximate causes of pore pressure exceedance and slope failure initiation in rainfall-triggered landslide systems [4,9]. This is an inherent constraint of the MCDM-based susceptibility paradigm, which is designed to characterize where failures are likely to occur rather than when. A related limitation concerns the completeness of the GSI-sourced inventory used for validation. Inventories compiled through government reporting are well documented to exhibit spatial bias toward road- and settlement-proximal events, while underrepresenting failures in remote high-relief terrain such as Munsyari and Dharchula [44]; the reported AUC of 0.887 should therefore be interpreted as reflecting discriminative capacity over the documented, accessible portions of the landscape rather than as a universal measure across the full topographic range of the study area. Future work should explore the integration of a dynamic triggering layer, derived from satellite soil moisture products, such as SMAP or Sentinel-1 SAR backscatter, or from CHIRPS-based rainfall anomaly indices, with the static FAHP-LSI surface generated herein. Such a combined predisposition–trigger framework would substantially advance the operational utility of this work for near-real-time early warning, particularly in the context of projected monsoon intensification across the Himalayan arc under future climate scenarios.