Most of the northern part of Pakistan is located in snow-covered mountains. The steep relief, snow, and glaciers, in the region, are exceptional but strong precipitation and a high seismicity contributes to the origin of widespread natural processes like debris flow, flash floods, earthquakes, rockfall or landslides  . In Karakoram Mountains eight various types of mass movements have been observed, rock falls, avalanche, rockslides, debris flow, flow slides, rotational slip, slumps and creep  . Among these debris flow and flow slides are the most prevailing and frequent type of the mass movement noted in the Karakoram. Debris flow is an abrupt mass movement which may cause by intense rainfall on unconsolidated steep hills  . According to  Hindu Kush Karakoram Himalaya (HKKH) region is facing increased flash flood and related hazards. The HKKH Mountains are especially prone to hydrogeological disasters, such as flash floods, landslides, and Glacier Lake Outburst Floods (GLOFs). Gilgit Baltistan (GB) is comprised of a rugged mountainous topography where mountains comprise 90% of the total area and susceptible to landslides, lakes formation and GLOFs  .
The Landslide is a consequence of multifarious interaction within various factors, such as meteorological, geomorphological and geological. The spatial information associated with mentioned factors is able to extract from remote sensing facts, land-based information, along with quite a lot of other data resources. Landslide susceptibility maps illustrate the comparative possibility of future landslide based exclusively on the fundamental properties of a setting or site. The AHP helps in generating the weight of every factor, which is then used in Geographic Information System (GIS) to create landslide susceptibility maps.
Even though Gilgit Baltistan is highly susceptible to landslides, merely a few studies have been carried out in the region  . Moreover no such studied has been conducted in district Ghizer. Besides this, landslide hazard maps or official landslide inventory maps are still lacking in Gilgit Baltistan This study aims to develop landslide susceptibility maps for district Ghizer by using freely available data through geospatial approaches. Hence, the study includes the assessment of main landslide causative factors in the district Ghizer and the temporal assessment of land cover change for past 16 years and its impacts on landslides in the study area.
1.1. Study Area
The district Ghizer lies in Hindu Kush region of Pakistan in the northern part of Gilgit-Baltistan, between latitude 36.0˚N to 37.0˚N and longitude 73.0˚E to 74.0˚E covering an area of 12,042 km2. The estimated terrain elevation above sea level is 3661 meters and the habitat is arid to semi-arid. The region has four Tehsils i.e. Gupis, Ishkoman, Punial and Yasin (Figure 1). The area is prone to different natural disasters although snow avalanches, landslides, and earthquakes
Figure 1. Location map of study area in Pakistan, created by first author (IR).
are frequent in the area  . Hindukush mountain region is considered as seismically active zone because of the occurrence of low-intensity earthquakes at frequent intervals. Various fault lines spread through entire Gilgit Baltistan region. Low to medium intensity earthquakes are frequent mainly responsible for the occurrence of landslides, Glacial Lake Outburst Flood (GLOF), avalanche, rock fall, and edge failure in the study area. The valleys are present in steeps hills and upper parts of the region get isolated for several months in the winters due to heavy snowfall, landslides and snow avalanches.
The road from district Ghizer to district Gilgit is highly susceptible to landslides because of erosion and rock fall as the slopes are made of muddy dust and loose sediments. There exists sixteen major landslide areas and thirteen small villages prone to rock fall, so as a whole 53 villages are considered to be at high risk to hazards in district Ghizer  .
1.2. Data Acquisition and Preparation of Causative Factors
In this study, twelve factors were used to generate the landslide susceptibility map. The factors were totally selected on the basis of their effectiveness and availability. According to Oh and Pradhan  , the assessment of the local landslide areas should be convenient and relevant and the factors should be illustrative and effectively available. The factors; slope, aspect, elevation, drainage network, SPI, and TWI were extracted from Digital Elevation Model (DEM) of 30 m resolution acquired from USGS Earth Explorer.
The slope angle ranges from 0˚ - 73.76˚ for the study area (Figure 2(a)). Slopes were reclassified into five classes i.e. <5˚, 5˚ - 15˚, 15˚ - 30˚, 30˚ - 45˚ and >45˚. According to  , highest landslide susceptible class is 30˚ - 40˚ slope angle which consists of steep slopes. The aspect map (Figure 2(b)) was divided in to nine classes based on dimensions flat (−1)˚, north (337.5˚ - 360˚, 0˚ - 22.5˚), north-east (22.5˚ - 67.5˚), east (67.5˚ - 112.5˚), south-east (112.5˚ - 157.5˚), south (157.5˚ - 202.5˚), south-west (202.5˚ - 247.5˚), west (247.5˚ - 292.5˚), and north-west (292.5˚ - 337.5˚). According to  the southwest and northwest facing slopes are highly susceptible to landslides.  suggest assigning higher
Figure 2. Landslide causative factor maps of study area: (a) Slope, (b) Aspect, (c) Elevation, (d) Drainage, (e) SPI, (f) TWI.
weights to southwest, west, and northwest facing slopes. The lowest point of district Ghizer is at elevation of 1662 m and the highest point of elevation is at 6789 m (Figure 2(c)). It was categorized into five classes 1162 - 2894 m, 2894 - 3668 m, 3668 - 4265 m, 4265 - 4829 m, 4829 - 6789 m. Landslide occurrence is linked to certain elevations  .  investigated that 64% of reported landslides were recorded at elevation of 2000 - 4000 m and 24% were observed at elevation of 1000 - 2000 m. Buffers were created around the drainage network and classified into 0 - 500 m, 500 - 1500 m, 1500 - 2500 m, 2500 - 5000 m and <5100 m (Figure 2(d)). Landslides increase if the distance to streams or rivers is decreased, due to slope instability which leads to erosion. Irregularities and fragmentations are caused by a river’s longitudinal profile due to slope failures  . The secondary attribute SPI (Figure 2(e)) tells the net erosion and net deposition in the areas of increased flow rate and decreased flow rate  . Furthermore, landslide susceptibility is higher with the higher SPI values. It was reclassified into four class values −13 - 0, 0 - 5, 5 - 10 and 10 - 14. The other secondary attribute TWI (Figure 2(f)) calculates the extent of water accumulation at a place; its higher values show higher accumulation causing more landslide susceptibility  . It was reclassified into three classes i.e. 3 - 9, 9 - 28 and <28.
From the geological map of scale 1:5,000,000 fourteen rock types were digitized, which was acquired from Geological Survey of Pakistan (Figure 3(a)). All the lithological units were categorized in view of their capability to trigger landslide Table 1. Each lithological unit has its own susceptibility towards landslides, so it needs to categorize lithological units accordingly  . Ranking of the rock types was based on their potential to cause landslide  . The rock formations in the study area are categorized into five classes from A to E; each
Table 1. Lithological units in the district Ghizer.
Figure 3. Landslide causative factor maps of study area: (a) Lithology, (b) Fault Lines, (c) Rainfall, (d) Roads, (e) Land Cover, (f) Soil.
rock formation is categorized into very highly stable, highly stable, moderately stable, less stable and least stable class; based on the geotechnical properties of the rock units present in the formation. Less stable rock units are highly prone to slope failures which cause landslides. Geological map of scale 1:5,000,000 was used to digitize the fault lines of the study area which were acquired from Geological Survey of Pakistan (Figure 3(b)). Buffers were created for distances of 0 - 3000 m, 3000 - 7000 m, 7000 - 11,000 m, 11,000 - 15,000 m and <16,000 m. Shearing causes rocks weak which are close to fault lines, consequently leading to landslide susceptibility  . Monthly average rainfall for different locations of district Ghizer was acquired from Tropical Rainfall Measuring Mission (TRMM) for the years 2006-2015 (Figure 3(c)). Rainfall is an important landslide triggering factor, but it is limited to the monsoon season  . Rainfall raster data map was prepared using the Inverse Distance Weighted (IDW) interpolation. The infrastructure of the study area is poor and no such complex road network exists (Figure 3(d)). Road network data was acquired from an online source http://www.mapcruzin.com. Buffers were created for roads in the study area at distance of 0 - 500 m, 500 - 1500 m, 1500 - 2500 m, 2500 - 5000 m and <5100 m. Cutting of slopes for road construction or road widening in the hilly regions could lead to slope failures causing landslide susceptibility  . Landsat 8 OLI image for the year 2015 and L4-5TM image for 1999 were acquired from USGS Earth Explorer. Land cover maps were prepared using maximum likelihood supervised classification techniques in ERDAS Imagine 14 (Figure 3(e)). The classes prepared were glacier, vegetation barren soil/ exposed rocks and water. Land cover images for 1999 and 2015 were compared to the land cover change and its impact on a landslide. Accuracy assessment of the classified images (1999 and 2015) was calculated to check the classification accuracy. The accuracy assessment was generated using 50 random points.
Soil sampling for soil texture was performed by taking twelve composite samples from each tehsil of the study area. The samples were air-dried and sieved through 2 mm size sieve. Forty ml of 1% sodium hexa meta-phosphate and 150 ml of distilled water was added to soil sample (40 g) and was kept overnight. The mixture was stirred for almost 10 minutes and was put in a graduated cylinder for readings, which was recorded with Boyoucos Hydrometer method  . IDW interpolation method was used to create the raster map of soil texture (Figure 3(f)).
2. Landslide Susceptibility Mapping
2.1. Analytical Hierarchy Process
The Analytical Hierarchy Process (AHP) is an adaptable tool which is created by  and it is used for various decisions makings such as suitability analysis and susceptibility analysis. It is considered to be a rational decision-making process for multi-criteria as well as for multi-target approach.
In the comparison matrix, the numerical value for each factor was between 1 and 9 (Table 2). The factors were organized hierarchically in the matrix and the Prioritized Factor Rating Value (PFRV) technique was used to assign a numerical value to the factors in the AHP on the basis of their importance as compared with other factors. The numerical value assigned to the factors was based on, expert knowledge, literature, observations, and experiences.
Table 2. Saaty’s proposed numerical scale.
Source: Saaty 1977.
The average of the hierarchically arranged factors was used to calculate the weights and rating value/eigenvalue along with the Consistency Ratio (CR), based on the prepositions of  .  expressed that the eigenvalue “λmax” and the total number of factors “n” are same for a consistent comparison matrix.
CI = Consistency Index which is as follow:
The consistency of the comparison matrix is checked through CR (Saaty 1977).
where RI = Random Consistency Index.
 have created RI by utilizing scales 1/9, 1/8, 1/7… 1… 8, 9. The average RI of 12 matrixes is given in (Table 3).
The calculated CR from the comparison matrix for the 12 factors was 0.028. This value demonstrates that the matrix of the factors is acceptable. The result of AHP showing weights of causative factors (Wj) and the factor rating values (wij) are given in the (Table 4).
Table 3. Random consistency index.
Source: Saaty 1977.
Table 4. Pair wise comparison matrix, factor weights and consistency ration of the data layers.
2.2. Weighted Linear Combination
Weighted Linear Combination (WLC) is comprised of both subjective and quantitative strategies and depends on the qualitative map combination approach (heuristic analysis)  . It is the last step in making the landslide susceptibility map in which all the weighted layers were combined using weighted overlay technique in ArcGIS 10.1. All the layers were reclassified to a typical scale and the vector layers were rasterized. The weights of the factors were linearly combined (WLC) to obtain the Landslide susceptible Index (LSI) according to the formula:
where LSI is Landslide susceptibility index, Wj is weight value for parameter j, wij is rating value or weight value of class I in parameter j and n is no. of classes.
3.1. Landslide Susceptibility Mapping
The weights of the factors; slope, aspect, elevation, drainage network, SPI, TWI, lithology, fault lines, rainfall, roads, land cover land use and soil were derived using AHP by Prioritized Factor Rating Value (PFRV) (Table 3). The final landslide susceptibility map was generated using these weights in the WLC. The resultant map showed that the pixel ranking value for landslide susceptibility varies from very low (1.53) to very high (4.43). The areas with high pixel values have more chance of landslide as compare to the low pixel values. The categorization of the pixel ranking values was obtained by natural breaks in GIS.
Based on the above categorization, the area and percentage of the five susceptibility classes were also determined (Figure 4). Very low susceptibility class covers an area of 8.66% while; low susceptibility class covers 16.96% of the area. In addition, a larger extent of the area lays in the moderate category i.e. 28.14%. Furthermore, the high susceptibility class is the one which covers a larger area in the district Ghizer i.e. 28.22%. The very high susceptibility class in the district falls over an area of 18.02%. Hence, in district Ghizer, a total of 74.38% of the surface area falls into the moderate to very high landslide susceptible zones whereas 25.62% of the area falls into low to very low landslide susceptible zones.
3.2. Susceptibility in Reaction to Land Cover Change in District Ghizer
The topographic, geologic, and hydrologic factors causing landslides are considered as stationary, while land cover is the factor that can change within a short time; therefore it is in a direct relation to landslide occurrence  . In this regard, Temporal assessment of land cover change was studied for the years 1999 and 2015, to analyze the difference in the land cover change over sixteen years in the district Ghizer and its impact on landslides. Hence, between the years 1999 till 2015 a number of landslide events have occurred and significant changes in
Figure 4. Landslide susceptibility map of district Ghizer derived through WLC model.
land cover have been observed. The changes are visible in the classified maps (Figure 5). It showed a major decline in the glaciers from 9.79% to 6.63%. As the district Ghizer is largely covered by barren soil/exposed rocks, it poses more vulnerability to landslides. Barren slopes have more chances of erosion as compared to areas with vegetation so they are more susceptible to landsliding  . However, the barren soil/ exposed rocks have reduced to 78.46% from 81.465. Vegetation cover has increased from 8.33% to 12.09%, while water class which was least area covering class in 1999, increased from 0.41% to 2.82%. The result of overall classification accuracies for the year 1999 and 2015 from the accuracy assessment were 80.0% and 80.01% respectively. In most of the studies overall classification accuracies target below of 85%  .
3.3. Validation of Susceptibility Map
There is number of methods to validate a susceptibility map. One such method is computing landslide frequency/density in the susceptibility classes  . In this study, landslide susceptibility map validation is prepared by computing landslide frequency in the susceptibility classes. For this, 34 observed landslide sites were considered (Figure 6).
The observed landslides in the very high susceptible zone were 38.2% with a landslide frequency of 0.0059, which was found to be the largest among other
Figure 5. Classified land cover land use change detection maps of district Ghizer.
Figure 6. Observed landslides in the study area are overlaid on landslide susceptibility map generated trough WLC model.
susceptibility classes. The high, moderate, low and very low classes showed frequencies of 0.0035, 0.001477, 0.001471 and 0.00096 respectively. The overall validation result shows that 88.1% of the landslides have occurred in the moderate to very high susceptibility zones (Figure 7).
Figure 7. Observed landslide frequencies in landslide susceptibility classes.
The weight values of each factor in AHP shows the level of in the landslide. Results showed that slope, distance from fault lines and lithology of the study area have the greatest impact on landslide hazard. It is evident from the results that most of the landslides occur in the gentle to moderate slopes. It has been observed that 20˚ to 40˚ slope angles are considered very susceptible to landslides  . From the literature, it was determined that slope angle was given highest value  . For this reason, the slope has been considered as an important factor in this study as well.  expresses that, according to the documented land and rock slides 44% of the slope instabilities are documented in the slope angles of 30˚ and 45˚. Hence gentle to moderate slopes are more susceptible to landslides. Moreover, the mountainous areas are more vulnerable to landslides with the presence of active fault lines. Main Karakorum Thrust and Trich Mir fault run across the district Ghizer. The two categories; high landslide susceptibility (28.22%) and very high landslide susceptibility (18.02%) are mostly present in the region where the slope is steep and the distance to fault lines is less. Thus, this shows that the slope angle and the fault lines are most important factors in landslide susceptibility.
Moreover, the finding demonstrated that the weaker rocks which are loosely held are more prone to falling. It is widely recognized that geology of an area, greatly influences the occurrence of landslides and rock falls in that particular area. Because every rock type has different composition and that leads to a difference in permeability  . The lithology of an area consists of different formations which are represented by the characteristics of rock type, which can cause landslides. The Kohistan Batholith Formation (KB) and Southern Karakoram Metamorphic Complex (Skm) were observed in high susceptibility classes, while low susceptibility classes were observed in rocks belonging to Eclogites (Ec), Shyok Suture Zone (Sv) and Hunza Plutonic Unit (HPU) Formations. Rocks belonging to KB and Skm Formation are highly deformed and lie in the most to medium sediment productivity class and inherently failure prone. Rainfall is taken into account in this respect, but it is almost same in all the parts of the study area and it receives 0 - 150 mm rainfall per year  . Therefore it is given a low weight. Rainfall is an important landslide triggering factor, but it is only limited to the monsoon season in the study area when the duration and intensity of rainfall are high. The drainage networks impact the weight of the soil only if a storm or substantial rain came. The streams can erode the slopes and cause landslide. In the study area, the drainage network only impacts the slopes during monsoon season  . The two factors soil and distance to drainage are associated with the rainfall in the study area, therefore, these are given a less value in the AHP. Aspect, TWI, and SPI are included in the study, but these are given less value according to literature.
Land cover has been considered an important factor in the study because barren slopes are widespread as the vegetation is mainly around the villages and few rangelands are present in the high mountains  . The landslide susceptibility map reveals that the areas covering vegetation were mostly observed in low landslide susceptibility zones. The land cover trend analysis of district Ghizer from the year 1999 to 2015 shows that glaciers are melting at a high pace and have reduced from 9.79% to 6.63%. The reason for this meltdown is global warming as the glaciers throughout the Himalayas are decreasing  . The debris material in these mountains is loosely held and is prone to flow or slide, which can cause flash floods, GLOFs, snow avalanches, and debris flows. The classified image of 2015 also shows a number of lakes and small water bodies exist near the areas where the glacier was present previously. And the water statistics shows that water has increased from 0.41% to 2.82%. Retreating glacier can frequently form glacial lakes near the glaciers  .
In the presented study, GIS techniques and AHP were applied to create landslide susceptibility map. Based on the achieved results, a large area in the district consists of moderate and high landslides prone zones. The produced susceptibility map was compared with randomly selected landslides for validation; landslide frequency/density was computed from observed landslides in the study area, which also indicated that highest frequency of landslides is in the very high susceptibility zone. Besides producing the landslide susceptibility map for the study area, temporal assessment of land cover change in the district Ghizer for the years 1999 and 2015 was investigated to study the impact of land cover on landslide susceptibility. Based on the results it can be stated that vegetation and water class has increased within the sixteen-year time span while the glaciers and barren soil/exposed rock classes have reduced. This approach can be applied to the landslide susceptibility mapping in other regions in the world. However, it is important to assign appropriate weights to the specific landslide-controlling factors, because it is mostly attributable to the nature of the terrain and type of landslide.
 Owen, L.A. (1996) Quaternary Lacustrine Deposits in a High Energy Semi-Arid Mountain Environment, Karakoram Mountains, Northern Pakistan. Journal of Quaternary Science, 11, 461-483.
 Chevalier, G.G., Medina, V., Hurlimann, M. and Bateman, A. (2013) Debris-Flow Susceptibility Analysis Using Fluvio Morphological Parameters and Data Mining: Application to the Central-Eastern Pyrenees. Natural Hazards, 67, 213-238.
 Calligaris, C., Poretti, G., Tariq, S. and Melis, M.T. (2013) First Steps towards a Landslide Inventory Map of the Central Karakoram National Park. European Journal of Remote Sensing, 46, 272-287.
 Kanwal, S., Atif, S. and Shafiq, M. (2016) GIS Based Landslide Susceptibility Mapping of Northern Areas of Pakistan, a Case Study of Shigar and Shyok Basins. Geomatics, Natural Hazards and Risk, 8, 348-366.
 Oh, H.J. and Pradhan, B. (2011) Application of a Neuro-Fuzzy Model to Landslide-Susceptibility Mapping for Shallow Landslides in a Tropical Hilly Area. Computers & Geosciences, 37, 1264-1276.
 Ahmed, M.F., Rogers, J.D. and Ismail, E.H. (2014) A Regional Level Preliminary Landslide Susceptibility Study of the Upper Indus River Basin. European Journal of Remote Sensing, 47, 343-373.
 Ahmed, M.F. and Rogers, J.D. (2014) Creating Reliable, First-Approximation Landslide Inventory Maps Using ASTER DEM Data and Geomorphic Indicators, an Example from the Upper Indus River in Northern Pakistan. Journal of Environmental & Engineering Geoscience, 20, 67-83.
 Pouydal, C.P., Chang, C., Oh, H.J. and Lee, S. (2010) Landslide Susceptibility Maps Comparing Frequency Ratio and Artificial Neural Networks: A Case Study from the Nepal Himalaya. Environmental Earth Sciences, 61, 1049-1064.
 Duman, T.Y., Can, T., Gokceoglu, C., Nefeslioglu, H.A. and Sonmez, H. (2006) Application of Logistic Regression for Landslide Susceptibility Zoning of Cekmece Area, Istanbul, Turkey. Environmental Geology, 51, 241-256.
 Leir, M., Michell, A. and Ramsay, S. (2004) Regional Landslide Hazard Susceptibility Mapping for Pipelines in British Columbia. Geo-Engineering for the Society and Its Environment. 57th Canadian Geotechnical Conference and the 5th Joint CGS-IAH Conference, Old Quebec, 24-27 October 2004, 1-9.
 Yalcin, A. (2008) GIS-Based Landslide Susceptibility Mapping using Analytical Hierarchy Process and Bivariate Statistics in Ardesen (Turkey): Comparisons of Results and Confirmations. Catena, 72, 1-12.
 Ayalew, L., Yamagishi, H. and Ugawa, N. (2004) Landslide Susceptibility Mapping Using GIS-Based Weighted Linear Combination, the Case in Tsugawa Area of Agano River, Niigata Prefecture, Japan. Landslides, 1, 73-81.
 Malek, Z., Zumpano, V., Schroter, D., Glade, T., Balteanu, D. and Micu, M. (2015) Scenarios of Land Cover Change and Landslide Susceptibility: An Example from the Buzau Subcarpathians, Romania. Engineering Geology for Society and Territory, 5, 743-746.
 Sarkar, S. and Kanungo, D.P. (2004) An Integrated Approach for Landslide Susceptibility Mapping Using Remote Sensing and GIS. Photogrammetric Engineering & Remote Sensing, 70, 617-625.
 DeGloria, S.D., Laba, M., Gregory, S.K., Braden, J., Ogurcak, D. and Hill, E. (2000) Conventional and Fuzzy Accuracy Assessment of Land Cover Maps at Regional Scale. Proceedings of the 4th International Symposium on Spatial Accuracy Assessment in Natural Resources and Environmental Sciences, Amsterdam, July 2000, 153-160.
 Kumar, R. and Anbalagan, R.J. (2016) Landslide Susceptibility Mapping using Analytical Hierarchy Process (AHP) in Tehri Reservoir Rim Region, Uttarakhand. Journal of the Geological Society of India, 87, 271-286.
 Ruff, M. and Czurda, K. (2008) Landslide Susceptibility Analysis with a Heuristic Approach in the Eastern Alps (Vorarlberg, Austria). Geomorphology, 94, 314-324.
 Kayastha, P., Dhital, M.R. and De Smedt, F. (2013) Application of the Analytical Hierarchy Process (AHP) for Landslide Susceptibility Mapping: A Case Study from the Tinau Watershed, West Nepal. Computers & Geosciences, 52, 398-408.
 Pradhan, B. and Lee, S. (2010) Delineation of Landslide Hazard Areas on Penang Island, Malaysia, by Using Frequency Ratio, Logistic Regression, and Artificial Neural Network Models. Environmental Earth Science, 60, 1037-1054.
 Roohi, R., Ashraf, R., Mustafa, N. and Mustafa, T. (2008) Preparatory Assessment Report on Community Based Survey for Assessment of Glacial Lake Outburst Flood Hazards (GLOFs) in Hunza River Basin. Water Resources Research Institute, National Agricultural Research Centre, Islamabad and UNDP Pakistan.