An integrated multi-hazard assessment using machine learning in the complex terrains of Northern Pakistan.

Skip to main content An official website of the United States government Here's how you know Here's how you know Official websites use .gov A .gov website belongs to an official government organization in the United States. Secure .gov websites use HTTPS A lock ( Lock Locked padlock icon ) or https:// means you've safely connected to the .gov website. Share sensitive information only on official, secure websites. Search Log in Dashboard Publications Account settings Log out Search… Search NCBI Primary site navigation Search Logged in as: Dashboard Publications Account settings Log in Search PMC Full-Text Archive Search in PMC Journal List User Guide PERMALINK Copy As a library, NLM provides access to scientific literature. Inclusion in an NLM database does not imply endorsement of, or agreement with, the contents by NLM or the National Institutes of Health. Learn more: PMC Disclaimer | PMC Copyright Notice Sci Rep . 2026 Mar 1;16:11499. doi: 10.1038/s41598-026-41029-w Search in PMC Search in PubMed View in NLM Catalog Add to search An integrated multi-hazard assessment using machine learning in the complex terrains of Northern Pakistan Muhammad Ismail Khan Muhammad Ismail Khan 1 National Centre of Excellence in Geology, University of Peshawar, Peshawar, Pakistan Find articles by Muhammad Ismail Khan 1, ✉ , Muhammad Shafique Muhammad Shafique 1 National Centre of Excellence in Geology, University of Peshawar, Peshawar, Pakistan 2 GIS and Space Application in Geosciences (GSAG) Lab, National Center of GIS and Space Application (NCGSA), Islamabad, 44000 Pakistan Find articles by Muhammad Shafique 1, 2 , Ghazanfar Ali Khattak Ghazanfar Ali Khattak 1 National Centre of Excellence in Geology, University of Peshawar, Peshawar, Pakistan Find articles by Ghazanfar Ali Khattak 1 , Abid Yahya Abid Yahya 3 Department of Electrical and Communications Systems Engineering, Botswana International University of Science and Technology, Palapye, Botswana Find articles by Abid Yahya 3 , Irfan Anjum Badruddin Irfan Anjum Badruddin 4 Department of Mechanical Engineering, College of Engineering, King Khalid University, Abha, 61421 Saudi Arabia Find articles by Irfan Anjum Badruddin 4 Author information Article notes Copyright and License information 1 National Centre of Excellence in Geology, University of Peshawar, Peshawar, Pakistan 2 GIS and Space Application in Geosciences (GSAG) Lab, National Center of GIS and Space Application (NCGSA), Islamabad, 44000 Pakistan 3 Department of Electrical and Communications Systems Engineering, Botswana International University of Science and Technology, Palapye, Botswana 4 Department of Mechanical Engineering, College of Engineering, King Khalid University, Abha, 61421 Saudi Arabia ✉ Corresponding author. Received 2025 Mar 22; Accepted 2026 Feb 17; Collection date 2026. © The Author(s) 2026, corrected publication 2026 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/ . PMC Copyright notice PMCID: PMC13056936  PMID: 41765999 Abstract The increasing frequency of natural hazards, intensified by climate change, poses substantial challenges to sustainable development worldwide. Northern Pakistan, particularly the Hunza district, is highly susceptible to multiple hazards, including landslides, earthquakes, glacier-induced floods, debris flows, and Glacier Lake Outburst Floods (GLOFs), driven by both climatic and tectonic factors. A multi-hazard assessment is essential to understand the complex interactions between these hazards, offering a comprehensive perspective on risk and facilitating more effective disaster preparedness and mitigation strategies. This study addresses the existing gap in multi-hazard assessments, which are often confined to single-hazard evaluations, by developing an integrated multi-hazard susceptibility map for the Hunza district in Northern Pakistan. The region’s complex topography, active tectonics, and accelerated glacier melting contribute to its high vulnerability to cascading and co-occurring hazards. The integrated assessment utilizes diverse data sources, including topographic attributes, geological, hydro-meteorological, environmental variables, and literature-derived hazard map for multi-hazard susceptibility analysis. A Machine Learning (ML) Forest-Based Classification and Regression (FBCR) model, Analytical Hierarchy Process (AHP), and Vs30-based site characterization was employed to classify and generate hazards individually and as integrated multi-hazard susceptibility map. The model incorporates eighteen geo-environmental variables for individual hazards assessment. The resulting multi-hazard susceptibility map indicates that 23.11% of the area is prone to landslides, 6.07% to flash floods, 4.66% to debris flows and flash floods, and 3.98% to a combination of flash floods, landslides, and debris flows. The highest multi-hazard zone, comprising seismic hazard, debris flows, landslides, and flash floods, covers 2.88% of the area, whereas low-hazard zones constitute 56.84% of the region. The landslide susceptibility model classifies 20% of the area as very high susceptible, while the flash flood, debris flow, and seismic hazard models indicate 5, 2, and 13% of the area, respectively, fall under very high susceptibility/hazard. This integrated multi-hazard approach provides a comprehensive risk assessment framework, supporting evidence-based disaster risk reduction policies and infrastructure planning in hazard-prone regions. The findings identify critical high-hazard zones, offering data-driven insights for targeted mitigation strategies and disaster risk reduction efforts. Keywords: Multi-hazard assessment, Machine learning, AHP, Northern Pakistan, Risk assessment Subject terms: Climate sciences, Environmental sciences, Environmental social sciences, Hydrology, Natural hazards Introduction The compounded effects of multiple hazards, intensified by their increasing frequency under climate change, pose substantial challenges to sustainable development worldwide. Climate change has altered environmental patterns, leading to erratic precipitation, temperature fluctuations, and accelerated glacier melting. These changes have intensified the frequency and magnitude of natural hazards such as landslides, debris flows, floods and Glacial Lake Outburst Floods (GLOF), resulting in widespread damage to communities, infrastructure, and economies 1 . Dilley et al. 2 reported that globally around 30 million km 2 of land and over 4.3 billion people are exposed to at least one natural hazard, 3.5 million km 2 area and 700 million people are exposed to two hazards and 0.5 million km 2 area and 105 million people are exposed to three or more hazards. Despite the significant risk from multiple hazards, hazard assessment studies have predominantly focused on single hazards, thereby underestimating the overall risk posed by multi-hazard scenarios, which may co-occur or trigger cascading effects 3 . Moreover, given the population growth, urbanization and deteriorating economic conditions particularly in developing countries, the vulnerability of infrastructure and already marginalized mountain communities is on the rise 4 . Consequently, multi-hazard risk assessments are critical and gaining momentum in hazard-prone countries to utilize for developing disaster risk reduction policies and measures 5 . Within this global context, Pakistan ranks among the most disaster-prone countries in South Asia, having incurred substantial social, physical, environmental, and financial losses over recent decades 6 , 7 . Given its complex tectonic, topographic and climatic conditions, the country has experienced and is still exposed to multiple hazards, including landslides, debris flows, floods, GLOFs, and earthquakes 8 . Intensifying climate variability, accelerating glacier retreat, and unplanned development within hazard-prone terrain have further amplified exposure, vulnerability and risk 7 , 9 – 11 . Between 2005 and 2025, more than 91,000 people lost their lives to disasters across Pakistan, resulting in economic damages exceeding $34 billion-of which approximately $10 billion stemmed from the 2010 floods and $5 billion from the 2005 earthquake and associated landslides 12 , 13 . In particular, northern Pakistan bears the most direct and severe impact of these climatic and tectonic stresses. During the monsoon season (July-September), intense precipitation and high temperatures frequently trigger multiple hazards throughout the Hindukush-Karakoram-Himalaya (HKH) mountain ranges, causing extensive human losses, infrastructure destruction, and socio-environmental degradation 14 , 15 . Rising temperatures-averaging an increase of about 0.85 °C since 1980-have shifted monsoon patterns northward and westward in Pakistan and accelerated glacier melt, amplifying the frequency and magnitude of extreme events. Topographically, deeply incised valleys, steep slopes, and extensive glaciofluvial and alluvial landforms enhance susceptibility to slope failures, while active fault systems generate frequent small- to large-magnitude earthquakes that often trigger cascading landslides and debris flows. A prominent example is the 2010 co-seismic Attabad rockslide, which dammed the Hunza River to form a 27 km-long lake, resulting in 19 deaths, the displacement of over 6,000 people, and the loss of vital infrastructure and agricultural land. Moreover, the region is highly exposed to flash and debris floods caused by accelerated melting of glaciers and GLOF 16 , 17 . Recent GLOF occurrences in June 2019, April 2020, May 2022 and August 2025, triggered by the Shisper Glacier in Hunza district, caused extensive damage to downstream communities, infrastructure, and the Karakoram Highway. These hydro-meteorological and geohazards pose ongoing challenges to multi-billion-dollar infrastructure projects under the China-Pakistan Economic Corridor (CPEC), leading to increased costs, delays, and disruptions to the sustainable supply chain along the Karakoram Highway—the only land route between China and Pakistan 18 . Given the intense concentration of hazards in Northern Pakistan and their cascading impacts on critical infrastructure, there is a growing need for comprehensive multi-hazard assessment frameworks capable of capturing these interconnected processes. However, the varying nature, impacts, causative and triggering factors, and spatio-temporal dynamics of different disasters make their integration into a single multi-hazard assessment highly complex and challenging 19 . Several methodologies have been developed for integrated multi-hazard mapping using multivariate statistical and decision-support approaches, such as logistic regression and the Analytical Hierarchy Process (AHP) 20 – 22 . More recently, machine-learning (ML) techniques, including tree-based algorithms such as Random Forest, have been increasingly adopted to capture complex and nonlinear relationships between conditioning and triggering factors in multi-hazard assessments 23 . These ML approaches have advanced susceptibility mapping by improving classification accuracy, handling large datasets efficiently, and facilitating more realistic and automated hazard assessments 24 . However, Northern Pakistan has witnessed and been exposed to multiple hazards; therefore, an integrated multi-hazard assessment is urgently required to capture the region’s compounded risk environment 15 . Therefore, this study aims to fill this crucial gap and develop an integrated multi-hazard map of the major hazards (flash flood, landslide, earthquake, debris flow, and GLOF) in the Hunza district of northern Pakistan, using a machine learning algorithm to contribute to evidence-based policymaking for effective disaster risk reduction. Materials and methods Description of study area The study area encompasses the Hunza district in northern Pakistan, a region characterized by its distinct topography, active tectonic settings, diverse geological features, and climatic conditions that contribute to the occurrence of multiple hazards, including earthquakes, landslides, debris flows, floods, and GLOFs. Geographically, the Hunza district is divided into three main regions: Upper Hunza (Gojal), Central Hunza, and Lower Hunza (Shinaki) (Valley and Dish, 2016). With a total area of 11,440 km 2 , elevation range of 1814 to 7749 m above sea level, and hosting some of the largest valley glaciers, including Batura (59 km), Passu (24 km), Gulkin (18 km), Hispar (61 km), and Ultar (10 km) (Fig. 1 , Panel 2). The Hunza River, which drains the region, is primarily fed by the melting of surrounding glaciers and serves as a major tributary of the Indus River. As a critical water resource, it plays a vital role in sustaining downstream water availability, making it essential for irrigation, hydropower, and overall water security in Pakistan 25 . However, the onset of summer accelerates the melting of glaciers and increases the water flow in streams and rivers, accelerating erosion, debris flows and landslides 26 , 27 . Fig. 1. Open in a new tab Administrative boundaries were obtained from open-source repositories Diva GIS ( https://diva-gis.org ) and manually adjusted for accuracy, (2) Location map of the district Hunza, generated in this study using ArcGIS Pro 2.8.3 software ( https://www.esri.com/en-us/home ). The study area is located in the region with the most active tectonics, a significant portion of Hunza Valley lies on the uplifted portion of the Main Karakorum Thrust (MKT) faults 28 , caused by the collision of the Indian and Eurasian tectonic plates 29 . Frequent earthquakes and active tectonic activity have extensively fractured the geological formations, making them prone to landslides. The geomorphology of the Hunza district consists of glacial-fluvial terraces, old glacier moraines, loose material on steep scree slopes, debris flow in the debris cone and colluvium deposits with split terrains, creating talus deposits at the bottom of towering cliffs 30 . The distinctive geological setting, characterized by thrust faults, fractured rock formations, and exposed geomorphology 30 , combined with human interference, including deforestation and infrastructure development 31 , have led to frequent and devastating disasters. To inform evidence-based policies for disaster risk reduction, a comprehensive multi-hazard analysis is essential for continuous monitoring and the implementation of proactive mitigation measures against natural hazards. Materials To acquire the topographic information, the ALOS-PALSAR digital elevation model (DEM) was utilized to compute the selected topographic attributes including terrain slope, aspect, plan curvature, profile curvature, Topographic Wetness Index (TWI), Relief Ratio (RR), Melton Ratio (MR), Stream Power Index (SPI), Sediment Transport Index (STI) and to extract streams and watershed boundary (Table 1 ). The distance from the river/stream and road was derived and classified into six buffer zones within ArcGIS Pro Environment 32 . Mean monthly rainfall data were acquired from the Global Precipitation Measurement (GPM) ( https://gpm.nasa.gov/data/directory ) at a spatial resolution of 0.1° (10 km) from 2014 to 2021. The rainfall data was resampled to 12.5 m in ArcGIS Pro-Platform to maintain spatial consistency with other predictor variables; however, this resampling does not increase the intrinsic precision of the rainfall data. Instead, it performs spatial averaging or interpolation to align datasets on a common grid, which is essential when integrating multiple variables. Rainfall data from 2014 to 2021, available from GPM provides adequate temporal coverage for analyzing variability and trends, given the lack of reliable meteorological stations in Hunza. A landcover map and Normalized Difference Snow Index (NDSI) were prepared based on the Sentinel-2 satellite images with 10 m spatial resolution, using the Support Vector Machine algorithm with a Kappa value of 0.85. The same image is used to compute the Normalized Difference Vegetation Index (NDVI) map. The lithology map was digitized from the geological map of northern Pakistan developed by Searle et al. 33 . GLOF hazard map was acquired from Ullah et al. 34 , with copyright clearance obtained through the publisher’s formal permissions process. Table 1. Details of geo-environmental factors used in multi-hazard susceptibility mapping, including their sources, scale, and the data from which they are derived. Data type Factors Data used Sources and scale Topographical 1. Elevation 2. Slope 3. Aspect 4. Plan Curvature 5. Profile Curvature 6. TWI 7. Streams 8. RR 9. MR 10. SPI 11. STI PALSAR DEM http://www.eorc.jaxa.jp/ALOS/en/aw3d30 (12.5 m*12.5 m resolution) Geological hydro-meteorological 12. Lithology 13. Faults 14. Rainfall 15. Roads Geological Map, GPM Data, Google Earth Pro Geological Survey of Pakistan https://gpm.nasa.gov/data/directory (0.1°) (2014–2021 Period) Environmental 16.Landcover 17. NDSI 18. NDVI Sentinel 2 A-2B (ESA) ESA Copernicus Open Access hub https://scihub.copernicus.eu/ GLOF hazard map 34 Open in a new tab Methods Multi-hazard modeling using machine learning In this study, we have integrated five hazards, including Landslide (LS), Debris Flow (DF), Flash Flood susceptibility (FFS), Seismic hazard (SH) and GLOF, using Forest-Based Classification and Regression Model (FBCR) (Fig. 3 ), an adaptation of Breiman’s random forest algorithm 35 . FBCR is a supervised machine-learning method that combines multiple decision trees, with each decision tree casting a unit vote for the most popular class 36 . Each tree in the forest is built using a bootstrap sample from the data, which may include duplicates and omit some original observations 37 . The FBCR is a natural, non-linear modeling tool that estimates variable importance, generalization error rates through the out-of-bag (OOB) error, and contributions to total risk 38 . The algorithm’s power lies in random feature selection and bagging 36 . As a result, FBCR can tackle both classification and regression problems, making it a popular choice for natural hazard forecasting 39 . The algorithm of FBCR is based on tree-structured classifiers. Each tree is trained on a bootstrap sample of the data, and the final prediction . Fig. 3. Open in a new tab Methodology for developing the multi-hazard susceptibility maps. is obtained by aggregating the outputs of all trees as follows: 1 where is the vector of hazard-related predictor variables, is the prediction of the tree, and n is the total number of trees in the FBCR model. Moreover, to integrate multiple hazards, each hazard map was classified (1 representing very low hazard, 2 for low hazard, 3 for moderate hazard, 4 for high hazard, and 5 for very high hazard (Fig. 16 ) and standardized to a common classification framework 40 . This ensured comparability between hazard maps derived from different modelling approaches (ML and AHP) and prevented any single hazard from dominating the integrated map due to differences in their native scales. Prior to integration, hazard weights derived from mean decrease Gini (MDG) analysis (details given in Sect. index importance degree) performed after model training, and further informed by field surveys, literature review and expert input, were applied to calibrate the combined map. All maps were classified and standardized to harmonize hazard-specific classification schemes, thereby ensuring consistency for direct comparison and subsequent analytical integration. The susceptibility modeling technique applied in this study requires samples of both hazards and non-hazards to generate multi-hazard susceptibility maps 41 . Therefore, an equal number of hazard and non-hazard locations was digitized in point format from each hazard map, and was assigned a value: (1) as hazard and (0) as non-hazard, which simplified the classification and integration process 15 , 24 . Following Pourghasemi, et al. 14 , Yousefi, et al. 15 , Pourghasemi, et al. 41 , the hazard locations were determined by identifying areas on each hazard map that showed the highest susceptibility levels, such as high and very high classes, while non-hazard locations were randomly sampled from moderate and low susceptibility classes. Random sampling also ensures the dataset’s spatial distribution throughout the study area, providing a balanced input for further processing. A total of 120 points were defined as LS points, while 120 as non-LS points, 90 points as SH, and 90 as non-SH points, 33 GLOF hazards points and 33 non-GLOF, 50 points for DF-50 non-DF, 60 FF points, while 60 represents non-FF points, and finally, the composite inventory of hazard locations was created (Fig. 4 ). Fig. 7. Open in a new tab Susceptibility maps showing the four hazards: LS map generated through a ML-FBCR model, SH map created by reclassifying slope and Vs 30 values, DF map developed using AHP, and FFS map generated through an ML- FBCR model. Fig. 4. Open in a new tab Composite multi-hazard inventory map developed for this study by the first author using ArcGIS Pro 2.8.0 software, showing sampled locations from each higher-hazard zones of individual hazard map. The number of hazard points varies due to the spatial distribution and intensities of each hazard map class. GLOF hazard is limited to a small region within particular glaciers, resulting in fewer points. LS is a frequent hazard across the district, leading to more hazard points. The valley and alluvial fans in the region are highly exposed and frequently impacted by debris flows, particularly during the rainy season or due to the accelerated triggering of glacial activity. Additionally, we divided the hazard and non-hazard point data using a random partitioning into two chunks of train and test by the ratio of 70 and 30%, where the train chunk was used for training, model tuning, and determining the variable’s importance, while the test chunk was used for validation 42 , 43 . Several rates, like 80/20 and 60/40 were tried to determine which provides the best model accuracy. However, 70/30 showed higher precision; thus, we selected 70/30 to construct the ML model. The point dataset was further utilized to extract the raster values of each hazard map as the input features (Predictor) for training the machine learning-FBCR model. These predictors enable the model to identify the spatial patterns between hazards and reflect the probability of hazard occurrence. Moreover, the predictor variables were exported from ArcGIS Pro 2.8.3 and processed in Python (scikit-learn), where the “RandomforestClassifier” was used to build the integrated multi-hazard map 44 . In addition, to reduce the bias caused by sampling in the data-sampling process, repeated sampling was performed over 100 runs, which is the standard value of the FBCR-model 45 , including the learning setting: number of trees = 500, leaf size = 5, and tree depth range = 47–50, and the number of predictors at each split was determined by the scikit-learn default (√number of predictors). However, the final map was imported and processed in ArcGIS Pro 2.8.3 software, including reclassification into hazard classes using the Natural Breaks (Jenks) method. Using out-of-bag estimates to monitor error We used the out-of-bag (OOB) method to estimate the generalization error of the model. By adjusting the number of decision trees, we observed a decrease in OOB error, which stabilized with minimal fluctuation beyond a sufficient tree count 46 . Lower OOB error indicates better model accuracy. Additionally, the OOB error was calculated as the ratio of misclassifications to the total samples 47 . Moreover, the OOB is constructed following the Eq. ( 2 ). 2 where is the average of all OOB predictions, n is the number of OOB observations in each tree, and MSE is the mean square error obtained during the construction of the classification trees 45 . Individual hazard assessment and influencing factors Flash flood and landslide inventory Hazard inventory is crucial in ML models and reveals the spatial interaction between the location of historical hazards and geo-environmental factors 48 . Therefore, in this study, an extensive flash flood and landslide inventory was developed (Fig. 5 ), comprising 56 flood-prone locations developed using Google Earth Pro, literature and field observations following Abbas, et al. 49 . Similarly, the landslide inventory consists of 344 locations, manually digitized through the visual interpretation of multi-temporal high-resolution Google Earth images (2010–2023), and verified during the field survey in September 2024 50 . Since we used the binary classification-FBCR model to develop the FFS and LS map, the model also requires the non-hazard locations (detail methodology is given in Sect. 2.3.1) 51 , which were sourced from areas that never experienced flooding or landslides, such as flat or gently sloping regions. Fig. 5. Open in a new tab Landslide-flood Inventory map for this study generated in ArcGIS Pro 2.8.3, the first author is depicted in (Panel 6), verifying landslide location with handheld GPS during the field survey, (Panels 6, 7, 9 and 10) magnified views of the landslide-flood areas, (Panel 8) Field photograph showing displaced boulders and damaged agricultural terraces opposite the recorded landslide area in Panel 7, (Panel 11) field photograph from Hassanabad along the Karakoram Highway, near the Excise and Taxation Department Hunza, showing a damaged structure as evidence of slope failure, (Panels 11 and 12) field photographs showing landslide-affected slopes along the Karakoram Highway, (Panel 13) destroyed houses and scattered rockfall or debris flow deposits near Ghulkin village, (Panel 10) Passu Gojal area showing flood impacts, (Panel 6) Karimabad areas frequently affected by landslides-floods, as reported through casual discussions of author with local residents during field visits, with field images as evidence. Influencing factors Based on the literature review 52 – 54 and the characteristics of the historical geohazards in the study area and multiple field observations, we selected eighteen geo-environmental variables for susceptibility mapping. Twelve variables, including altitude, slope, aspect, plan curvature, profile curvature, distance from river, distance from road, land use, NDVI, lithology, rainfall, and TWI were selected for FFS, out of which ten variables, such as slope, rainfall, elevation, aspect, profile curvature, distance from river, distance from road, land use, lithology, and fault, were chosen for the LS, and eight variables, slope, rainfall, TWI, Melton ratio, relief ratio, SPI, NDSI, and STI, were selected and mapped within ArcGIS Pro 2.8.3 platform for DF (Fig. 14 ). Fig. 5. Open in a new tab Multi-hazard susceptibility factors, computed using open-source datasets and generated in ArcGIS Pro 2.8.3 for this study. The terrain aspect is an important factor in determining flood occurrence because it is directly associated with the convergence and direction of water flow. Additionally, altitude impacts the natural movement of water, and thus influences flood susceptibility 55 . Terrain slope determines the flow and velocity of runoff and significantly influences damage distribution and intensity, and its angle impacts DF, FF, and LS, with flooding common in flat areas and steep slopes prone to destabilization 53 , 56 . The magnitude and intensity of rainfall are the primary triggers for all three hazard occurrences 57 , 58 . Topographic characteristics of an area were understood by plan and profile curvature 26 . Plan and profile curvature influence FFS by affecting water flow and runoff patterns, and LS by indicating slope steepness and change rate through profile curvature 59 . Landcover and NDVI affect the distribution and impacts of flash floods, with vegetation density impacting the velocity and runoff. Concurrently, landcover modifies hydrological and geomorphological dynamics, influencing landslide occurrence by altering soil stability and erosion processes 60 , 61 . The lithology determines the permeability and infiltration and hence influences the surface runoff and destructiveness of flash floods 62 . Similarly, lithology and distance to faults collectively influence LS by affecting geological stability and soil strength 32 . TWI, which illustrates soil wetness, water depth, and saturation, informs topographical controls and hydrological processes for FFS, while also indicating the water-saturated zones affecting infiltration and runoff, and susceptibility to LS and DF 63 . Distance from rivers affects the probability and magnitude of flooding and LS by influencing terrestrial water storage, flood events, and destabilizing slopes through river erosion and water accumulation 57 , 60 . Roads built on high-altitude areas alter water flow patterns, inducing high-intensity FFS that typically destroy road infrastructure and compromise slope stability, increasing LS through changed drainage patterns and erosion 61 . The MR, a ruggedness index, indicates topographic complexity, which can facilitate DF initiation by increasing runoff concentration and erosion potential 64 , 65 . The RR influences DF susceptibility by increasing the likelihood of slope failure due to its direct impact on gravitational stress 65 . STI exacerbates DF potential by optimizing sediment transport and deposition, whereas the NDSI influences DF by quantifying snowmelt-driven runoff variability, thereby enhancing the likelihood of debris flow events during summer months 54 , 66 . Debris flow hazard (DF) The DF hazard assessment methodology follows Shah, et al. 54 and Ullah, et al. 67 , whose study of in the same region topographically complex terrain aligns closely with the geomorphological and climatic conditions of the current study area. Given the absence of historical DF inventory data, a prerequisite for data-driven methods like ML (utilized for LS and FF with available inventories), the AHP was selected. AHP, validated in analogous data-scarce mountainous environments 68 , derives factor weights through expert-guided pairwise comparisons, with consistency ratios computed to ensure reliable judgment matrices (Fig. 5 , Panel 6). The selected DF influencing factors are consistent with Shah, et al. 54 and Ullah, et al. 67 , underwent rigorous multi-collinearity evaluation to ensure factor independence (Table 2 ). For a detailed description of the methodology, the reader is referred to Shah, et al. 54 Ullah, et al. 67 . Table 2. Multi-collinearity analysis for selected independent variables. Variables Collinearity test VIF Tol Altitude 1.87 0.48 Slope 2.60 0.41 Aspect 0.90 0.93 Plan curvature 1.57 0.64 Profile curvature 1.44 0.69 Distance from river 1.46 0.68 Distance from road 1.38 0.72 Rainfall 1.47 0.68 Land use 1.54 0.65 Lithology 2.03 0.49 TWI 1.67 0.60 NDVI 1. 25 0.97 NDSI 1.39 0.67 MR 1.82 0.53 RR 2.06 0.50 SPI 2.21 0.43 STI 1.92 0.53 Open in a new tab Seismic hazard (SH) Shear wave velocity of the top 30 meters surface (V s 30 ) is one of the important indicators of the seismic response amplification and is frequently utilized for seismic site characterization around the globe 69 – 71 . In regions with limited field measures, the topographic slope and geological information are effectively utilized as a proxy for the V s 30 -based seismic site characterization 72 , 73 . In this study, a V s 30 -based seismic site characterization map is developed by reclassifying the terrain slope map using the thresholds after Allen and Wald 74 . Index importance degree The FBCR model computes the relative importance of each predictor (referred to here as the ‘index importance degree’) to quantify its contribution to the modelled hazard susceptibility. In this study, variable importance was evaluated using the mean decrease in Gini impurity (MDG) of the Random Forest (RF) algorithm for the integrated multi-hazard mapping, FFS, and LS analyses. MDG measures how well a split divides samples in a node, considering feature importance and splitting frequency and is mathematically represented as: 3 where is the Gini importance, which specifies how often a specific feature is chosen to split and how important that feature is. Assigning the weight defines any imbalance in the distribution of classes during the learning Algorithm. while (, ) represents the change in Gini impurity resulting from splitting the data based on feature at a specific threshold τ and node M . Generally, natural hazards are threshold-dependent and are influenced by various effective factors 75 . Therefore, determining the relative importance of each hazard is key to evaluating the multi-hazard risk patterns 76 . Python (scikit-learn) was used to compute the Mean MDG from the RF model 36 , 77 , and validated the relative importance of each hazard in the multi-hazard susceptibility map through field surveys, literature review, expert opinions, and decision-maker engagement 78 . This validation approach ensured alignment between the calculated hazard weights and real-world conditions. Performance metrics ROC curves Validation is an essential component in the development of susceptibility mapping and the determination of its quality 79 . The FFS and multi-hazard susceptibility maps were validated by the receiver operating characteristic (ROC) curve. The Area Under Curve (AUC) value for this was estimated considering the true positive (TP) and false positive (FP) values of susceptibility modeling. TPs were pixels that were correctly estimated to be susceptible to hazards, while FPs were pixels that were incorrectly estimated to be susceptible to hazards 61 . Further, the ROC curve with AUC values has a range of 0.5–1.0. A higher AUC indicates a better model performance 80 , 81 . Results Individual hazard assessment Multi-collinearity analysis To assess the multi-collinearity between the causative factors selected for the multi-hazard susceptibility assessment, we conducted a multi-collinearity test using VIF and TOL limits (Table 2 ) and analysis affirmed that all eighteen variables were suitable for evaluating the multi-hazard susceptibility in the study area. Index importance degree analysis The importance of the geo-environmental variables for multi-hazard susceptibility mapping was estimated (Fig. 15 ). The most important variables for predicting FFS were rainfall (100%), altitude (81.89%), and distance from the river/streams (75.03%). For the LS the most important variables were slope (23.5%), lithology (18.7%), and elevation (18.3%). While the DF most important variables were Rainfall (25%), MR (21%), and RR (18%). Other variables were associated with moderate to lower importance. Fig. 6. Open in a new tab Relative importance of geo-environmental variables for individual hazard models, LS and FFS importance weights were derived from the mean decrease in Gini impurity of the RF, while DF variable importance was obtained from the weights assigned in the AHP model. Landslide susceptibility The areal coverage of LS zones covers an area of 20% (2335.84 km 2 ) Very High LS, 25% (2826.42 km 2 ) High LS, 14% (1614.55 km 2 ) Moderate LS, 16% (1831.93 km 2 ) low LS, and 24% (2776.02 km 2 ) Very Low LS (Fig. 16 ). The very high and high susceptibility zones are located in places with steep slopes (> 30°) in the northern, northwestern, northeastern, and central parts. Proximity to fault lines (< 200 m) and rivers (< 300 m) also increases susceptibility. Further, factors contributing to elevated landslide risk include low vegetation cover, deforestation, and anthropogenic activities such as road construction and agricultural terracing. These factors, combined with the district’s rugged topography, make these areas highly prone to slope failures. Flash flood susceptibility The flood and non-flood points containing the values of flash flood conditioning factors were used as input datasets in Python (scikit-learn) to build FFS models. The areal coverage of very high to high FFS areas was 5% (536.306 km 2 ) and 1% (131.43 km 2 ), while the moderate, low, and very low susceptible zones and the areal coverage of these zones was 18% (1983.02 km 2 ), 48% (5492.67 km 2 ), and 28% (3157 km 2 ), respectively (Fig. 16 ). The FFS map shows the regions associated with lower altitudes (1778–3553 m) and gentle slopes, specifically in the lower Hunza, central Hunza, and northwest and southwest regions along riverbanks are highly susceptible to frequent flash floods. This is due to intense rainfall events, glacier melt, and water accumulation. Land cover types like waterbodies, sparse vegetation, and built-up areas further increase flood risk in these regions. Conversely, regions with higher altitude (4311–7729 m) and steeper slopes reflect lesser susceptibility, as water moves away easily. However, these results show the combined impact of terrain features and hydrological factors on flood susceptibility in District Hunza. Seismic hazard The SH map classifies the district into: 11% (1288.63 km 2 ) very high SH, 23% (2689.86 km 2 ) high SH, 28% (3078.39 km 2 ) moderate SH, 20% (2263.38km 2 ) low SH, and 18% (2118 km 2 ) very low SH (Fig. 16 ). The reclassification of terrain slope and Vs30 reflects that areas with higher Vs30 values, typically found in steeper slopes and harder rocks, have lower seismic hazards. In contrast, areas with lower Vs 30 values, found in softer sediments and flatter terrain, are more prone to seismic hazards. These areas tend to amplify ground shaking during an earthquake. Notably, the central and northeastern regions, with softer sediments, show higher seismic hazards, making them more vulnerable. Regions near fault lines also face increased seismic risk due to lower Vs 30 values, highlighting the impact of geology and topography on SH. Debris flow hazard The DF hazard map categorizes catchments into five hazard levels using the natural breaks (Jenks) method (Fig. 16 ). The lower part of the district contains catchments with very high to high hazard levels, characterized by steep terrain, high rainfall, and loose sediment. This combination increases susceptibility to debris flow initiation and propagation. The results indicate ten very high-hazard catchments, covering 2% of the district (252.26 km²), and nine high-hazard catchments, covering 5% of the district (598.48 km²). Additionally, ten moderate hazard catchments, spanning 9% of the district (1035.25 km²), are identified. Driven by terrain conditions and intense rainfall, these zones pose significant debris flow risks. Multi-hazard analysis results Index importance degree analysis The MDG index, extensive field survey, available literature, local experts’ opinions, and engagement with decision-makers were considered to determine the importance of each hazard: SH was given (5.93) importance, LS (8.06), GLOF hazard (8.36), FFS (5.88), and DF (9.69). A synthesized multi-hazard susceptibility map All the selected hazards were integrated into a synthesized multi-hazard map. Pixels of the final integrated multi-hazard map were reclassified through ArcGIS Pro 2.8.3 natural break method into eight classes: Low Hazard, LS, FF, DF + FF, FF + LS+DF, SH + DF+LS + FF, DF + LS+GLOF + SH, and FF + SH+GLOF + LS+DF (Fig. 17 ), and the area were calculated for different classes of various hazards in percent and in km 2 (Table 3 ). Fig. 8. Open in a new tab Multi-hazard susceptibility map produced using FBCR model, synthesizing the five hazard maps for district Hunza (LS Landslide, FF flash flood, DF debris flow, SH seismic hazard and GLOF), and visualized in ArcGIS Pro 2.8.3. Table 3. Showing the areas of multi-hazard classes. Multi-hazard Area (km) 2 % Low hazard 6,469 56.84 LS 2630 23.11 FF 690.34 6.07 DF + FF 529.91 4.66 FF + LS+DF 453.11 3.98 SH + DF+LS + FF 327.66 2.88 DF + LS+GLOF + SH 190.07 1.67 FF + SH+GLOF + LS+ DF 91.29 0.80 Total 11380.95 100.00 Open in a new tab Areas of only LS 2630 km 2 (23.11%) or FF 690.34 km 2 (6.07%), areas with DF + FF is 529.91 km 2 (4.66%) and FF + LS+DF 453.11km 2 (3.98%), while the areas of SH + DF+LS + FF 327.66 km 2 (2.88%), DF + LS+GLOF + SH 190.07 km 2 (1.67%), FF + SH+GLOF + LS+DF 91.29 km 2 (0.80%), and Low hazard 6469 km 2 (56.84%). Significant variations can be seen in the integrated map in the northern part of Khunjerab Pass, especially in the NW and around Bara Khun Village and Pidakkesh Village (Fig. 17 , Panels 23 and 24), where there is an overlap between multiple hazards such as debris flows, landslides, and flash floods. In the NW region and close to Bara Khun Village of Khunjerab Pass, seasonal livestock overgrazing has resulted in decreased vegetation cover, increased erosion, and degradation of the land. In addition, the lower part of District Hunza shows substantial variations and overlap in all five hazards including the villages of Murtazabad (Fig. 17 , Panel 18), Hassanabad (Fig. 17 , Panel 19), Karimabad and Aliabad (Fig. 17 , Panel 20), and Hussaini (Fig. 17 , Panel 21), while, the map shows that Shimshal Valley (Fig. 17 , Panel 25) is mainly susceptible to landslides and flash floods. This is due to steep terrain, prevailing climatic conditions, intense rainfall, melting of glaciers, proximity to riverbanks, and human activities like infrastructure construction and population growth. The result indicates that no area is free from hazards even low-hazard areas are mainly found at higher altitudes and are primarily blanketed in snow/glaciers, showing that even these regions are susceptible to hazards. Using out-of-bag estimates to monitor error The results showed that the number of cycles increased from 250 to 500, and the number of trees in all MSE tests showed a decrease in every set of hazards from 29.44 to 28.09, respectively, as shown in Fig. 26 . This accounts for a variation explaining percentage ranging from 3.65 to 1.86 for DF, 32.55 to 29.50 for FF, GLOF 9.64 to 7.52, LS 27.51 to 25.93, and the SH is 33.44 to 32.84. These percentages indicate how well the model captures the variability in each hazard, with higher values of 33.44 to 32.84 indicating a better fit of the model to the data. Fig. 9. Open in a new tab Showing the FBCR out-of-bag error with the number of decision trees. Model performance validation The FBCR-model displayed the highest success-rate curve AUC value (0.866) for the synthesized multi-hazard map, the AUC value (0.822) for FFS map, while 0.843 value for LS and the AHP results with AUC value (0.791) for the DF. Discussion Northern Pakistan is susceptible to hydro-meteorological and geological hazards making multi-hazard studies essential. The present study takes a crucial step forward by developing the synthesized multi-hazard map for district Hunza Northern Pakistan using an ML model 36 . This novel approach, integrated multiple-hazards-FFS, DF, LS, GLOF, and SH (Fig. 16 ) to advance disaster risk reduction efforts in Hunza district, northern Pakistan. Natural hazard research has traditionally concentrated on individual threats in isolation, ignoring the critical connections between them, this oversight may lead to miscalculations of risk 82 . By examining how multiple hazards intersect and overlap in a specific area, we can gain a deeper understanding of their collective risk than is yielded by simply combining the results of single-hazard 19 . For this purpose, procedures have been developed worldwide based on varying situations and with varying amounts of available data, but they have advanced the modeling process and have revealed the spatial distributions of the natural hazards in many study areas. In this study, FBCR-supervised ML model utilized to develop an integrated multi-hazard map in a mountainous region of northern Pakistan. Besides determining variable importance of each hazard, we employed an integrated approach including the mean decrease (Gini Impurity) of FBCR model, extensive field surveys, available literature, local experts’ opinions, and engagement with decision-makers. This approach yielded the following importance values for each hazard: Seismic hazard (5.938), Landslide Susceptibility (8.064), GLOF hazard (8.364), Flash Flood (5.888), and Debris Flow (9.694). Given the region’s geological and climatic characteristics, no area is completely hazard-free, with most regions susceptible to at least one hazard. A substantial portion of the district faces the threat of multiple hazards, with spatial overlap between flash flood, landslides, seismic, GLOFs, and debris flow hazards. Particularly, our study finds that, while some areas are susceptible to a single hazard, numerous other areas might face two or even all five-hazards simultaneously, underscoring the complexity and severity of the hazard in Hunza (Fig. 17 , Panel 18). Furthermore, the results also highlight the evidence of spatial overlap of multiple hazards in low-lying locations, where the built-up area is primarily dispersed along lower elevations with scattered settlements 83 . The lower and central regions of the district appeared to be slightly more susceptible to multiple hazards due to the combination of increased constructed developments and climate events. Moreover, the Karakorum highway (KKH) and other essential road network is located in the Hunza district. The study confirms a multi-hazard susceptibility along the entire Karakorum highway. The FBCR model ROC-AUC value (0.866) indicates high performance and no overfitting, indicating it as an effective ML technique for multi-hazard susceptibility assessment 84 . These results are consistent with the study of other authors, Pourghasemi, et al. 14 created an integrated map considering-landslides, foods, and earthquakes using the ensemble model called SWARA-ANFIS-GWO. They found 17.14% of the region is influenced by no hazard, and most parts were prone to landslide and flood hazards together (33.70%). They also asses’ accuracies of 84% and 80% for flood and landslide maps. Ahmad, et al. 18 developed hazard inventories for upper Indus basin using remote sensing, field observations, and historical data to measure geohazards (such as landslides and debris flows) and then created a multi-hazard map using statistical methods, including SE, FR, WoE, and LR, the final susceptibility map shows the regions with overlapping hazards. Ullah, et al. 24 developed a multi-hazard susceptibility map by combining flash flood, landslide, and debris flow in district Shangla using the deep learning algorithm coupled with-historical hazards data, topography, hydrology, and environmental data where the final multi-hazard map demarcates that 62.43% of region is exposed to hazards while 37.57% of the study area are hazard-free. Rehman, et al. 85 have developed the database inventories for the geo-hazards (landslides and floods) and applied the AHP, frequency ratio and geospatial technique to produce a compound multi-hazard susceptibility map, and the map showed that majority of these geo-hazards were found in the high and very high susceptibility zones in Northwest Himalayas Pakistan. Further, integrated multi-hazard map will empower sustainable development and adaptive management in Hunza District by providing a comprehensive view of environmental hazards 86 . It will enable planners to implement integrated strategies, initiatives, and resource allocation, accounting all hazards simultaneously, for a resilient and sustainable future that it seems more susceptible to climate change, which is projected to intensify the frequency and harshness of extreme weather events (IPCC, 2014; Ministry of Climate Change, Pakistan, 2018). The synthesized multi-hazard map will also support prevention and mitigation programs that lower risk and improve resilience through policy and decision-making at the local, national, and international levels. Communities can identify intersections and cascade effects of natural disasters, develop comprehensive risk reduction strategies, and may help to raise awareness by identifying high-risk places. This strategy enhances emergency response planning, encourages involvement of stakeholders, and supports sustainable development. Additionally, we also believe that our findings will be useful for international readership because the developed maps and the approach will be used by relevant individuals or organizations involved in major construction projects in areas susceptible to multiple-hazards. However, this study also has several limitations such as lacks of glacier inventory data, reliant on past statistics, and lacks access to real-time data on hydro-meteorological hazards. The study’s spatial resolution might not be able to capture local-scale diversity, and the FBCR-model has a probability of overfitting 87 . Furthermore, multi-hazard susceptibility mapping was not as effective as it may have been since ML ignored the spatial structures of hazards while accounting for their spatial attributes. To address these problems, research on the integration of multi-hazard mapping using ML models combined with deep learning algorithms is being suggested 88 , which can make full use of both the spatial structure and attribute information of spatial objects. Secondly, integrated multi-hazard mapping has further limitations, with only one available GLOF map for the entire Hunza region, which restricts the accuracy of risk assessment. Therefore, this study recommends that further GLOF events should be generated to enhance the estimations of integrated hazard risks. In addition, the integrated multi-hazard map needs to be used to evaluate the exposure of elements that at risk, vulnerability of the region, and the total risk from the combination of integrated multi-hazard. Looking ahead, it is crucial to consider future exposure of elements at risk and develop strategies for a more resilient tomorrow. Conclusion In this study, integrated multi-hazard modeling was conducted for Hunza District, Northern Pakistan, by generating a composite hazard inventory to construct an integrated multi-hazard susceptibility map. Five hazards were considered: FFS, DF, SH, (GLOF), and LS. By integrating these hazards, we developed an effective Integrated Multi-Hazard map which exposes the reality that a substantial portion of the district faces the threat of multiple hazards. Particularly, the central and lower parts of the district, areas along transportation corridors, demonstrate relatively higher multi-hazard susceptibility. These findings highlight the importance of adopting integrated hazard assessment approaches in mountainous environments where hazard interactions are common. Although this study focuses on susceptibility rather than direct risk quantification, the generated multi-hazard maps provide a useful spatial framework to support land-use planning, infrastructure development, and disaster risk reduction strategies. Incorporating multi-hazard perspectives into planning processes may contribute to more informed decision-making and improved resilience in regions exposed to complex hazard dynamics. Acknowledgements The authors are grateful to the Higher Education Commission of Pakistan for supporting the study through the CRG-CPEC-130 project. The authors also extend their appreciation to the Deanship of Research and Graduate Studies at King Khalid University for funding this work through the Large Research Project under grant number RGP.2/487/46. Author contributions Muhammad Ismail khan: conceptualization, methodology, visualization, writing—original draft, field validation. Muhammad Shafique: conceptualization, supervision, methodology, field validation, writing and editing. Ghazanfar Ali Khattak.: field validation, writing and editing. Abid Yahya.: data curation, writing and editing. Irfan Anjum Badruddin: data curation, writing and editing. Data availability All datasets generated for this research are provided in the manuscript, and external datasets from literature were included with proper permissions under the relevant licenses. Declarations Competing interests The authors declare no competing interests. Ethical approval All ethical guidelines were followed; no formal human experiments were conducted, and informal discussions with local residents during field surveys required no formal consent or ethical approval. Consent for publication This manuscript is not under consideration for publication elsewhere. Footnotes Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Change history 4/30/2026 The original online version of this Article was revised: In the original version of this Article the Acknowledgements section was incorrect. The Acknowledgements section now reads: “The authors are grateful to the Higher Education Commission of Pakistan for supporting the study through the CRG-CPEC-130 project. The authors also extend their appreciation to the Deanship of Research and Graduate Studies at King Khalid University for funding this work through the Large Research Project under grant number RGP.2/487/46.” The original article has been corrected. References 1. Faramarzi et al. Using machine learning techniques in multi-hazards assessment of Golestan National Park, Iran. Nat. Hazards 1–25 (2023). 2. Dilley et al. Natural Disaster Hotspots: A Global Risk Analysis 5 (World Bank, 2005). 3. Aksha, S. K., Resler, L. M., Juran, L. & CarstensenJr L. W. A geospatial analysis of multi-hazard risk in Dharan, Nepal. Geomatics Nat. Hazards Risk . 11 , 88–111 (2020). [ Google Scholar ] 4. Eshrati et al. Multi hazards risk assessment, a new methodology. Int. J. Health Syst. Disaster Manage. 3 , 79 (2015). [ Google Scholar ] 5. Ward, P. J. et al. Natural hazard risk assessments at the global scale. Nat. Hazards Earth Syst. Sci. 20 , 1069–1096 (2020). [ Google Scholar ] 6. Ullah, F., Saqib, S. E., Ahmad, M. M. & Fadlallah, M. A. Flood risk perception and its determinants among rural households in two communities in Khyber Pakhtunkhwa, Pakistan. Nat. Hazards . 104 , 225–247 (2020). [ Google Scholar ] 7. Ali, S., Haider, R., Abbas, W., Basharat, M. & Reicherter, K. Empirical assessment of rockfall and debris flow risk along the Karakoram Highway, Pakistan: S. Ali et al. Nat. Hazards . 106 , 2437–2460 (2021). [ Google Scholar ] 8. Khan, I., Ullah, A., Zaidi, A. Z. & Panhwar, V. Assessing glacial lake outburst flood potential using geospatial techniques: a case study of western part of Gilgit-Baltistan, Pakistan. Arab. J. Geosci. 16 , 49 (2023). [ Google Scholar ] 9. Wasson, R. A debris flow at Reshūn, Pakistan Hindu Kush. Geogr. Annaler: Ser. Phys. Geogr. 60 , 151–159 (1978). [ Google Scholar ] 10. Uitto, J. I. & Shaw, R. In Sustainable Development and Disaster Risk Reduction 1–12 (Springer, 2016). 11. Hussain, M. et al. GIS-based multi-criteria approach for flood vulnerability assessment and mapping in district Shangla: Khyber Pakhtunkhwa, Pakistan. Sustainability 13 , 3126 (2021). [ Google Scholar ] 12. Downton, M. W., Miller, J. Z. B. & Pielke, R. A. Jr Reanalysis of US National Weather Service flood loss database. Nat. Hazards Rev. 6 , 13–22 (2005). [ Google Scholar ] 13. Bronkhorst, V. B. et al. Disaster risk management in South Asia: regional overview. The World Bank (2012). 14. Pourghasemi, H. R., Gayen, A., Panahi, M., Rezaie, F. & Blaschke, T. Multi-hazard probability assessment and mapping in Iran. Sci. Total Environ. 692 , 556–571. 10.1016/j.scitotenv.2019.07.203 (2019). https://doi.org:. [ DOI ] [ PubMed ] [ Google Scholar ] 15. Yousefi, S. et al. P. A machine learning framework for multi-hazards modeling and mapping in a mountainous area. Sci. Rep. 10 , 12144 (2020). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 16. Awan, S. A. The climate and flood risk potential of northern areas of Pakistan. Sci. Vis. 7 , 100–109 (2002). [ Google Scholar ] 17. Baig et al. District-level disaster risk and vulnerability in the Northern mountains of Pakistan. Geomatics Nat. Hazards Risk . 12 , 2002–2022 (2021). [ Google Scholar ] 18. Ahmad, H. et al. Geohazards susceptibility assessment along the upper indus basin using four machine learning and statistical models. ISPRS Int. J. Geo-Information . 10 , 315 (2021). [ Google Scholar ] 19. Kappes, M. S., Keiler, M., von Elverfeldt, K. & Glade, T. Challenges of analyzing multi-hazard risk: a review. Nat. Hazards . 64 , 1925–1958 (2012). [ Google Scholar ] 20. Conoscenti, C. et al. Gully erosion susceptibility assessment by means of GIS-based logistic regression: A case of Sicily (Italy). Geomorphology 204 , 399–411 (2014). [ Google Scholar ] 21. Yalcin, A. GIS-based landslide susceptibility mapping using analytical hierarchy process and bivariate statistics in Ardesen (Turkey): comparisons of results and confirmations. catena 72 , 1–12 (2008). [ Google Scholar ] 22. Mohseni, U., Jat, P. K. & Siriteja, V. Multi-criteria analysis-based mapping of the cyclone-induced pluvial flooding in coastal areas of India. DYSONA-Applied Sci. 6 , 309–321 (2025). [ Google Scholar ] 23. Youssef, A. M., Mahdi, A. M., Al-Katheri, M. M., Pouyan, S. & Pourghasemi, H. R. Multi-hazards modeling using machine learning algorithms in Southwestern Saudi Arabia. (2022). 24. Ullah, K., Wang, Y., Fang, Z., Wang, L. & Rahman, M. Multi-hazard susceptibility mapping based on Convolutional Neural Networks. Geosci. Front. 13 , 101425 (2022). [ Google Scholar ] 25. Tahir et al. Snow cover dynamics and hydrological regime of the Hunza River basin, Karakoram Range, Northern Pakistan. Hydrol. Earth Syst. Sci. 15 , 2275–2290 (2011). [ Google Scholar ] 26. Mukhtar, M. A. et al. Integrated flood risk assessment in Hunza-Nagar, Pakistan: unifying big climate data analytics and multi-criteria decision-making with GIS. Front. Environ. Sci. 12 , 1337081 (2024). [ Google Scholar ] 27. Owen, L. A. Mass movement deposits in the Karakoram Mountains: their sedimentary characteristics, recognition and role in Karakoram landform evolution. Z. für Geomorphologie . 35 , 401–424 (1991). [ Google Scholar ] 28. Ahmad, H. et al. Landslide risk assessment integrating susceptibility, hazard, and vulnerability analysis in Northern Pakistan. Discover Appl. Sci. 6 , 7 (2024). [ Google Scholar ] 29. Searle, M., Khan, M. A., Fraser, J., Gough, S. & Jan, M. Q. The tectonic evolution of the Kohistan-Karakoram collision belt along the Karakoram Highway transect, north Pakistan. Tectonics 18 , 929–949 (1999). [ Google Scholar ] 30. Hewitt, K. Catastrophic landslides and their effects on the Upper Indus streams, Karakoram Himalaya, northern Pakistan. Geomorphology 26 , 47–80 (1998). [ Google Scholar ] 31. Ali et al. The road to deforestation: An assessment of forest loss and its causes in Basho Valley, Northern Pakistan. Glob. Environ. Change . 15 , 370–380 (2005). [ Google Scholar ] 32. Pourghasemi, H. R. & Kerle, N. Random forests and evidential belief function-based landslide susceptibility assessment in Western Mazandaran Province, Iran. Environ. Earth Sci. 75 , 1–17 (2016). [ Google Scholar ] 33. Searle et al. Geological Map of North Pakistan and Adjacent Areas of Northern Ladakh and Western Tibet.(Western Himalaya, Salt Ranges, Kohistan, Karakoram, Hindu Kush), 1: 650 000 . (1996). 34. Ullah et al. An integrated approach for GLOF hazard, vulnerability and risk assessment in the Karakoram Mountain Range of northern Pakistan. J. Mountain Sci. 1–14 (2025). 35. Arabameri, A. et al. Novel ensembles of COPRAS multi-criteria decision-making with logistic regression, boosted regression tree, and random forest for spatial prediction of gully erosion susceptibility. Sci. Total Environ. 688 , 903–916 (2019). [ DOI ] [ PubMed ] [ Google Scholar ] 36. Breiman, L. Random forests. Mach. Learn. 45 , 5–32 (2001). [ Google Scholar ] 37. Yeh et al. Going-concern prediction using hybrid random forests and rough set approach. Inf. Sci. 254 , 98–110 (2014). [ Google Scholar ] 38. Krkač et al. Method for prediction of landslide movements based on random forests. Landslides 14 , 947–960 (2017). [ Google Scholar ] 39. Lee, S., Kim, J. C., Jung, H. S., Lee, M. J. & Lee, S. Spatial prediction of flood susceptibility using random-forest and boosted-tree models in Seoul metropolitan city, Korea. Geomatics Nat. Hazards Risk . 8 , 1185–1203 (2017). [ Google Scholar ] 40. Disha, R. A. & Waheed, S. Performance analysis of machine learning models for intrusion detection system using Gini Impurity-based Weighted Random Forest (GIWRF) feature selection technique. Cybersecurity 5 , 1 (2022). [ Google Scholar ] 41. Pourghasemi, H. R. et al. Assessing and mapping multi-hazard risk susceptibility using a machine learning technique. Sci. Rep. 10 , 3203 (2020). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 42. Cao, J., Zhang, Z., Wang, C., Liu, J. & Zhang, L. Susceptibility assessment of landslides triggered by earthquakes in the Western Sichuan Plateau. Catena 175 , 63–76 (2019). [ Google Scholar ] 43. Remondo, J. et al. J. F. Validation of landslide susceptibility maps; examples and applications from a case study in Northern Spain. Nat. Hazards . 30 , 437–449 (2003). [ Google Scholar ] 44. Pedregosa, F. et al. Scikit-learn: Machine learning in Python. J. Mach. Learn. Res. 12 , 2825–2830 (2011). [ Google Scholar ] 45. Pumhirunroj, B. et al. Machine-Learning-Based Forest Classification and Regression (FCR) for Spatial Prediction of Liver Fluke Opisthorchis viverrini (OV) Infection in Small Sub-Watersheds. ISPRS Int. J. Geo-Information . 12 , 503 (2023). [ Google Scholar ] 46. Fang, L., Wang, Q., Yue, J. & Xing, Y. Analysis of Optimal Buffer Distance for Linear Hazard Factors in Landslide Susceptibility Prediction. Sustainability 15 , 10180 (2023). [ Google Scholar ] 47. Fang wt al. Analysis of Optimal Buffer Distance for Linear Hazard Factors in Landslide Susceptibility Prediction. Sustainability 15 , 10180 (2023). [ Google Scholar ] 48. Nachappa, T. G. et al. Flood susceptibility mapping with machine learning, multi-criteria decision analysis and ensemble using Dempster Shafer Theory. J. Hydrol. 590 , 125275 (2020). [ Google Scholar ] 49. Abbas, N. et al. Environmental disaster assessment using geospatial techniques for Hunza-Nagar district, Gilgit-Baltistan, Pakistan. Sci. Int. 28 , 5319–5328 (2016). [ Google Scholar ] 50. Pham, B. T., Pradhan, B., Bui, D. T., Prakash, I. & Dholakia, M. A comparative study of different machine learning methods for landslide susceptibility assessment: A case study of Uttarakhand area (India). Environ. Model. Softw. 84 , 240–250 (2016). [ Google Scholar ] 51. Duy et al. Flood risk assessment using machine learning, hydrodynamic modelling, and the analytic hierarchy process. J. Hydroinformatics . 26 , 1852–1882 (2024). [ Google Scholar ] 52. Wang et al. Urban flood prediction under heavy precipitation. J. Hydrol. 577 , 123984 (2019). [ Google Scholar ] 53. Chapi, K. et al. A novel hybrid artificial intelligence approach for flood susceptibility assessment. Environ. Model. Softw. 95 , 229–245 (2017). [ Google Scholar ] 54. Shah, N. A., Shafique, M., Owen, L. A., AI-Mulla, Y. & Ullah, Y. Morphometric analysis of debris flow hazard and risk assessment in the mountain terrains of northern Pakistan using remote sensing and field data. Earth Sci. Inf. 18 , 295 (2025). [ Google Scholar ] 55. Das, S. Geospatial mapping of flood susceptibility and hydro-geomorphic response to the floods in Ulhas basin, India. Remote Sens. Applications: Soc. Environ. 14 , 60–74 (2019). [ Google Scholar ] 56. Wang et al. Spatial prediction of landslide susceptibility based on GIS and discriminant functions. ISPRS Int. J. Geo-Information . 9 , 144 (2020). [ Google Scholar ] 57. Tehrany, M. S., Kumar, L., Jebur, M. N. & Shabani, F. Evaluating the application of the statistical index method in flood susceptibility mapping and its comparison with frequency ratio and logistic regression methods. Geomatics Nat. Hazards Risk (2018). 58. Ullah, K. & Zhang, J. GIS-based flood hazard mapping using relative frequency ratio method: A case study of Panjkora River Basin, eastern Hindu Kush, Pakistan. Plos one . 15 , e0229153 (2020). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 59. Basharat, M., Shah, H. R. & Hameed, N. Landslide susceptibility mapping using GIS and weighted overlay method: a case study from NW Himalayas, Pakistan. Arab. J. Geosci. 9 , 1–19 (2016). [ Google Scholar ] 60. Azarafza, M., Azarafza, M., Akgün, H., Atkinson, P. M. & Derakhshani, R. Deep learning-based landslide susceptibility mapping. Sci. Rep. 11 , 24112 (2021). [ DOI ] [ PMC free article ] [ PubMed ] [ Google Scholar ] 61. Band, S. S. et al. Flash flood susceptibility modeling using new approaches of hybrid and ensemble tree-based machine learning algorithms. Remote Sens. 12 , 3568 (2020). [ Google Scholar ] 62. Bui, D. T., Tsangaratos, P., Ngo, P. T. T., Pham, T. D. & Pham, B. T. Flash flood susceptibility modeling using an optimized fuzzy rule based feature selection technique and tree based ensemble methods. Sci. Total Environ. 668 , 1038–1054 (2019). [ DOI ] [ PubMed ] [ Google Scholar ] 63. Mandal, S., Mondal, S., Mandal, S. & Mondal, S. Machine learning models and spatial distribution of landslide susceptibility. Geoinformatics and Modelling of Landslide Susceptibility and Risk: An RS & GIS-based Model Building Approach in the Eastern Himalaya 165–175 (2019). 64. Melton, M. A. Geometric properties of mature drainage systems and their representation in an E4 phase space. J. Geol. 66 , 35–54 (1958). [ Google Scholar ] 65. Nikolova et al. Morphometric analysis of debris flows basins in the Eastern Rhodopes (Bulgaria) using geospatial technologies. Nat. Hazards . 105 , 159–175 (2021). [ Google Scholar ] 66. Ahmad, I. et al. Application of hydrological indices for erosion hazard mapping using Spatial Analyst tool. Environ. Monit. Assess. 191 , 482 (2019). [ DOI ] [ PubMed ] [ Google Scholar ] 67. Ullah, Y., Shafique, M., Shah, N. A., Khattak, G. A. & Ullah, S. Multi-scale debris flow risk assessment using morphometrics in high-altitude catchments of Northern Pakistan. Nat. Hazards . 122 , 7. 10.1007/s11069-025-07833-9 (2025). [ Google Scholar ] 68. Yin, G., Luo, J., Niu, F., Lin, Z. & Liu, M. Machine learning-based thermokarst landslide susceptibility modeling across the permafrost region on the Qinghai-Tibet Plateau. Landslides 18 , 2639–2649 (2021). [ Google Scholar ] 69. Ahmed, T., Rehman, K., Shafique, M., Ali, W. & Shah, N. A. Proxy-based Vs30 modelling of the Muzaffarabad, northern Pakistan. Environ. Earth Sci. 83 , 188. 10.1007/s12665-024-11502-8 (2024). [ Google Scholar ] 70. Shafique, M., Hussain, M. L., Asif Khan, M., van der Meijde, M. & Khan, S. Geology as a proxy for Vs30-based seismic site characterization, a case study of northern Pakistan. Arab. J. Geosci. 11 , 298. 10.1007/s12517-018-3648-6 (2018). [ Google Scholar ] 71. Castellaro, S. & Mulargia, F. VS30 estimates using constrained H/V measurements. Bull. Seismol. Soc. Am. 99 , 761–773 (2009a). [ Google Scholar ] 72. Wills, C. J., Gutierrez, C. I., Perez, F. G. & Branum, D. M. A Next Generation VS30 Map for California Based on Geology and Topography. Bull. Seismol. Soc. Am. 105 , 3083–3091. 10.1785/0120150105 (2015). [ Google Scholar ] 73. Iwahashi, J., Kamiya, I. & Matsuoka, M. Regression analysis of Vs30 using topographic attributes from a 50-m DEM. Geomorphology 117 , 202–205. 10.1016/j.geomorph.2009.11.004 (2010). [ Google Scholar ] 74. Allen, T. I. & Wald, D. J. Topographic slope as a proxy for seismic site-conditions (VS30) and amplification around the globe (Report No. 2331 – 1258). (Geological Survey (US), 2007). 75. Conforti, M., Aucelli, P. P., Robustelli, G. & Scarciglia, F. Geomorphology and GIS analysis for mapping gully erosion susceptibility in the Turbolo stream catchment (Northern Calabria, Italy). Nat. Hazards . 56 , 881–898 (2011). [ Google Scholar ] 76. Conoscenti, C., Di Maggio, C. & Rotigliano, E. GIS analysis to assess landslide susceptibility in a fluvial basin of NW Sicily (Italy). Geomorphology 94 , 325–339 (2008). [ Google Scholar ] 77. Ai, F. et al. Application of random forests to select premium quality vegetable oils by their fatty acid composition. Food Chem. 143 , 472–478 (2014). [ DOI ] [ PubMed ] [ Google Scholar ] 78. Gallina, V. et al. A. A multi-risk methodology for the assessment of climate change impacts in coastal zones. Sustainability 12 , 3697 (2020). [ Google Scholar ] 79. Pourghasemi et al. Landslide susceptibility mapping using a spatial multi criteria evaluation model at Haraz Watershed, Iran. Terrigenous Mass Movements: Detection, Modelling, Early Warning and Mitigation using Geoinformation Technology 23–49 (2012). 80. Jaafari, A. et al. Meta optimization of an adaptive neuro-fuzzy inference system with grey wolf optimizer and biogeography-based optimization algorithms for spatial prediction of landslide susceptibility. Catena 175 , 430–445 (2019). [ Google Scholar ] 81. Jaafari, A. LiDAR-supported prediction of slope failures using an integrated ensemble weights-of-evidence and analytical hierarchy process. Environ. Earth Sci. 77 , 42 (2018). [ Google Scholar ] 82. Gill, J. C. & Malamud, B. D. Reviewing and visualizing the interactions of natural hazards. Rev. Geophys. 52 , 680–722 (2014). [ Google Scholar ] 83. Akbar et al. Multi-hazard susceptibility mapping for disaster risk reduction in Kargil-Ladakh Region of Trans-Himalayan India. Environ. Earth Sci. 82 , 68 (2023). [ Google Scholar ] 84. Golkarian, A. & Rahmati, O. Use of a maximum entropy model to identify the key factors that influence groundwater availability on the Gonabad Plain, Iran. Environ. Earth Sci. 77 , 1–20 (2018). [ Google Scholar ] 85. Rehman, A. et al. Multi-hazard susceptibility assessment using the analytical hierarchy process and frequency ratio techniques in the Northwest Himalayas, Pakistan. Remote Sens. 14 , 554 (2022). [ Google Scholar ] 86. Tiepolo et al. Multihazard risk assessment for planning with climate in the Dosso region. Niger Climate . 6 , 67 (2018). [ Google Scholar ] 87. Cheng, C. H. Evaluating naval tactical missile systems by fuzzy AHP based on the grade value of membership function. Eur. J. Oper. Res. 96 , 343–350 (1997). [ Google Scholar ] 88. Yang, L. & Cervone, G. Analysis of remote sensing imagery for disaster assessment using deep learning: a case study of flooding event. Soft. Comput. 23 , 13393–13408 (2019). [ Google Scholar ] Associated Data This section collects any data citations, data availability statements, or supplementary materials included in this article. Data Availability Statement All datasets generated for this research are provided in the manuscript, and external datasets from literature were included with proper permissions under the relevant licenses. Articles from Scientific Reports are provided here courtesy of Nature Publishing Group ACTIONS View on publisher site PDF (13.9 MB) Cite Collections Permalink PERMALINK Copy RESOURCES Similar articles Cited by other articles Links to NCBI Databases Cite Copy Download .nbib .nbib Format: AMA APA MLA NLM Add to Collections Create a new collection Add to an existing collection Name your collection * Choose a collection Unable to load your collection due to an error Please try again Add Cancel Follow NCBI NCBI on X (formerly known as Twitter) NCBI on Facebook NCBI on LinkedIn NCBI on GitHub NCBI RSS feed Connect with NLM NLM on X (formerly known as Twitter) NLM on Facebook NLM on YouTube National Library of Medicine 8600 Rockville Pike Bethesda, MD 20894 Web Policies FOIA HHS Vulnerability Disclosure Help Accessibility Careers NLM NIH HHS USA.gov Back to Top