District-level disaster risk and vulnerability in the Northern mountains of Pakistan

Abstract Mountain systems of Northern Pakistan are highly vulnerable due to their spatial locality, low natural capacity, insufficient or no alternatives to cope and to reduce the adverse impacts posed by natural stresses. This case study aims to provide a holistic approach to estimate district-level disaster risk and vulnerability rankings in the northern mountains of Pakistan (∼72,000 km2) by integrating inferences from the hazard and vulnerability assessments. Based on eight hazard and nine coping capacity indicators, districts are scaled by using Analytic Hierarchy Process and overlay weight analysis. The result shows that moderately exposed mountain communities with lower coping capacity is likely to cause calamities due to mountain-specific constraints. Specifically, the level of risk as the likelihood of destruction and severity is high in Shigar, Diamer and Ghanche. Moreover, Gilgit, Shigar and Ghizer lie in the Karakoram portion of Northern Pakistan are highly exposed to natural hazards. The ability of mountain communities residing in central Karakoram (Baltistan region) is lowest to withstand and cope with adverse effects of hazardous phenomena. This approach is expected to find possibilities for learning the relationship between hazards, likelihood of disasters and coping capacities of mountain communities.


Introduction
Pakistan is amongst the most disaster-prone countries in the South Asia, that has undergone an approximate 18 billion US$ losses and damages during the last decade (Khan 2017). Especially, the Climate Vulnerability Index of 2019 has ranked Pakistan 8th amongst the 10 most affected countries, because of extreme weather conditions between 1998 and 2017 (Eckstein et al. 2019). Mountain systems of Northern Pakistan are highly vulnerable due to their spatial or socio-political locality, low inherent capacity to cope or are not resilient or there are insufficient or no alternatives to reduce the adverse impacts posed by natural stresses. Especially, due to its the geographic location, occupying portions of three mountain ranges (Hindu Kush, Karakoram and Himalaya), high seismic processes, steep slopes and a combination of a variety of glaciers, glacier lakes and high-pitched snow-covered peaks (Baig et al. 2020), this region is under the threat of an extensive diversity of typological disasters (e.g. floods, landslides, mudflows and avalanches). Continued glacier surging (Hewitt 2014;Bhambri et al. 2017) and declining has progressively generated many unstable glacier lakes that are controlled by moraines (Gilany et al. 2020). Furthermore, as a result of atmospheric hazards, seasonal low depressions emerging over the Bay of Bengal or the Arabian Sea and regular flooding in the Indus river basin trigger secondary hazards like building collapse, glacial lake outburst floods (GLOFs) (Baig et al. 2020) and mudslides. Potential risk to downstream mountain communities is increasing due to the absence of hazard and disaster preparation plans to counter adverse impacts of geological hazards (Shah 2013).
This study focuses on ten districts of Gilgit-Baltistan, the northern mountains of Pakistan within the Upper Indus Basin (UIB) and the junction of three mountain ranges (Hindu-Kush, Himalaya and Karakoram). Many studies focused on perception-based (Khan, Qureshi et al. 2020;Khan, Rana, et al. 2020) or comprehensive examination of individual geo-hazards in northern Pakistan (Rehman et al. 2014;Ali et al. 2015;Khan et al. 2016;Calligaris et al. 2017;Rahim et al. 2018). But ten districts of this mountainous region occupying $72,000 km 2 area exhibit varieties of exposures to different risk scenarios. Impacts of these scenarios on each district's specific mountain systems and the ability of their systems to cope with the adverse impacts are different as well. In Diamer district, tectonic movements at frequent intervals (Rafiq and Blaschke 2012), about 4 ± 1 mm yr À1 slip rate (world's greatest strike-slip) along the Karakoram fault (Brown et al. 2002) between the colliding boundaries of the Eurasian and Indian tectonic plates regularly generate geological hazards. The unique geomorphological characteristics (Derbyshire et al. 2001) of high pitched mountains peaks (e.g. K2, Gashaburum I, II, III) of districts (e.g. Shigar), lies in the Karakoram, cause extensive variety of geological hazards (e.g. avalanches, GLOFs and landslides). Hydro-metrological hazards particularly glaciers-induced flash floods are permanent phenomena in Hunza, Shigar and Skardu districts. Especially, Astore district of Gilgit-Baltistan is extremely predisposed to flash and river floods during spring, summer and winters (Easterling et al. 2000). Eight GLOFs pose a great threat for the downstream population of Hunza district due to glacier advancing (Hewitt 2014;Bhambri et al. 2017) and retreat (Baig et al. 2018). The residents of central Hunza residing at a substantial distance from the volatile lakes (e.g. Khurdupin in Shimshal valley) and Attabad-lake (Cook and Butz, 2015) are exposed to danger. Therefore, for a better disaster preparedness strategy, it is imperative to scale the contributing environmental and social factors of Gilgit-Baltistan districts which have collectively increased the vulnerability of the mountain communities. Therefore, this article aims to provide a holistic approach to estimate risk (adverse interaction of humans and the environment) rankings of these mountainous districts of UIB by integrating inferences from the hazard (UNISDR, 2009) and vulnerability (referring to socio-economic conditions especially 'poverty' in social sciences and referring to climate change in environmental sciences) (Alwang et al. 2001) assessments in the form of an approximation of the potential threat to downstream mountain communities and the environment.
Natural disasters occur from the natural processes or origins that affect humans adversely. In order to establish the relationship between potential affected downstream communities as vulnerability and potential natural hazards, data pertaining to two sets of indicators are employed: (I) the sources or origins of hazards and (II) mountain communities' response mechanisms for sustainable livelihood indicators. The input data for indicator I are topography, land cover, environment and seismicity, which may trigger a particular hazard that causes destruction and/or damage. For example, debris-flow (geology) due to glacier movement from higher to lower elevation (topography) (Muhammad et al. 2020) can damage agricultural land (land cover). Similarly, 'soil' plays a vital role in erosional processes, therefore, is employed as a contributing factor. Indicator II reflects the ability of a community or a region to respond to or prepare for a hazard (Birkmann 2007).
Proximity, population density, scientific understanding, early warning systems and communication lines, building types, cultural factors, public education and awareness, are major factors which plays vital role in assessments of vulnerability to hazards and disasters. To assess the indicator II, data pertaining to the environment, Human Development Index (HDI), socio-economic conditions of mountain communities are analyzed. For district-wise comparison, an index value is maintained for both hazards and vulnerability. Risk factor of each district is calculated by multiplying potential severity scaling obtained from hazard and vulnerability assessments. An equal ratio Figure 1. Map of the study area (ten districts), six districts lie in the Karakoram, and three in the Himalayan region along with a portion of Kharmang district and a significant portion of Diamer district in the Hindukush region of Northern Pakistan. Elevation is overlaid on district boundaries. from vulnerability and hazard assessment is used to calculate total risk of each district. Geographic Information System (GIS)-based Analytic Hierarchy Process (AHP) (Saaty 1990) and overlay weight methods are employed to obtain the results. Moreover, probability of population at risk is calculated as the total risk divided by the population at risk.

Topography and spatial extent of study area
The mountainous region of Northern Pakistan comprises the western part of the Hindu Kush-Karakoram-Himalayas (HKH) region which stretches from Afghanistan to Myanmar. The study area is spread over $725,000 km 2 and falls within the Pakistan-administered territory of Gilgit-Baltistan (GB), bordering with Afghanistan, China and Indian-administered Kashmir ( Figure 1). The Karakoram begins in the Wakhan Corridor (Afghanistan) in the west and encompasses the majority of Gilgit-Baltistan (Pakistan) and extending into Ladakh (India) and the Aksai Chin region of China. Six districts namely Ghizer (11,700 km 2 ), Gilgit ($ 4000 km 2 ), Hunza (11,200 km 2 ), Nagar ($3000 km 2 ), Shigar ($9000 km 2 ) and Ghanche (8200 km 2 ) and a small portion of Kharmang ($2800 km 2 ) district lie within the Karakoram range. Astore ($5200 km 2 ) and Skardu ($7000 km 2 ) districts entirely lie in the Himalayan range. Almost half of the Diamer districts ($7000 km 2 ) lie in the Karakorum and half in the Hindu-Kush range. Ghizer district is the largest district of Gilgit-Baltistan followed by Hunza and Shigar.
Much of Gilgit-Baltistan is a cold desert with the northern and north-western valleys receiving only between 100-200 mm of precipitation annually (Archer and Fowler 2004) but higher elevations intercept about 1000-2000 mm water equivalent, resulting in snowfall and glacier formation (Hewitt 2014). Majority of these districts host many giant peaks >7000 metres above sea level (m) and roughly 700 peaks above 6000 m (Wester et al. 2019;Ahmad et al. 2020). Snow and glacial melt water in summer is the source of many streams, springs, lakes, wet meadows, marshes and peat lands in which a variety of fauna and flora flourish.
Altitudinal zones along with major rivers of Northern Pakistan, hosts many giant peaks >7000 meters above sea level (m) and roughly 700 peaks above 6000 m (Wester et al. 2019;Ahmad et al. 2020). These peaks have made this region prone to glacier related hazards as well.

Hazard and vulnerability analysis flowchart
Methodology flowchart ( Figure 2) comprises names of hazard and vulnerability indicators. Details pertaining to datasets used for this analysis are presented in Table 1. Relative importance of one indicator over another (in %) is presented in the flowchart by assigning weightages to each indicator as influencing factor or contribution to risk estimation.
The input data for hazard mapping is seismicity, fault lines, earthquake incidents, geology (Searle 2011), soil types and land cover, elevation and slope. Environmental indicator comprised of land cover classes (Figure 3(c)). Although, the employed land  (25%) is assigned to geological formation as it is the main cause of triggering earthquakes, mud and landslides in Northern Pakistan followed by seismic zone (12%), fault lines (12%) and earthquake incidences (11%). 'Vulnerability' in this study refers the way a hazard will affect human lives and properties. For vulnerability mapping, this study employs environmental and economic indicators of the districts (Table 2). Although, the employed forest and agriculture taxonomy as contributing indicators plays vital role in mitigation processes but highest weightage (20%) is assigned to degree of child dependency followed by poverty (15%). The levels of hazard depend on the natural settings, sensitivity or behaviour of origins or sources that has the potential to cause destruction or harm the environment or population. A relationship between hazard likelihood and its sensitivity presented in Table 3 is established.
Vulnerability of districts to hazards and disasters as a key component in the risk calculation is linked with the positive capacities to cope, survive and recover from the adverse impacts of hazards. The capacities could be hazard-specific managerial or operational, to reduce the range of hazard impacts and the extent of vulnerability. A relationship between coping capacity (1 ¼ high to 10 ¼ low) and severity of exposure of population presented in Table 4 is established.

Hazard and vulnerability mapping
3.2.1. Hazard mapping Methodology to reclassify hazard related datasets based on a common measurement scale by using hazard specific methods to generate a hazard map is presented in Figure 4.
For hazard mapping, topographic factors (elevation and slope) (Figure 3(a)) are extracted from a world-wide compiled, arranged and filtered ASTER Global Digital Elevation Model (GDEM2) (Fujisada et al. 2012) obtained from https://earthexplorer. usgs.gov/. DEM-based analysis by using ArcGIS 10.1 (ESRI) enabled us to generate range of slopes from topographic factor (elevation) (Figure 3(a)). A preliminary seismo-tectonic map (scale: 1:2,000,000) comprised of seismic zones, earthquake incidence locations and fault lines, is digitized and extracted the data specific to Northern Pakistan (Figure 3 Four levels-of-hazards (e.g. minor, insignificant, moderate and major) are assigned to the reclassified classes obtained from weighted overlay method.

Vulnerability to hazard mapping
Proximity of settlements to hazards and population density are key elements of vulnerability to hazard mapping. Therefore, geographic locations of settlements (villages) are used as a baseline for vulnerability to hazard map generation. The communities living in these villages have different HDI, surrounding environment and socio-economic conditions. Therefore, these points are populated with the tabular data obtained from many sources (Table 1) related to vulnerability indicators (Table 2) and produced vulnerability to hazard maps.

Risk mapping
In order to produce risk map, out of the total hazard map of Gilgit-Baltistan (including likelihood of potential damage to populated and non-populated areas), hazardous region within 5 km buffer zone around villages is extracted and reclassified based on common measurement scale. The resulting hazard map is multiplied with reclassified vulnerability map.

Standardization of hazard and vulnerable maps
Reclassification of raster dataset is the process of re-assigning a value or a range of values, or a list of values in a raster to new output values. The 'reclassify' tool of ArcGIS replaces values of raster-based datasets to alternative values as a common measurement scale by using a variety of methods based on a certain criteria. The values of all classes of each hazard indicator are replaced with hazard contribution or weightage (w) between w ¼ 1 (low) and w ¼ 10 (high) ( Table 5). For example, moderate slopes create favourable settings for flooding (Skilodimou et al. 2019), Therefore, the highest hazard weightage (w ¼ 10) was replaced with index values of slope class range between 32 -38 . Similarly, areas lying in 'major' category of seismic zone create favourable settings for earthquakes, therefore, we replaced with w ¼ 10 or 'high' weightage to the class values lying within 'major' category of seismic zone. The values of class 'maximum likely earthquake magnitudes (5-7)' are replaced with highest weightage as well. Fault-lines that crosscut the district area could have the hardest impact, therefore, we replaced with the highest weightage to this class. Soil related dataset comprise of 8 classes (e.g. Ao72-2b, glaciers (GL), I-Y-2c, etc.). In mountain environment, highly glaciated (GL) regions of Northern Pakistan are an indication of high degree of potential hazard damage to downstream population, therefore, we replaced highest weightage with class representing GL or glaciers. Weightages to other soil types are provided based on the qualitative ordinal scale famously known as Mohs's Hardness Scale (Tabor 1954). Finally, all reclassified hazard maps are combined to produce one single hazard map by using 'weighted overlay' function of ArcGIS. Similarly, for vulnerability analysis, pixels of each class of vulnerability maps are replaced with a contribution weightage between w ¼ 1 (high) and w ¼ 10 (low) ( Table 6). Proximity of villages to hazards and population density are given utmost importance as mountain communities relatively poor often reside on higher elevations close to risky mountain environments and hazard prone glaciers and GLOFs. Child dependency, poverty and topography are considered to be the root causes of vulnerability of mountain communities, therefore, given highest vulnerability weightage followed by those having little agricultural land, lower income opportunities suffer more adversely from disasters like landslides, avalanches and mud-sliding. Therefore, highest weightage is given to child dependency and poverty. Majority of private or public health facilities are available on low elevation areas. Communities living on higher elevations cannot afford transportation charges and highly cost private health facilities due to limited resources. Rugged topography of mountains places them in vulnerable groups prone to disasters. Individual coping capacities maps are combined to produce one single vulnerability map based on common measurement scale are by using the same method. The scale of ability of mountain communities to withstand and cope with the effects of hazardous phenomena depend on the conditions, resources and possession of sources that has the potential to decrease destruction or harm downstream population. Levels of coping capacity-based vulnerability (e.g. highly vulnerable, vulnerable, between safe and vulnerable, safe and extremely safe) are assigned to the reclassified classes obtained from weighted overlay method.

Hazard, vulnerability to hazards and risk estimation
Again, we applied 'weighted overlay' function to combine resulting hazard (Table 5) and vulnerability (Table 6) raster datasets (maps) by using a common measurement scale and equal weights or equal importance. In countries like Pakistan, lack of broad chronological data on disastrous events makes hazard analysis extremely difficult (Freeman and Kunreuther 2002). Secondly, the choice of a method for calculating risk or vulnerability is just as subjective as well as the variables used in these methods. Following risk estimating notations are few examples: This notation revolves around 'conventional risk' (UNISDR, 2004) referring particularly to the physical characteristics of vulnerability ('Degree of loss (from 0% to 100%) resulting from a potentially damaging phenomenon'.) (UNDHA (United Nations Department of Humanitarian Affairs) 1992).
Where 'coping' is the capacity to recover quickly from difficulties (Kelman 2018).
(UNISDR, 2009) Where 'Exposure' is 'People, property, systems, or other elements present in hazard zones that are thereby subject to potential losses'.  Risk ¼ ProbabilityxSeverity Where Severity ¼ Magnitude (potential injury or death) -Mitigation (internal or external strengths or resources).
In this study, 'risk' is referred to the 'relationship between "sensitivity of hazard" and its "likelihood" and between "coping capacity" and "level of exposure to origins of natural hazards"' (see Tables 3 and 4) (see Section 3.1). Levels of hazard, disasters and risk are calculated based on these tables by using reclassification method (see Section 3.2.4). Subsequently, districts are ranked between (0 ¼ low level of exposure or safe) to (10 ¼ high degree of exposure) based on total risk estimates. The weighted overlay table of ArcGIS allowed multiple-criteria analysis between Raster sets of these factors. Raster form of hazard and vulnerability maps are overlaid and pixel values are multiplied to estimate total risk of each district. 3.2.6. Estimation of probability of population at risk Probability of population at risk is calculated as the total risk divided by the population at risk (Van Westen 2014). For example, Gilgit district have an area of $ 4000 km 2 with a total population of 0.33 million people. About 10% population is under 'high' or 'very high' risk of deaths or injuries from natural hazards, the individual risk of being killed or injuries by a natural hazard in Gilgit district is:

Result and analysis
4.1. Risk ranking of districts Figure 5(a,b) shows the graphical output of the risk analysis of ten districts of Gilgit-Baltistan under all potential hazard and vulnerability indicators. Potential severity   index or ranking (10-0-10) is represented by the colourful bar lines. They infer that vulnerability and hazard ranking represented to the middle (5-0-5) of the ranking could be less damaging or favourable and on the extreme right and left (10-6-10) are more damaging or unfavourable. Two separate rankings are provided in each districtwise plot to demarcate the upper and lower limits of hazard and vulnerability estimation. The resulting gap between the decision rankings expands out towards the left and right corners of the graph. Probability of population at risk for all districts is presented in Table 7 and Figure 6. Values of relative likelihood of occurrence of a disaster event for districts are represented between 0.0 (impossibility) and 1.0 (certainty). Overall, districts of Gilgit-Baltistan prone to risks are very diverse and strongly associated to their unique geography, topography, HDI, environment and socio-economic conditions. At regional scale, the possible consequences of mountain communities being exposed to natural hazards and the likelihood of occurring are between 'medium' and 'high risk'. At the same scale, major portion of the study area falls within 'major' category of potential hazard severity scale followed by 'moderate' and 'insignificant'. Moreover, susceptibility of mountain systems or communities is 'vulnerable' or 'highly vulnerable' in terms of coping capacity due to mountainspecific socio-economic constrains and settings.
At district-level, significant portion of population of Shigar district falls within 'very high risk' category of severity 'risk' scale. Specifically, the three districts having the highest estimated risk are Shigar ($9000 km 2 ) with risk index (44) and probability (6.8 Â 10 À4 ), Diamer ($7000 km 2 ) with risk index (43) and probability (1.2 Â 10 À3 ) and Ghanche ($8000 km 2 ) with risk index (36) and probability (2.3 Â 10 À3 ). Shigar Figure 5. Vulnerability and hazard ranking of districts based on potential severity or ranking (10 ¼ highly damaging or unfavourable, 1 ¼ lowest damaging or favourable and 0 ¼ No hazard). Size of each bar represents the severity level of hazard or vulnerability of each district. district lies within the seismic zone and fault-line potential to be impacted by a hazard while Diamer district has lowest natural capacity, insufficient or no alternatives to cope and to reduce the adverse impacts posed by natural stresses. Similarly, based on risk probability index, the three districts having the lowest estimated probability of risks are: Hunza ($11,000 km 2 ) with risk index (16) and probability (2.8 Â 10 À3 ), Kharmang ($2800 km 2 ) with risk index (20) and probability (8.8 Â 10 À4 ) and Ghizer ($11,000 km 2 ) with risk index (26) and probability (4.3 Â 10 À4 ). This trajectory reflects the basic principle that a higher risk increases the robustness of the hazard to vulnerability. Table 7. Estimation of population at risk taking total potential hazard, vulnerability and risk into account.

District
Total hazard (see Figure 7) Total vulnerability (see Figure 8) Total risk (see Figure 6) Total population 1 Estimated population at risk (see Figure 6)  Figure 6. Risk map representing population at risk lies within 5 km buffer zone: population at risk is scaled between 'low' to 'major' based on vulnerability and hazard assessments.

Multi-hazard map
The multi-hazard analysis that we implemented in Section 3.2.1 highlighted the potential hazard severity of districts (Figure 7 and Table 5). The continuous weightages of the five hazard severity assessment maps were classified into eight categories. Three districts namely Gilgit, Shigar and Ghizer lies is the Karakoram are prone to 'very high' degree of hazardous events. Kharmang district lies within the Karakoram portion of Northern Pakistan is 'extremely safe' followed by Astore lies in Himalayan portion and Ghanche in the Karakoram portion. Out of the eight hazard indicators, potential adverse impact cause by seismic zone factor is high in the Gilgit-Baltistan of Northern Pakistan followed by geological, soil and fault-line factors. Furthermore, although the slope gradient of Hunza district is high but the average elevation of Shigar district is the highest among all districts. Diamer district is prone to potential geological impact while the soil composition makes Hunza hazard-prone due to presence of $1300 glaciers. Three fault lines pass through Ghizer district and two Shigar district which make them disaster prone risky regions of Gilgit-Baltistan. Highest numbers of earthquake incidents were recorded in Ghizer district followed by Diamer and Hunza. Potential adverse impact of land cover is high in Shigar district possibly due to highest glaciated region in Upper Indus followed by Hunza and Sheyok.

Multi-vulnerability map
The vulnerability assessment that we implemented in Section 3.2.3 highlighted the potential vulnerability ranking of all districts (Figure 8 and Table 6). The continuous Figure 7. Coping-capacity-based vulnerability map representing exposure of community settlements (black dots), coping capacity between 1-4, 5-6, 7-8, 9-10 and level of vulnerability between 'extremely safe' to 'highly vulnerable'. Five km buffer-zone around settlements or villages (black dots) represents population, forest and agricultural land, health facilities, income range and poverty details.
weightages of the vulnerability are scaled into two categories: 'vulnerable' (red) and 'safe' (green). Multi-vulnerability is scaled into five categories, 'highly vulnerable', 'vulnerable', 'between safe and vulnerable', 'safe' and 'extremely safe'. Overall, the result shows that Ghanche and Kharmang districts of Gilgit-Baltistan lie within the Karakoram, are most vulnerable followed by Shigar and Diamer. The mountain communities' of these three districts have low level of capacities or coping strengths to deal with adverse effects of potential hazardous events. Due to unfavourable circumstances, people are at very high risk. On the other hand, Hunza district of Gilgit division of Northern Pakistan is extremely safe followed by Ghizer and Gilgit due to high level of capacities or coping strengths to deal with adverse effects of potential hazardous events. Overall, out of the nine vulnerability indicators, child dependency is reported to be the highest factor effecting mountain communities of Gilgit-Baltistan followed by poverty. Lowest literacy rate of Ghanche make it most vulnerable district among all districts of Gilgit-Baltistan. Except a single Basin Health Unit, secondary and tertiary health facilities are not available in Kharmang district while lowest life expectancy makes Ghizer disaster vulnerable as well. Highest numbers of earthquake incidents were recorded in Ghizer district followed by Diamer and Hunza districts. Hunza district has lowest agricultural and forest land as compared to other districts of Gilgit-Baltistan. Comparatively fragile but vast agricultural land of Ghizer district makes it highly safe place. The under-utilized land can be cultivated for agricultural production and food security purposes. Figure 8. Multi-hazard map representing a relationship between hazard sensitivity and likelihood of ten districts of Northern Pakistan: severity scale of potential hazard is represented with four classes between 'minor (rare/unlikely)' to 'major/certain to occur)'.

Multi-risk map
The risk assessment that we implemented in Section 3.2.4 highlighted the population at risk in all districts of Northern Pakistan ( Figure 6, Table 7). Final risk map is a representation of multiplication of pixel values obtained from hazard and vulnerability maps or datasets. Population at risk lies within 5 km buffer zone from the population centres are scaled between 'low risk' to 'major risk' based on vulnerability and hazard assessment. The mountain community of Shigar district is at 'major risk' while Kharmang and Ghanche districts are at 'moderate risk'. Community residing in the remote or rural areas (e.g. Darel and Tangir) of Diamer district are at 'major risk' compared to population in the urban areas (e.g. Chilas). Similarly, community residing in the far flung areas (e.g. Ishkoman, Budswat, Chatorkhand) of Ghizer district are at 'major risk' compared to population in the urban areas (e.g. Punial or Yasin).

Discussion
The terrible nature of the losses from the July 5, 2005 earthquake that centred in Muzaffarabad area of Pakistan administered Azad Kashmir forced the Government of Pakistan to take steps like establishment of Earthquake Rehabilitation and Reconstruction Authority (ERRA) in 2005 and National Disaster Management Authority (NDMA) in 2007 to deal with disasters and their management. There is a need for developing a long term disaster preparation and mitigation programs at district or community level in Pakistan generally and mountainous Gilgit-Baltistan particularly to devise disaster emergency response mechanisms to mitigate adverse effects of potential hazards. This study is a first step in this endeavour. Although, the Gilgit-Baltistan Disaster Management Authority (GBDMA) provides emergency response to mountain communities but this study provides a basis for hazard, risk and vulnerability assessments for initiating disaster risk reduction programs in risky districts (e.g. Shigar, Diamer and Ghanche). Especially, identification of geographically rugged areas and most susceptible and vulnerable mountain communities under the threat of potential hazards, will greatly contribute to assess the ability of those sectors to withstand and cope with the effects of hazardous phenomena. Moreover, aggregation and analysis of scattered spatio-temporal data to scale the districts based on their geographic location, geological formation, land cover and other parameters is a new addition to the hazard, vulnerability and risk related knowledge of this region. But data scarcity is one of the major hindrances. Especially, lack of scientific data available with all the organizations responsible for disaster management is one of the major issues to initiate mitigation plans and preparation studies.
Specification, based on our result, we assume that potential risks especially seismicity and poor socio-economic conditions of Shigar, Diamer and Ghanche districts are likely to worsen existing challenges and increase some prevailing risks, and to result in the generation of new secondary risks like unemployment, displacement and migration. If we project current vulnerability status of the districts of Gilgit-Baltistan across the mountain communities of HKH, along with existing and earlier disaster events, new natural disasters are likely to occur in coming years. We opined that socio-economic conditions of mountain communities especially increased density of vulnerable population pockets in risk prone regions, will contribute to rising indirect losses. Communities living in these mountain regions have a long history of adjusting themselves to climatic-induced hazards. However, recent socio-economic inflation and mountain environmental degradation (e.g. deforestation), construction of giant dams (e.g. Dassu and Diamer-Basha) and constructions of alternatives routes in Diamer, are likely to aggressively effect mountain eco-systems, and could increase the challenges of downstream mountain communities to protect their livelihood opportunities.
Currently employed mechanism for vulnerability assessment is adequate for scaling of districts based on their mitigation strengths (e.g. income, agricultural land, literacy rate, health facilities, etc.) and weaknesses (e.g. child dependency, population density, poverty, etc.) but inadequate to precisely find disaster prone spots within the districts. To compensate limited number of hazard-specific indicators, a number of assumptions were made in this study, setting certain limitations on the results. However, they are sufficient for a comparison and estimation of differences between multi-hazard and risk scaling of districts and risk prone regions. Secondly, the multi-hazard risk ranking was established at district level. However, the criteria (weightage) allocated to various indicators is a subjective process, therefore, the judgments and observations involved might vary from other studies. We conclude that this study identify only the most susceptible districts and population at risk within 5 km buffer zone encompassing from human settlements/villages, devoid of sufficient spatial resolution for the identification of vulnerable settlements or villages within the districts. An integration of environmental, socio-economic and physical dimensions of vulnerability is used for the estimation of the risk ranking and coping capacity like Habib et al. (2017). However, due to the versatile nature of vulnerability as well as data limitation, an all-inclusive approach is mandatory.

Conclusion
The applied multi-hazard multi-vulnerability assessment within the scope of this study achieved relatively reliable results regarding the ranking of districts that are highly vulnerable or unfavourable and least vulnerable or most favourable with respect to seismicity, geological formation, land cover, topographic features and socio-economic status of the mountain communities. We conclude that the mountain communities of the Karakoram in the Northern Pakistan are vulnerable due to their spatial locality, low socio-economic capacity, insufficient or no alternatives to cope and to reduce the adverse impacts posed by natural stresses. Specifically, the multihazard map indicates that the high or very high unfavourable areas (prone to hazards) are located primarily in the central part (Gilgit, Ghizer and Diamer) of the study area, while the entire Baltistan lies in the Central Karakoram have lower level of exposure to natural hazards but higher level of exposure to socio-economic vulnerabilities. Spatial distribution and exact demarcation of the seismic-zones remained uncertain due to small-scale (1:12,000,000) of the seismic-tectonic maps provided by the Geological Survey of Pakistan. For large scales, the approach could identify the smaller areas within each district that might be prone to natural hazards. Standard multi-hazard multi-vulnerability risk assessment procedures are desirable to study the HKH region as a multi-hazard and multi-vulnerable mountainous environment. Successful disaster management plans are desirable in not only affected regions solely by natural hazards (Shigar, Gilgit, Diamer) or vulnerability (Kharmang, Diamer, Ghanche), but across the entire Northern Pakistan. To execute such plans, the planners need disaster management assessment tools for mountain regions as well. Lastly, the reasons for relying on datasets from other sources are due to the lack of up-todate vulnerability-related datasets (socio-economic, local coping strategies, etc.), which not only constrained our approach to correctly identify population at risk at the village levels, but also limited the number of possible indicators. We conclude that response mechanism, preparedness plans, coordination among line departments of GB, information management, establishment of early warning systems, resource mobilization, public training and education and community-based disaster preparedness, are required for disaster preparedness in Northern mountain regions.

Data availability statement
Authors have attached the zip file of data and materials that supported the results or analyses presented in their paper. Additionally, data obtained from other sources are cited and referenced in our manuscript. Some data that support the findings of this study related to vulnerability will be made available on request from the corresponding author.

Disclosure statement
We wish to confirm that there are no known conflicts of interest associated with this publication and there has been no significant financial support for this work that could have influenced its outcome.