Analysis of Landslide Movements Using Interferometric Synthetic Aperture Radar: A Case Study in Hunza-Nagar Valley, Pakistan

Abstract

From a geological standpoint, northern Pakistan is one of the most active and unstable areas in the world. As a consequence, many massive landslides have occurred in the area in historical times that have destroyed infrastructure, blocked the Hunza River, and damaged the Karakoram Highway repeatedly. However, despite the high frequency of large magnitude landslide events, and the consequent damages, the entire area is largely understudied, mainly due to the difficult logistics and the large distances involved. This work is aimed at applying the potential use of Interferometric Synthetic Aperture Radar (InSAR) for landslide identification and investigation for the Hunza-Nagar Region. Sentinel-1 images covering a period of more than two years (February 2017–August 2019) were used and processed by adopting the small baseline subset (SBAS) method. The obtained deformation rate measured along the line of sight (V_LOS) varies from −114 to 20 mm/year. The downslope velocity deformation rates (V_slope) range from 0 to −300 mm/year. The V_slope stability threshold for our study area was calculated to be −14 mm/year from the V_slope standard deviation. Four active landslides with V_slope exceeding 14 mm/year were recognizable and have been confirmed by field inspection. The identified landslides listed from the most active to least active are the Humarri, Mayoon, Khai, and Ghulmet landslides, respectively. V_LOS exceeding 114 mm/year was observed in the Humarri landslide, which posed a threat of damming a lake on the Hispar River and was also a risk to the Humarri Village located below the landslide. The maximum mean deformation detected in the Ghulmet, and Mayoon landslide was in the order of 30 mm/year and 20 mm/year, respectively. More importantly, it was found that in places, the slope deformation time series showed a patchy correlation with precipitation and seismic events in the area. This may indicate a complex, and possibly uncoupled, relationship between the two controlling agents promoting the deformation. However, the collective impact of the two factors is evident in the form of a continuously descending deformation curve and clearly indicates the ground distortion. The results indicate a potentially critical situation related to the high deformation rates measured at the Humarri landslide. On this specific slope, conditions leading to a possible catastrophic failure cannot be ruled out and should be a priority for the application of mitigation measures.

1. Introduction

Landslides are the most common and frequent hazards in the mountainous regions of the world which, directly or indirectly, affect the lives of the native dwellers [1]. In the period between 1998 and 2017, landslides, and their associated hazards, affected an estimated 4.8 million people and caused more than 18,000 deaths globally [2]. By improving the overall understanding of the science of this hazard, the safety of the communities could be ensured, as well as the annual damage and death toll could be significantly reduced. The northern rugged mountainous regions of Pakistan are prone to landslide hazards [3]. A higher number of recorded landslide events is linked to population growth and development in these areas [4].

The Hunza-Nagar region in northern Pakistan is well known for its unique geologic and geomorphologic features. Topographically, the Hunza-Nagar valley is characterized by deeply incised valleys and rugged mountains, with elevations ranging from 1746 m to 7315 m above sea level (Figure 1d). The topography has a strong influence on local climate in the region. The high relative relief causes rainfall to concentrate in the high-altitude areas (≥1000 mm/year), with much less precipitation in the valleys (250–500 mm/year) [5]. Northern Pakistan is located in a region which is seismically very active and has experienced many disastrous earthquakes in the past. The seismic activity in the region is demonstrated by about 300 M > 5 and seven M > 7 earthquakes events since 1980 (Kashmir Earthquake October, 2005: M = 7.6; Afghanistan October, 2015: M = 7.5) according to the USGS Earthquake Catalog, [6] (Figure 1b). Undercutting of landslide toes by Hunza River below the rock/debris fall upslope are major concerns in this area [7]. All the conditions mentioned above, coupled with the active tectonics, makes this region prone to geohazards. The authors of [8] have generated a landslide inventory map for the Hunza-Nagar valley by using the visual interpretation of the SPOT-5 satellite imagery and mapped a total of 172 landslides. Numerous landslide scars can be seen on the valley walls along with mass movement deposits commonly found on the valley floors and the lower part of the valley walls (Figure 2). The old landslide deposits and associated layers of very fine-grained late sediment of variable thickness in an otherwise steep gradient channel point towards the frequent damming of the Hunza River valley in the past [9].

In 1858, a landslide event occurred near the Salmanabad village, 3 km southwest of Attabad village, blocked the Hunza valley (Figure 2a). The lake formed because of the landslide and was sustained for six months [3]. The landslide dam was ultimately breached, causing a catastrophic flood that damaged substantial life and property along a >400 km long valley [9]. In 1992 (7–9 September), intense rainfall triggered the slope collapse at Sarat village, which affected the Karakoram Highway (KKH) and deposited enormous volumes of boulders several meters in diameter downstream [3,10]. In January 2010, a sizeable Attabad landslide dammed the Hunza River raising the water level to 100 m from the normal flow. At least 20 people died in the event that blocked the Hunza River, creating a landslide dam and lake that progressively extended 23 km upstream, submerging four villages. The landslide also blocked the KKH, a vital link to China, cutting off 26,000 people in Upper Hunza. The debris obstructed approximately three kilometers of the Hunza River and a long stretch of the KKH [11]. Even though the area is familiar with such events, this was the most significant landslide events in the last two decades.

The increasing frequency of such geohazard events can be attributed to the increasing population pressure in this region, which is rapidly pushing the human settlements closer to the areas which have a known history of recurring natural hazards. It thus has strong socio-economic effects on the overall development of the region because it directly threatens the infrastructure and human settlements [13]. The recent investments in the Hunza Valley under the initiative of the China Pakistan Economic Corridor (CPEC) has profoundly increased the installation of infrastructure including dams, road, tunnels, bridges, and ports, which have also increased the importance of hazard monitoring in the valley [14]. As CPEC will be the major direct trade route from China to Europe, the Middle East, and Africa, the probable blockage in trade due to landslide hazard will have global effects.

To ensure the safety of the communities, the investments in the area, and the continuity of trade in the region, the detection, and monitoring of such devastating effects of landslides is essential [13,15]. However, to detect and monitor the slope movements of such a complex area with a large number of slopes is a challenge [16].

Remote sensing (RS) is a powerful tool for landslide displacement detection [17,18]. It provides a solution in the form of a survey or early warning in those unstable areas where massive landslides occur frequently and accelerate suddenly [13,19]. The RS images offer a regional understanding of a region both in terms of its spatial spread as well as its chronological evolution. However, compared to traditional RS techniques such as Global Navigation Satellite System (GNSS), geophysical surveys, and optical remote sensing, the Interferometric Synthetic Aperture Radar (InSAR) technique has been widely used in landslide research due to its broad coverage, high spatial resolution, and ability to operate under all weather conditions [20,21,22]. InSAR is a remote sensing technique used to detect and quantify the superficial geomorphological changes on the Earth, which in ordinary cases might happen unnoticed in an inaccessible area. It has been widely used for 30 years for Earth remote sensing and reporting of many landslide events [23]. The InSAR technique can be applied at regional scales for slope movements monitoring, ground subsidence, glacier monitoring, and for studying the movements of the tectonic plate [24]. Many studies [13,25,26,27,28,29,30] have used the InSAR technique because of its wide-ranging applications which are achievable through the data with a high spatial and temporal resolution. Mainly, a permanent scatters InSAR (PS-InSAR) method introduced by [31] Ferretti et al. and a small baseline subset InSAR (SBAS-InSAR) method proposed by Berardino et al. have been widely used in related fields with millimetric accuracy [32,33].

SBAS method is based on the selection of multi-master interferograms having a short spatial and temporal baseline. Short spatial and temporal baselines overpower the spatio-temporal decorrelation effects [34], and therefore, it is valuable for ground deformation detection. The SBAS-SAR interferometry technique can detect broad areas of ground deformations with high accuracy, making this method appropriate to investigate surface deformation of both regional and local scale [22,35].

The previous work in the region [5,8,36,37,38] is based on deterministic and statistical relationships and regression investigation of landslides and their contributing factors. This paper, for the first time, uses the SBAS-InSAR technique to investigate the surface ground deformation in a highly active collision zone between two continents. This makes it a novel approach to detect landslide movements in the northern Pakistan in specific and in southeast Asia in general. The aim is to detect and monitor landslides in areas with rugged topography, a typical challenge for InSAR analysis, to highlight the presence of potentially hazardous landslides which were previously unknown to the attention of local authorities, and to ultimately prompt dedicated field surveys and confirmation studies of the priority cases. This technique would also offer a basis for the reduction and possible prevention of the landslide risk over a large area which is otherwise difficult to monitor by using field methods, due to the difficult logistics and long distances. Incidentally, this area is also home to several important local communities and hosts strategic water resources and transport networks of paramount importance for the connections between Pakistan and its north-border countries.

2. Geological and Geomorphological setting of the Study Area

The research area lies in the strategically important Hunza-Nagar Districts, northern Pakistan. The section of the Karakorum Highway (KKH) along the narrow Hunza Nagar valley connecting Pakistan and China can potentially impede the trade route in case of a major landslide event. This area hosts some of the highest peaks and active glaciers of northern Pakistan. The general structural trend in the Karakoram ranges along the Hunza-Nagar valley in the NW-SE direction (Figure 1c). The topography of the area is defined by steep-walled glacially scoured valleys dominated by arete and horn geometries. The valley shows enormous relative relief, and the topographic slopes of 60° and above adjacent to the glaciers are significant in contributing large masses of snow and rocks to glacier surfaces. The Passu, Hussaini, Gulmit, and Hoper are important glaciers of the area while Hunza and Haspir rivers control the drainage system of the area [7].

Tectonically the area is in the active collision zone of the Indian and Asian plates (Figure 1b) [12]. The crustal shortening, active faulting, and subduction are continuing with convergence and uplift rates of ~4–5 cm/year [12] and ~7mm/year, respectively [39,40]. The significant tectonic features in this area are Main Karakoram Thrust (MKT) and Karakorum Fault (KF), which are responsible for brittle deformation (Figure 1c) [41]. Due to this brittle deformation, the rock mass is highly fractured and jointed. The MKT is the collision zone of the southern margin of the Eurasian plate and extends into Baltistan through the Hashupa, Shigar, and Shyok valleys, respectively. MKT is a seismically active thrust fault that has a high angle along which many earthquakes occurred [42].

The study area includes granites, meta-volcanic, and meta-sedimentary rocks, as well as quaternary deposits (glacial moraines, debris deposits, flood plain deposits, and lake deposits) (Figure 1c) [43]. Sediments of various genesis like alluvial, debris-flow, morainic, lacustrine, and aeolian are irregularly dispersed along the Hunza-Nagar valley [44]. These deposits having a thickness of more than hundred meters unconformably overlaying the bedrock through the whole valley [43]. Most of these deposits are unconsolidated and lack cementation. Rugged topography, massive amounts of unconsolidated scree bounded by erosional river valleys, accompanied by steep slopes, define the geomorphological character of the study area, making it vulnerable to large-scale landslides. A combination of geomorphology with the sharp expected ground accelerations from seismic activity amplifies the landslide triggers.

3. Data and Methodology

SAR data is one of the most significant and available remote sensing resources which can be used in the investigation of the Earth’s surface deformation [45]. It can be applied at various scales from regional to local, for example, for faults monitoring, landslide detection, and analysis of the movement of a landslide. In this study, 37 scenes of Sentinel-1 multi-temporal dataset IW (Interferometric Wide-swath) mode SLC (Single Look Complex) acquired from February 2017 to August 2019 were used to obtain land deformation information. The data were acquired from the Copernicus Open Hub. Spatially, this series of images covers part of northern Pakistan and contains the area of interest for this research. To eliminate the residual atmospheric and topographic noise, digital elevation model 1 arc-seconds SRTM was used.

The SBAS-InSAR technique is adopted to process Sentinel-1 data in the ENVI SARScape module to get the surface deformation map of the study area. This technique is an advanced InSAR method that could improve the monitoring precision [46]. The SBAS-InSAR method is based on a suitable combination of differential interferograms within the thresholds of temporal and spatial baselines, so the geometric decorrelation is insignificant [33,47,48,49]. Moreover, the SBAS-InSAR method is more appropriate to regions, which lack persistent enough scatters [50].

The processing of the SBAS-InSAR technique starts from the connection graph in which pairs within small spatial and temporal baseline were chosen (Figure 3). The image acquired on 2019-03-15 is selected as the super master image (SMI). The SMI is the reference image of the entire operation, and all the pairs co-registered on SMI reference geometry. The connection graphs are drawn with segments, consistent each to an interferometric couple to be processed. Each related point links to an input acquisition. The chosen acquisitions are illustrated as green points and the SMI as the yellow point. In this study, the maximum normal baseline was set to 45% of the critical baseline, and a maximum temporal baseline was set to 120 days. One hundred and forty-one differential interferograms were generated in this step (Figure 4).

In the next step of the interferometric process, the Goldstein filter method was implemented to improve the signal/noise ratio of the interferograms. The minimum cost flow (MCF) method and 3D-unwrapping method was chosen to finish phase unwrapping, and a coherence threshold of 0.35 was selected [51]. Refinement and re-flattening were performed to analyze and remove the residual phase coefficients and phase ramps from the unwrapped phase stack using 34 ground control points (GCPs). The GCPs were selected from the better-unwrapped phase area with the coherence value >0.7 [52]. Some interferometric pairs with poor unwrapping and low coherence were removed. As for the estimation of the extra height and the displacement related information through steps of refinement and re-flattening and inversion, preliminary displacements were estimated by the cubic model that is robust and commonly used [52].

The residual atmospheric phase was removed from the displacement by high-pass temporal filtering and low-pass spatial filtering with window size of 365 days and 1000 m respectively. Subsequently, geocoding in the line of sight (LOS) direction with a resolution of 10 m is employed. To further study the characteristics of the unstable slope, the velocity in the line of sight (V_LOS) is converted into velocity in the slope direction (V_slope) using the equation adopted from [13].

$${\mathrm{V}}_{\mathrm{slope}}=\frac{{\mathrm{V}}_{\mathrm{los}}}{\mathrm{Index}}$$

(1)

the Index can be formulated as below.

$$\mathrm{Index}={\mathrm{n}}_{\mathrm{los}}\times {\mathrm{n}}_{\mathrm{slope}}$$

(2)

$${\mathrm{n}}_{\mathrm{los}}=\left(-{\sin \mathsf{\theta }\cos \mathsf{\alpha }}_{\mathrm{s}},{\sin \mathsf{\theta }\sin \mathsf{\alpha }}_{\mathrm{s}},\cos \mathsf{\theta }\right)$$

(3)

$${\mathrm{n}}_{\mathrm{slope}}=(-\sin \mathsf{\alpha }\cos \mathsf{\phi },-\cos \mathsf{\alpha }\cos \mathsf{\phi },\sin \mathsf{\phi })$$

(4)

In this equation, ${\mathrm{V}}_{\mathrm{slope}}$ is the deformation rate along the slope angle, ${\mathrm{V}}_{\mathrm{los}}$ is the deformation rate of the LOS direction, α is the aspect of slope, θ is the radar incident angle, φ is the gradient degree of slope, ${\mathsf{\alpha }}_{\mathrm{s}}$is the angle between the direction of the satellite orbit and true north (which is the radar satellite flight direction), ascending data are positive, and descending data are negative. Finally, the time-series deformation was obtained and mapped to the study area.

4. Results and Analysis:

4.1. Spatial Distribution of Surface Deformation

The results of the regional deformation time series from February 2017 to August 2019 are shown in Figure 5 and Figure 6. The surface displacement rate was extracted for points with temporal coherence threshold >0.5. The total numbers of derived coherent targets (CTs) was 520,750 with average density of 612 CTs per Km². The deformation rate in the LOS varies from −114 to 20 mm/year (Figure 5). Positive values indicate that CTs move towards the satellite along LOS, whereas negative values indicate that CTs move away from the satellite along LOS. The CTs cover most of the hillslopes in the valley with the exception those around Attabad lake (Figure 5). The reasons of the low density around Attabad lake are the layover and the foreshortening caused by the westward and eastward hillslopes orientations with a high slope angle over 40°. CTs with high V_LOS (>30 mm/year) were detected around Attabad, Nagar, and Ghulmet villages. Moreover, large surface deformation was also detected on Passu glacier deposits with V_LOS >50 mm/year (Figure 5). Since, depending on the slope orientation, even a positive V_LOS could indicate a downslope movement [13,53], the V_LOS was converted into displacement rates along the steepest slope direction (V_slope) to highlight the actual downslope component of the slope deformation. The SBAS-derived downslope velocity deformation rates (V_slope) range from 0 to −300 mm/year (Figure 6), and the location of CTs with higher V_slope is consistent with V_LOS. Considering that a positive value of V_slope indicates uphill movement, CTs with positive value were not included in the final results which induced a lower density of CTs in V_slope map (Figure 6). Several previous studies identify the unstable areas based on the displacement rate overcoming a constant threshold calculated as the standard deviation of V_slope when the slope cell contains a sufficient number of CTs with V_slope over the stability threshold [21,22,54,55,56]. The V_slope stability threshold for our study area was calculated to be −14 mm/year from the V_slope standard deviation. And the active landslides and unstable slopes were identified as a slope cell containing at least three CTs with V_slope exceeding 14 mm/year.

A field survey of the anomalous areas was conducted to check the reliability of InSAR results. Among the coherent locations characterized by high downslope deformation rates and field surveys, we selected the landslides which either show strong ground movement or pose a direct threat to the nearby population. The detail results obtained on such landslides are as follows.

4.1.1. Mayoon Landslide

Mayoon landslide is located in the north of the Mayoon village, which rests on the Hunza River’s right bank and is 15 km northwest of the Rakaposhi peak. The landslide covers a surface area of around 1500 m². The general direction of the slide is towards the south, and its slope gradient varies from 20° to 45°. The landslide body is consisting of Yasin group sediments overlain by quaternary moraine deposits (Figure 1c). Deformation with mean V_LOS >20 mm/year and V_slope >30 mm/year was detected on the landslide (Figure 7a). Due to continuous slip along the slide, a series of parallel cracks on the main scarp face initially widen before the major failure, while new cracks emerge on the yet undisturbed scree slope up the hill. The cracks’ width varies with the continuous sliding movement; however, a rough estimate of the opening along the cracks range between 0.1 and 5 m (Figure 7 b,d). According to local knowledge, the landslide event was first triggered in 1976, damaging the agricultural land. Another slope failure occurred on the eastern end of the escarpment in 2011, damaging the infrastructure. Another event of relatively smaller magnitude occurred in 2012, which prompted the evacuation of 20 families from the area.

4.1.2. Ghulmet Landslide

The Ghulmet landslide covers about 0.7 km² area, located on the right bank of the Hunza River near Ghulmet village. It is an active rotational landslide with a slope direction towards the south (Figure 8). Tectonically the landslide is bounded by MKT and Chalt faults (Figure 1c). The landslide body is consisting of meta-sediments overlain by quaternary moraine deposits (Figure 1c). The slope gradient varies from 30° to 55°. A maximum mean deformation rate of >30 mm/year was detected in the LOS direction (Figure 8a). The upper sector of the landslide exhibits high V_LOS as compared to the lower section of the landslide (Figure 8a). Deformations with V_slope >70 mm/year were found in the slope (Figure 6). Another significant slope deformation was located in the downslope area of morphologically distinct scarps, which shows the existence of relict landslide masses in this zone. A giant crack is visible on the optical satellite imagery, and it has been observed during the field survey as well (Figure 8b). Based on the sliding direction of the slope and the gradient in this area, the landslide is moving towards the Hunza River and the KKH, posing a potential threat as a landslide dam for the Hunza River and the highway in case of sudden mass failure.

4.1.3. Khai Landslide

The Khai landslide is a small landslide of about 500 m in length. The landslide covers a 150,000 m² area and ranges in elevation from 2100 m to 2370 m from sea level. The scar of the landslide is visible on satellite imagery acquired on 05 May, 2017 (Figure 9a). High surface deformation rates are detected in the middle part of the landslide. V_LOS ranges from −10 to −30 mm/y (Figure 9a).

Despite the notable active displacement of the previous three cases, the most pronounced surface deformation identified in this study is situated shortly south of Khai, near Humarri village (Figure 5). Here, a significant slope deformation has been identified, which has not been detected previously. The deformation results from the Humarri village are of great importance and are addressed in detail in Section 5.

4.2. Temporal Change in Displacement Pattern in Relationship to the Precipitation and Seismicity

In addition to analyzing the mean deformation rates, we investigated the temporal correlation of the identified surface deformations, in the period between February 2017 and August 2019, with the rainfall pattern and the seismicity in the same time frame. To this end, we plotted the values of the displacement times series against the monthly rainfall data and the temporal occurrence of significant seismic events in the area. The data were obtained from the Meteorological Department and Water and Power Development Authority (WAPDA) Pakistan (

Table 1. The significant earthquake recorded from February 2017 to August 2019 in the northern Pakistan.

Date	Magnitude	Depth (km)	Location	Distance of the Research Area from Epicenter
2019-06-12	5.2	12.90	34.681°N, 72.857°E	230 km
2019-02-05	5.2	20.67	35.103°N, 74.169°E	120 km
2018-07-07	4.3	83.27	35.692°N, 74.446°E	55 km
2018-05-09	6.2	116.00	36.994°N, 71.382°E	270 km
2017-08-23	4.6	49.61	35.407°N, 74.578°E	85 km
2017-05-06	4.2	77.30	36.771°N,73.982°E	70 km
2017-03-06	4.4	90.32	37.051°N, 73.504°E	130 km

).

The relationship of the total deformation and the monthly average rainfall over February 2017 to August 2019 can be found in Figure 7c, Figure 8b, Figure 9c, and Figure 10c. The study area is typically humid from May to October and becomes very dry from November to April, revealing an apparent periodicity. Time series analysis reveals that the displacement in Ghulmet, Mayoon, and Khai landslides may be partially correlated to the precipitation data; however, in general, the displacement curve does not truly follow the precipitation or the seismic events. Therefore, it is possible to assume that there is no clear causal relationship between rainfall and seismic stress in the selected cases. In some cases, however, some weak correlation is visible, as highlighted in the following. The time series of displacements shows that the Ghulmet landslide has reactivated over the time covered by this study (Figure 8c). Analysis of the deformation pattern indicates that the movement of the landslide has some significant relationships to the changes in rainfall seasonal regime. For example, the deceleration in the beginning of November 2017 ends after two months into the new rainfall season, in October 2018. Again, the short stationarity observed from February to April 2019 follows a short dry period in January 2019. However, the rainfall peaks at the end of the recorded period do not seem to force any significant response into the displacement pattern, implying that different factors are controlling the slope deformation at this stage. The Mayoon landslide showed a sudden increase in the deformation corresponding to an increase in rainfall intensity in August and September 2017 with a one-month delay (Figure 7c). Other rainfall events show only a scanty correlation with rainfall, though. For example, the dry period from January to April 2018 does not seem to influence the displacement velocity in the main landslide body (Figure 7a, points b and c) and it is only capable of inducing a deceleration upslope of the main scrap (Figure 7a, point a). A slightly stronger positive correlation is found between deformation time series and rainfall data in the Khai landslide, which indicates that the part of landslide might possibly be triggered by rainfall (Figure 9b). However, the response to rainfall variations is much more visible for the central part of the landslide body (Figure 9a, point a) than for the right flank (Figure 9a, point b).

Seismic activities in and around Hunza valley are also thought to perform an essential part in predisposing and triggering landslides [5]. In addition to the rainfall, the correlation between deformation time series and earthquake events during the investigated periods shows that the deformations may be, in part, significantly corresponding to the earthquake events in the area. The temporal evolution of the deformation along the selected points in Ghulmet and Mayoon landslides show a distinct acceleration in deformation in March and August 2017 (Figure 7c and Figure 8c). These accelerations in the deformation are in functional correlations with the earthquake events in these months. Higher mean deformation rates were recorded in the Humarri landslide following the March 2017 earthquake (Figure 10c).

By associating the time series of deformation with the precipitation and earthquakes events, it can be concluded that the movement of the landslides partly corresponds to the rainfall and earthquake events around the study area. This partial correspondence can be attributed to the complex interference of triggering agents and their overlapping role in controlling a specific landslide event. The interplay between seismic forcing, monthly rainfall changes and possibly existing, dynamic factors in the area produce the spurious single-variable correlations seen in the above cases. A deeper understanding of such relationships may be possibly gained by focusing on local slope scale studies supported by high-frequency rainfall measurement and seismic ground response analysis.

5. Humarri Landslide

The Humarri landslide is situated on the hanging wall of the MKT in the south-west of Humarri village in district Nagar (Figure 5). This village is located on the valley slope of Glacio-fluvial deposits underlain by high grade metamorphic rocks of Karakoram block (Figure 11a) on the left bank of the Hispar River. A local fault named the Humarri fault trending in the direction of NE-SW passes through the Landslide (Figure 8a). It is an active landslide with steep slopes ranging 30° to 70°. The Humarri landslide’s crest is highly fractured, and the fractures in its crown area are even visible on the satellite imagery (Figure 10a and Figure 11a). This region is located between the Hunza River and the Hoper glacier in the MKT zone (Figure 11a). The mass movements’ quantity in the hanging wall block of the MKT has been observed previously to be higher as compared to the footwall block [38]. In order to demarcate the geometry of the landslide material, subsurface model is constructed using surface geology. The landslide area is estimated to be 3.4 km², with a volume of 4.1 × 107 m³. This is a massive landslide in size, with a length of 2.3 km, a width of 1.3 km, and a maximum thickness of about 250 m (Figure 11a,b).

The area with significant rates of deformation was observed to be concentrated around and close to the Humarri fault (Figure 10a). The slope is orientated in the north-east direction. The sliding body consists of loose glacio-fluvial deposits, while the bedrock is mainly composed of ortho and para-gneisses. The deformation map of the hillslope that shows average velocities within the time series for the selected CTs is given in Figure 10. V_LOS ranges between −30 and −114 mm/year, and locally even higher values, have been observed inside the landslide compared to stable points beyond its periphery (Figure 10a). Large deformation rates (>100 mm/year) are detected in the middle and lower part of the landslide. The significant slope deformation has never been detected previously in any study related to this area, and the considerable deformation is a clear indication of a potential major landslide. The different sectors of the slide show different velocity patterns. The deformation time series of the landslide reveal different deformation rates corresponding to their different gradients. These different deformation rates may be related to the nature of the materials involved in the displacement, which are mainly loose glacio-fluvial deposits and thus have an intrinsic heterogeneity.

According to the slope orientation and gradient in the region, this landslide is moving towards the Humarri village, Hispar River, and Nagar Khas road (Figure 11a). The valley is only about 80 m wide at this place. In the event of Humarri landslide triggering, which is much larger than Attabad landslide, the potential of a huge lake damming the Hispar River is possible (Figure 11c). This situation is threatening in the form of a possible landslide dam that would also, most probably, block the road connecting Nagar district to KKH and would be devastating to the community living in the near vicinity of the landslide. The large active deformation rate in the area is evident from the large fractures in the walls of many houses of the Humarri area (Figure 11c,d). The Humarri village alone comprises 80 households and 700 individuals approximately. This inhabited area may be in danger if any event triggers and accelerates the landslide movement.

6. Analysis of Potential triggers to Humarri Landslide

Karakoram-Himalayas is one of the most active mountain belts of the world that originated after the continent–continent collision of the Indian Plate with the Eurasian Plate. Being an ongoing collision, the major faults, along with their active splays, frequently release seismic energy. Major tectonic faults, including the MKT and the Karakoram Fault, play a significant role in triggering seismic events. Seismic activities in and around the region of Hunza and Nagar are thought to play a critical role in predisposing and triggering landslides. Earthquakes in 2002 (Mw 6.3), 2004, and 2005 (Mw 7.8) are the most evident among many tremors and jolts [57]. The more recent Mw 7.5 earthquake in the region of Yamgan, Afghanistan, on October 26, 2015, with its epicenter being around 300 km away from the Hunza-Nagar region, also triggered numerous landslides in the study area [6]. Interviews of the local people were conducted during the field survey, and it was deduced that the cracks and damages in the houses started to appear after the 2015 earthquake. However, no indicative satellite images were available for the detailed examination of the impact of the mentioned earthquake event. Figure 12a,b represents that the scar and cracks were already present before 2013, and no sudden accelerations in the scars can be seen from the satellite acquired after the 2015 earthquake. The slow-moving landslides typically respond to precipitation that results in an upsurge in pore-water pressure and a reduction in the effective normal stress along basal shear zones [58]. The temporal correlation analysis between the precipitation and the movement of the landslide were examined through a time series analysis. The association of the total deformation and the monthly average precipitation over February 2017 to August 2019 can be found in Figure 8c. Due to the unavailability of meteorological station information for the investigated landslide, rainfall data from Hunza station was used to evaluate the results. However, a correlation between the precipitation and the deformation was not evident. The deformation rates from February 2017 to August 2019 are quite linear and uncorrelated to precipitation. That could mean that other factors are at play, such as irrigation canals in the area, river incision, and tectonics that produce slope dynamics leading to a fast (but linear) displacement trend that may typical be of loose materials [59] (Figure 10c).

Irrigation is one of the factors which can trigger landslide by increasing groundwater level [60,61]. In north Pakistan, the irrigation channels are diverted from the glaciers at a higher elevation, and then by maintaining a gentle slope, these irrigation channels run for kilometers at higher elevation before they reach the target population cluster. Because these irrigation channels are openly dredged into mud, debris, and sand, the water easily percolates through the channel base and weakens the slope. In the upper part of Humarri landslide near the fresh escarpments, there are two irrigation channels that supply water for agricultural purposes (Figure 12d). These channels are derived from the glacial valley in the east and seem active as vegetation could be traced along with them. The irrigation system may leak a large amount of melt water in the predicted landslide area, and the slope is being affected continuously. The springs along the slope are indicative of soil saturation in the terrace area (Figure 11a).

The continuous erosion and destabilization of the toe of the Humarri landslide by the Hispar river adds to the slope instability (Figure 12c). When the soil at the toe crumbles, slumps, or liquefies, large areas of the bluff toe support can then disastrously fail in the form of a landslide. The different sectors of the slide show different velocity patterns. The deformation time series of the landslide reveals different deformation rates corresponding to their different gradients. These different deformation rates may be related to the nature of the materials involved in the displacement, which are mainly loose glacio-fluvial deposits.

Interactions between tectonic uplift, landslide erosion on hillslopes, and valley incision lead to the steep topography of the mountainous landscapes [62,63]. Hillslopes in rapidly elevating landscapes are believed to respond to river incision into bedrock by steepening to a threshold [64]. The removal of soil material by river erosion causes variations in the stress distribution in the adjacent soil and thus can make the slope unstable. These changes might contribute to disturbing the static equilibrium that makes a slope stable and consequently results in slope failure, and thus, a landslide occurs [65]. The seismic activity, toe cutting, and irrigation system leakage are probably the most prominent triggering agents for the landslides in the research area followed by precipitation and the unconsolidated nature of the scree slope. However, all these triggering factors together contribute to the ultimate slip pattern creating a net result.

7. Conclusions

The focus of the current study was to verify the capacity of the InSAR technique in rugged terrain using Sentinel-1 data for the detection of landslide deformation. The SBAS time series InSAR technique, applied to process Sentinel-1 images, successfully identified the landslides in the Hunza-Nagar area. This ultimately helped to highlight the rates of movement of these landslides, which in some of the cases raise concerns on the possible development of accelerations. The landslides identified through the InSAR technique, and the displacement velocities measured, remarkably correspond to the actual ground situation surveyed on the field, where scars, cracks on the ground and on buildings, and infrastructures were clearly observed. Comparing ground conditions at the studied sites to the known recent events of landslides in the broader North Pakistan region, prompts the conclusion that these unstable slopes may progress towards potential catastrophic failures in the near future. This is especially true for the Humarri landslide, previously undetected, that shows V_LOS velocities overcoming 100 mm/year in the surveyed period. The successful application of InSAR analysis based on freely available radar images such as ESA’s Sentinel-1 is a potential breakthrough for the application of the same methods over much larger areas. Such areas may include as the whole northern Pakistan including the CPEC route for ensuring the identification of all possible locations for the landslides which have the potential to either directly impact the communities or could create ground deformation related geohazards.