Integrating InSAR and Channel Steepness for AI-Based Coseismic Landslide Modeling in the Nepal Himalaya

1. Introduction

The 25 April 2015 M_w 7.8 Gorkha earthquake and its major aftershocks, including the M_w 7.3 event of 12 May 2015, triggered tens of thousands of landslides across central Nepal, causing widespread destruction, blocking critical road networks, damming rivers, and contributing to significant casualties and economic losses [1,2,3]. In steep, rapidly uplifting mountain belts such as the Nepal Himalaya, earthquake-induced landslides amplify shaking impacts and prolong disaster cascades by disrupting access, damaging infrastructure, and mobilizing large volumes of sediment. Coseismic landslide susceptibility mapping (LSM) and probability modeling are therefore essential for post-event response, reconstruction planning, and corridor-scale risk management [4,5]. Recent case studies in Nepal have demonstrated that landslide occurrence is strongly associated with geological structures, lithological contrasts, and antecedent rainfall, underscoring the need for multi-factor conditioning frameworks that integrate tectonic, hydrological, and terrain controls [6].

The Nepal Himalaya is characterized by a complex tectonic architecture, with the Main Himalayan Thrust (MHT) accommodating convergence between the Indian and Eurasian plates and driving active uplift and erosion across the orogen [7,8]. Lithology varies systematically from south to north: Lesser Himalayan metasediments and granites, Higher Himalayan crystalline rocks, and Tethyan sedimentary sequences, each with distinct mechanical properties and susceptibility to slope failure [1,8]. Structural discontinuities including thrust faults, shear zones, and lithological contacts concentrate stress and weaken hillslopes, making them preferential sites for earthquake-triggered landslides [1,2]. Coseismic landslide inventories for the 2015 Gorkha event show that landslide density correlates with slope gradient, peak ground acceleration, surface deformation, and proximity to mechanically weak lithologies and large plutonic intrusions [1,2]. Antecedent rainfall further modulates slope stability by elevating pore-water pressure and reducing effective stress, particularly in weathered and fractured rock [9]. Monsoon-driven erosion and tectonic uplift sustain high relief and over-steepened hillslopes that are prone to failure during strong ground motion [7]. Integrating geological, geomorphic, hydrologic, and seismic conditioning factors is therefore essential for robust coseismic landslide probability modeling in the Nepal Himalaya [6,10].

Over the past two decades, statistically based and machine learning (ML) approaches have become standard tools for LSM [4,5,6,11,12,13,14]. Classical methods such as logistic regression, generalized additive models, and Bayesian approaches have been increasingly complemented or superseded by ensemble tree methods (e.g., Random Forest, gradient boosting, XGBoost) and, more recently, deep neural networks [11,12,13,14,15,16,17,18]. These models can flexibly exploit multi-source geospatial data, capture complex nonlinear relationships among conditioning factors, and produce spatially continuous probability fields suitable for operational decision-making under limited data availability and strong class imbalance [5,11,12,13,14,15,16,17,18,19]. Recent work has also emphasized the role of deep learning architectures, including convolutional and recurrent networks, in improving landslide detection, susceptibility mapping, and inventory compilation when combined with high-resolution remote sensing data [16,17,18,19,20,21].

In parallel, advances in satellite and airborne remote sensing have dramatically expanded the availability of high-resolution topography, land cover, and seismic shaking products for landslide modeling. Near-global DEMs, multispectral optical imagery, and long time series of SAR acquisitions now allow for systematic characterization of terrain, vegetation, hydrology, and ground deformation. In the context of event-based coseismic LSM, this enables models that not only rely on static conditioning factors but also incorporate transient signals related to deformation, coherence loss, and hydrometeorological forcing [5,22,23,24,25,26,27]. For example, recent studies have integrated multi-temporal InSAR deformation velocities or persistent scatterer (PS) deformation rates into ML-based susceptibility frameworks, demonstrating clear gains in model skill compared with purely static terrain-based approaches [20,25,26,27,28]. Such dynamic, InSAR-enhanced LSM frameworks are particularly powerful in alpine gorge and corridor settings where continuous monitoring of slope response is required but in situ observations are sparse [25,26,28].

Most LSM studies still rely heavily on static conditioning factors such as slope, aspect, curvature, elevation, lithology, distance to faults or drainage, PGA, rainfall, and land use/land cover (LULC) [4,5,29,30,31]. However, static terrain and environmental variables only partially capture the transient deformation and damage associated with large earthquakes. Spaceborne interferometric synthetic aperture radar (InSAR) offers unique capabilities to measure ground deformation and coherence loss at regional scales, even under cloud cover and during the monsoon [22,23,24,31,32]. InSAR has been successfully used to map coseismic crustal deformation for the Gorkha earthquake [31,32] and to investigate landslide processes and unstable slopes in other regions, including along major transportation corridors and in rapidly urbanizing basins [22,23,25,26,27]. Integrating InSAR-derived products such as line-of-sight (LOS) displacement and coherence-based damage proxy maps (DPM) into coseismic landslide models therefore has the potential to highlight zones of deformation, coherence loss, and surface disturbance that are not readily inferred from static terrain metrics alone. Yet, systematic comparisons of ML and deep learning (DL) models with and without InSAR inputs remain limited for the Nepal Himalaya.

Topographic characterization for landslide modeling and geomorphic analysis increasingly relies on near-global digital elevation models (DEMs) such as the Shuttle Radar Topography Mission (SRTM) 1 arc-second product, which provides approximately 30 m resolution elevation data for most land surfaces between 60°N and 56°S [33]. The SRTM 1 arc-second DEM has become a standard input for terrain derivatives (slope, aspect, curvature), stream network extraction, and channel steepness calculations, as well as for coregistering and geocoding interferometric synthetic aperture radar (InSAR) products in mountain regions where high-resolution local DEMs are scarce [33].

At the same time, geomorphic indices derived from DEMs, notably the normalized channel steepness index (Kₛₙ), provide insight into long-term tectonic forcing, base-level change, and hillslope stability [34,35]. High Kₛₙ values are often associated with steep, actively incising channels, knickzones, and over-steepened hillslopes that are prone to slope failure [34,35,36,37]. Recent work has highlighted the potential of combining Kₛₙ with remote sensing data and statistical or ML methods to improve detection and mapping of landslides and other mass movements [36,37,38]. Because Kₛₙ can be derived systematically from widely available DEMs such as SRTM, it offers a promising, physically interpretable proxy for longer-term tectonic–geomorphic predisposition to landsliding that can be integrated with more transient InSAR-based indicators of ongoing deformation [33].

In parallel with these physical and geomorphic advances, there has been a rapid growth in the use of deep learning and geospatial artificial intelligence for landslide studies. Recent reviews document how convolutional, recurrent, and hybrid neural networks, when combined with optical, LiDAR, and synthetic aperture radar (SAR) data, can improve landslide inventory mapping, susceptibility forecasting, and early warning, while also highlighting key challenges such as data imbalance, limited regional transferability, and model interpretability [11,15,16,17,18,19,20,21]. InSAR coherence and phase time series have similarly been exploited using deep learning for damage and deformation mapping, where recurrent networks act as probabilistic anomaly detectors to separate true damage from background surface changes [19,39,40]. A growing body of literature demonstrates that integrating InSAR-derived deformation velocities or active deformation with machine learning and deep learning algorithms, including graph convolutional networks, gradient-boosted decision trees, and window-based atmospheric correction methods, significantly enhance landslide susceptibility mapping in alpine gorge regions and geotectonically complex settings [26,28]. These developments motivate a closer integration of deep segmentation architectures (e.g., U-Net, DeepLab-style models) with physically meaningful InSAR and geomorphic predictors for coseismic landslide probability mapping.

This study addresses these gaps by jointly integrating InSAR-derived and geomorphic predictors into an event-based coseismic landslide probability modeling framework for the 2015 Gorkha earthquake. We derive LOS displacement from ALOS-2 PALSAR-2 ScanSAR interferometry and compute a normalized channel steepness index (Kₛₙ) from DEM data, then combine these predictors with a suite of 14 conditioning factors that capture topography, geomorphology, hydrology, seismic shaking, and land cover. We train and compare five representative models: Random Forest (RF), XGBoost, a fully convolutional CNN, U-Net, and DeepLabV3.

The specific objectives of this study are to:

(1) describe the derivation of LOS displacement and DPM from ALOS-2 ScanSAR and of Kₛₙ from a DEM, and integrate these into a multi-factor coseismic landslide probability framework for the Gorkha 2015 region;

(2) compare the performance of RF, XGBoost, CNN, U-Net, and DeepLabV3 using discrimination, detection, and calibration metrics (AUC-ROC, AUC-PR, CSI, Brier score);

(3) quantify the contribution of InSAR-derived layers by contrasting models with and without LOS and DPM inputs; and

(4) identify dominant conditioning factors, especially Kₛₙ, and discuss their implications for earthquake-induced landslide occurrence and geomorphic controls across the Nepal Himalaya.

2. Materials and Methods

2.1. Study Area

The study area encompasses the region strongly affected by the 2015 Gorkha earthquake sequence in central Nepal, including districts such as Gorkha, Sindhupalchok, Rasuwa, Nuwakot, Dhading, and parts of Kathmandu and surrounding valleys. The area lies within the Nepal Himalaya, characterized by high relief, steep slopes, active thrust faulting (e.g., Main Himalayan Thrust, MHT), and monsoon-driven erosion [7,8]. Lithology varies from Lesser Himalayan metasediments and granites to Higher Himalayan crystallines and Tethyan sediments, with structural and lithological controls on slope stability [9]. The M_w 7.8 mainshock and M_w 7.3 aftershock produced strong ground shaking and widespread slope failures, with landslide inventories documenting thousands of triggered landslides [1,2]. The spatial extent of the analysis is defined by the coverage of the ALOS-2 ScanSAR interferometric products and the available landslide inventory, yielding a raster domain of approximately 3,467 × 5,573 pixels at the chosen resolution (Section 2.2).

Figure 1. Study area map showing the Central Nepal Himalaya, the epicenter of the 2015 Gorkha earthquake, and the spatial extent of the landslide inventory. Landslide data are compiled from Roback et al. [2,41] and the U.S. Geological Survey (USGS) [3].

Figure 1. Study area map showing the Central Nepal Himalaya, the epicenter of the 2015 Gorkha earthquake, and the spatial extent of the landslide inventory. Landslide data are compiled from Roback et al. [2,41] and the U.S. Geological Survey (USGS) [3].

2.2. Data

2.2.1. Landslide Inventory

The landslide inventory used in this study is a binary raster (landslide vs. non-landslide) derived from post-event satellite imagery and the landslide map data of Roback et al. [41] for the 2015 Gorkha earthquake, available via U.S. Geological Survey ScienceBase. Landslide pixels were assigned a value of 1 and non-landslide pixels 0; areas with no data (e.g., clouds, water bodies) were masked. The inventory exhibits severe class imbalance typical of event-based LSM: landslide pixels represent a small fraction (~0.47%) of the total valid pixels. This imbalance was explicitly addressed in both ML and DL pipelines (Section 2.4).

2.2.2. InSAR-Derived Products (ALOS-2 ScanSAR): Derivation of LOS Displacement and DPM

ALOS-2 L-band ScanSAR (WBDR1.1) single-look complex (SLC) imagery was acquired from the Japan Aerospace Exploration Agency (JAXA) on 22 February, 5 April, and 17 May 2015, covering frame 3050. The M_w 7.8 mainshock of 25 April 2015 occurred between the April and May acquisitions. Two interferometric pairs were formed: a pre-seismic pair (22 February–5 April 2015) to provide pre-event coherence for the damage proxy map, and a coseismic pair (5 April–17 May 2015) spanning the earthquake for line-of-sight (LOS) displacement and post-event coherence. The acquisitions used for LOS displacement and the DPM are summarized in Table 1 and Table 2, respectively. All InSAR processing was performed using the InSAR Scientific Computing Environment (ISCE) version 2.6.3 [42,43,44].

The processing workflow followed the ALOS-2 ScanSAR burst-mode InSAR chain implemented in ISCE 2 [42]. Burst synchronization was applied to remove unsynchronized signals, followed by SLC coregistration using water-body masking and DEM-based refinement with the SRTM 1 arc-second global product [33]. Interferograms were formed with 7 range and 2 azimuth looks, filtered using the Goldstein algorithm (strength 0.8, window 64, step 4), and unwrapped with water body masking. Ionospheric phase estimation and correction were applied to mitigate L-band ionospheric effects [45]. Final products were geocoded to geographic coordinates using bicubic interpolation. The InSAR processing flow diagram is shown in Figure 2.

The LOS displacement d_LOS (cm) was computed from the unwrapped interferometric phase φ (radians) as:

$$d_{LOS}=-{\displaystyle \frac{\left(\phi \times \lambda \right)}{4\pi }}\times 100$$

(1)

where λ = 23.6 cm is the ALOS-2 L-band wavelength. The negative sign conforms to the convention that positive phase corresponds to increasing range (motion away from the sensor). The geocoded LOS displacement from the coseismic pair is shown in Figure 3(a).

The Damage Proxy Map (DPM) was derived from the difference between pre-event and post-event coherence. Pre-event coherence γ_pre was obtained from the 22 February–5 April interferogram, and post-event coherence γ_post from the 5 April–17 May interferogram. The two coherence maps were registered to a common grid, and the DPM was computed as:

DPM = γ_pre − γ_post

(2)

DPM values were clipped to [No damage, High damage] for visualization, with higher values indicating greater coherence loss and thus stronger damage proxy. The DPM is shown in Figure 3(b).

2.2.3. Tectonic Geomorphology and the Channel Steepness Index (Kₛₙ)

The normalized channel steepness index (Kₛₙ) was computed to assess spatial variations in river-profile steepness and infer relative tectonic activity across the study area. The Kₛₙ is derived from the stream power incision model, which relates channel slope S to drainage area A through a power-law relationship [46,47]:

$$S=kₛA^{-\theta }$$

(3)

where _ks is the channel steepness index and θ is the concavity index. Under steady-state conditions, steepness scales with rock uplift rate relative to erodibility [35,47]. To facilitate comparison across drainage basins with varying concavities, a reference concavity θ_ref is applied to define the normalized steepness index [34]:

$$kₛₙ=S{A}^{\theta \_ref}$$

(4)

where S is channel slope (dimensionless, or in radians), A is drainage area (m²), and θ_ref is the reference concavity. A commonly used value is θ_ref = 0.45, representative of bedrock rivers in equilibrium [34].

Data and Processing

A digital elevation model (DEM) was used to extract the stream network and compute Kₛₙ. Processing was performed in MATLAB using TopoToolbox [48]. The workflow comprised: (1) filling sinks in the DEM; (2) computing flow direction and flow accumulation to obtain drainage area; (3) extracting the stream network above a drainage area threshold of 1 km²; (4) computing local channel slope using an 8-direction gradient (in radians); (5) sampling drainage area and slope at each stream node; and (6) calculating Kₛₙ at each node via Equation (4). The mean Kₛₙ per stream segment was then computed and exported as a polyline shapefile for mapping and overlay with other datasets.

The resulting Kₛₙ map highlights channels that are steeper than expected for their drainage area. Elevated Kₛₙ values are interpreted as indicating higher relative rock uplift or lower erodibility and may coincide with active structures or transient landscape response to tectonic forcing (Figure 4).

2.2.4. Other Conditioning Factors

The study utilizes fourteen conditioning factors representing topography, geomorphology, hydrology, land cover, and seismic intensity. These include CosAspect, SinAspect, DEM, Distance to River, DPM, Drainage Density, Kₛₙ, LULC, LogSPI, PGA, Rain, Slope, and the absolute values of Curvature and LOS displacement. Peak ground acceleration (PGA) was obtained from the USGS ShakeMap [49] and earthquake database [50]. Land use/land cover (LULC) was derived from the ICIMOD Regional Land Cover Monitoring System (RLCMS) for the Hindu Kush Himalaya (HKH) region (2015) [51]. The DEM used for Kₛₙ derivation and InSAR processing was the SRTM 1 arc-second global product [33]. The sine and cosine transformations of aspect were used to appropriately represent its circular nature and to avoid artificial discontinuities at the 0°/360° boundary, thereby enabling the machine-learning models to capture directional effects more effectively.

In addition to the gridded predictor layers, two auxiliary data visualizations were prepared to support interpretation of the physical controls on earthquake-triggered landslides and to justify the inclusion of selected conditioning factors. Figure 5a presents a circular (rose) diagram of landslide counts by slope aspect derived from the mapped inventory. This exploratory analysis was used to evaluate directional bias in landslide occurrence and to motivate the explicit encoding of aspect using its cosine and sine components within the modeling framework, consistent with previous earthquake-induced landslide studies [1,2,30,31].

Figure 5b shows the monthly rainfall time series from September 2014 to August 2015 for the major earthquake-affected districts. This hydrometeorological summary was compiled to characterize antecedent moisture conditions prior to the April 2015 coseismic period and to support the inclusion of rainfall as a conditioning factor. The district-wise rainfall data used to generate Figure 5b are provided in Appendix Table A1. Prior Himalayan research indicates that elevated antecedent precipitation can increase pore-water pressure and reduce slope stability [9].

Normalization and Absolute Value Transformations

To ensure effective training, all continuous factors were resampled to a common 30-meter grid and normalized using the Min-Max formula. All spatial analysis, raster processing, and map preparation were performed using ArcGIS Pro 3.6.0 (Esri, Redlands, CA, USA).

$$X_{\left\{norm\right\}}={\displaystyle \frac{X–X_{\left\{min\right\}}}{X_{\left\{max\right\}}-X_{\left\{min\right\}}}}$$

(5)

where:

X = original pixel value

X_min = minimum value of the factor

X_max = maximum value of the factor

X_norm = normalized value scaled between 0 and 1

Calculating the absolute values for LOS and curvature was a critical refinement, as it collapsed bi-directional risks (toward/away or convex/concave) into a single magnitude-based susceptibility indicator.

Figure 6. Spatial distribution of the fourteen conditioning factors used for landslide modeling: (a) absolute LOS displacement, (b) Kₛₙ, (c) DPM, (d) DEM, (e) distance to river, (f) PGA, (g) cos(aspect), (h) sin(aspect), (i) absolute curvature, (j) LULC, (k) annual rainfall (2014–2015), (l) slope, (m) logSPI, and (n) drainage density.

Figure 6. Spatial distribution of the fourteen conditioning factors used for landslide modeling: (a) absolute LOS displacement, (b) Kₛₙ, (c) DPM, (d) DEM, (e) distance to river, (f) PGA, (g) cos(aspect), (h) sin(aspect), (i) absolute curvature, (j) LULC, (k) annual rainfall (2014–2015), (l) slope, (m) logSPI, and (n) drainage density.

2.3. ML and DL Model Architectures

In this study, we evaluated a suite of machine learning (ML) and deep learning (DL) architectures for coseismic landslide probability mapping. Two complementary modeling tracks were used: (i) pixel-wise ML baselines (Random Forest and XGBoost) and (ii) patch-based DL segmentation models (CNN, U-Net, and DeepLabV3 with a ResNet-50 backbone). All models were trained using a multi-band raster stack of 14 conditioning factors, including topographic (slope, aspect, curvature, DEM), seismic (PGA), hydrological (drainage density, log-transformed SPI, rainfall), geomorphic (Kₛₙ), InSAR-derived (LOS displacement and DPM), distance to river, and LULC.

2.3.1. Machine Learning Models

Random Forest (RF):

RF was implemented as a per-pixel classifier using the 14-band predictor stack. Pixels containing NaN values in any predictor or in the landslide label were excluded. Because the dataset was highly imbalanced (≈19.2 million non-landslide pixels, 99.6%; 89,416 landslide pixels, 0.4%), majority-class undersampling was applied to enforce a 20:1 negative-to-positive ratio. This yielded 89,048 positive and 1,780,960 negative samples. The balanced data were split 80/20 for training and validation. RF was trained with 200 trees (n_estimators = 200) using inverse-frequency class weighting (class_weight = "balanced"). The model produced per-pixel landslide probabilities in the range (0, 1).

XGBoost:

XGBoost was implemented as a gradient-boosted tree benchmark using the same feature matrix and sampling strategy as RF. Majority undersampling maintained the 20:1 class ratio, followed by an 80/20 train–validation split. The model was trained with 200 boosting rounds (n_estimators = 200). Class imbalance was handled using scale_pos_weight = n(neg)/n(pos) computed from the training data. The final model outputs per-pixel landslide probabilities in (0, 1).

2.3.2. Deep Learning Models

Deep learning models were designed to generate pixel-wise landslide probability maps using spatial context. All networks accept input patches of size (B, 14, 128, 128) and produce full-resolution outputs. The first convolution in each architecture was modified to accommodate the 14-band predictor stack. Because landslide pixels remain sparse within patches, training used a combined weighted loss (binary cross-entropy + Dice, 0.5/0.5) with BCE pos_weight = 15.

CNN:

A lightweight fully convolutional CNN was implemented with four sequential convolutional stages that preserve the original 128×128 resolution (no pooling or downsampling). Each stage uses a ConvBlock with two 3×3 convolutions followed by batch normalization and ReLU activation. Channel depth increases progressively (14 → 64 → 128 → 128). A final 1×1 convolution produces the single-channel logit map. The network contains approximately 857,857 trainable parameters.

U-Net:

A standard encoder–decoder U-Net was implemented for pixel-wise segmentation. The architecture includes four encoder levels, a bottleneck, and a symmetric decoder with skip connections. Using a base width f = 64, the encoder increases feature depth through f, 2f, 4f, and 8f via 2×2 max pooling, with the bottleneck expanding to 16f. The decoder restores spatial resolution using 2×2 transposed convolutions and concatenates encoder features at each scale before applying ConvBlocks. A final 1×1 convolution generates the per-pixel logits. The resulting U-Net contains approximately 31.04 million trainable parameters.

DeepLabV3:

DeepLabV3 was implemented using the torchvision DeepLabV3-ResNet50 model (PyTorch) [52,53]. The pretrained architecture was adapted by modifying the first convolutional layer of the ResNet-50 backbone to accept the 14-band input stack (Conv2d(14, 64, 7×7, stride=2)). The model retains the standard ASPP head and produces a single-channel segmentation map.

Figure 7. Architectures of the deep learning models used for coseismic landslide mapping: (a) lightweight CNN, (b) U-Net encoder–decoder with skip connections, and (c) DeepLabV3 with a modified ResNet-50 backbone and ASPP module, all producing single-channel probability maps.

Figure 7. Architectures of the deep learning models used for coseismic landslide mapping: (a) lightweight CNN, (b) U-Net encoder–decoder with skip connections, and (c) DeepLabV3 with a modified ResNet-50 backbone and ASPP module, all producing single-channel probability maps.

2.4. Training and Evaluation

2.4.1. Data Split and Sampling

Machine Learning Models (Random Forest, XGBoost)

Training samples were derived from the pixel-level predictor stacks and the binary landslide inventory. Because landslide pixels represent approximately 0.4% of the study area, the majority (non-landslide) class was undersampled to reduce bias toward the negative class while maintaining computational tractability. Only negative pixels were undersampled, whereas all positive pixels were retained. The target ratio was 20 negative pixels per positive pixel (e.g., 89,048 positive pixels and 1,780,960 negative pixels in the balanced dataset).

The balanced dataset was divided into training (80%) and validation (20%) subsets using stratified random sampling so that both subsets preserved the 20:1 negative-to-positive ratio. The split was performed with a fixed random seed (random_state = 42) to ensure reproducibility. Performance metrics were computed on the held-out validation subset that was not used during model fitting.

For the pixel-based ML models, stratified random sampling was adopted to enable stable model comparison under severe class imbalance. While spatial blocking is often desirable in landslide studies, the present work focuses on comparative model behavior under identical sampling conditions.

Deep learning segmentation models:

For the convolutional models, image patches were extracted from the predictor stacks using a sliding window approach (patch size = 128 × 128 pixels; stride = 64 pixels, i.e., 50% overlap), yielding dense spatial coverage of the study area. Each patch was assigned a binary label: positive if it contained at least one landslide pixel and negative otherwise.

To mitigate patch-level class imbalance, positive (landslide-containing) patches were oversampled by a factor of two (pos_patch_oversample = 2). The number of negative patches was optionally capped to keep the training set computationally manageable while still ensuring sufficient exposure to positive examples.

Training and validation patch sets were kept geographically independent through a spatial partition. Specifically, patches whose top-left row coordinate fell within the upper 80% of raster rows were assigned to the training set, whereas patches located within the lower 20% were reserved for validation. This north–south hold-out reduces spatial overlap between training and validation data and helps limit optimistic bias due to spatial autocorrelation in convolutional models. If this rule yielded fewer than 50 validation patches, the pipeline reverted to a stratified random patch split to maintain statistical reliability.

Within each patch, pixel-level imbalance was further addressed using a composite loss consisting of weighted binary cross-entropy (BCE) and Dice loss. The positive class weight in BCE was set to 15, and the Dice component was weighted at 0.5 (total loss = 0.5 × Dice + 0.5 × BCE). This formulation encourages sensitivity to rare landslide pixels while improving spatial overlap with the reference mask.

Together, the pixel-based balanced sampling for ML models and the spatially partitioned patch sampling for DL models allow a consistent comparison between paradigms while respecting their different learning units (pixels vs. patches).

2.4.2. Training Configuration

Traditional Machine Learning Models:

Random Forest was trained using bootstrap aggregation with class weighting (class_weight = “balanced”) to compensate for residual imbalance in the 20:1 training data. Standard hyperparameters (e.g., number of trees, depth, and minimum leaf size) followed the implementation settings; no exhaustive grid search was performed in the present experiments.

XGBoost was trained using gradient boosting with class weighting via scale_pos_weight set to the negative-to-positive ratio in the balanced training subset. To improve probabilistic reliability, probability calibration was applied using isotonic regression within a 5-fold cross-validation framework (CalibratedClassifierCV).

Both ML models were trained on the same stratified 80% training subset and evaluated on the identical held-out 20% validation subset.

Deep Learning Segmentation Models:

All three architectures (CNN, U-Net, DeepLabV3) were optimized using the Adam optimizer with an initial learning rate of 1 × 10⁻³ and weight decay of 1 × 10⁻⁴. Training progress was monitored using validation loss (composite BCE–Dice loss on the spatial validation patches). The learning rate was reduced on plateau (ReduceLROnPlateau), and early stopping with a patience of 12 epochs was applied to limit overfitting. The model checkpoint with the lowest validation loss was retained for final evaluation.

2.4.3. Evaluation Metrics

Model performance was assessed using complementary discrimination, detection, and calibration metrics:

• AUC-ROC: Measures overall discrimination across thresholds and is relatively insensitive to class imbalance.
• AUC-PR (Average Precision): More informative for rare-event detection, emphasizing performance on the landslide class.
• Critical Success Index (CSI): Defined as TP / (TP + FP + FN), evaluating the balance between correct detections and false alarms.
• Brier score: Quantifies probabilistic calibration; lower values indicate better-calibrated predictions.
• Confusion matrix (row-normalized, %): Reports True No-LS predicted as No-LS, False alarms (No-LS → LS), Misses (LS → No-LS), and Correct detections (LS → LS).

2.5. Variable and Feature Importance

To identify the most influential conditioning factors controlling landslide occurrence, permutation-based importance analysis was performed for both ML and DL models using the 14 predictor bands (slope, aspect cos/sin, absolute curvature, DEM, Kₛₙ, drainage density, logSPI, PGA, distance to river, rainfall, LULC, LOS_abs, and DPM).

For the ML models (Random Forest and XGBoost), permutation importance was computed using the scikit-learn permutation_importance function with AUC-PR as the scoring metric. Each feature was randomly permuted in the validation set, and the resulting drop in AUC-PR was averaged over five repetitions to reduce variance.

For the DL models (CNN, U-Net, and DeepLabV3), a band-wise (block) permutation approach was used. Each input band was sequentially zeroed in the validation patches, and the corresponding decrease in AUC-PR was recorded. This ablation-based analysis quantifies the relative contribution of each conditioning factor to model performance.

3. Results

3.1. Model Performance Comparison

Model performance was evaluated using AUC-ROC, AUC-PR, CSI, and Brier score. Results are summarized in Table 3 and Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12.

3.1.1. Discrimination and Ranking

ROC curves indicate strong discrimination for all models. The tree-based classifiers achieved the highest AUC-ROC (XGBoost = 0.972; Random Forest = 0.969), followed by DeepLabV3 (0.937), CNN (0.903), and U-Net (0.892). However, because ROC-AUC can appear optimistic under severe class imbalance, AUC-PR provides a more informative measure of landslide detection performance.

3.1.2. Precision–Recall and Threshold Performance

Based on AUC-PR, DeepLabV3 performed best (0.768), followed by U-Net (0.667), CNN (0.664), XGBoost (0.654), and Random Forest (0.643). At the fixed probability threshold of 0.5, DeepLabV3 also achieved the highest CSI (0.49), indicating the most favorable balance between detection and false alarms. Random Forest (0.40) and U-Net/XGBoost (~0.38) showed moderate performance, while CNN had the lowest CSI (0.30) due to higher false positives.

In terms of calibration, the tree-based models produced the lowest Brier scores (XGBoost = 0.025; RF = 0.043), whereas CNN showed the highest Brier score (0.276), indicating less reliable probability estimates.

3.1.3. Confusion Matrix Interpretation

Row-normalized confusion matrices highlight distinct model behaviors:

• Random Forest: Strong detection with moderate false positives.
• XGBoost: Conservative predictions at the 0.5 threshold, resulting in many missed landslides.
• CNN: Very high landslide recall (~94%) but excessive false alarms.
• U-Net: More balanced than CNN but with lower recall.
• DeepLabV3: Best overall trade-off, combining high landslide recall (~90%) with strong non-landslide accuracy (~84%).

3.2. Spatial Probability Maps

The trained models were applied to the full study area to generate pixel-wise landslide probability maps for Random Forest, XGBoost, CNN, U-Net, and DeepLabV3. Each map represents the predicted likelihood of coseismic landslide occurrence (0–1), enabling visual comparison of spatial susceptibility patterns among models. Areas consistently identified as high probability indicate robust agreement on landslide-prone terrain, whereas differences reflect the distinct learning behavior of each architecture. The ML models (RF and XGBoost) produce full-coverage maps at native pixel resolution, while the DL outputs were tiled and mosaicked to obtain continuous probability surfaces.

Figure 11 presents a visual comparison between the zoomed in landslide probability maps generated by XGBoost, Random Forest, U-Net, CNN, and DeepLabV3 and the reference landslide inventory. All models successfully delineate the primary landslide concentration along the valley corridor; however, deep learning models, particularly CNN and DeepLabV3, better capture fine-scale spatial patterns and exhibit closer agreement with the ground truth. Tree-based models produce smoother susceptibility patterns and show relatively higher spatial overprediction in surrounding areas.

3.3. Predictor Importance, InSAR Role, and Kₛₙ Dominance

Permutation analyses (Section 2.5) show that multiple terrain and seismic factors contribute to landslide prediction, including DEM, PGA, slope, and drainage density, although their relative influence varies by model.

3.3.1. Role of InSAR

InSAR-derived variables (LOS displacement and DPM) consistently ranked among the most important predictors, particularly in the DL models. Their inclusion improved model discrimination and calibration, highlighting the value of integrating spaceborne deformation and coherence products for coseismic landslide mapping in high-relief, data-scarce regions such as the Nepal Himalaya.

In the illustrative example, we compare three panels: (a) a probability map generated without InSAR inputs, (b) a probability map generated with InSAR (LOS and DPM) included, and (c) the mapped landslide inventory (ground truth). In the no-InSAR case, the model largely fails to assign elevated probabilities to a known landslide zone, effectively missing an important cluster of failures. When LOS displacement and DPM are added, the same area is highlighted as a high-probability zone, and this pattern closely matches the observed landslides in the ground-truth panel. This contrast demonstrates that, without InSAR products, some coseismic landslide zones remain undetected in the probability maps, whereas incorporating InSAR substantially improves the model’s ability to capture real failures.

Figure 13. Effect of InSAR inputs on landslide prediction: (a) model output without InSAR variables, (b) model output with LOS displacement and DPM included, and (c) reference landslide inventory (ground truth).

Figure 13. Effect of InSAR inputs on landslide prediction: (a) model output without InSAR variables, (b) model output with LOS displacement and DPM included, and (c) reference landslide inventory (ground truth).

3.3.2. Dominance of Kₛₙ

The normalized channel steepness index (Kₛₙ) emerged as the single most influential predictor across all ML and DL models. This result underscores the strong control of fluvial incision and tectonically driven hillslope steepening on earthquake-triggered landslides. High Kₛₙ values reliably delineate over-steepened, failure-prone terrain, supporting its use as a robust and transferable geomorphic indicator in active mountain belts.

To further quantify this relationship, we derived a Kₛₙ-density surface and classified it into five zones (very low, low, medium, high, and very high). We then extracted four south–north profile lines across the study area to examine how mapped landslides are distributed relative to these Kₛₙ zones as shown in Figure 14. The majority of landslides fall within the high-Kₛₙ zone along these profiles, providing additional, spatially coherent evidence that steep, high-Kₛₙ reaches concentrate earthquake-induced slope failures and reinforcing the dominant role of Kₛₙ in controlling landslide probability.

4. Discussion

This study demonstrates that integrating machine learning and deep learning with InSAR-derived products and geomorphic indices provides a powerful framework for event-based coseismic landslide probability modeling in the Nepal Himalaya. The comparative evaluation of Random Forest, XGBoost, CNN, U-Net, and DeepLabV3 shows that all models are capable of reproducing the broad spatial pattern of landslides triggered by the 2015 Gorkha earthquake, but they differ in discrimination, calibration, and spatial detail. In particular, the combination of InSAR line-of-sight (LOS) displacement, coherence-based damage proxy map (DPM), and normalized channel steepness index (Kₛₙ) emerges as a key driver of model performance.

4.1. Model Performance

Tree-based models (Random Forest and XGBoost) achieve the highest AUC-ROC and provide the best-calibrated probabilities. However, when class imbalance is considered, DeepLabV3 performs best based on AUC-PR and CSI, while CNN and U-Net show comparable skill. Spatially, the deep learning models, especially CNN and DeepLabV3, better capture fine-scale landslide clusters and boundaries. In contrast, tree-based models produce smoother maps and slightly overpredict surrounding areas. This reflects their different strengths: tree models excel in feature discrimination, while convolutional networks better exploit spatial context.

4.2. Role of InSAR-Derived Products

LOS displacement and DPM significantly improve model results. When these layers are removed, some known landslide zones are underestimated or missed. Including InSAR data improves both detection and calibration, particularly for DeepLabV3 and U-Net. This confirms that deformation and coherence information provide critical signals that static terrain variables alone cannot capture, especially in steep, data-scarce regions like the central Nepal Himalaya.

4.3. Dominance of Kₛₙ and Geomorphic Controls

Across all models, Kₛₙ consistently emerges as the most influential predictor, highlighting the strong geomorphic control on coseismic landsliding. High Kₛₙ values correspond to steep, actively incising channels and over-steepened hillslopes that are more susceptible to failure during strong shaking.

Density and profile analyses further support this relationship: most mapped landslides occur within high-Kₛₙ zones, and peaks in Kₛₙ density align with landslide clusters. This confirms that Kₛₙ is a physically meaningful proxy for long-term tectonic–geomorphic predisposition.

4.4. Implications, Limitations, and Future Work

The integration of InSAR metrics and Kₛₙ offers a practical pathway for regional hazard screening in data-poor mountain belts using only SAR and DEM data. Such event-based models can support rapid post-earthquake prioritization and corridor-scale risk assessment.

Nevertheless, several limitations remain. The models are trained on a single event (2015 Gorkha), and temporal transferability is untested. Uncertainties in the landslide inventory, InSAR processing, and DEM-derived Kₛₙ propagate into predictions. Deep learning models also require greater computational resources and careful handling of class imbalance.

Future research should therefore focus on (i) validating the models in other segments of the Nepal Himalaya and in different tectonic and climatic settings, (ii) assessing temporal transfer by applying the trained models to landslides from subsequent earthquakes, (iii) exploring additional InSAR time-series metrics and multi-orbit combinations, and (iv) systematically evaluating alternative geomorphic indices and Kₛₙ calculation schemes. Integrating physics-based ground motion and slope stability models with data-driven probability mapping also represents a promising direction for improving interpretability and robustness.

5. Conclusions

This study compared five ML and DL models for coseismic landslide probability mapping in the 2015 Gorkha earthquake using a 14-band predictor stack that includes InSAR LOS displacement, DPM, and DEM-based geomorphic indices such as Kₛₙ. All models capture the overall spatial pattern of landslides, but their performance characteristics differ. Tree-based models (RF and XGBoost) exhibit the highest AUC-ROC and best calibration (lowest Brier scores), while DeepLabV3 achieves the highest AUC-PR and CSI, reflecting superior performance under severe class imbalance and better spatial delineation of landslide clusters.

Deep learning segmentation models produce sharper and more spatially coherent probability maps, whereas tree-based models generate smoother but slightly overpredicted susceptibility fields. InSAR-derived LOS and DPM substantially improve discrimination and calibration across all architectures and are especially valuable for deep learning approaches that exploit spatial context in multi-band predictor stacks. Across models, the normalized channel steepness index (Kₛₙ) consistently emerges as the dominant predictor, confirming the strong geomorphic and tectonic control on earthquake-triggered landslides.

Overall, our results indicate that integrating InSAR products with Kₛₙ and related DEM-derived metrics provides a transferable, remotely sensed framework for coseismic landslide probability mapping in high-relief, data-scarce regions. Ensemble tree models offer efficiency, robustness, and interpretability for regional screening, while deep learning architectures add fine-scale spatial detail at higher computational cost. Further validation across different earthquakes and regions is needed to confirm the broader applicability of this framework and to refine best practices for combining InSAR, geomorphic metrics, and ML/DL methods in landslide hazard assessment.

Appendix A.1.

Table A1. Monthly rainfall data (mm) used to generate the graph in Figure 5b.

Monthly Rainfall in mm
Month	Gorkha	Dhading	Rasuwa	Nuwakot	Kathmandu	Lalitpur	Sindhupalchok	Kavre	Dolakha	Ramechhap
Sep-14	198.72	205.7	132.2	524.6	325.6	118.9	437	150.9	82.05	126
Oct-14	53.25	100.2	0.2	104.6	87.2	121.3	89.625	64.45	28.35	45.5
Nov-14	0	0	0	0	0	0	0	0	0.5	0
Dec-14	40.825	82.7	39.5	28.2	30.2	22.7	44.725	18.4	13.65	8.5
Jan-15	82.88	13.5	16.4	5	15.8	8.6	22.375	18.05	0.05	1
Feb-15	66.2	62.8	4.5	51.8	44.8	31.5	39.45	24.6	8.85	20
Mar-15	93.3	31.8	103.8	91.4	90.1	78.9	129.625	65.35	35.15	45.5
Apr-15	30.56	49.3	31	68	8.8	49.7	54.675	49.25	58.1	52.3
May-15	36.425	61.2	8.8	128.6	10	36	133.7	47.8	21.35	26.3
Jun-15	181.8	246.7	118.6	198.6	317.8	180.7	319.1	47.5	62.15	63.5
Jul-15	416.8	359.8	133.3	778	556.4	299.8	480.275	294.45	212.55	267.5
Aug-15	312.26	181.2	281.7	775.8	692.8	280.3	612.55	257.5	159.3	158

Table A1. Monthly rainfall data (mm) used to generate the graph in Figure 5b.

Monthly Rainfall in mm
Month	Gorkha	Dhading	Rasuwa	Nuwakot	Kathmandu	Lalitpur	Sindhupalchok	Kavre	Dolakha	Ramechhap
Sep-14	198.72	205.7	132.2	524.6	325.6	118.9	437	150.9	82.05	126
Oct-14	53.25	100.2	0.2	104.6	87.2	121.3	89.625	64.45	28.35	45.5
Nov-14	0	0	0	0	0	0	0	0	0.5	0
Dec-14	40.825	82.7	39.5	28.2	30.2	22.7	44.725	18.4	13.65	8.5
Jan-15	82.88	13.5	16.4	5	15.8	8.6	22.375	18.05	0.05	1
Feb-15	66.2	62.8	4.5	51.8	44.8	31.5	39.45	24.6	8.85	20
Mar-15	93.3	31.8	103.8	91.4	90.1	78.9	129.625	65.35	35.15	45.5
Apr-15	30.56	49.3	31	68	8.8	49.7	54.675	49.25	58.1	52.3
May-15	36.425	61.2	8.8	128.6	10	36	133.7	47.8	21.35	26.3
Jun-15	181.8	246.7	118.6	198.6	317.8	180.7	319.1	47.5	62.15	63.5
Jul-15	416.8	359.8	133.3	778	556.4	299.8	480.275	294.45	212.55	267.5
Aug-15	312.26	181.2	281.7	775.8	692.8	280.3	612.55	257.5	159.3	158