Integrated Approach for Landslide Risk Assessment Using Geoinformation Tools and Field Data in Hindukush Mountain Ranges, Northern Pakistan

: Landslides are one of the most recurring and damaging natural hazards worldwide, with rising impacts on communities, infrastructure, and the environment. Landslide hazard, vulnerability, and risk assessments are critical for landslide mitigation, land use and developmental planning. They are, however, often lacking in complex and data-poor regions. This study proposes an integrated approach to evaluate landslide hazard, vulnerability, and risk using a range of freely available geospatial data and semi-quantitative techniques for one of the most landslide-prone areas in the Hindukush mountain ranges of northern Pakistan. Very high-resolution satellite images and their spectral characteristics are utilized to develop a comprehensive landslide inventory and predisposing factors using bi-variate models to develop a landslide susceptibility map. This is subsequently integrated with landslide-triggering factors to derive a Landslide Hazard Index map. A geospatial database of the element-at-risk data is developed from the acquired remote sensing data and extensive ﬁeld surveys. It contains the building’s footprints, accompanied by typological data, road network, population, and land cover. Subsequently, it is analyzed using a spatial multi-criteria evaluation technique for vulnerability assessment and further applied in a semi-quantitative technique for risk assessment in relative risk classes. The landslide risk assessment map is classiﬁed into ﬁve classes, i.e., very low (18%), low (39.4%), moderate (26.3%), high (13.3%), and very high (3%). The developed landslide risk index map shall assist in highlighting the landslide risk hotspots and their subsequent mitigation and risk reduction.


Introduction
Around the globe, landslides are becoming more frequent and devastating, with 200,000 casualties in the twentieth century and an estimated exposure of around 3.7 million km 2 area and 300 million population [1][2][3][4].Given changes in the global climate and microclimate, increases in population, urbanization, ever-expanding communication networks, and other infrastructural development on fragile slopes, the landslide-induced risk to society and the economy is growing [5][6][7][8].The assessment of landslide risk is critical to demarcate landslide-prone areas, infrastructure, and communities and accordingly plan and implement risk reduction strategies.Space and airborne remote sensing offer a variety of geospatial data for Geographic Information Systems (GIS) to perform landslide hazard, vulnerability, and risk assessments from a local to regional scale [9][10][11].
Landslide susceptibility is the probability of the spatial distribution of landslides driven by the topography, geology, tectonics, climate, land cover, and human activities [12,13].A landslide susceptibility analysis begins by developing a comprehensive landslide inventory, which is often developed from space and airborne remote sensing data using a range of image classification techniques [14,15].Techniques for landslide susceptibility assessment can be categorized in the evaluation of landslide inventories and include heuristic techniques, geomorphological mapping, and physical and statistical modeling approaches [16][17][18].Usually, the landslide susceptibility assessment is further analyzed with the temporal distribution, landslide frequency-magnitude relationship and major triggers, including rainfall and earthquakes, for the landslide hazard assessment [19][20][21].The population, buildings, economic activities, public utilities, infrastructure, and environmental features are often impacted by landslides and thus considered as the elements at risk [22][23][24][25][26][27].The degree of loss for elements at risk is determined by the landslide vulnerability assessment, which is influenced by the landslide type and magnitude and the nature of the exposed elements [22,28,29].The typological information of the buildings, including their foundation type, building material, design etc., strongly affects the vulnerability of the buildings to landslides and other hazards [26].Landslide risk refers to the expected damage to the people, property, and economy in a region and a given period [30].
Landslide risk can be assessed through qualitative and quantitative techniques, where the qualitative methods mostly rely on subjective evaluation, and the quantitative methods aim to quantify the landslide-induced losses [31].Contrary to the qualitative risk assessment that leads to relative rankings (high, moderate, low), quantitative risk assessment (QRA) leads to quantifying the landslide-induced risk in terms of integrated costs of the exposed element at risk [31].QRA is effective in comparing the landslide-prone sites and administrative units, and therefore, it is important for community awareness, resource prioritization, landslide mitigation planning, and policy development for landslide risk reduction [26,32,33].However, a comprehensive QRA requires a range of data on the spatial and temporal distribution of landslides with associated characteristics, historical records of the potential landslide triggers, calculated costs of the elements at risk, and scenarios of human activities with time, which are rarely available in the developing countries [34][35][36].Alternatively, semi-quantitative risk assessment techniques have been developed that apply weighting and rating of the input parameters, which are eventually integrated to compute the risk [35,37].A semi-quantitative risk assessment approach is encouraging in landslide studies [38,39].Crawford et al., (2022) [40] proposed that a practical level exists between qualitative and quantitative risk assessments for a significant landslide risk assessment at a regional scale with limitations regarding hazard behavior and asset data.
Given the rough terrain, active tectonic and denudation, glaciation, climatic conditions, and human factors, landslides are one of the most damaging natural disasters in the Hindukush mountains in northern Pakistan [41][42][43].Landslides in the region are investigated by various researchers using field surveys, geological and tectonic records, and geological settings to derive the landslide susceptibility map [42, [44][45][46][47][48].However, despite the increasing landslide frequency and associated damages that are exacerbated by climate change, an integrated risk assessment is lacking for the region.The aim of the present study is to utilize the freely available geospatial and extensive field data and integrate these data in a semi-quantitative approach for landslide risk assessment.The proposed study is initiated with quantitative data for hazard and vulnerability assessments, and the derived risk is classified into relative classes from very high to very low risk levels.The developed methodology can be replicated in other landslide-prone regions to assist in risk reduction.

Study Area
The study area is located in the Hindukush mountain ranges in the Ghizer district of northern Pakistan, with a total area of 7040 km 2 and elevation ranging from 1672 to 6239 m above sea level (ASL) (Figure 1).The study area is located in the Hindukush mountain ranges in the Ghizer district of northern Pakistan, with a total area of 7040 km 2 and elevation ranging from 1672 to 6239 m above sea level (ASL) (Figure 1).The area is located in an active tectonics zone attributed to the collision between the Eurasian plate and the Indian tectonic plate.The area lies at the Kohistan Island arc and Karakorum block, with the Main Karakorum Thrust between them.The lithology of the area is comprised of the Kohistan Batholiths, Chalt Volcanic group, Greenstone complex, Quaternary deposits, Darkot group, and Yasin group [49].The topography of the area is characterized by steep to gentle slopes favoring widespread landslides, rock falls, scree slopes and debris flows [50][51][52] (Figure 2).The area is located in an active tectonics zone attributed to the collision between the Eurasian plate and the Indian tectonic plate.The area lies at the Kohistan Island arc and Karakorum block, with the Main Karakorum Thrust between them.The lithology of the area is comprised of the Kohistan Batholiths, Chalt Volcanic group, Greenstone complex, Quaternary deposits, Darkot group, and Yasin group [49].The topography of the area is characterized by steep to gentle slopes favoring widespread landslides, rock falls, scree slopes and debris flows [50][51][52] (Figure 2).

Materials
The landslide boundaries and road network were visually digitized from Google Earth Pro and subsequently verified in the field.The resolution of the data has a signifi-

Materials
The landslide boundaries and road network were visually digitized from Google Earth Pro and subsequently verified in the field.The resolution of the data has a significant impact on the parameters for deriving the causative factors of landslides, landslide hazards, presentation of the element-at-risk database, and terrain representation [53,54] The ALOS PALSAR Digital Elevation Model (DEM) with a 12.5 m spatial resolution was exploited to compute the topographic features, including the terrain elevation, aspect, slope, curvature, and stream network (Table 1).The geology and fault lines were digitized from the geological map after [49].The landcover map was developed from the Sentinel-2 data (https:// sentinel.esa.int/web/sentinel/missions/sentinel-2,accessed on 25 July 2022).For landslide hazard assessment, the precipitation record was derived from the global precipitation measurements (GPM), and the peak ground acceleration (PGA) values were adopted from [55].For the landslide exposure and vulnerability analysis, the population data were acquired from the 2017 census records.The building footprints were visually demarcated from the base map in the ArcGIS software, and the typological information, including building use (houses, schools, and hospitals), roof and wall materials, physical condition, and age of the buildings, was collected for each of the buildings through an extensive field survey and subsequently utilized for the landslide vulnerability and risk assessments.A landslide inventory depicts the type, volume, magnitude, date, and place of occurrence of landslides in an area.A landslide inventory of the study area was developed through the visual interpretation of Google Earth images following [56].The landslides were demarcated considering their morphology, changes in vegetation, and the presence of eroded deposits and subsequently updated in the field and classified after [57].Among the mapped landslides, 70% of the landslide data were used for susceptibility modeling, and 30% were used for validation.

Landslide Causative Factors and Susceptibility Assessment
Considering the influence on the spatial distribution of landslides in the area, nine landslide causative parameters were selected for the susceptibility mapping.Elevation has an influence on the weather, drainage, and other factors influencing the spatial distribution of landslides [58,59] (Figure 3a).Terrain slope determines the shear forces acting on hill slopes and thus impacting the occurrence of landslides [60][61][62] (Figure 3b).The terrain aspect (Figure 3c) is the measurement of the direction of the slope, impacting the vegetation cover, moisture retention, and strength of the soil [63,64].Terrain curvature (Figure 3d) is the change in slope along a small arc of the curve and is often applied as a controlling factor for landslides [65].Water runoff contributes to slope stability through infiltration, saturation and slope toe-cutting [66,67] (Figure 3e).Excavation for road construction and repair in mountainous areas often destabilizes the slope and triggers landslides [68] (Figure 3f).The Euclidean distance from the selected parameters is classified into five classes following literature and field observations.Variation in land cover affects the hydrological conditions, shear strength and slope instability [41].A land cover map of the area was developed from the Sentinel-2 data using the maximum livelihood classification algorithm and classified into six categories, i.e., water bodies, settlements, barren land, snow, forest, and agriculture (Figure 3g).Proximity to a fault (Figure 3h) was reclassified following [61].Lithology and faults influence the rock's strength and largely drive the spatial concentration of landslides [69].The area of study is comprised of eight lithological units (Figure 3i).For the landslide susceptibility assessment, the frequency ratio (FR) approach evaluating the relationship between the landslide inventory and the selected causative factors is frequently and effectively applied [70][71][72][73][74].The frequency ratio (FR) is the ratio of the landslide's occurrence area to the total area and leads to higher accuracy compared to other models [46,[75][76][77].

Landslide Triggers and Hazard Mapping
The intensity and duration of rainfall are among the major triggers for widespread landslides and are often used for landslide hazard assessments [41].The area-specific threshold of rainfall intensity and duration are often used to predict landslides [78].Global precipitation measurement (GPM) data for seven years (1st January 2015 to 31st December 2021) were utilized with mean rainfall ranges from 1002 to 2873 mm/year [79].The earthquakeinduced shaking depicted through the peak ground acceleration (PGA) leads to widespread co-seismic landslides.
For the landslide hazard assessment, the temporal distribution of landslides and their frequency-magnitude records are important; however, they are rarely available in data-poor regions, such as the study area [35].Due to the lack of temporal and intensity records of landslides, a static hazard map was developed in this research work.Following [35,[80][81][82][83], weights of 0.3 to PGA and 0.7 to rainfall are assigned and integrated into the landslide susceptibility to derive the landslide hazard map.
where TI is the triggering indicator, PGA is the peak ground acceleration, and RF is the rainfall precipitation.The exposure of buildings and roads to landslides was assessed through an overlay analysis [31].

Landslide Vulnerability and Risk Assessment
Vulnerability indicates the susceptibility of the exposed elements to the impacts of landslides [22,82,[84][85][86].Physical vulnerability refers to damage to buildings, utilities, and infrastructure, social vulnerability is the proneness of people and relates to fatalities or injuries [26], and environmental vulnerability refers to the impact on land cover, i.e., stream environments, forest cover, agriculture, etc.In this study, physical, social, and environmental vulnerabilities were derived and, using the multi-criteria analysis approach, integrated into a vulnerability index map (Figure 4).

Physical Vulnerability
An extensive database of 23456 building footprints was developed through onscreen digitization on the ArcGIS base map and subsequently verified in the field.For each building, the building use and associated typological data comprising the structural system, geometry, material properties, state of maintenance, levels of design codes, foundation

Physical Vulnerability
An extensive database of 23,456 building footprints was developed through onscreen digitization on the ArcGIS base map and subsequently verified in the field.For each building, the building use and associated typological data comprising the structural system, geometry, material properties, state of maintenance, levels of design codes, foundation and superstructure details, number of floors, and other factors, were recorded during field visits [26].Applying a multi-criteria analysis, weights were assigned to each parameter using expert opinion and AHP after [87] (Table 2).Recently constructed buildings with no structural damage are generally more resistant to hazards compared to older and damaged buildings [88].Therefore, buildings with no structural damage were assigned a weight of 0.15, and buildings with minor structural damage were assigned a weight of 0.35; those with major structural damage were assigned a weight of 0.5.Dry rubble stone without mortar was given a weight of 0.4, dressed stone masonry with mortar was given a weight of 0.2, cement block masonry was given a weight of 0.3, and reinforced confined masonry was assigned a weight of 0.1.Horizontal metallic sheets supported from underneath were given a weight of 0.4, metallic sheets with pitched sloping were given a weight of 0.3, wood planks covered by mud were given a weight of 0.2, and reinforced concrete cement was assigned a weight of 0.1.Considering the use of the building, educational institutes, i.e., schools/colleges were assigned a weight of 0.35, hospitals/dispensaries/clinics were assigned a weight of 0.25, houses were assigned a weight of 0.2, commercial buildings/hotels/offices were assigned a weight of 0.13, and worship places were assigned a weight of 0.07.Mapped roads are classified as main roads traversing through the main valley and link roads in the side valleys.The main roads were assigned a weight of 0.7, and the linking roads were assigned a weight of 0.3.The vulnerability of roads was integrated into the buildings to generate a physical vulnerability map.

Social and Environmental Vulnerability
Social vulnerability refers to the likelihood of human losses determined by the nature, magnitude, and intensity of the landslide and the coping capacity of the community [89].Using a dissymmetric mapping approach, the average population per building was computed.Following [83], the population data were divided into four classes ranging from low to very high population, and accordingly, weights were assigned.
Using the developed land cover map, the highest weightage was assigned to agriculture and forest, i.e., 0.5 and 0.4, respectively, due to their environmental and economic significance.Water bodies were given a weight of 0.1.The physical, social, and environmental vulnerabilities were integrated to develop a vulnerability index map.

Landslide Risk Assessment
A landslide risk map was developed by integrating the developed landslide hazard and vulnerability index maps [33,90].For landslide risk assessment, the landslide hazard and vulnerability indexes were classified into five classes ranging from very low to very high.The landslide risk index was developed by integrating the normalized landslide vulnerability map with the landslide hazard map through multi-criteria analysis following Equation (2) and after [83].
where R is the landslide risk index, H is the landslide hazard, and V is the landslide vulnerability index.The derived risk values are normalized and classified into 5 classes, including very high, high, moderate, low, and very low.The classes are the relative ranking of the risk classes.The very high risk is the riskiest area with very high hazard and very high vulnerability.The elements at risk overlay these very high risks and are very prone to damage due to landslides.

Landslide Inventory
In the studied area, 1748 landslides were mapped with a total area of 344 km 2 (Figure 5).The types of landslides distributed in the area included rock falls (102), rockslides (62), scree (1382), debris flows (78), and debris cones (124).Landslides are largely concentrated along streams.A few large rock falls were identified in the study area along the Gilgit river, leading to the blockage of the river and natural dams, which were later drained due to overflow.The rock avalanches are composed of boulders and scree deposits of the Cretaceous to Jurassic age volcanic and volcanoclastic sediments, containing greenschist and slates.scree (1382), debris flows (78), and debris cones (124).Landslides are largely concentrated along streams.A few large rock falls were identified in the study area along the Gilgit river, leading to the blockage of the river and natural dams, which were later drained due to overflow.The rock avalanches are composed of boulders and scree deposits of the Cretaceous to Jurassic age volcanic and volcanoclastic sediments, containing greenschist and slates.

Landslide Causative Factors
Landslides in the area are largely concentrated along the slope ranging from 20-40°.The frequency ratio (FR) for the landslides reaches up to 1.27 for the slope class of 20-30° and 1.15 for the 30-40° class.Landslides in class < 10° and >60° have the lowest FR.The slope aspects toward the north, northeast, west, south, and southwest show the highest FR and lowest for flat areas.The FR is at its maximum for the elevation classes of <2500, 2500-3000, and 3000-3800, and it decreases in higher elevations, which are usually covered with glaciers.The FR for concave and flat terrain is 1.03, and the FR for convex terrain is 0.68.Landslides have an inverse relation with distance from streams.In summers, glaciers melt, increasing water in the streams, and monsoon rainfalls also add more water, triggering slope-toe erosion that often leads to landslides.The stream class of 100-200 m has the highest frequency ratio of 1.82, and the ratio decreases to 0.80 as the distance from the stream decreases (Figure 6).The road distance buffer zones along the roads have the highest FR values, ranging from 0.96-2.81due to slope cut-offs, uncontrolled blasting during construction, and vibration due to heavy traffic.The distance from fault lines shows that the highest FR is computed for <100 m, i.e., 1.49, reflecting increased landslides along the faults.The geological units of the Shyok Suture Zone (Sv) and Kohistan Batholiths (KB) are highly susceptible to landslides, with the highest FR ratio, i.e., 1.64 and 1.13, respectively.The barren land landcover class has an FR value of 1.8 (Figure 6).

Landslide Causative Factors
Landslides in the area are largely concentrated along the slope ranging from 20-40 • .The frequency ratio (FR) for the landslides reaches up to 1.27 for the slope class of 20-30 • and 1.15 for the 30-40 • class.Landslides in class < 10 • and >60 • have the lowest FR.The slope aspects toward the north, northeast, west, south, and southwest show the highest FR and lowest for flat areas.The FR is at its maximum for the elevation classes of <2500, 2500-3000, and 3000-3800, and it decreases in higher elevations, which are usually covered with glaciers.The FR for concave and flat terrain is 1.03, and the FR for convex terrain is 0.68.Landslides have an inverse relation with distance from streams.In summers, glaciers melt, increasing water in the streams, and monsoon rainfalls also add more water, triggering slope-toe erosion that often leads to landslides.The stream class of 100-200 m has the highest frequency ratio of 1.82, and the ratio decreases to 0.80 as the distance from the stream decreases (Figure 6).The road distance buffer zones along the roads have the highest FR values, ranging from 0.96-2.81due to slope cut-offs, uncontrolled blasting during construction, and vibration due to heavy traffic.The distance from fault lines shows that the highest FR is computed for <100 m, i.e., 1.49, reflecting increased landslides along the faults.The geological units of the Shyok Suture Zone (Sv) and Kohistan Batholiths (KB) are highly susceptible to landslides, with the highest FR ratio, i.e., 1.64 and 1.13, respectively.The barren land landcover class has an FR value of 1.8 (Figure 6).

Landslide Susceptibility and Hazard Mapping
The FR weights were integrated to develop a Landslide Susceptibility Index map (LSI) and classified into five classes (Figure 7a).The assessment of the accuracy of the susceptibility map was performed using the 30% dataset of landslides.The area under the curve (AUC) is a useful indicator to validate the prediction performance of the model [91].The susceptibility map was classified into 100 classes with a 1% cumulative interval.Subsequently, the validation landslide set was correlated with the classified susceptibility map to identify the landslide area in each of the 100 susceptible classes.Qualitative analysis of the area under the curve of the success rate curve shows 88% accuracy for the susceptibility map (Figure 7b).The landslide susceptibility is further integrated with the weights of the landslide triggers (precipitation and PGA) to compute the landslide hazard map (Figure 8).The landslide hazard map indicates that 14% of the area lies in the very low hazard class, while 29.7% lies in the low hazard class, 30.3% lies in the medium hazard class, 19.2% lies in the high hazard class, and 6.6% is in the very high hazard class.

Landslide Susceptibility and Hazard Mapping
The FR weights were integrated to develop a Landslide Susceptibility Index map (LSI) and classified into five classes (Figure 7a).The assessment of the accuracy of the susceptibility map was performed using the 30% dataset of landslides.The area under the curve (AUC) is a useful indicator to validate the prediction performance of the model [91].The susceptibility map was classified into 100 classes with a 1% cumulative interval.Subsequently, the validation landslide set was correlated with the classified susceptibility map to identify the landslide area in each of the 100 susceptible classes.Qualitative analysis of the area under the curve of the success rate curve shows 88% accuracy for the susceptibility map (Figure 7b).The landslide susceptibility is further integrated with the weights of the landslide triggers (precipitation and PGA) to compute the landslide hazard map (Figure 8).The landslide hazard map indicates that 14% of the area lies in the very low hazard class, while 29.7% lies in the low hazard class, 30.3% lies in the medium hazard class, 19.2% lies in the high hazard class, and 6.6% is in the very high hazard class.
The spatial comparison of the buildings and roads with the hazard map shows that 118 buildings (0.38%) are located in the very low hazard zone, 2838 buildings (12.1%) are located in the low hazard zone, 7037 buildings (30%) are located in the moderate hazard zone, 8233 (35.1%) are located in the high hazard zone, and 5230 buildings (22.3%) are located in the high hazard class.The exposure assessment of roads indicates that 33% of roads are in the very high hazard zone, 27% of roads are in the high hazard zone, 20% are in the moderate hazard zone, 13% are in the low hazard zone, and 7% are in a very low hazard zone.The spatial comparison of the buildings and roads with the hazard map shows that 118 buildings (0.38%) are located in the very low hazard zone, 2838 buildings (12.1%) are located in the low hazard zone, 7037 buildings (30%) are located in the moderate hazard zone, 8233 (35.1%) are located in the high hazard zone, and 5230 buildings (22.3%) are located in the high hazard class.The exposure assessment of roads indicates that 33% of roads are in the very high hazard zone, 27% of roads are in the high hazard zone, 20% are in the moderate hazard zone, 13% are in the low hazard zone, and 7% are in a very low hazard zone.

Physical Vulnerability
The typological information of 23456 buildings mapped in the area (Table 3) indicates that the wall materials were dry rubble stone without mortar in 9645 buildings (41%), dressed stone masonry with mortar in 5741 buildings (24%), cement block masonry in 7867 buildings (34%), and reinforced concrete masonry in 203 buildings (1%).The roof materials were horizontal metallic sheets supported from underneath by wood/steel studs in 6310 buildings (27%), metallic sheets with pitched sloping in 9437 buildings (40%), wood planks covered by mud in 5709 buildings (24%) and reinforced concrete cement (RCC) in 2055 buildings (9%).There were 16996 buildings with different physical conditions showing no structural damage (72%), 4383 with minor cracks in the walls/plaster and broken seepage (19%), and 2077 with cracks in the beams, reinforcement exposed, and major cracks in the walls and roofs (9%).The ages of the buildings show that 5537 (24%) buildings are less than 10 years, 3861 (16%) buildings are 11-20 years of age, 5931 (25%) buildings are 21-30 years old, 4929 (21%) buildings are 31-40 years of age, and 3891 (14%) buildings are older than 40 years.Although a large part of the study area lies in the low and intermediate hazard classes, a large portion of the population is situated in high and very high hazard zones.Buildings with weak physical conditions, older age, and walls and roofs with building materials having less resistance occur in the high vulnerability class.All the educational institutes (schools and colleges), hospitals, and dispensaries are located in very highly vulnerable classes.Following the multi-criteria analysis technique and the weights derived from Table 2, a physical vulnerability map was  The typological information of 23456 buildings mapped in the area (Table 3) indicates that the wall materials were dry rubble stone without mortar in 9645 buildings (41%), dressed stone masonry with mortar in 5741 buildings (24%), cement block masonry in 7867 buildings (34%), and reinforced concrete masonry in 203 buildings (1%).The roof materials were horizontal metallic sheets supported from underneath by wood/steel studs in 6310 buildings (27%), metallic sheets with pitched sloping in 9437 buildings (40%), wood planks covered by mud in 5709 buildings (24%) and reinforced concrete cement (RCC) in 2055 buildings (9%).There were 16996 buildings with different physical conditions showing no structural damage (72%), 4383 with minor cracks in the walls/plaster and broken seepage (19%), and 2077 with cracks in the beams, reinforcement exposed, and major cracks in the walls and roofs (9%).The ages of the buildings show that 5537 (24%) buildings are less than 10 years, 3861 (16%) buildings are 11-20 years of age, 5931 (25%) buildings are 21-30 years old, 4929 (21%) buildings are 31-40 years of age, and 3891 (14%) buildings are older than 40 years.Although a large part of the study area lies in the low and intermediate hazard classes, a large portion of the population is situated in high and very high hazard zones.Buildings with weak physical conditions, older age, and walls and roofs with building materials having less resistance occur in the high vulnerability class.All the educational institutes (schools and colleges), hospitals, and dispensaries are located in very highly vulnerable classes.Following the multi-criteria analysis technique and the weights derived from Table 2, a physical vulnerability map was developed and subsequently classified into five classes ranging from very low to very highly vulnerable (Figure 9).The vulnerability of roads indicates that the main road is highly vulnerable to landslides because of the highly hazardous zone along the main road and the high traffic passing through it.The linking roads are less vulnerable due to less traffic and are largely located in the low to medium vulnerability class.

Social Vulnerability
Population density is directly related to residential buildings.In the study area, locations such as the Gahkuch, Sherkilla, Gupis, and Teru are more populated, leading to high social vulnerability (Figure 10a).The social vulnerability map was classified into five classes (very low, low, medium, high, and very high).Of the total population of 117,280, 14,425 (12.3 %) are in the very highly vulnerable class, 29,785 (25.4 %) people are in the highly vulnerable class, 34,011 (29 %) people are in the medium class, and the remaining 39,288 (33.5 %) are in very low and low vulnerable classes.

Social Vulnerability
Population density is directly related to residential buildings.In the study area, locations such as the Gahkuch, Sherkilla, Gupis, and Teru are more populated, leading to high social vulnerability (Figure 10a).The social vulnerability map was classified into five classes (very low, low, medium, high, and very high).Of the total population of 117,280, 14,425 (12.3 %)

Social Vulnerability
Population density is directly related to residential buildings.In the study area, locations such as the Gahkuch, Sherkilla, Gupis, and Teru are more populated, leading to high social vulnerability (Figure 10a).The social vulnerability map was classified into five classes (very low, low, medium, high, and very high).Of the total population of 117,280, 14,425 (12.3 %)

Environmental Vulnerability
Forest and agriculture have high environmental and economic impacts, and therefore, the highest weights were assigned to these areas.Barren land, snow, and built-up land are in the medium to low vulnerable classes in the environmental vulnerability map (Figure 10b).
The physical, social, and environmental vulnerability maps were combined into an integrated vulnerability map (Figure 11).Forest and agriculture have high environmental and economic impacts, and therefore, the highest weights were assigned to these areas.Barren land, snow, and built-up land are in the medium to low vulnerable classes in the environmental vulnerability map (Figure 10b).
The physical, social, and environmental vulnerability maps were combined into an integrated vulnerability map (Figure 11).The landslide vulnerability map and landslide hazard map were integrated by assigning 0.6 weightage to the hazard map and 0.4 weightage to the vulnerability map to compute the landslide risk assessment map (Figure 12).The landslide risk map was re-

Landslide Risk Map
The landslide vulnerability map and landslide hazard map were integrated by assigning 0.6 weightage to the hazard map and 0.4 weightage to the vulnerability map to compute the landslide risk assessment map (Figure 12).The landslide risk map was reclassified into five classes according to the natural break method, namely, very low risk, low risk, moderate risk, high risk, and very high risk.In the study area, 18% of the area is in the very low-risk class, 39.4% is in the low-risk class, 26.3% is in the moderate-risk class, 13.3% is in the high-risk class, and 3% of the area is in the very high-risk class.

Discussion
This is one of the first studies on landslide risk assessment in the complex topographic, climatic, and tectonic environment of northern Pakistan.Landslides are one of the most frequent and damaging hazards in the region; however, the available literature on the area is limited to landslide susceptibility assessments [47,[92][93][94][95][96][97].From the global perspective, most of the landslide risk papers employing semi-quantitative risk assessments utilize the elements of at-risk data from secondary sources and often map the settlements collectively as large polygons without the demarcation of individual houses.In this study, 23456 buildings were mapped separately, and subsequently, their typological information was collected from extensive fieldwork to assess the landslide vulnerability In the study area, 18% of the area is in the very low-risk class, 39.4% is in the low-risk class, 26.3% is in the moderate-risk class, 13.3% is in the high-risk class, and 3% of the area is in the very high-risk class.

Discussion
This is one of the first studies on landslide risk assessment in the complex topographic, climatic, and tectonic environment of northern Pakistan.Landslides are one of the most frequent and damaging hazards in the region; however, the available literature on the area is limited to landslide susceptibility assessments [47,[92][93][94][95][96][97].From the global perspective, most of the landslide risk papers employing semi-quantitative risk assessments utilize the elements of at-risk data from secondary sources and often map the settlements collectively as large polygons without the demarcation of individual houses.In this study, 23456 buildings were mapped separately, and subsequently, their typological information was collected from extensive fieldwork to assess the landslide vulnerability and risk to each of the buildings considering their spatial location in the landslide hazard class and typological information determining its vulnerability to landslides.In this paper, specific buildings are identified that are high risk, as well as the underlying reasons for these classifications.
The number of landslides in this study varies from other studies in the area, including [46,98,99].A comprehensive landslide inventory was developed from the VHR images, which were subsequently verified and updated in the field.The developed landslide inventory is largely concentrated at the lower elevations.In this area, landslides rarely exist at higher elevations due to the presence of snow cover and glaciers.The majority of the landslides in the study area are scree slopes, rock falls, and debris flows.Debris flows are very fast to extremely fast-flowing slides of saturated debris that occur in steep channels, with a significant impact on the downstream settlements and road networks [57,100,101].Considering the fine resolution of images from Google Earth Pro, a detailed landslide inventory can record smaller landslides, which are often difficult to map from the relatively coarser resolution satellite images provided by Landsat, ASTER and MODIS (Moderate Imaging Spectroradiometer) [56,99].Moreover, it is often convenient to visualize landslides in 3D through Google Earth Pro to precisely demarcate their boundaries.
The frequency ratio model used for the susceptibility assessment in this study lead to higher accuracy than the other frequently applied techniques, including the weight of evidence, logistic regression, and artificial neural networks [94].The critical issue for the landslide hazard assessment is largely the lack of temporal records of landslides in developing countries, and therefore, the standard procedure for the landslide hazard assessment using the spatial, temporal, and magnitude/size probabilities is often difficult to apply in the study area [102].Eventually, following [90,103], we integrated the landslide triggers, including the rainfall and PGA, to acquire the landslide hazard map.However, the inclusion of the temporal probability of landslides shall improve the reliability of the landslide hazard map.
A vulnerability assessment methodology relies heavily on the typology of exposed elements.The structural system, geometry, material properties, state of maintenance, levels of design codes, foundation and superstructure details, number of floors, and other factors are among the typological parameters that determine a building's ability to withstand landslide actions [26].Researchers have used different approaches for vulnerability assessment.Ref. [104] proposed the back analysis of real event damage data to acquire the correlations between the intensity of landslides and building vulnerability through regression analysis.
Ref. [105] calculated vulnerability curves as a function of the differential settlements of a reinforced concrete frame building.Ref. [106] developed a debris flow intensity index that considers the flow height and velocity to calculate the probability of damage.Ref. [88] generated the vulnerability of buildings according to the type of building structure by the visual interpretation of high-resolution imagery, while current research work interprets individual buildings by performing building typology.
For landslide risk assessments, many researchers have focused on the quantitative assessment of vulnerability [107][108][109].Refs.[32,34,110] performed qualitative landslide risk assessments at a regional scale.In this study, a semi-quantitative approach was adopted for risk assessment through the analysis of the collected extensive field-based data of the element-at-risk and applied methods.The landslide hazard and vulnerability maps were derived using quantitative approaches.However, to integrate the hazard and vulnerability maps into a risk map, both maps were normalized and then integrated into the risk map.Field visits verified that Gahkuch, the most populated zone in the study area with a high hazard probability of landslide, mostly lies in the moderate-to high-risk classes.Similarly, Phander and Teru (Figure 12) are the other populated areas that lie in the high-risk class.The quantitative validation of a risk index is often ignored due to the lack of high-quality data [111] and, therefore, could not be included in this study; however, it shall be accomplished in future quantitative landslide risk assessment studies in the region.This study is an intermediate methodology between purely qualitative and quantitative approaches, including quanitative data, and providing relative risk levels.As the performance of risk mapping is largely dependent on data availability, the lack of temporal probability and more detailed data for elements at risk lead to some limitations in the research work.

Conclusions
In this research work, we conducted a semi-quantitative risk assessment for landslides at a regional scale in the Ghizer district in northern Pakistan based on the integration of statistical analysis, field surveys, and spatial (GIS) analysis.The study concentrated on the potential damage to buildings, loss of life for the populations in the buildings, roads, and environmental elements at risk.Due to the lack of a temporal record of landslides, the landslide susceptibility map was integrated with landslide triggers (rainfall and PGA) to compute the landslide hazard map.The results show that the major causative factors for landslides in the study area are roads and stream networks; other parameters, such as the slope, elevation, and faults, also dominantly influence the occurrence of landslides, and the main triggering indicator is rainfall.The extensive database of the mapped building footprints associated with the field-based typology data assists in the realistic vulnerability and risk assessment.The derived vulnerability and risk maps identify and quantify the number of buildings, roads, schools, and hospitals prone to different levels of landslide hazards.The derived landslide risk map shall assist the concerned organizations in developing and implementing landslide mitigation strategies and prioritizing resource allocations and land use planning.

Figure 1 .
Figure 1.Location map of the study area; (a) Pakistan boundary, (b) study area boundary, and (c) topography, locations, roads, and drainage network.

Figure 1 .
Figure 1.Location map of the study area; (a) Pakistan boundary, (b) study area boundary, and (c) topography, locations, roads, and drainage network.

Sustainability 2023 , 23 Figure 2 .
Figure 2. Field photographs of landslides in steep slopes.(a) Catchment area of debris flow in Khalti village; (b) debris slide; (c) scree captured from the toe; (d) rockfall along the main road in Teru village.

Figure 2 .
Figure 2. Field photographs of landslides in steep slopes.(a) Catchment area of debris flow in Khalti village; (b) debris slide; (c) scree captured from the toe; (d) rockfall along the main road in Teru village.

Figure 5 .
Figure 5. Spatial distribution of landslides in the area.

Figure 5 .
Figure 5. Spatial distribution of landslides in the area.

Figure 6 .
Figure 6.Graph showing the influence of classes of causative factors on the landslides using the frequency ratio modal.The blue lines show the values of frequency ratio in each class of causative factors.

Figure 6 .
Figure 6.Graph showing the influence of classes of causative factors on the landslides using the frequency ratio modal.The blue lines show the values of frequency ratio in each class of causative factors.

Figure 7 .
Figure 7. Landslide susceptibility index map and its validation.(a) Susceptibility map showing the zones prone to landslides, developed using frequency ratio model and (b) accuracy of the model by success rate curve, shown on the x-axis.

Figure 7 .
Figure 7. Landslide susceptibility index map and its validation.(a) Susceptibility map showing the zones prone to landslides, developed using frequency ratio model and (b) accuracy of the model by success rate curve, shown on the x-axis.

Figure 8 .
Figure 8. Landslide hazard index map of the study area, showing the landslide hazard zones from very low to very high hazard.

Figure 8 .
Figure 8. Landslide hazard index map of the study area, showing the landslide hazard zones from very low to very high hazard.

Figure 9 .
Figure 9. Physical vulnerability map of the study area.The two insets show the zoom view of some buildings.

Figure 9 .
Figure 9. Physical vulnerability map of the study area.The two insets show the zoom view of some buildings.
are in the very highly vulnerable class, 29,785 (25.4 %) people are in the highly vulnerable class, 34,011 (29 %) people are in the medium class, and the remaining 39,288 (33.5 %) are in very low and low vulnerable classes.

Figure 9 .
Figure 9. Physical vulnerability map of the study area.The two insets show the zoom view of some buildings.
are in the very highly vulnerable class, 29,785 (25.4 %) people are in the highly vulnerable class, 34,011 (29 %) people are in the medium class, and the remaining 39,288 (33.5 %) are in very low and low vulnerable classes.

Figure 10 .
Figure 10.Social and environmental vulnerability maps of the area: (a) social vulnerability and (b) environmental vulnerability.The black rectangles show a zoomed-in view of some villages.Figure 10.Social and environmental vulnerability maps of the area: (a) social vulnerability and (b) environmental vulnerability.The black rectangles show a zoomed-in view of some villages.

Figure 10 .
Figure 10.Social and environmental vulnerability maps of the area: (a) social vulnerability and (b) environmental vulnerability.The black rectangles show a zoomed-in view of some villages.Figure 10.Social and environmental vulnerability maps of the area: (a) social vulnerability and (b) environmental vulnerability.The black rectangles show a zoomed-in view of some villages.

Sustainability 2023 , 23 Figure 12 .
Figure 12.Landslide risk map of the study area.The black rectangles show a zoomed-in view of some villages.

Figure 12 .
Figure 12.Landslide risk map of the study area.The black rectangles show a zoomed-in view of some villages.

Table 1 .
Data and source used in this study.

Table 2 .
Weights assigned to buildings' typology parameters and roads for physical vulnerability.

Table 3 .
Information on the typology data of buildings.