GIS based landslide susceptibility mapping of northern areas of Pakistan, a case study of Shigar and Shyok Basins

ABSTRACT This study presents the use of geographical information system (GIS) datasets and methods to investigate landslide susceptibility in the rugged mountainous terrain of Shigar and Shyok Basins in northern areas of Pakistan. Study area is situated in Karakorum mountainous region where catastrophic landslides occur frequently and pose a serious threat to local living conditions. Landslide susceptibility index maps were prepared by combining four main indicators groups: (1) human induced parameters, which contain Landsat 8 imagery and distance from roads; (2) topographical parameters, which include slope, aspect and plan curvature; (3) hydrological parameters, which contain stream power index, topographic wetness index and distance from drainage and (4) geological parameters, which consist of lithology and distance from major thrusts and faults. These layers were prepared in GIS environment, and analytical hierarchy process (AHP) based heuristic approach was adopted to generate final landslide susceptibility map for this preliminary regional level landslide hazard study. Three different susceptibility zones have been identified in the region. AHP weights with 0.04 consistency ratio were used in the preparation of final map. The final landslide susceptibility map of the study area indicates that the low, moderate, and high landslide susceptibility classes are respectively covering 23.2% (4745 km2), 45% (6737 km2), 23.2% (3444 km2) when distance from roads, distance from drainage and distance from major thrusts and faults parameters were not incorporated in the analysis. Including distance from roads, distance from drainage and distance from major thrusts and faults parameters layers to the final landslide susceptibility analysis shows that moderate susceptibility class is the predominant landslide susceptibility category covering 88% (13,412 km2) of the study area characterized by steep slopes (30°–45°), low positive plan curvature (0–1) and weak lithology.


Introduction
Construction of new roads and highways, expansion of urban built-environment and population growth have resulted in reduction of natural forest covers and raised concerns about environmental safety (Turner & Jayaprakash 1996) because disturbance in the place's environmental setting may allow various factors to contribute towards occurrence of natural hazards like slope failure, soil erosion, floods, etc. Different natural hazards have tendency to combine, superimpose and chain-up

Study area
Study area includes Shigar Valley and Shyok Basins, located in the central Karakoram Range and covers some eastern part of Upper Indus Basin (UIB), as shown in Figure 1. This area exhibiting high topographic relief and climatic variability lies in semi-arid climatic zone receiving 0À50 mm/yr rainfall and experiences glacial erosion adding to the susceptibility of the study area to landslides. Located along the right bank of Indus River falls in the northern areas of Pakistan, the Shigar Valley is a gateway to the high mountains of Karakorum. Its estimated population is more than 70,000 (2009 estimation). It lies between latitude from 35 15 to 35 45 N and longitude from 75 30 to 75 45 E covering about 400 km 2 . Extensive alluvial and debris fans make up the Shigar Valley. While Shyok Basin lies between latitude from 34 39 N to 35 42 N and longitude from 75 56 to 77 27 E, bounded by Jammu and Kashmir Disputed Territory in south, Indus River basin in west and China in north-east (Gilany 2014). The elevation ranges from 2500 m asl to more than 7700 m asl.

Geological setting
A part of study area lies partly on the northern end of the KLIA (Kohistan-Ladakh Island Arc) and partly on southern edge of the Asian plate. The northern suture or the MKT (Main Karakoram Thrust) passing through the centre of the Shigar Valley separates the meta-sediments of the Asian plate from the volcano-clastic rocks of the KLIA. Shroder et al. (2011) mapped the following four landslide complexes throughout the Shigar and Braldu valleys and along Balturo glaciers: the Ghoro Choh rock avalanche; the Buspersackung slope failure; the Gomboro slope failure; the Urdokas rock slide. The Ghoro Choh rock avalanche is adjacent to the place where Shigar Fault strand of the MKT suture and the Karakoram metamorphic complex intersect each other (Maheo et al. 2004). The presence of two major slope-failure complexes in such proximity to each other as well as to where the MKT crosses the upper Shigar Valley from one side to the other, further suggests possible genetic interrelationships. Proximity of fractured slopes to rivers, faults and roads is a major factor initiating slope failures, thereby, causing landslides and destruction (Kamp et al. 2008).

Data used
In landslide susceptibility assessment studies, determination of the causative factors for the landslides is the first step. In fact, the regional landslide assessments should be practical and applicable for that particular area to be investigated; therefore, the input parameters should be representative, reliable and easily accessible (Oh & Pradhan 2011). Factors such as rainfall, drainage and proximity to roads may contribute to generate slope failure. However, not all of these can be acquired, mapped and used cost-effectively in landslide hazard assessment. The basic landslide conditioning factors such as slope angle, slope aspect, lithology, distance from streams, distance from faults and thrusts, topographic wetness index (TWI), stream power index (SPI), major land covers and plan curvature were employed in the landslide susceptibility analysis. These input factors were divided into four main groups. First group consists of human induced factors, i.e. land cover map and distance from roads. Second group includes slope, aspect and plan curvature factors. The third group includes factors representing hydrological conditions in the study area and contains distance from stream network, SPI and TWI parameters. Fourth group consists of lithology and distance from faults and thrust which were prepared from geological maps of the study area.

First group À human induced parameters
This group includes Landsat 8 images and roads dataset. Landsat 8 images were processed to map major land covers in the study area. Landsat satellite remote sensing programme started in 1972 with its first satellite Landsat-1, and Landsat-8 being the latest satellite launched in 2013. The high spatial resolution and presence of variety of bands make Landsat suitable for this study.
Landsat 8 images for year 2013 with less than 10% cloud cover were selected and downloaded from USGS Earth Explorer. For convenience in image handling, bands 2À7 were stacked together to generate a final image containing 6 bands. All images were calibrated into reflectance measures and then visually analyzed for cloud cover. Since the amount of cloud cover was less, clouds were visually detected and manually removed.
Land cover mapping was performed using per-pixel classification approach. For each land cover class, 40À60 sample plots were chosen as training areas. Supervised classification was performed using Erdas Imagine 2011 with 80% overall accuracy. The following categories were identified: water bodies, rocks/barren land, vegetation and glaciers/snow, as shown in Figure 2. The choice of the classes is based on the land cover classification published for mountainous regions.
Theoretically, barren land is more susceptible to landslide than any other land cover. It is because there may be no deep root which could hold the soils firm in place. Contrarily, vegetative areas tend to curb the landslide occurrences by providing anchorage through their roots. It also provides both hydrological and mechanical effects that generally are beneficial for the stability of slopes. Therefore, vegetative areas were considered less prone to landslide due to erosion hazard. Removing vegetation cover for cultivation or any other purpose and even an improper way of farming may play a significant role in initiating landslide.
The presence of road is a triggering factor for landslide mainly due to change in stress and slope stability resulting from the undercutting of the slope, excavation, additional loads, changes in hydrological conditions and drainage (Pourghasemi et al. 2012). When the distance from the roads increases, the tremors caused by the vehicle goes on decreasing, hence, the chances of slope instability hazard resulting from roads proximity also decrease. Six buffer zones with 500 m interval, as shown in Figure 3, were generated to account the influence of distance from the road network on landslide.

Second group À topographical parameters
Advanced Space-borne Thermal Emission and Reflected Radiometer (ASTER) Global Digital Elevation Model version 2 (GDEM2) dataset was used to derive this set of parameters. It is better than older SRTM (90 m, 3 arc second) because of its higher resolution, better tonal definition of surface topography and improvement in accurate measurements in steep mountainous areas due to its radar characteristics (Ahmed et al. 2014) as well as offers improved absolute vertical accuracy compared to its ASTER GDEM version 1 (Tachikawa et al. 2011;Li et al. 2012). ASTER GDEM2 tiles for the study area were mosaicked together using ArcGIS 10, co-registered together with respect to WGS 1984 UTM Zone 43 system and further processed to generate sink free DEM.
Slope map was derived in degree values ranging from 0 to 87.5 . In the reclassified map, slope values were subdivided into five classes, as shown in Figure 4: (1) very gentle slopes, <5 , (2) gentle slopes, 5 À15 , (3) moderately steep slopes, 15 À30 , (4) steep slopes, 30 À45 and (5) escarpments, >45 . Generally, landslides are not expected to occur on gentle slopes due to lower sheer stress. It is expected that landslide susceptibility increases along with the increase in slope steepness reaching the maximum in the fourth class.
Compared to convex slope, soil cover on a concave slope tends to contain more water after precipitation and may retain it for a longer period. Therefore, concave slopes increase the probability of landslides occurrence. On the other hand, at many places, convex slopes mark the outcrop of strong bedrock among looser rocks. For the study area, the plan curvature map was prepared from depression corrected ASTER GDEM2, as shown in Figure 6. Negative values indicate a higher probability of landslide occurrence and vice versa.

Third group À hydrology
This group contains SPI, TWI and distance from stream network parameters.
The TWI is considered as a very important terrain attribute for landslide susceptibility studies. Its values are positive non-zero and increase as the catchment area increases and as the slope angle decreases (Hengl & Reuter 2008). For study area, TWI value ranges from 8 to 28 and 63% area fall in TWI class 8À12, as shown in Figure 7. An increase in TWI increases the potential for landslide (Lee & Min 2001;Gorsevski et al. 2006;Hengl & Reuter 2008).
The SPI measures stream erosion power and contributes toward stability within an area. Higher SPI values generally represent steep, straight, scoured reaches and bedrock gorges while low SPI values correspond to broad alluvial flats, floodplains and slowly subsiding areas, where the valley fill is usually intact and deepening. The SPI map of the study area is shown in Figure 8. Approximately, 26% of study area falls in SPI category from ¡6 to 0 and 72% lie in 0%À5% category.
In the study area, terrain hill slopes are disturbed and modified by rills, gullies and streams through erosion. Slope instability particularly increases in the area where drainage follows the cutting, bank and head-ward erosion (Miller & Sias 1998). Drainage was extracted from depression corrected ASTER GDEM 30 m using HEC-HMS Module in ArcGIS environment. Five buffer zones   with 500 m intervals were generated to account the influence of distance from the drainage network on landslide, as shown in Figure 9. Increase in distance from the streams, extracted in the form of drainage network here, reduces landslide hazard.

Fourth group À geology
This group is comprised of lithology and distance from major thrusts and faults. These parameters affect rock's strengths and soil permeability, thereby, controlling landslide vulnerability. Based on geological map at scale 1:600,000, various rock formations in the study area, as shown in Figure 10, were grouped into different classes to prepare lithology map layer using ArcGIS 10. All lithological units were grouped into six categories and reclassified based on their potential to trigger landslide, given in Table 1.
Details of each of the lithological formation unit are given in Table A1 (see Appendix). The Karakorum Fault and two main thrusts named Main Karakorum Suture (MKS) Zone and MKT were identified being responsible for increased seismicity in the region. Landslide occurrences increase with proximity to geological and tectonic structures. Four buffers around faults and thrusts were generated in meter scale with 4000 m interval.

Methodology
Landslide susceptibility mapping is an important step prior to landslide assessment planning, management and disaster mitigation (Fell et al. 2008). The choice of method for susceptibility mapping not only depends on the size of the study area, but it is also conditioned by the data availability. Due to the limitation of landslide occurrences and soil data availability for the whole study area, a qualitative index based heuristic approach in GIS environment was applied to assess the landslide susceptibility in the region which assumes that the causative factors triggering future landslides will be the same factors that have influenced landslides in the past. Ten thematic layers were selected for spatial mapping of landslide susceptibility in the study area. Analytical hierarchy process (AHP) was used to determine the weight of all parameters on the basis of the previous studies, nature and physical properties of each parameter and subjective expert's opinions (Pachauri et al. 1998; Each parameter was subdivided into number of classes and rating value between 1 and 9, 1 being least susceptible and 9 most susceptible to landslides, is assigned to each class based on their potential to cause slope instability. Again, rating values were assigned to each class based on previous studies and subjective expert's opinion (El Khattabi & Carlier 2004;Calligaris et al. 2013;Ahmed et al. 2014). The AHP method is used in this study to systematically assign preferences based on numerical scale proposed by Saaty, given in Table 2.
All these elements were compared against each other in a pair wise comparison matrix (PCM) which allows for an independent evaluation of each factor contribution and redundancy thereby reduces the measurement error as well as produces a measure of consistency of the judgments comparison, called consistency index (CI), is defined as follows: where λ Max D maximum eigenvalue of the matrix; n D number of participating parameters. Saaty randomly generated reciprocal matrixes using scales 1/9, 1/8,…, 1,…, 8, 9 to evaluate a socalled random consistency index (RI) (Saaty 2000), as given in Table 3.
Consistency ratio (CR) is determined as follows: If the value of the CR is smaller or equal to 0.1, the inconsistency is acceptable, but if the CR is greater than 0.1, the subjective judgment needs to be revised (Saaty 1977). In this way, based on AHP and through normalization of AHP-PCM, different weights were assigned to each input parameter with respect to their anticipated landslide susceptibility level. This is a systematic way to generate weights for heuristic weighted overlay method and avoid inconsistencies in the weights. Table 4 indicates that slope has the highest influence. Then, the order of degree of influence for landsliding is distance from thrusts and faults, lithology, land covers, distance from roads, plan curvature, SPI, TWI, distance from drainage and aspect. This is in accordance with the subjective expert's judgments, given in the literature. The values given in Table 4 show that CR values for each factor and their classes are less than 0.1 and consequently prove that criteria are consistent. Therefore, relative weights given in Table 4 were found appropriate to be subsequently used in the final analysis.

Landslide susceptibility index (LSI)
Once the criterion scores were standardized and weights were computed, final landslide susceptibility index (LSI) map was produced using following equation (Voogd 1983): where W j D weight value for parameter j; w ij D rating value or weight value of class I in parameter j; n D number of parameters.

Results
In the first run, buffered layers were not used in the analysis and the output susceptibility map was divided into three main groups according to their respective level of landslide hazard. Later, a shape file of extent of glaciated areas in the study area was overlapped on the final preliminary susceptibility map to exclude these areas, as shown in Figure 11. Details of each zone are described in Table 5.
The landslide susceptibility mapping was evaluated qualitatively to ensure and improve the prediction accuracy of the landslide susceptibility map.
In the second run, all the 10 layers were incorporated in WLC analysis and final susceptibility map was divided into three susceptibility classes same as in the previous run, as shown in Figure 12. Statistics of each zone is given in Table 6.
Moderate susceptibility zone corresponds to the areas where the slope gradient is elevated and the plan curvature is low positive (0À1), reflecting the significance of these two parameters. A shape file representing extent of glaciated areas was overlapped on the final susceptibility map to exclude the glaciated areas, as shown in Figure 13. Weaker lithological units are also contributing to its vulnerability to slides. While, flat areas, some covered by glaciers, are relatively stable. Figure 11. Landslide susceptibility map of the study area (Run 1) giving somewhat equal weights to all the participating layers. A layer containing point location of major settlements in the study area was overlaid on the LSI map to produce the landslide risk map of the study area, as shown in Figure 13. A large number of settlements fall in the category labelled as moderate susceptibility zone. However, there are few settlements which fall in the susceptibility zone of high landslide potential. These results can be useful as basic data to help slope management, land use planning and avoiding the landslide risk.

Validation of the susceptibility maps
Susceptibility map validation techniques normally compare susceptibility map with the landslide occurrences and assume that future landslides will occur in the same places as existing ones (Kamp et al. 2010). In our study, landslideÀsusceptibility validation was done using two separate sets of input data: previously generated landslide susceptibility map and landslide events data of the study area, previously not employed in susceptibility analysis. A total of 46 landslide events were available for the study area (Hewitt et al. 2011). All of these events occurred in only 51% of the study area only, as shown in Figure 14. This is the same area where most of the population in this region is residing. No information on landslide event occurrence was found in rest of the area.
The ratio between the distribution of the landslides over the susceptibility classes and the distribution of susceptibility class in study area (Lee & Min 2001) was used for evaluating landslide susceptibility maps, as shown in Table 7. Availability of landslides occurrence data for some parts of the study area limited the success of validation. Figure 12. Landslide susceptibility map of the study area (Run 2) when all of the causative factor layers were incorporated using weights obtained from AHP-PCM analysis.

Discussion
The size of the study area brings along some challenges in this study regarding the availability of input data in terms of spatial coverage and map scale. Generally, a full landslide inventory covering the study area is desirable which is by far the most basic requirement for landslide hazard mapping. This may help to establish a temporal analysis of landslide hazard in the region and also to check if any past slides are still fresh and/or reactivated. Unfortunately, most places in the world lack the historical database of landslide inventory including Pakistan. This presented us a big challenge in the selection of suitable method for landslide susceptibility mapping in the area. The choice of a method for regional landslide susceptibility mapping was conditioned by the size of the study area and availability of data at suitable scale. This includes some obvious limits due to the scale of the available geological map (1:600,000), the resolution of the ASTER DEM (30 m £ 30 m), the lack of some fundamental information such as the landslide occurrences for whole study area and top soil cover. Therefore, a qualitative index based heuristic approach is used in this study which is an optimal approach to map the landslide susceptibility on regional scale especially in the absence of landslide inventory data. It assumes that landslide predisposing factors are similar to those observed in past. It allows identification of the areas which have an intrinsic susceptibility to landslides. In this study, more focus has been on landslide triggering due to slope instability as slope is by far the most critical factor. Moreover, presence of active thrust and fault system in this rugged mountainous region makes it more vulnerable. After geology, land cover has been considered as the next very important factor. Closeness to road features also poses certain risk mainly due to undercutting of the slope and may result in slope failure leading to landslides. Precipitation is not taken into account because, here, it is relatively homogenous and low (approximately 0À50 mm/y) (Calligaris et al. 2013;Ahmed et al. 2014). However, other datasets like top soil cover, climate data may also be included to understand the role of these factors in landslides in the region and identify other slide type like rainfall triggered.

Conclusions
In the present study, a qualitative index based heuristic approach was used to map landslide susceptibility in the study area. The empirical results showed that the slope, geology and land cover were the most important and highly weighted factors. The investigated area was divided into three classes of susceptibility, low, medium and high, which showed that about 88% of the area has moderate and 3% has high susceptibility for landslides. The quality of these results largely depends on the quality of the input data, but the approach used is so simple that results can easily be updated on demand and it provides relatively quick results. Validation of the landslide susceptibility results using the available 46 landslide events data showed that these events occurred mostly in moderate susceptibility zone which was the same area where most of the population in this region was living. No information on landslide event occurrence was found in rest of the area.
Approach used in this study provides relatively quick results and is so simple that results can easily be updated on demand. These maps can be considered as base map for landslide hazard evaluation when intending to avoid or reduce future hazard's impacts and improve decision-making prior to any development in this region requiring detailed investigation. Results are useful for explaining the conditioning factors for triggering landslides, thereby, supporting the efforts to mitigate future