The Hindu Kush Himalaya (HKH) is one of the most disaster-prone regions, where active seismic zones, extreme topographic gradients, seasonal monsoon-driven precipitation, and fragile geological formations intersect, creating a highly dynamic risk landscape [1]. Nepal, which lies in the HKH region, frequently experiences cascading natural hazards, particularly floods and landslides, that often interact and trigger chain reactions from upstream to downstream. The causes of these events include underlying with young mountainous terrain, unplanned urbanization, poor watershed management, and high frequencies of earthquakes, which create conditions for natural disaster.

The cascading effects, however, can be triggered by several factors, which might be time-dependent or occur simultaneously. For instance, intense monsoon rainfall can trigger flooding and landslides concurrently. On the other hand, time-dependent causes include wildfires that destabilize slopes, often resulting in debris flow [2], and earthquakes that can directly trigger landslides or Glacier Lake Outburst Flood [3]. The Melamchi flood (2021) mainly due to a cascading chain reaction[4] devastating the impact. These examples highlight that disasters are no longer defined by natural extremes but by the complex interactions among multiple hazards, social vulnerability and systematics governance gaps [5]. Since floods and landslides can occur in the same region due to the same triggering factor, this study seeks to identify the overlay and flow zones. Decision makers today not only face challenges of mitigating single hazards but they need to respond to the chains of cascading events which require understanding of systematic interaction and interrelationship. The cascading effect concepts are relatively new in natural risk governance, with limited methodologies and field applications, yet it is critical in understanding cumulative damage, where multiple hazards affect the same element of risk (people, infrastructure and economy). In a simple term, cascading events can be understood as a snowball process, in which an initial trigger evolves into a series of adverse outcomes, each branch of the event tree represents a cause-and-effect chain that amplifies vulnerability [6]

Despite this evolving risk landscape, Nepal lacks an integrated national scale understanding of where floods and landslide intersect as cascading hazards, limiting long-term planning for infrastructure, settlements, agriculture and emergency preparedness. Mountain settlement in Nepal is particularly vulnerable due to their location on steep slopes, proximity to river and glacier lakes, and dependence on limited infrastructure. Identifying zones where flood and landslide interact is therefore critical for ensuring the safety and resilience of these communities, guiding sustainable infrastructure development and informing risk informed land use strategies. This will also directly support global frameworks including the Sendai Framework for Disaster Risk Reduction and Sustainable development Goals (SDG11), which emphasizes understanding risk, mitigating risk, evidence-based planning, early warning systems and resilient infrastructure development[7, 8] required for sustainable human settlement. Although Nepal has advance national risk management strategies[9], the cascading hazard dimension remains largely unexplored.

Previous studies in Nepal have primarily employed Geographic Information System (GIS)-based multi-criteria analysis (MCA) to assess individual hazards, relying on expert-defined weights for variables such as slope, land cover, elevation, and other parameters [10–12]. Although these methods suffer from subjectivity and limited transferability, this was only the solution for mapping susceptibility. With the advancement in machine learning techniques, there are several methodologies focusing on the linear and non-linear relationships of the parameters. This study applied the Random Forest model, which offers a more data-driven weightage, which will overcome the problem of human bias by automatically ranking variable importances[13, 14] to generate national-scale susceptibility maps for both floods and landslide using recognized drivers including drivers include topographical, hydrological, land characteristics and vegetation characteristics[15, 16] and also applied unsupervised machine learning to identify the cluster for both flood and landslide. To our knowledge, this is the first national-level assessment in Nepal that combines machine learning based susceptibility modeling for floods and landslides to identify their cascading interaction zones to assess risks to mountain settlement.

The objective of this study is therefore to:

Nepal, with its diverse topography, complex geography, and highly varying climate, provides an ideal study ground for risk susceptibility mapping. Globally, Nepal is ranked as the 4th most climate-vulnerable nation and among the top 20 countries that face multiple hazards, both natural and human-induced [17]. The figure of loss and damage reflects that an average of two people died, and more than 300 families have been affected by all kinds of disasters in the past 50 years [9]. The past semicentennial analysis reflects that floods and landslides are major risk concerns in the Nepalese context. The figure for the loss of life due to landslides is an average of 113 deaths each year [18]. Floods cause 9.33% of total disaster-related deaths and affect 61.60% of families [9].

The methodology is divided into two sections. First, flood susceptibility maps were generated using a Random Forest regression model integrated with various flood related variables. Similarly, landslide susceptibility was modeled using dependent variables and validated against the ground truth samples and existing inventories. Second, an unsupervised machine learning model was employed to cluster the flood and landslide susceptibility data. Finally, a detailed impact assessment was conducted for single hazards, cascading overflow area as shown in Fig. 2.

After compiling all the required data, a stratified random sampling method was used to ensure a more accurate representation of the entire population by dividing it into distinct subgroups using Eq. 1, thereby reducing overfitting and spatial autocorrelation [29]. This sampling strategy was designed to provide a 99.9% confidence level, using Eq. 2 [30].

eq1

Here is the sample size, is the confidence level, and the population of the partition to be used. is the margin of error. The adjusted sample size formula is used when the population size is known:

eq2

Here is the sample size calculated from Eq. 1, N is the population size and

The total population size represented by the image pixels, which cover the whole of Nepal at a resolution of 30 meters, is approximately thirty-seven million. Using Eq. 1 and Eq. 2 a total number of samples required was found to be 1600 samples, of which 70% of the data was used for training purposes and 30% for testing the models.

The flood and landslide inventory were then transformed into point datasets. In the case of stratified random sampling, we created a grid with a resolution of 400m x 400 m, which spread the dataset all over the study area at least 400 meters apart. It was used to minimize spatial autocorrelation.

This study used a random forest regression model for susceptibility mapping of floods and landslides using GEE. The small sample size was used for hyperparameter tuning between the number of trees for prediction from 400–600, and it was found that the optimal number of trees for best model accuracy was 600. The model deals with binary data, which is flooded and non-flooded points or landslide and non-landslide points and is implemented using the bootstrapping technique to develop multiple decision trees [13]. The predictions from individual trees are averaged, enabling the model to effectively capture the relationship between the input variables and the target phenomenon. The model produces continuous outputs ranging from 0 to 1, representing susceptibility scores using Eq. 3.

Mathematically, it is expressed as:

eq3

Here, is the estimated probability distribution

is the number of trees

is the prediction of b^th individual trees.

We integrated various environmental parameters (features) for flood and landslide susceptibility mapping. Each feature was used to train the model, where multiple decision trees were constructed from bootstrap samples of the training data. At each node of the tree, a random subset of features was selected, and among them, the best split features were used. Randomization reduces the correlation between trees and decreases the model’s variance.

The importance of each parameter in the model was determined using the permutation importance method as shown in Eq. 4. This method determines the importance of each parameter based on the model’s ability to predict output in terms of an increase in error metric when the feature values are randomly shuffled [13]

Mathematically, the importance of a feature is computed as:

eq4

Where:

is the importance of score of feature j

is the mean square error of the model

is the mean square error after shuffling the values of the feature

Accuracy Assessment

The accuracy of the model was evaluated using the Receiver Operating Characteristics (ROC) curve, which is a graphical representation that illustrates the model’s ability to distinguish between actual events and non-events[31]. The curve is plotted between the true positive rate (TPR) and the false positive rate (FPR). The TPR, also known as sensitivity, is calculated as:

Eq. 5

Eq. 6

In these formulas, TP represents True positives, FN represents False Negatives, FP represents False Positives, and TN represents True Negatives [32]. The Area Under the Curve (AUC) is a single number summary of the overall performance of the ROC curve, ranging from 0 to 1, where 1 denotes a perfect classifier and a value of 0.5 denotes a worthless classifier.

Published data was used for validation. This involved several tests between the two results, which include Root Mean Square Error (RMSE), Bias, and Index of Agreement (see Supplementary 2).

The RMSE measure provides the standard deviation of the residuals, indicating how this product varies.

² Eq. 7

The bias measures the average error of the two results, indicating the accuracy of the results. It helps to identify which results are overestimated.

Bias b= Eq. 8

The index of agreement (D) is used to compare the similarities between two sets of data, indicating reliability.

Index of agreement

Eq. 9

Here is the total number of observations. is the data from susceptibility is the validation data from a different source and is the mean of the observed data.

Unsupervised machine learning involves algorithms finding patterns in unlabeled data without a target outcome. The goal is to discover inherent data structure through methods like clustering or dimensionality reduction. The flood susceptibility above 0.6 and landslide above 0.8 was applied to identify the clustered where flood and landslide events might occur separately and in clustered. In this study K-means clustering was applied which relies on minimizing the object function which represents the within-cluster sum of squares:

Eq. 10

This equation sums the squared Euclidean distance between each grid cell, and it’s assigned the clustered centroid optimally grouping areas with similar hazard susceptibility [33]. This highlights three zones safe zone, flood risk zone and landslide risk zone. The landslide risk zone was further classified into two landslides occurring alone and flood and landslide interaction zone using the equation below:

Let, is Flood susceptibility above 60% and landslide susceptibility area above 80%

Eq. 11

indicates that the landslide susceptibility area does not overlap with any flood susceptibility area A

means the boundaries of the landslide susceptibility area and the flood susceptibility area A intersects with that which is spatially adjacent.

The susceptibility thresholds for defining high risk flood and landslide classes were selected based on the empirical comparison with historical hazard event datasets. The multiple threshold values were tested and overlaid with past events. A flood susceptibility value of > 0.60 showed the highest spatial correspondence with known flood with cascading events. Similarly, landslide susceptibility threshold of > 0.80 provided the best agreement with mapped landslide scars available in supplementary (S15).

Flood risk susceptibility model was assessed using a Random Forest model incorporating eleven key conditioning factors. The parameters were ranked by importance accessibility > flow accumulation > landform > landcover > slope > hydrological conductivity > flow direction > terrain roughness index > elevation > terrain wetness index > HAND as shown in supplementary S9 a). The overall accuracy of flood susceptibility model with validation dataset is 0.84 as shown in ROC curve supplementary S8 a). The susceptibility data was further validated with hydrological model dataset which was developed by METROR’s 100- year return period flood data which have index of agreement of 74% with RMSE of 0.48 and positive bias of 0.18 as shown in supplementary S 10 a). The positive bias indicate that susceptibility data outperform than METROR dataset. This might be due to flooded and non-flooded training dataset capture the latest incidents where the METROR dataset only used DEM, Rainfall and landcover as primary input parameter.

Approximately 19% of Nepal falls into medium to very high susceptible zones, predominantly in the southern low land regions Fig. 3a). Exposure analysis reveals significant vulnerability: nearly 900 thousand people and over 3.4 million infrastructures are in the medium to very high-risk zones as shown in supplementary ST1.

Landslide susceptibility was assessed using a Random Forest model incorporating twelve key conditioning factors. These variables are ordered by their importance for prediction (slope > aspect > NDVI > DEM > sand > silt > distance from roads and rivers > TWI > landform > precipitation > land cover), are recognized indicators of slope stability as shown in supplementary S9 b). The model achieved robust performance with an Area Under the Curve (AUC) value of 0.85 as shown in ROC curve supplementary S8 b). The susceptibility data strong performance with manually digitized 2015 earthquake induced landslide inventory with index of agreement of 81%, a small negative bias of (-0.18) and an RMSE of 0.43 as shown in supplementary S 10 b). The negative bias may be attributed to the inclusion of NDVI as a predictor, which reflect vegetation cover. Nevertheless, the overall agreement with observed landslide and robust discrimination ability indicate that model provides reliable and conservative estimates for landslide susceptibility.

Approximately 40% of the total area fall in medium to very high-risk zones predominantly in the northern part of Nepal Fig. 3b). Exposure analysis indicates that two hundred thousand individuals and six hundred thousand infrastructures fall within medium to very high landslide susceptibility zones, revealing localized but critical vulnerabilities, particularly in mountainous and hilly regions of Nepal as shown in supplementary ST2.

Both flood and landslide susceptibility data were used as a primary parameter to develop the clustered area which highlights the low hazard zone, flood zone, and landslide zone. The flood risk zones were observed mostly in southern part of Nepal as shown in Fig. 4D), however landslide zone predominately observed in the northern part of Nepal Fig. 4C). The cluster of both flood and landslide was observed in the middle part of Nepal as shown in Fig. 4A) and B). In a country scale, the total area of low hazard zones includes 81%, flood risk zone 9%, and landslide risk zone is 10%.

To understand the complex interactions between the flood and landslide, we classified clustered analysis results into four hazard zones using Eq. 11. 81% of the area was found to be low hazard zone, 9% high flood only risk zone, 5% both cascading zones and high landslide only risk zone (supplementary S18). The low hazard zone exhibits minimal susceptibility to both hazards (flood threshold < 0.6, landslide threshold > 0.8) and spans most of Nepal. High flood zones occur primarily in the south regions and river valleys, while landslide-only zones are spatially isolated but often parallel to flooded prone area. Finally, cascading zones are defined as areas of high landslide probability immediately adjacent to flood-prone areas, highlighting locations where landslides could trigger secondary flooding by blocking river channels.

To delineate cascading landslide hazard impact, we derived the cascading flood hazard zone defined as flood-prone area that intersect with high probability landslide zones. The runoff from these cascading landslide zones was modeled to extend downstream until the nearest settlement, capturing potential impact to communities. If the cascading flood encounters another landslide prone area along the way, the flow is further propagated to the next downstream settlement, iteratively accounting for compounded hazard interactions. With this approach it was found at cascading flood zone covers 7588 km² impacting 88 km² of built-up areas (12% of the total) and 1722 km² of cropland, underscoring the exposure of settlement and agricultural land to compounding hazards. The spatial patterns show minimal cascading effects in the central region Fig. 5E), where stable geomorphic conditions and gentle slopes limit simultaneous hazard occurrence. In contrast, steep terrain with dense river networks and fragile landcover ( Fig. 5B, C,D,F) exhibits higher likelihood of compounding events.

This study provides first national-scale, machine-learning based cascading flood-landslide zonation in Nepal. By integrating Random Forest model with large geospatial dataset and unsupervised clustering, the analysis presents a more robust, objective and reproducible understanding of the country’s multi-hazard landscape compared to traditional, expert-driven methods.

The RF model demonstrated high predictive accuracy, benefiting from the hyperparameter tuning and large-scale computation offered by Google Earth Engine (GEE). Compared to traditional approaches such as Multi-Criteria Decision Analysis (MCDA) or Analytical Hierarchy Process (AHP), which are limited in handling numerous parameters and complex pairwise comparisons [34]. The machine learning approach allowed the incorporation of 11 variables for floods and 12 for landslides. This increased the robustness, objectivity, and generalizability of hazard predictions, a trend corroborated by other geospatial applications like wetland mapping and temperature estimation [35]. The RF models automated variable importance ranking reduces human bias and improves transferability to other regions within Nepal or the border Himalayan.

The susceptibility maps reveal distinct spatial patterns of flood and landslide risks in Nepal. High flow accumulation areas significantly increase flood risk in the southern part of Nepal by concentrating runoff [36], whereas high hydrological conductivity exhibits lower flood susceptibility as soils and vegetation improve water absorption increasing threats to both population and infrastructure. Steep slopes, loose terrain, and weathering processes make northern hilly and mountainous regions prone to landslides.

Anthropogenic activities like road construction, land use and land cover change, slope modification, riverbank encroachment and deforestation further intensify susceptibility[37, 38]. Expanding settlements into flood plains and agricultural intensification in the southern part of Nepal, exposure, while in the central region hill cutting and road expansion in northern part mountainous regions weaken slope stability. These anthropogenic drivers highlight the urgency of integrating geospatial evidence into land use decisions.

The susceptible data was further validated with hydrological model dataset and earthquake induced landslide which contain positive and negative bias. This highlighted a crucial implication for disaster management models that should be updated with recent inventories to provide more accurate and timely risk assessments,

The use of unsupervised machine learning like kmeans clustering algorithms identifies natural groupings in this multidimensional data space. The landslide susceptibility and flood susceptibility were combined based on its characteristics resulting in distinct zones flood susceptible, landslide susceptible zones and low hazard zones. This provides clear, data driven visualization of risk distribution. This visualization is a powerful tool for sustainable mountain settlement planning. Decision makers can use the low hazard zones or no interaction zones for developing new settlements and infrastructure, as these areas inherently possess combinations of factors that correlate with low risk. The flood risk zone can be targeted for specific mitigation strategies such as implementing structural defenses, enforcing building codes, or focusing on environmental measures. Similarly, for the landslide risk zone can be targeted for specific mitigation strategies such as planting the forest in slope area.

A key advancement of this study is the identification of cascading hazard zones, where floods and landslides susceptibility co-occur or are spatially adjacent. Such areas are susceptible to multi-hazard chains, including landslide dammed lakes, sudden outburst floods, and slope failures triggered by prolonged inundation. While previous studies in the region typically addressed these hazards in isolation [12, 39–41], our methodology provides a data driven framework for spatially explicit interaction mapping. These areas represent potential “snowball” processes [42], in which an initial hazard triggers subsequent events. For example, a landslide may block a river, forming a temporary dam that can fail suddenly and generate an outburst [43]. Such a cluster was identified in the northeast -to-northwest part of Nepal. Conversely, intense flooding can satisfy slopes, increasing landslide probability. Such interactions are particularly relevant in Nepal’s monsoon season, when tectonics, land-use changes, prolonged droughts, and forest fires can amplify cascading risks [44, 45].

The cascading hazard maps have direct relevance for Nepal development planning, climate adaptation and disaster risk reduction agendas. Low-hazard zones offer opportunities for safer settlement expansion and infrastructure placement, high flood susceptible zones require stricter zoning, flood resistant designs and improved watershed and wetland management. Landslide dominant areas necessitate slope stabilization, bioengineering and controlled infrastructure development. Cascading hotspots are particularly important for hydropower planning, road alignment and community based early warning systems. The cascading hotspots area needs to redetermine the threshold for early warning considering the damming of river. Incorporating this map into local government planning process avoid high-risk areas, enhance preparedness and support evidence-based investment aligned with Sendai Framework and Sustainable cities and communities. In flood-prone southern regions, adherence to flood-resistant building codes, such as FEMA’s BRIC principles, can enhance infrastructure resilience[46]. In northern areas, it is vulnerable to landslides, policies promoting reforestation and restricting hill cutting can stabilize slopes.

Furthermore, the maps guide risk-sensitive land use planning. Sustainable practices such as terracing, contour farming, wetland restoration, and controlled levees (Kasumi-tei) can mitigate cascading hazard impacts [47–50]. Restricting new construction in high-risk zones or relocating settlements transform hazard-prone areas into safer land uses, such as forests or terrace fields. Integration of susceptibility maps with policy tools allows evidence-based decision-making to enhance resilience and reduce socio-economic losses.

The methodology developed in this study is easily transferable to other multi-hazard contexts because it relies primarily on freely available datasets and is implemented using open-source tools. This makes the framework adaptable for researchers aiming to investigate additional hazard interactions such as wildfire–landslide, drought–flood, earthquake–landslide, or avalanche–GLOF cascading processes. However, applying this workflow to new regions or hazard types requires careful calibration of susceptibility thresholds based on local environmental conditions, hazard characteristics, and inventory quality. Developing region-specific and process-specific thresholds will be essential to ensure accurate, context-driven multi-hazard assessments.

Although this methodology provides valuable national-scale cascading zone mapping and is suitable for high-level strategic planning. This study does not incorporate temporal dynamics which are essential for developing scenario based multi-hazard modeling. The current methodology based on susceptibility mapping identifies where interactions are likely to occur but does not predict when or how they will evolve over time. Therefore, while the national-scale data guides broad policy, future studies should focus on integrating temporal dynamics and developing specific landslide hazard scenarios. This work could involve incorporating time-series rainfall data, citizen reported inventories, geomorphological field validation and dynamic simulation model to generate scenario-based cascading hazard forecasts, including potential one-day to multi-day river blockage breach sequences.

For the first time, this study presents a high-resolution, national-scale mapping of cascading flood and landslide risk susceptibility in Nepal using machine learning and the cloud-based Google Earth Engine platform using open datasets. By integrating diverse geospatial, climatic, soil, and geological datasets along with historical hazard inventories, the Random Forest model demonstrated strong predictive performance while efficiently processing large and complex data. This study identified distinct spatial patterns, with flood risk primarily concentrated in the southern plains and landslide susceptibility highest in the northern hilly and mountainous region, reflecting differential impacts across Nepal’s varied terrain. Importantly, this study identifies significant cascading hazard zones where flood and landslide risks intersect and interact, highlighting the critical need for multi-hazard approaches in disaster risk management. The cascading high risk zones provide valuable inputs for advanced hydrodynamic and dam breach modeling, early warning systems, and disaster preparedness planning to design resilient infrastructure and human settlements. Despite these advances, the challenges remain in estimating landslide volumes and modeling secondary impacts such as sediment transport and flood wave dynamics. These aspects need further research. Future directions include dynamic multi-hazard simulation frameworks, integrating real-time hydrometeorological data, and establishing open-access databases of cascading hazard events. This novel framework could go beyond isolated hazard assessments toward integrated risk management, directly supporting climate adaptation, disaster resilience, and planned settlement development within Nepal, and scale out to other Hindu Kush Himalaya countries and beyond.