LANDSLIDE SUSCEPTIBILITY MAPPING BY USING GIS ALONG THE CHINA PAKISTAN ECONOMIC CORRIDOR ( KARAKORAM HIGHWAY ) , PAKISTAN

LANDSLIDE SUSCEPTIBILITY MAPPING BY USING GIS ALONG THE CHINA PAKISTAN ECONOMIC CORRIDOR (KARAKORAM HIGHWAY), PAKISTAN. Sajid Ali, Peter Biermanns, Rashid Haider and Klaus Reicherter 1Neotectonics and Natural Hazards, RWTH Aachen University, Lochnerstr. 4-20, 52056 Aachen, Germany. 5 2Department of Earth Sciences, COMSATS Information Technology, Abbottabad, Pakistan. 3Geological Survey of Pakistan, Islamabad, Pakistan.


Introduction
The Karakoram Highway (KKH), a part of the China-Pakistan Economic Corridor (CPEC), connects northern Pakistan with western China (Fig. 1).It passes through rapidly rising mountain ranges of Himalaya, Karakoram, Hindu Kush forming the junction between the Indian and Eurasian plates including the Kohistan Island Arc (Derbyshire et al., 2001).The area is characterized by fractured and weathered rockmass, diverse lithologies (igneous, metamorphic, and sedimentary), high seismicity, deep gorges, high relief, arid to Monsoon climate and locally high rates of tectonic activity.These conditions make the study area a unique geohazards laboratory.Starting with its construction in 1979, KKH's stability has been endangered by a variety of geohazards.Landslides constitute an appreciable threat, having blocked KKH for several times.
Landslides are a result of different geodynamic processes and represent a potentially momentous type of geohazard, causing economic and social loss by damaging infrastructure and buildings (Vallejo and Ferrer, 2011).Landslides are mainly caused by conditioning and triggering factors.Conditioning factors include relief, lithology, geological structure, geomechanical properties and weathering, whereas precipitation, seismicity, change in temperature and static or dynamic loads are triggering factors.Variations in these factors affect the occurrence of landslides.Assessment of risk related to landslides remained a longtime challenge for geologists but was significantly facilitated with the eventual availability of remote sensing data (Shahabi and Hashim, 2015).Landslide susceptibility mapping is the spatial prediction of landslide occurrence by considering causes of previous events (Guzzetti et al., 1999).It largely depends upon the knowledge of slope movement and controlling factors (Yalcin, 2008).It has hitherto been carried out by many researchers in order to denominate potential landslide hazard zones through evaluation of responsible factors (Basharat et al., 2016;Komac, 2006;Lee et al., 2002Lee et al., , 2004;;Shahabi and Hashim, 2015;Süzen and Doyuran, 2004).Preparation of their maps was based broadly on qualitative and quantitative approaches.
Early research work (Nilsen et al., 1979) was largely quantitative, utilizing deterministic and statistical correlations and regression analysis of landslides and their controlling factors.In this context, safety factors, calculated on the base of engineering parameters are adduced to imply deterministic methods.In more recent works, statistical methods are favoured, attempting correlation between spatial distribution of landslides and their controlling factors.Among these are e.g.analytical hierarchy process, bivariate, multi variate, logistic regression neural network, fuzzy logic etc. (Basharat et al., 2016;Guzzetti et al., 1999;Komac, 2006;Lee et al., 2002Lee et al., , 2004;;Shahabi and Hashim, 2015;Süzen and Doyuran, 2004).These techniques were proven to be better options for comparatively large and complex areas (Cardinali et al., 2000).Expert opinion and landslide inventories are the decisive components of qualitative approaches (Yalcin, 2008).In most cases, landslide inventories were adduced to estimate failure susceptibility based on previous hazards in locations with similar geological, geomorphological and hydrological setups.Some geoscientists incorporated statistical techniques (analytical hierarchy approach and weighted linear combination) in qualitative methods to provide the identified factors with a numerical weightage.
Analytical hierarchy process (AHP) is a simple and flexible method to analyze and solve complex problems (Saaty, 1987(Saaty, , 1990)).It facilitates the estimation of influence that different factors might have on landslide development by comparing them in possible pairs in a matrix.As it entirely depends upon the knowledge of experts, thus different interpretations may change the results.The combination of AHP with qualitative approaches resulted in semi-qualitative or semi-quantitative approaches that are considered a better option for regional studies (Soeters and van Westen, 1996).Furthermore, the use of Geographic Information Systems (GIS) facilitated the extraction of geomorphic and hydrological parameters required for susceptibility assessment.Digital Elevation Models (DEM) are commonly processed in GIS to extract crucial parameters for susceptibility assessment such as elevation, slope, aspect, curvature, watershed etc. (Ayalew et al., 2004(Ayalew et al., , 2005;;Basharat et al., 2016;Ohlmacher and Davis, 2003;Rozos et al., 2011;Shahabi and Hashim, 2015;Süzen and Doyuran, 2004).

General situation of the study area
The study area is the 840 km long Highway (KKH), N35, located in the Karakoram Mountains, Himalaya.The area hosts some of the highest reliefs and highest peaks (Nanga Parbat: 8126 m, Rakaposhi: 7788 m) in the world (Hewitt, 1998).Goudie et al. (1984) termed the study area the steepest place on the earth where elevation drops from 7788 m to 2000 m over a horizontal distance of 11 km (Fig. 2b).From Abbotabad, the Highway leads northwards through the Sub-Himalayas entering the Indus valley at Thakot, and the Hunza valley at Gilgit, running parallel to the eponymous rivers.From Thakot onwards, it passes through deeply incised valleys and gorges.
Weather conditions along KKH are not uniform and are characterized by a wide range of annual mean temperatures (-5 °C to 46 °C) and precipitation (15 mm to 1500 mm).The distribution of precipitation is additionally strongly fluctuating throughout the year.During the westerlies (January, February and March) and the monsoon period (July, August), the study area receives heavy rainfall.According to meteorological data, average annual precipitation between Abbottabad and Dassu is 1444 mm.
However, north of Dassu, an abrupt change from monsoon to semi-arid to arid conditions is recorded which is owed to a change in valley orientation from north-south to east-west (Fig. 1 and Fig. 3).Furthermore, vertical climatic zonation exists in the Hunza valley along KKH.The surrounding peaks and slopes higher than 5000 m receive precipitation greater than 1000 mm, whereas the valley floor below is characterized by a semi-arid to arid climate (Hewitt, 1998).
The Karakoram and Himalaya host some of the world's longest continental glaciers with the steepest gradient and highest glacial erosion rates (Goudie et al., 1984).The snouts of some glaciers (Batura, Ghulkin, Pasu, Gulmit, Gulkin) are close to KKH and partly cross it (Fig. 2a).Relatively warm temperature in the valleys results in sudden melting of ice, frequent surges, catastrophic debris flows and blockage of rivers.

Geology along the KKH
Tectonically, the area is characterized by orogenic features that started forming with the onset of the Indo-Eurasian collision 50 Myr ago.Crustal shortening, subduction, and active faulting are still ongoing with convergence rates of ~4-5 cm/y (Jade, 2004) and uplift rates of ~7mm/y (Zeitler, 1985).Main Mantle Thrust (MMT), Kamila Jal Shear zone (KJSZ), Raikot Fault, Main Karakoram Thrust (MKT) and Karakorum Fault are important tectonic features responsible for brittle deformation along the highway (Fig. 4) (Bishop et al., 2002;Burg et al., 2006;DiPietro et al., 2000;Goudie et al., 1984).Due to this brittle deformation, the rockmass is highly jointed and fractured.The general geology along the KKH consists of sedimentary, igneous and metamorphic rocks.Highly active landslide zones were identified from the distribution of existing landslides (Fig. The geology of the Jijal-Dassu section is composed of ultramafic and low to high grade metamorphic rocks.The Mansehra granite, the Besham group, the Jijal complex, the Kamila amphibolite and the Chilas complex are important lithological units in this section (Fig. 5a).The Besham group comprises biotitic gneisses, cataclastic gneisses and quartzite which were metamorphosed during the Himalayan orogeny ∼65 Ma ago (Ding et al., 2016).It shares a faulted contact with the Jijal complex (Kohistan Island arc) along the northward dipping MMT (Williams, 1989).The ultramafic rocks with garnet granulites and Alpine type metamorphic rocks between Jijal and Pattan are collectively termed as Jijal complex (Tahirkheli and Jan, 1979).
The Besham group shows signs of crushing of individual minerals and staining of quartz whereas the Jijal complex is massive and sheared (Khan et al., 1987).Owed to the contacts of the Jijal Complex with the MMT in the north and the Pattan Fault in the south, it is highly tectonised and deformed.The Kamila amphibolite consists of sheared basic lavas and intrusive plutons (Treloar et al., 1996).It is classified into two types: garnet bearing and garnet free amphibolites.The former is massive and sheared due to the presence of Pattan and Kamila Jal shear zones (KJS), whereas the latter is banded.The garnet free amphibolite shares a sheared contact with the Chilas complex, mafic intrusions of predominantly gabbro-norites, sheared gabbro-norites and diorites (Searle et al., 1999).
The Raikot Bridge section (Fig. 5b) is a highly active landslide zone being located at the seismically active western limb of the MMT known as Raikot fault, a strike-slip fault with right lateral movement.It is marked by concentration of hot springs and a large shear zone.Granitic gneisses, quartzites, gabronorite, schists and Quaternary sediments are the main lithological units in this section.Continuous erosion by Indus River and highly deformed rocks are responsible for landslide events.
MKT is a part of the Hunza valley section and is responsible for many landslide events.The pre-dominant local lithologies of the Baltit Group, Chalt schists, Karakoram batholith as well as Quaternary sediments are highly tectonised and deformed (Fig. 5c).
The highly deformed Misgar slates, along with the Gujhal dolomite and Kilk formation, are the main components of the Sost section (Fig. 5d).The Karakoram fault is an important tectonic feature in this section.Highly fissile and closely jointed slates are important sources of scree on steep slopes along the Highway.Intense weather conditions aggravate the situation in this section.

Seismology
The highway passes through one of the seismically most active areas in the world.The presence of active thrusts and strike slip faults gives rise to earthquakes, anon triggering numerous landslides.The seismic activity along KKH is demonstrated by 317 M>5 and 10 M>7 recorded earthquake events (Muzaffarabad Oct, 2005: M=7.6, Afghanistan Oct, 2015: M=7.5) since 1904 (Zhiquan and Yingyan, 2016).The highway passes through important seismic zones: the Indus Kohistan seismic zone Sazin section of the highway is a part of HSZ, an active seismic zone hosting recent events with magnitudes of 3 to 6.2.The active Raikot fault traverses RSSZ and is responsible for shallow seismicity of magnitudes 3 to 6.3.Both the fault and the KKH run, in direct vicinity, on the western banks of Indus River.The YSZ encompasses the region surrounding the small town of Yasin.It is characterized by earthquakes with magnitudes between 3 and 5 and focal depths of less than 50 km.The MKT is suspected to be the main source of seismicity for this seismically active region.

Causative Factors and Spatial Distribution analysis
Geological, morphological, seismo-tectonic, topographic and climatic factors are generally considered as landslide-controlling parameters (Kamp et al., 2008).The following ten causative factors were considered for preparation of the map: Lithology, distance from faults, seismicity, elevation, slope, aspect, curvature, land cover, rainfall intensity and distance from drainage.
Size and distribution of the landslides varies locally, depending on the presence of the parameters mentioned above.Thus, the creation of an accurate and precise GIS-based landslide susceptibility map is entirely dependent on the availability of data related to controlling factors (Ayalew et al., 2004).Rock slides, debris slides, rock avalanches, rock fall, toppling, wedging, mud flows and debris flows are important landslide processes along the KKH.

a. Lithology
Time and type of slope failure is determined by the slope-building lithology.Each lithology is unique in terms of response to stresses and therefore features a particular susceptibility to potential slope failure (Vallejo and Ferrer, 2011).The KKH traverses a great variety of lithologies comprising sedimentary, igneous and metamorphic rocks.According to spatial analysis, Quaternary deposits, the Jijal complex and the Misgar slates exhibit the highest numbers of mass movement events (Fig. 6f).

b. Distance from faults
The Main Mantle Thrust (MMT), the Kamila Jal Shear zone (KJSZ), the Kamila Strike-slip Fault (KSF), the Raikot Fault, the Main Karakoram Thrust (MKT) and the Karakoram Fault are important structural features in close proximity of the Highway (Fig. 4).Landslides are concentrated along these active faults where rockmass is highly deformed (Ali et al., 2017).The fact that 54% of mapped landslides were found within a maximum distance of 1 km from these faults, while 69% were found within a 2km range, impressively substantiates the postulated strong control of structural features (Fig. 6e).

c. Geomorphologic Factors
Slope angle is an important geomorphic factor responsible for initiation of slope movements (Lee et al., 2004), to be considered for preparation of landslide susceptibility maps.Steep slopes are more susceptible to failure as compared to gentle ones.The study area demonstrates variation in topography ranging from steep to gentle slopes, high plains to narrow gorges and high cliffs.Slope steepness in the area has been divided into five classes (Fig. 7b).More than 50% of landslides occurred in class III (30°-45°) areas, whereas least landslide events (2%) occurred in class I and class V (0 -15° and >65°) areas (Fig. 6a).In addition to slope and elevation, aspect and curvature were also considered important factors for preparation of the landslide susceptibility map.However, in our area these parameters seem to have a reduced influence on landslide occurrence (Fig. 6b, 6c, 6d).

d. Hydrology
The proximity of streams has been considered for the preparation of landslide susceptibility maps by many researchers (Akgun et al., 2008;Basharat et al., 2016;Kanwal et al., 2016;Wang et al., 2015).In our study area, many small tributaries feed main rivers, Indus and Hunza.During Monsoon season and after heavy precipitation events, these streams exhibit high-energy flows and large discharge, and are the main source of mud and debris flows.Heavy precipitation during monsoon and the westerlies triggers many landslides by increasing pore water pressure in unconsolidated sediments.A strong association between precipitation and mass movements along the highway has been found (Ali et al., 2017).

e. Land Cover
Variations in land cover control the spatial distribution of landslides along with other conditioning parameters (lithology, seismology, slope geometry) (Malek et al., 2015).Changes in land cover influence the hydrological condition of the slopes, leading to slope instability.Generally, vegetation tends to resist the erosion process whereas bare rock or soil is more susceptible to slope failure (Reichenbach et al., 2014).Restrepo and Alvarez (2006) found a strong relationship between land use and landslide events.Due to variations in the mountain ecosystem from Hassan Abdal to Khunjab pass, the KKH is surrounded by vegetation, bare rock/soil, water bodies and ice/snow covered slopes.The section of the highway between Hassan Abdal and Thakot is heavily vegetated due to considerable mean annual rainfall.From Thakot to Sazin, slopes are sparsely vegetated.From Sazin onward, barren rock slopes (Fig. 8) characterize the area.

Methodology
The flow chart (Fig. 9) describes the steps and techniques involved in preparation of the susceptibility map, involving multiple techniques, literature review, field reconnaissance and remote sensing.catalogues (Khan et al., 2000).FWO is responsible for clearance and maintenance of the KKH after its potential blockage.A landslide inventory map was prepared using ArcGIS 10.3, based on Geological Survey of Pakistan publications as well as realtime data acquired from FWO.A comprising geological and seismo-tectonic map (Fig. 4 and 5) of the area was prepared by digitizing and compiling various pre-existing maps (Khan et al. 2000).Data related to lithology, faults and seismicity has been extracted from these maps.Rainfall data from six weather stations along the KKH has been acquired from Pakistan Meteorological department to prepare a precipitation map.

b) Field Reconnaissance
Locations of landslides, lithological contacts and faults were validated and supplemented during a field visit.In addition to locations, types, size, failure mechanisms and structural control of landslides were determined.Acquired data was further used to prepare landslide inventory map within 2 km 2 around the Highway.

c) Remote Sensing
Geomorphological parameters (elevation, slope, aspect, curvature) and drainage were extracted from Shuttle Radar Topography Mission (SRTM) based DEM (30 m × 30 m).Thematic layers were prepared and classified using the natural break method in ArcGIS 10.3 (Fig. 7).Satellite images of Landsat 8 (19-21 November 2017) were acquired from USGS web portal and then pre-processed by using QGIS 2.18's Semi-Automatic Classification Plugin (SCP), followed by supervised classification of composite images in Arc GIS 10.3.The map was divided into four classes: vegetation, water bodies, glacial ice/snow and bare rock/soil (Fig. 8).The final product of the map was assessed using the confusion matrix method.Randomly distributed test pixels for each class were taken on the same image that had previously been classified.The accuracy of the land cover map added up to 87%.

d) Analytical hierarchy process (AHP)
Analytical hierarchy process (AHP) is a multi-criteria decision making approach to prepare landslide susceptibility maps.It has been used by previous researchers to assign weightage to landslide-controlling factors (Basharat et al., 2016;Kanwal et al., 2016;Shahabi and Hashim, 2015;Yalcin, 2008).It is based on the user's decision to weigh factors through their pairwise comparisons.Each factor is assigned a score (1-9) depending upon its relative importance, with increasing impact from 1 to 9 (Table 1, Saaty, 1990).The values assigned are based on spatial analysis of data, field observations and experience of the user.
If the parameter on the x-axis is more important than the one on y-axis, the value ranges between 1 and 9. Conversely, when the factor on y-axis is more important, the values are in reciprocals (1/2-1/9).Subsequently, the quality of the process is evaluated by consistency ratio (CR): Where CR is consistency ratio, CI is consistency index and RI is random index.
The value of CR indicates the inconsistency in the expert's decision during weighing the parameters.Values of CR lower than 0.1 prove the decision consistent, while values greater than 0.1 indicate inconsistency and suggest a revision of judgement.
The value of CR=0.07 in our study proves the consistency in procedure so that input parameters are further used for weighted overlay.

e) Weighted Overlay Method
For the weighted overlay technique, we used an overlay of raster layers of all controlling factors to prepare a susceptibility map.Raster layers of each controlling factor were reclassified and weighted according to their importance determined by AHP (Table 2).The cumulative weight of all input layers was maintained at 100.The completion of this process resulted into the ultimate production of a landslide susceptibility map of the highway.

Results
The produced landslide susceptibility map was classified in four classes: low susceptibility, moderate susceptibility, high susceptibility and very high susceptibility (Fig. 10 and 11).
Areas of 49.9% and 10.4% of the classified map respectively belong to the high susceptibility and very high susceptibility classes, particularly owed to the presence of active faults, seismic zones and steep slopes (Table 3).34.1% and 5.4% of the highway fall into intermediate and low susceptibility areas, respectively.
Owed to lucidity and space reasons, the final version of the map was divided into two parts: Sections (1) Hassan Abdal-Chilas section and (2) Chilas-Khunjrab Pass section (Fig. 10 and Fig. 11).The highway section 1 until Thakot, is characterized by broad valleys and gentle slopes covered by vegetation and therefore falls into the low to intermediate susceptibility zones (Fig 9).Contrastingly, the following section north of Thakot, particularly close to Jijal, lies in high and very high susceptibility zones (Fig. 10a).Threads arise from the presence of the southern suture (MMT), poor rockmass quality, the active IKSZ and steep slopes.In the section from Pattan to Sazin, more than half of the highway is at high susceptibility and very high susceptibility (Fig. 10b and 10c).This is because multiple shear zones (KJS) cross the highway between Pattan and Dassu and the surroundings of Sazin fall into the reach of Kamila strike slip fault (KSF) and the active HSZ (Fig. 10C).Two locations close to drainage features (Samar and Harbon Nala) near Sazin (Fig. 10c) also fall in the very high susceptibility zone.
The second section starts in Chilas and ends at the Khunjrab Pass (Chinese border).Some parts of the highway were found at very high susceptibility (Fig. 11).The Raikot bridge section is the most dangerous part of the highway as it lies directly over aggravating the situation (Fig. 11a).Due to presence of the Northern suture (MKT), loose glacial deposits and steep slopes in the Hunza section, some locations of the highway are declared very high susceptibility zones (Fig. 11b).Also north of Sost, two locations (Kafir Pahar and Notorious Killing zone) were found in very high susceptibility zones (Fig. 11c).

a) Accuracy assessment
Accuracy assessment of the map is an essential component of the whole process.In the past, different statistical techniques have been employed to check the predictive ability of a landslide susceptibility map: Predication Rate Curve (PRC), Landslide Density Analysis (LDA), Receiver Operating Curve (ROC) and Area Under Curve (AUC) (Deng et al., 2017).In this study, we also used ROC and LDA to estimate the predictive accuracy of the map.
In the first step, map classes were compared with landslide densities in their respective classes.Spatial analysis of the map and landslide events was performed on ArcGIS 10.3 using the tabulate area tool.According to the obtained results, most of the landslide events were found in high and very high susceptibility areas and very few landslides were present in moderate and low susceptibility zones respectively (Table 4).These statistics confirm a strong connection between susceptibility zonation and landslide events.Thus, this assessment indicates an adequate accuracy of the map.
ROC is the product of a graphical plot opposing true positive rate (TPR) and false positive rate (FPR) (Fig. 12).TPR indicates the correctly predicted events and is plotted on the y-axis while FPR indicates falsely predicated events, and is plotted against the x-axis.Area under curve (AUC) in a graphical plot explains the efficiency of the model.AUC may range from 0.5 to 1 in different cases depending on the accuracy of model.A value close to 0.5 indicates random results while values close to 1 indicate a perfect model.In this study, we used 72 landslide locations to validate the final version of the map.AUC was found 0.72 indicating a reputable accuracy (72%) of the map (Basharat et al., 2016;Deng et al., 2017).

Conclusions
A set of Landslide Susceptibility Maps of the KKH (CPEC) was prepared using ArcGIS 10.3 involving multiple techniques: literature review, remote sensing, field surveys.Ten controlling parameters (lithology, seismicity, rainfall intensity, distance from fault, elevation, slope angle, aspect, curvature, land cover and hydrology) were considered of which each one was assigned a numerical weight using AHP.Thematic layers of these parameters were then overplayed using the WOL tool of Although the part of the highway between Hassan Abdal and Thakot receives heavy precipitation in Monsoon season, the area is stable owed to mature geomorphology.According to results, active faults, slope gradient, seismicity and lithology have a strong influence on landslide events along the highway.In the final step, the predictive accuracy of the map was determined by using LDA and ROC.The accuracy of the map was rated to a satisfactory 72%, which is suitable for mitigation planning.
ArcGIS 10.3.Four different classes of landslide susceptibility were then applied to the final map: low susceptibility, moderate susceptibility, high susceptibility and very high Susceptibility.10 % and 50% of the highway were found in very high and high susceptibility zones, respectively.Active faults (MMT, KJS, KSF, RF, MKT, KF), seismic zones (IKSZ, HSZ, RSSZ) and steep slopes are responsible for the associated susceptibility in these areas.The highway is characterized by a variety of mass movements: rockfall, debris fall, rockslide, debris slide, debris flows and mudflows.The threatened sections are Nat.Hazards Earth Syst.Sci.Discuss., https://doi.org/10.5194/nhess-2018-39Manuscript under review for journal Nat.Hazards Earth Syst.Sci. Discussion started: 5 March 2018 c Author(s) 2018.CC BY 4.0 License.especially unstable in Monsoon and Westerlies seasons every year.Altogether, a detailed investigation is inevitable to enable hazard free and safe traveling.About 40% of the highway lies in low and moderate susceptibility zones, which remain almost stable throughout the year.The highway sections from Hassan Abdal to Thakot, near Chilas, Gilgit and Sost fall into these moderate and low susceptibility zones and are quite stable due to their course in broad U-shaped valleys with gentle slopes.

Figure 2 :
Figure 2: a) Glaciers along the highway (b) Passu Glacier's snout approaching the Highway (c) Relief along the highway (d) Profile drawn along axis drawn in 2c.

Figure 3 :
Figure 3: Precipitation map of the study area (Pakistan Meteorological Department)

Figure 8 :
Figure 8: Land Cover of the study area

Figure 9 :
Figure 9: Flow Chart showing multiple steps involved in the preparation of the susceptibility map.