Multi-Hazard Susceptibility Assessment Using the Analytical Hierarchy Process and Frequency Ratio Techniques in the Northwest Himalayas, Pakistan

: Globally, natural hazards have become more destructive in recent times because of rapid urban development and exposure. Consequently, signiﬁcant human life loss, the damage to property and infrastructure, and the collapse of the environment directed the attention of geoscientists to control the consequences and risk management in relation to geo-hazards. In this research, an effort was made to produce a compound map, geo-visualizing the susceptibility of multi-hazards, to select suitable sites for sustainable future development and other economic activities in the region. Muzaffarabad District was chosen as a case research area due to the high magnitude of hydro-meteorological and geological hazards. On the one hand, both selected geo-hazard inventories were developed using the ﬁeld survey and remote sensing data. The subjective and objective weight of all the causative factors and their classes were calculated using the assembled geospatial techniques, such as the Analytical Hierarchy Process (AHP) and Frequency Ratio (FR) in the Geographic Information System (GIS). The results reveal that the most suitable areas are distributed in the southern and northwestern parts, which can be used for future sustainable development and other economic activities. In contrast, the eastern and western regions, including Muzaffarabad City, are within high and very susceptibility zones. Finally, more than 50% of the land area is located in very low and low susceptibility zones. The validation of the proposed model was checked by using three different techniques: the Receiver Operative Characteristic (ROC) curve, Seed Cell Area Index (SCAI), and Frequency Ratio (FR). Both ROCs, the Success Rate Curve (SRC) and the Predictive Rate Curve (PRC), showed the goodness of ﬁt for both the selected geo-hazards: landslides (81.3%) and ﬂoods (93.2%), at 80.1% and 91.7%, respectively. All the validation techniques showed good ﬁtness for both the individual and multi-hazard maps. The proposed model sets a baseline for policy implementation for all the stakeholders to minimize the risk and sustainable future development in areas of high frequent geo-hazards.


Introduction
The world has been facing various natural and human-made disasters resulting in life losses, and the damage to property, environment, economy, infrastructure, and other aspects of human life throughout history. In recent decades, the frequency of natural

Study Area
Muzaffarabad District was chosen as our case research area in which to test the multi-hazard susceptibility model. Geographically, the study area is situated between 34 • 04 00" N to 34 • 35 00" N latitudes and 73 • 30 45" E to 73 • 59 00" E longitudes, having an area of approximately 1300 square kilometers ( Figure 1). Politically, the research area is the capital of the Azad State of Jammu and Kashmir, part of Pakistani Administrated Kashmir (PAK). The Muzaffarabad District contains two tehsils: the first is Muzaffarabad itself, and the second is the Pattika, also known as the Naseerabad. Muzaffarabad City is located on the confluence of the River Neelum and the River Jhelum. The study area's relief is low in the south (Jhelum Valley), increasing towards the northeast and southeast (Neelum Valley) as high as 4375 m due to the collision zone of the Indian and Eurasian tectonic plates. The research area lies in the humid-to-high subtropical land receiving maximum precipitation (1400 mm) during the summer monsoons (late June to early September). During the summer monsoon, torrential and intense rainfall spells often flood the major rivers and their tributaries. The mean minimum temperature for January and the mean maximum temperature for June are −2.65 • C and 36.75 • C, respectively.
Tectonically, the research area is highly seismically active and part of the northwestern Himalayas. Geologically, the rock formations present in the research area range from Precambrian to Quaternary. Lithologically, the area mainly consists of slate, phyllite, metasediments, metavolcanics, marble, limestone, dolomite, sandstone, shales, unconsolidated silt, sand, and clay deposits.
The research area, including Muzaffarabad City and part of Neelum valley and Jhelum, has a population of 726,000 with an annual growth rate of 2.80% [31]. Being a capital city, all the government, non-government, and nonprofit organization (NGO) offices are located here. All roads that are linked to the major cities are paved and connected with bridges. From an economic point of view, the study area has great importance due to its natural resources and beauty that attracts many tourists locally and internationally.

Methodology
Four steps accompanied this research work: (1) geo-hazard inventories and all the causative factors were prepared and their spatial database generated in ArcGIS 10.5 software; (2) the subjective weights (SWs) and objective weights (OWs) of both geo-hazards were calculated using the assembled FR-AHP geospatial techniques, (3) the weights of the individually produced geo-hazard maps were synthesized for the multi-hazard susceptibility assessment; and (4) the accuracy assessment and model validation was assessed to check the performance of the produced multi-hazard susceptibility assessment by using the ROC, SCAI, and frequency ratio methods. A complete methodology flowchart is given in Figure 2.

Data Collection and Selection of the Causative Factors
Data was collected from primary and secondary data sources to accomplish the research objectives. Intensive fieldwork was carried out for the primary data collection. Geological and topographical maps were obtained from the Geological Survey of Pakistan [32]. In addition, satellite images were downloaded from the United States Geological Survey (USGS), Earthexplorer, and the Alaska Satellite Facility (ASF) websites. The complete detail of the data and their sources are shown in Table 1.

Field Survey and Preparation of the Inventories
Field visits were carried out to better understand the actual ground conditions and real-time data collection. During the field visits, the spatial distribution of the geo-hazards and their causative factors were well examined. These visits were beneficial for preparing the database inventories for the landslides and floods. The landslide inventory map was prepared from the SPOT-5, QuickBird, Google Earth, and Sentinel 2 satellite images (from 2005 to 2018) using visual interpretation and the supervised image classification tool in the ArcGIS 10.5 software. A total of 670 landslides were identified ( Figure 1). The interpretation of all the landslides was identified at a scale of 1:50,000. During the field visits, 68 landslides (app. 5% of the total known/interpreted landslides) were confirmed using the Global Positioning System (GPS) to check the accuracy of the landslide inventory map. A total of 55 landsides out of 68 (80.7% accuracy assessment) was confirmed and accurately located. All the identified landslides were included with equal importance during the analysis, despite their type, size, and shape. The 2010 and 2014 flood events inundation data and flood plain locations were confirmed using high-resolution satellite images and field visits for the flood inventory map. A total of 180 locations out of 208 (86.5% accuracy assessment) were confirmed and accurately located in Figure 1. Later, both fields' varied landslide and flood locations were used for the model validation and accuracy assessment rather than dividing the data into training and validating data sets, as was suggested by many researchers [26,33,34]. The study area had a very rugged topography and steep slopes. Some remote areas were inaccessible due to a lack of road links. During the field visits, a digital photograph was taken and verified, and the accessible locations were chosen to confirm the accuracy of the inventory maps.

Geospatial Techniques
Multicriteria decision analysis (MCDA) is a complex phenomenon used to evaluate several causative factors simultaneously. Sometimes, MCDA, also known as Multicriteria Decision Making (MCDM), is a systematic method used to address complex problems by dividing the complex problem into small parts, analyzing them, and integrating them in a meaningful and logical manner. The resultant solution was then obtained in a single composite [35]. MCDA assembled with GIS (GIS-MCDA) effectively provides a wide range of geographically spatial data manipulation costs and time. A GIS-MCDA offers procedures and techniques to evaluate, analyze, design, and prioritize value judgments. Many studies show the importance of GIS-MCDA for hazard assessment [19,[36][37][38][39][40]. Many qualitative and quantitative GIS-MCDA techniques and procedures were developed to analyze and manipulate spatial data in recent years. Most widely used geospatial techniques are bivariate and multivariate, such as Frequency Ratio (FR) [12,41,42], AHP [20,22,43,44], logistic regression [9,33,42,45], fuzzy logic [46], the weight of evidence [34,47], and their assembled techniques, such as FR-AHP, FR-WoE, FR, LR-WoE, and LR-AHP [33,47,48]. In this research, the FR-AHP assembled geospatial technique was used to calculate the weights of the selected causative factors and their subclasses.

Analytical Hierarchy Process (AHP)
The AHP was introduced by Saaty [49]. This method is also known as a pairwise comparison method because, in this method, a comparison matrix is used to calculate the weight of the selected factors by calculating the matrix's eigenvectors [50]. The AHP method is a decision-making process through which a set makes complex decisions regarding the choices or alternatives. Sometimes, these choices/alternatives are qualitative, quantitative, and occasionally conflicting [51][52][53]. In this method, two selected factors are simultaneously evaluated, and preference is given to one factor over another based on experienced judgments. These judgments are measured using a preference scale of 1-9, in which a lower number (1) indicates a lesser preference and a higher number (9) shows more importance of that particular factor over another. The preference score values are given in Table 2, which was used in this research. The AHP method also provides the accuracy assessment and consistency of the weights and rank scores of the priorities (set of alternatives), based on experienced judgments [37,42,55]. The consistency (Consistency Ratio (CR)) can be calculated from the following Equation (1): where CR is the consistency ratio of a particular AHP matrix, the Consistency Index (CI) can be calculated from Equation (2), and RI is the random index, which depends on the number of factors (n).
where λ max is the consistency vector (average value of each calculated eigenvectors), and "n" represents the number of selected factors. The RI was estimated from 500 different matrices by Saaty [56], and their values are given in Table 3. According to Lootsma [57] and Arabameri, and Rezaei [58], the CR should be within the limits of 10% (CR < 0.01), otherwise the assigned score values should be reconsidered. The Frequency Ratio is the most simple and widely bivariate statistical technique, often used in many research types, such as the earth sciences, environmental sciences, hazard assessments, and physical sciences. FR enables the evaluation of the impact of a phenomenon (such as floods and landslides) by the attributes (subclasses) of a causative factor (such as slope, lithology, and curvature). Moreover, the FR is defined by Lee and Pradhan [59] and Lee and Sambath [60], which is the ratio of the probability of a phenomenon occurrence to that of no phenomenon occurrence for any given attribute. According to Yilmaz [61], the FR is the most understandable and straightforward statistical technique. The FR of each causative factor class can be expressed as the ratio of the percent of the class and the percentage of the total hazard (e.g., landslides and floods) in that class [48]. All the causative factors for both selected geo-hazards were classified to calculate the OW. Subsequently, the inventory maps were overlaid for each causative factor, and the FR was calculated within each causative factor class (Figures 3 and 4). A ratio value greater than 1 indicates a substantial correlation with a geo-hazard occurrence and/within a causative factor class, and a ratio value of <1 indicates the weaker relationship [48,62]. For the final rank (weights of a class) calculation, the FR was further integrated with AHP (FR-AHP). The AHP matrix assigned the score values according to their FR value, instead of an expert opinion value judgment.

Assembled Geospatial Techniques
To avoid the biases and uncertainties of assigning the weights of all causative factors and their classes, the assembled geospatial technique was used in this research, namely as the Frequency Ratio assembled to the Analytical Hierarchy Process (FR-AHP). The aim of selecting the assembled geospatial techniques was to avoid assigning class weights based on so-called experienced judgments or expert opinions. These techniques have proven to be the best fit, together [48,58,63]. The SWs of the geo-hazards, as mentioned above, were evaluated by an intensive literature review, site visits, and experienced judgments. Additionally, a survey was conducted by 10 relevant experts/geoscientists. Their expert opinions and value judgments were then averaged, and an appropriate score value was assigned to each causative factor. The SWs were then determined with the AHP application using the Expert Choice 11 software. First, the FR was calculated for each causative factor class for the determination of the OWs. After that, the final OWs were calculated using the AHP application according to their relevant FR values.

Multi-Hazard Susceptibility Assessment
On the one hand, geo-hazards, as mentioned above, were produced individually by the application of Weighted Linear Combination (WLC) in ArcGIS 10.5 as follows: where LSI is the Landslide Susceptibility Index, n is the total number of causative factors, CF i represents the SWs of the causative factors, i and SSCW j i are the OWs that are subjectively determined by the jth class within causative factor i using the AHP [48].
where FSI (Equation (4)) is the flood susceptibility index, n is the total number of causative factors, CF i represents the SWs, and i and SSCW j i are the OWs that are subjectively calculated. The individual geo-hazard susceptibility maps that we produced in the previous sections were classified into an equal number of classes by using the same classifying scheme, SD. However, the simple summation of individuals producing susceptibility maps is not an appropriate way to create the single multi-hazard map [27]. Since every geohazard has different intensity and interaction/consequences with other hazards, sometimes a single hazard triggers multiple hazards [23]. To overcome this problem, [26,27] used the qualitative approach for the relative weight calculation of individual geo-hazards. However, in this paper, we used the same qualitative approach to assess the relevant weight of both geo-hazards based on the historical record and field knowledge.
Furthermore, a questionnaire survey was also conducted by 10 experts. The assigned weights were then averaged and used for the final multi-hazard map. The first preference, with a 70% relative weight, was assigned to LSA because the research area was situated in the northwestern Himalayas, which is highly seismically active. The frequency of landslides can be seen from a single event (the 2005 Kashmir earthquake) that triggered thousands of landslides in the research area [4,10,64,65]. The 30% relative weight was assigned to the FSA. The overall relative weights for the final multi-hazard map were then synthesized by using the following Equation (5): where MSI is the multi-hazard susceptibility assessment, n represents the number of geohazards, H i is the geo-hazard, and Wi is the relative weight of the geo-hazard i.

Accuracy Assessment
The accuracy assessment and validation of the models are essential for susceptibility/predicting analysis [22,66]. In the current research work, we used three different techniques, the Receiver Operating Characteristics (ROC) curve, the Seed Cell Area Index (SCAI) method, and the frequency Ratio model, to measure the models' accuracy and future prediction power.
ROC is a graphical plot that reflects the effectiveness and predictive power of the model in a qualitative hierarchy by calculating the Area Under the Curve (AUC). A higher percentage of the area below the curve indicates a higher accuracy of the model, and the lower area below the curve indicates a lower accuracy of the model to predict the future occurrences of the phenomenon. On the other hand, ROC also quantitatively measures the model accuracy as AUC <0.6 (<60%), 0.6-0.7 (60-70%), 0.7-0.8 (70-80%), and >0.9 (>90%) indicate the modeling accuracy as weak, moderate, good, and very good, respectively [67]. The Success Rate Curve (SRC) and Prediction Rate Curve (PRC) were used for both individual geo-hazards (LSA and FSA) and the multi-hazard map. For calculating the AUC for SRC and PRC, the dependent variables (such as landslides and floods) must be in their proportion of training and validating data sets. The different schemes are available in the literature to divide the data sets (dependent variable) into their training and validated data sets, as 50:50, 60:40, and 70:30, respectively. In this research, to validate the selected model's actual field data sets were used that were collected during the field visit and the accuracy assessments of the inventories. Many researchers chose the latter classification scheme for its popularity and recommendation [34,47].
The SCAI method represents the area coverage of each hazard susceptibility class in percentage over the seed cells (hazard density in percentage) in each susceptibility class [68,69] as: SCAI = % area coverage/% seed cell (6) It is accepted that the SCAI value should be very high for the very low and low susceptibility classes, and the SCAI value should be very low for the high and very high susceptibility classes, accordingly [19,33].
The frequency ratio model is considered as quite popular in the scientific community to evaluate the model accuracy [9,70]. This method is quite similar to the SCAI method but in a reverse manner. The susceptibility maps were first classified by the same classification schemes with an equal number of susceptibility classes using the later models. Subsequently, the inventories maps were superimposed onto each susceptibility class, and the SCAI and Frequency Ratio were calculated.

Preparation of the Geo-Hazard Maps
As mentioned above, the study area was regularly affected by the geo-hazards due to the unfavorable geological conditions, seismicity, rugged topography, and many anthropogenic activities. Based on the extensive literature review, site characteristics and data availability, as well as eleven causative factors, namely the slope, elevation, hydrology, aspect, curvature, road, lithology, structure, NDVI, LULC, and rainfall, were selected for the LSA, out of which eight causative factors, such as the slope, elevation, hydrology, lithology, structure, NDVI, LULC, and rainfall, were chosen for the FSA. The complete detail of each causative factor and individually produced geo-hazard map is presented in the following sections.

Slope
A slope is an important causative factor when considering landslides and floods [30,55]. Usually, steep slope areas are more susceptible to landslides [71] and a low slope gradient, usually flat areas, are more susceptible to flooding [55]. The slope gradient was generated from the DEM image. The slope gradient ranged from 0 to 86 degrees, and was then classified into five and seven numbers of classes for the LSA and FSA. After the classification, the results showed that most of the research area (29.86%) was found in the 28-35 • slope class ( Figure 3A), and Muzaffarabad City and the other residential areas along the River Jhelum were found in the lowest slope class range ( Figure 3B). While comparing these results with the inventory maps, the highest landslide frequency (2.08) was found in the higher slope range (46-86 • ) and the lowest slope class was found to have the highest flood frequency at 41.88 (Table 4).

Elevation
The higher altitude areas are important for slope instability analysis, and the low lying areas, such as flat areas, foothills, and valley floors, are more susceptible to flooding when interlinked with discharge channels [14,66,72]. The elevation map was generated by using the DEM image. The altitude difference range, from 521 to 4435 m, was accordingly classified into seven and nine numbers of classes for the LSA and FSA, respectively, as shown in Figure 4A,B. The relationship between the location distribution (inventory) and elevation classes revealed that the lowest elevation class had the highest landslide and flood frequency ratios, as 2.28 and 101.2. The calculated FRs and the OWs of the elevation classes are given in Table 5.

Hydrology
The hydrology of an area directly controls the groundwater [73], flooding [24,40,72], and landslides [12]. The watershed of the research area was delineated from the DEM image by using the Arc hydro tool 10.5 extensions according to the Strahler classification scheme [74]. After delineating the hydrology, the Euclidean distance tool generated the buffer around the streams. The distance buffer was classified into five classes: <50, 50-100, 100-150, 150-200, and >200 m ( Figure 5). The calculated buffer distance was considered around all the delineated streams with equal importance because the landslides were found in every stream buffer class. While calculating the FR, all the stream buffer classes were found to have a substantial relationship (FR > 1) for the landslide distribution, except for the stream buffer class range with a distance of >200 m. Additionally, most of the flooded locations were found near the streams and are presented in Table 6.

Aspect
The slope facing (aspect) is also suggested by [45,75,76], which incorporates sunlight exposure, evapotranspiration, precipitation towards soil erosion, and the vegetation type [77], and was calculated from the DEM image. Most of the research area (16.36%) was found on a southwestern facing slope, and the lowest area (0.05%) was found on a flat surface ( Figure 6). The southern facing slopes (SE, S, SW) were found, FR > 1, indicating a positive correlation with the landslides. In contrast, the northern facing slopes (NE, N, NW) showed a weaker relationship with the landslides, occurring as FR < 1 ( Table 7). As expected, the flat areas showed no occurrence of landslides (FR = 0).  The curvature is the slope curve derivative that indicates the slope's steepness at a certain point by the line tangent, indicating the slope change rate at that particular point. The slope curvature is an essential topographic factor in landslide susceptibility [22,78]. The slope curvature map was generated from the DEM image. The plan curvature type suggested by the above researchers was used. The calculated plane curvature type was classified into three classes: the negative value indicates a concave slope, the values from −0.1 to 1.0 indicate the flat plane, and >0.1 indicates a convex slope, respectively, as shown in Figure 7 [66,76]. After the classifications, the results showed that most of the area lay on both concave and convex classes with an approximate equal-area percentage (47.9 and 48.0), and only 4.1% was found on the flat surface. The relationship between the curvature and landslide distribution indicated that both curvature classes had quite similar landslide distributions, but the concave class shows a slightly higher FR = 1.09. For this reason, the concave class attained a higher OW than the convex class (Table 8). Road construction, especially in the hilly areas, led to mass movement and slope instability due to the undercutting of the steep slope and the removal of material near the foothill [34,59]. The areas under and sides of the road also experienced the heavy load and vibration that caused slope failure. The road networks of the research area were digitized from eight mosaicked topographical maps of 1:50,000 and a high resolution (0.6 × 0.6 m) QuickBird image. The buffer was classified with distances of 0-50, 50-100, 100-150, and >150 m (Figure 8). The highest FR 2.27 was found within a distance of 0-50 m, followed by the 50-100 and 100-150 m as 1.63 and 1.37, respectively ( Table 9). The higher FR near the roadsides indicated that the road construction and road cuts positively correlated with the landslides and slope instability in the research area.

Lithology
An area's geology (lithological units) is an essential and classic parameter for the geo-hazard susceptibility assessment [5,47,79,80]. The eight geological maps covering the research area at the scale of 1:50,000, were acquired from the Geological Survey of Pakistan (GSP) and published maps [14,[81][82][83][84]. The geological map of the research area was constructed and ten different lithological rock formations were identified (Figure 9).
The landslide distribution results found that the Muzaffarabad Formation has the highest FR, 5.1, while the Recent Alluvium class showed that the highest flood distribution was 6.85 (Table 10).

Structure
The structural features, such as the lineaments, contact boundaries between different lithological units, and faults lines, produced weaker fracture zones, leading to slope instability [14,66], and control the ground and surface water [85]. The structural map of the research area was digitized from geological [32] and published maps [34]. Four major fault lines were identified as the Panjal Thrust (PT), Main Boundary Thrust (MBT), Muzaffarabad Fault (MF), and Jhelum Fault (JF). A buffer was generated with distances of 0-250, 250-500, 500-750, 750-1000, and >1000 m, presented in Figure 10. The developed buffer also revealed that most of the settlements were near the fault lines, especially Muzaffarabad City, sandwiched between MBT and MF. The landslide distribution results (Table 11) found that the FR possesses a negative correlation with increases in the distance from fault lines. For testing the structure control on flooding, this causative factor was tested for the first time in the research area. However, the smaller lineament features and the contact boundaries between the different lithological units were not included, and only the major faults lines were considered. Table 11 illustrates the OWs and the correlations between the flood locations and the fault buffer classes, with a substantial correlation of FR > 1.    Many anthropogenic factors, such as tunnels, roads, bridges construction, and deforestation, often lead to slope failure [19,22], and these artificial surfaces lead to low or no infiltration and more runoff. Sensitive areas (near the stream channel) are becoming more susceptible to flooding [86,87]. A medium-to-high resolution (10 × 10 m) Sentinel-2 satellite image was used to generate the LULC map of the research area. The "Maximum Likelihood Classification" command was used for supervised image classification [3]. A total of six LULC classes were identified as barren land, built-up land, forest, grassland, snow cover, and water bodies (Figure 11.) The highest OWs were assigned to the barren land (FR = 2.30) and built-up land (FR = 1.27), due to the substantial correlation with landslide occurrences ( Table 12). The LULC with the flood locations is given in Figure 11. Three LULC classes, barren land, built-up land, and water bodies, were found to have a strong positive correlation with the flooding, with a value of FR > 1. Other land use classes were found to have a weaker correlation with FR < 1 (Table 12).  The NDVI represents the vegetation/biomass density of an area. The vegetation cover shows a negative correlation with the slope instability [76] and flooding [88]. The NDVI was calculated from Sentinel-2 images by using the following equation. The derived NDVI ranges from −0.32344 to 0.848607, and was then classified into three vegetation types, densely vegetated areas (forest), sparsely vegetated areas (grassland and rangeland), and no vegetated areas (built-up, barren, and water bodies), which are shown in Figure 12. Table 13 represents the FRs and OWs for the distribution of the geo-hazards, and it was found that the no vegetation class had the highest landslide and flood distributions.

Rainfall
Rainfall is a critical factor responsible for triggering landslides and slope instability [12,20]. The research area lies in the monsoon region and receives most of its rainfall during the summer, which often causes flooding in major and minor streams. Despite its great importance, rainfall is not included in most previous research of the study area. For the construction of the rainfall map, the rainfall data from 20 metrological stations for the last 30 years were used. The rainfall data were averaged, and a database was constructed in ArcGIS. For generating the prediction surface map, a geo-statistical technique Inverse Distance Weighting (IDW) method was used in ArcGIS. The rainfall map was then classified into three classes: 1330-1350, 1350-1370, and 1370-1390 mm ( Figure 13). While calculating the FR and OW, it was observed that the medium rainfall class (1350-1370 mm) positively correlated with landslide and flood occurrence (Table 14).   After determining all the selected causative factor classes, the relevant OWs and SWs were calculated using the FR-AHP method. The complete procedure was discussed earlier in the Methodology Section, and the SWs for the LSA and FSA are given in Tables 15 and 16, accordingly. The assigned weights influenced the final resultant maps; therefore, great care was taken, and the Consistency Ratio (CR) was firstly prompted, before the generation of the final maps, and both Tables 15 and 16 show good CRs, as 9% and 7%. Subsequently, all the OWs and SWs were combined, and the Landslide Susceptibility Index (LSI) and Flood Susceptibility Index (FSI) maps were calculated. For the susceptibility classes, the index maps were classified using the standard deviation method into the equal number of classes, very low, low, medium, high, and very high, as suggested by [24,27]. The final LSA and FSA maps are shown in Figure 14A,B.

Multi-Hazard Susceptibility Assessment (MSA)
The individual geo-hazard maps were synthesized using the relative score values of 70% (LSA) and 30 % (FSA) via the weighted sum tool in ArcGIS. The complete procedure is given in the Methodology Section. For the generation of the multi-hazard susceptibility classes, the synthesized index map was classified into five classes, described in the previous section. The multi-hazard map represents both geo-hazards. To verify the impact, both historical events in the form of inventories were overlaid with a multi-hazard map, and the highest densities of both geo-hazards were found in the very high and high susceptibility zones ( Figure 14C).

Landslide Susceptibility Assessment (LSA) Map
The LSA map was derived by calculating the assembled FR-AHP and GIS methods. Accordingly, the derived map was classified into five classes: very low, low, medium, high, and very high susceptibility zones ( Figure 14A). The landslide inventory map was then superimposed over the LSA map to check the landslide impact on the research area. The calculated result within each landslide susceptibility zone is given in Table 17. More than 50% of the land area was found within the very low and low susceptibility zones, and the lowest land area (app 16%) was found within the high and very high susceptibility zones. In contrast, the lowest landslide density was found in the first two zones and the latter zones showed the highest landslide occurrence. The results also indicate that the northwestern and southeastern areas, such as Muzaffarabad City, including the vicinity areas of Dhani, Nausari, and Langarpura, are the most landslide-affected areas due to the active seismicity, unfavorable lithological conditions, and many anthropogenic activities, and these results can be seen from the SWs in Table 15, and correlated with the previous research [5,11,64,84,89].

Flood Susceptibility Assessment (FSA) Map
According to Figure 14B, the high and very high flood susceptibility zones are near the major streams. Muzaffarabad City was found within a very high susceptibility zone, in which major streams converge. When comparing this with the historical events (flood inventory), it is also evident that 83% of the flood locations were situated within very high and high flood susceptibility zones. The areal extent and spatial distribution of flood events within each flood susceptibility zone are given in Table 17. The SW also indicates that the major streams and tributaries are responsible for triggering and enhancing the flood susceptibility, followed by the elevation, slope, and lithology (Table 16). The distance from fault lines, which we considered the first time in the research area, showed less importance for the SW calculation of 0.083 (8.3%). Still, objectively, the area near the fault zones shows a substantial correlation with the flood events (Table 11 and Figure 14B).

Multi-Hazard Susceptibility Assessment (MSA) Map
The result of the MSA map is shown in Figure 14C, which directed the distribution of both geo-hazards in the research area. The majority of these geo-hazards were found in the high and very high susceptibility zones, indicating that the research area was severely affected in the past. The percentage of the area coverage of the multi-hazard susceptibility zones and geo-hazard densities are given in Table 17. The calculated results from the relative weights that we assigned during the production of the MSA map also proved the best fit for the MSA model, for predicting both geo-hazards in the research area. Regarding the high and very high susceptibility zones, most of the settlements already fall into these zones, indicating the great potential for harm during the next geo-hazard event. In addition to this, the model also predicted the very low and low susceptibility zones in the southern, western, and further northeastern parts suitable for future planning. However, these parts lie in remote areas.

Accuracy Assessment and Validation of the Selected Model
The accuracy assessment and validation of the models are important for susceptibility analysis [22,66]. We used three different techniques, Receiver Operating Characteristics (ROC) curve, Seed Cell Area Index (SCAI), and Frequency Ratio model, to measure the models' accuracy and future prediction power.
To observe the smooth results of the generated multi-hazard maps and their ability to predict geo-hazards, both the Success Rate Curve (SRC) and Prediction Rate Curve (PRC) were used.
The AUC was calculated by plotting the cumulative susceptibility area percentage on the x-axis and the cumulative landslide and flood occurrence percentages on the y-axis on a scatterplot. The success rates (SRCs) for the LSM and FSM were calculated as 0.743 (74.3%) and 0.958 (95.8%), respectively. For the MSA, the SRC was calculated as 0.743 (74.3%) and 0.906 (90.6%) for both geo-hazards, as shown in Figure 15A and Table 17. Only the SRC cannot represent the efficiency of predicting the power of the models. To measure the real predicting power of the models, the same procedure was repeated for the PRC; however, instead of using training data-set, the validating data sets of known landslides and flood locations were used, which were collected during the field visits for the construction inventory maps. After calculating the AUC, the findings showed results that were relatively similar to those that we calculated for the SRC. This indicated the best fit of the models to predict the occurrence of future geo-hazards (Table 17). The calculated PRC and their details are given in Figure 15B. The second approach for testing the accuracy of the used models and the efficiency of the predicting power was carried out using the SCAI method, suggested by [68,69].
This method represents the area coverage of each susceptibility class in percentage, over the seed cells (hazard density in percentage) in each susceptibility class. Using the SCAI, it is accepted that the SCAI value should be very high for very low and low susceptibility classes, and should be very low for high and very high susceptibility classes, respectively [19,33]. The calculated SCAI shows all the models' goodness of fit because all the classified maps show the very high SCAI values for the very low and low susceptibility classes and vice versa (Table 17).
Reclassified susceptibility maps were also validated using the frequency ratio suggested by [9,70]. This method is quite similar to the SCAI method, but in a reverse manner. The susceptibility maps were first classified by the same classification schemes with an equal number of susceptibility classes using the Frequency Ratio model. Following this, the inventory maps were superimposed onto each susceptibility class, and the frequency ratio was calculated.
Theoretically, it is accepted that the frequency of the phenomenon will gradually increase from very low to very high susceptibility classes. The calculated results are given in Table 17, indicating that the Frequency Ratios of LSA, FSA, and MSA gradually increase from very low to high susceptibility classes ( Figure 15C).

Conclusions
Landslides are more common and disastrous in a high terrain that has a fragile slope. The rugged topography, active seismicity, unfavorable geological conditions, heavy monsoon rainfall, and rapid urban growth increase the susceptibility of Muzaffarabad to landslides. After the interpretation of the results, it can be concluded that the following geo-environmental conditions play a vital role in triggering landslides: the slope ranges from 46 to 86 degrees; elevation variations between 521 to 1048 m; stream buffers between 50-100 m; rainfall ranges between 1350-1370 mm in a southern facing slope; positive plane curvature values; unconsolidated and highly fractured rock units; fault distances of 0-250 m; and variations in the LULC. Multicriteria analysis integrated with geospatial techniques was used for the LSA. Additionally, an effort was made to develop a suitable model for the LSA and geo-visualizations by utilizing the historical records. The validations of the model show quite satisfactory goodness of fit and future predicting power as 74.3% and 73.8%, respectively. The selected model produced the satisfactory goodness of fit (74.3% and 90.6%) and predictive power (73.1% and 90.1%) to identify both selected geo-hazards. It can also be concluded that the LSA model is cost-effective and can be used in other geographical locations with the same geo-environmental conditions.
Flooding is a major threat to Pakistan. During the summer monsoon, major and minor streams flooded and often inundated nearby agricultural and populated areas. The Muzaffarabad District was selected to test the susceptibility model due to its ever-growing population and many economic activities near the stream channels without proper planning and site investigations. After the interpretation of all the results, it can be concluded that flood distribution is primarily governed by geo-environmental conditions, such as the slope degrees <5 • ; elevation <621 m; distance from both minor and major streams (stream buffer) 0-50 m; unconsolidated silt, clay, and gravel deposits; distance from major faults (fault buffer) 500-700 m; no vegetation of NDVI; water bodies and settlements of the LULC; and rainfall class between 1350-1370. In this research, multicriteria geospatial techniques were used for FSA, based on the flood inventory map after the 2010 and 2014 flood events. The implemented susceptibility model showed a goodness of fit of 90.6% and a predicting power of 90.1%, by using the ROC area under-curve. We can also concluded that the FSA model is a cost-and time-effective model that can be used in other geographical locations with the same geo-environmental conditions.
The multi-hazard susceptibility model further concluded that there is a positive correlation between the floods and landslides that affect the research area. The multi-hazard analysis used in the current study yielded reliable results in the delineating regions appropriate for urban development, based on the location of previous landslides and flood events.
The current research provides the initial baseline for multi-hazard susceptibility assessment in the northwest Himalayas. The current study's implications will help to understand the most susceptible areas to minimize the risks for all stakeholders, even though the current study consisted of two geo-hazards due to the tight time frame and resource allocation.
Nevertheless, the proposed model shows the goodness of fit; therefore, it is strongly recommended that more geo-hazards be added for future detailed susceptibility assessments.