Inferring potential landslide damming using slope stability, geomorphic constraints and run-out analysis; case study from the NW Himalaya

1 Prediction of potential landslide damming has been a difficult process owing to the 2 uncertainties related to the landslide volume, resultant dam volume, entrainment, valley 3 configuration, river discharge, material composition, friction, and turbulence associated with 4 material. In this study, instability pattern of landslides, geomorphic indices, post failure run-5 out predictions, and spatio-temporal pattern of rainfall and earthquake are explored to predict 6 the potential landslide damming sites. The Satluj valley, NW Himalaya is chosen as a case 7 study area. The study area has witnessed landslide damming in the past and incurred $ ~30M 8 loss and 350 lives in the last four decades due to such processes. Forty-four active landslides 9 that cover a total ~4.81 ± 0.05 x 10 6 m 2 area and ~34.1 ± 9.2 x 10 6 m 3 volume are evaluated to 10 identify those landslides that may result in the potential landslide damming. Out of forty-four, 11 five landslides covering a total volume of ~26.3 ± 6.7 x 10 6 m 3 are noted to form the potential 12 landslide dams. Spatio-temporal varying pattern of the rainfall in the recent years enhanced 13 the possibility of landslide triggering and hence of the potential damming. These five 14 landslides also revealed 24.8 ± 2.7m to 39.8± 4.0m high debris flow in the run-out 15 predictions. 16


INTRODUCTION
Landslide damming is a normal geomorphic process in narrow river valleys and poses significant natural hazard (Dai et al. 2005;Gupta and Sah 2008;Delaney and Evans 2015;Fan et al. 2020).Many studies have explored damming characteristics (Li et al. 1986;Costa and Schuster 1988;Takahashi and Nakawaga 1993;Ermini and Casagli 2003;Fujisawa et al. 2009;Stefanelli et al. 2016;Kumar et al. 2019a).However, studies concerning the prediction of potential landslide dams and their stability at regional scale have been relatively rare, particularly in Himalaya despite a history of landslide damming and flash floods (Gupta and Sah 2008;Ruiz-Villanueva et al. 2016;Kumar et al. 2019a).In order to identify the landslides that have potential to form dams, the following factors have been main requisites; (i) pre-and post-failure behaviour of landslide slopes, and (ii) landslide volume, stream power, and morphological setting of the valley (Kumar et al. 2019a).
To understand the pre-failure pattern, Finite Element Method (FEM) based slope stability evaluation has been among the most widely used approaches for complex slope geometry (Griffiths and Lane 1999;Jing 2003;Jamir et al. 2017;Kumar et al. 2018).However, the selection of input parameters in FEM analysis and the set of assumptions used (material model, failure criteria, and convergence) may also result in uncertainty in the final output (Wong 1984;Cho 2007;Li et al. 2016).Uncertainty from input parameters can be resolved by performing parametric analysis, whereas the utilization of most appropriate criteria can minimize the uncertainty caused by assumptions.Post-failure behavior of landslides can be understood using run-out analysis (Hungr et al. 1984;Hutter et al. 1994;Rickenmann and Scheidl 2013).These methods could be classified into empirical/statistical and dynamical categories (Rickenmann 2005).Owing to the flexibility in rheology, solution approach, reference frame, and entrainment, dynamic models have been relatively more realistic for site-specific problems (Corominas and Mavrouli 2011).Though the different numerical models have different advantages and limitations, Voellmy rheology (friction and turbulence) (Voellmy 1955;Salm 1993) based Rapid Mass Movement Simulation (RAMMS) (Christen et al. 2010) has been used widely owing to the inclusion of rheological and entrainment rate flexibility.
Apart from the pre and post-failure pattern, landslide volume, stream power and morphological setting of the valley are crucial to infer the potential landslide damming.The Morphological Obstruction Index (MOI) and Hydro-morphological Dam Stability Index (HDSI) have been widely used geomorphic indices to infer the potential of landslide dam formation and their temporal stability (Costa and Schuster 1988;Ermini and Casagli 2003;Stefanelli et al. 2016).
The NW Himalaya has been one of most affected terrains by the landslides owing to the active tectonics and multiple precipitation sources i.e., Indian Summer Monsoon (ISM) and Western Disturbance (Dimri et al. 2015).The NW Himalaya has accommodated ~51 % of all the landslides in India during yrs.1800-2011 (Parkash 2011).The Satluj River valley, NW Himalaya is one such region where landslides and associated floods have claimed ~350 lives and resulted in the loss of minimum 30 million USD in the last four decades.This region holds a high potential for future landslide damming and resultant floods (Ruiz-Villanueva et al. 2016;Kumar et al. 2019a).Therefore, the Satluj valley is taken as a case study area, and 44 active landslides belonging to the different litho-tectonic regimes are modeled using the FEM technique.Multiple slope sections and a range of values of different input parameters are used to perform the parametric study.In order to determine the human population that might be affected by these landslides, census statistics are also used.The MOI and HDSI are used to determine the potential of landslide dam formation and their stability, respectively.In view of the role of rainfall and earthquakes as main landslide triggering factors, the spatio-temporal regime of these two factors is also discussed.Run-out prediction of certain landslides is also performed to understand the role of run-out in the potential landslide damming.This study provides detailed insight into the regional instability pattern, associated uncertainty, and potential landslide damming sites and hence it can be replicated in other hilly terrain witnessing frequent landslides and damming.
The Main Central Thrust (MCT) fault separates the HHC from the underlying schist/gneissic rockmass of the LHC.The LHC comprises mica schist, carbonaceous schist, quartzite, and amphibolite.A thick zone of gneiss i.e., Wangtu Gneissic Complex (WGC) is exposed in the LHC, which comprises augen gneiss and porphyritic granitoids.The LHC is delimited at the base by the Munsiari Thrust (MT) fault that is thrusted over the Lesser Himalaya Sequence (LHS) rockmass.The MT contains breccia, cataclastic, and fault gouge (Sharma 1977;Vannay et al. 2004;Kumar et al. 2019b).The LHS in the study area consists of quartz-arenite (Rampur Quartzite) with bands of phyllite, meta-volcanics, and paragneiss (Sharma 1977).
The present study covers forty-four active landslides (20 debris slides, 13 rock falls, and 11 rock avalanches) along the study area (Table 1) that have been mapped recently by Kumar et al. (2019b).Field photographs of some of these landslides are presented in Fig. 2. The TS and LHS in the study area have been subjected to relative tectonic tranquility with exhumation rates as low as 0.5 -1.0 mm/yr, whereas the HHC and LHC region have undergone 1.0 -4.5 mm/yr rate of exhumation (Thiede et al. 2009).The MCT fault region and the WGC are noted to have maximum exhumation rate (i.e., ~4.5 mm/yr) that is evident from the deep gorges in these regions (Fig. 2c,2e).A majority of the earthquake events in the study area in the last 7 decades have been related to the N-S oriented Kaurik -Chango Fault (KCF) (Kundu et al. 2014;Hazarika et al. 2017;International Seismological Centre Catalogue 2019).The climate in the study area shows a spatial variation from humid (~800 mm/yr mean annual precipitation) in the LHS to the semi-arid (~200 mm/yr) in the TS (Kumar et al. 2019b).The HHC acts as a transition zone where climate varies from semi-humid to semi-arid in the SW-NE direction.
This transition has been attributed to the 'orographic barrier' nature of the HHC that marks the region in its north as 'orographic interior' and the region to its south as the 'orographic front' (Wulf et al. 2012;Kumar et al. 2019b).
Landslides in the study area have been a consistent threat to the socio-economic condition of the nearby human population (Gupta and Sah 2008;Ruiz-Villanueva et al. 2016;Kumar et al. 2019a).Therefore, the human population in the vicinity of each landslide was also determined by considering the nearby villages/town.Notably, a total of 25,822 people reside within 500 m extent of the 44 landslide slopes, and about 70 % of this population is residing in the reach of debris slide type landslides.Since the Govt. of India follows a 10 year gap in census statistics, the human population data was based on last official data i.e., Census-2011.The next official census is due in year 2021.The population density in the Indian Himalayan region was estimated to be 181/km² in the year 2011 that might grow to 212/km² in 2021 with a decadal growth rate of 17.3% (https://censusindia.gov.in,retrieved on 02 Sep 2020; http://gbpihedenvis.nic.in,retrieved on 02 Sep 2020).

METHODOLOGY
The methodology involved field data collection, satellite imagery analysis, laboratory analyses, slope stability modelling, geomorphic indices, rainfall/earthquake pattern and run-out modelling.Details are as follows;

Field data, satellite imagery processing, and laboratory analyses
The field work involved rock/soil sample collection from each landslide location, rockmass joint mapping, and N-type Schmidt Hammer Rebound (SHR) measurement.Joints were included in the slope models for the FEM based slope stability analysis.The dataset involving the joint details is available in the data repository (Kumar et al. 2021).The SHR values were obtained as per International Society of Rock Mechanics (ISRM) standard (Aydin 2008).
Cartosat-1 satellite imagery and field assessment were used to finalize the location of slope sections (2D) of the landslides.Cartosat-1 imagery has been used widely for the landslide related studies (Martha et al. 2010).The Cartosat-1 Digital Elevation Model (DEM) having 10m spatial resolution, prepared using the Cartosat-1 stereo imagery, was used to extract the slope sections of the landslides using the Arc GIS-10.2 software.Details of the satellite imagery are mentioned in Table 2.
The rock/soil samples were analyzed in the National Geotechnical Facility (NGF) and Wadia Institute of Himalayan Geology (WIHG) laboratory, India.The rock samples were drilled and smoothed for Unconfined Compressive Strength (UCS) (IS:    -1986).If the soil samples contained < 5% fines (< 75 mm), the hydrometer test was not performed for the remaining fine material.In the direct shear test, soil samples were sheared under the constant normal stress of 50, 100 and 150 kN/m 2 .The UCS test of soil was performed under three different rates of movements i.e., 1.25 mm/min, 1.50 mm/min and 2.5 mm/min.

2 Slope stability modelling
The Finite Element Method (FEM) was used along with the Shear Strength Reduction (SSR) technique to infer the critical Strength Reduction Factor (SRF), Shear Strain (SS), and Total Displacement (TD) in the 44 landslide slopes using the RS2 software.The SRF has been observed to be similar in nature to the Factor of Safety (FS) of the slope (Zienkiewicz et al. 1975;Griffiths and Lane 1999).To define the failure in the SSR approach, non-convergence criteria were used (Nian et al. 2011).The boundary condition with the restraining movement was applied to the base and back, whereas the front face was kept free for the movement (Fig. 3).In-situ field stress was adjusted in view of dominant stress i.e., extension or compression, by changing the value of the coefficient of earth pressure (k).A value of k = σh/σv = 0.5 was used in extensional regime, whereas k = σh/σv = 1.5 was used in compressional regime.The Tethyan Sequence has been observed to possess the NW-SE directed extensional regime.The structures in the upper part of the HHC are influenced by the east directed extension along the SD fault.The lower part, however, is characterized by the SW directed compression along the Main Central Thrust.In contrast to the HHC, structures in the Lesser Himalaya Crystalline and Munsiari Thrust region are influenced by the compressional regime.In the Lesser Himalaya Sequence region, the SW directed compressional regime has been observed on the basis of the SW verging folds, crenulation cleavage, and other features (Vannay et al. 2004).
The soil and rock mass were used in the models through the Mohr-Coulomb (M-C) failure criterion (Coulomb 1776;Mohr 1914) and Generalized Hoek-Brown (GHB) criterion (Hoek et al. 1995), respectively.The parallel-statistical distribution of the joints with normal-distributed joint spacing in the rock mass was applied through the Barton-Bandis (B-B) slip criterion (Barton and Choubey 1977;Barton and Bandis 1990).Plane strain triangular elements having 6 nodes were used through the graded mesh in the models.Details of the criteria used in the FEM analysis are mentioned in Table 3.The dataset of input parameters used in the FEM analysis is available in the data repository (Kumar et al. 2021).It is to note that the FEM analysis was performed under the static load i.e., field stress and body force.The dynamic analysis was not performed, at present, in absence of any major seismic events in the region in the last 4 decades (sec.4.3) and lack of reliable dynamic load data of nearby major seismic events.
To understand the uncertainty caused by the selection of 2D slope section, multiple slope sections were taken, wherever possible.More than one slope section was modeled for each debris slide, whereas for the rock falls/ rock avalanche only one slope section was chosen due to the limited width of the rock falls/rock avalanche in the study area.To find out the relative influence of different input parameters on the final output, a parametric study was performed.
The selection of these landslides for the parametric study was based on the following two factors; (1) to choose the landslides from different litho-tectonic regime, and (2) to represent varying stress regime i.e., extensional, compressional, and relatively stagnant.The Parametric study of the debris slide models involved following 9 parameters; field stress coefficient, stiffness ratio, cohesion and angle of friction of soil, elastic modulus and Poisson's ratio of soil, rockmass modulus, Poisson's ratio and uniaxial compressive strength of rock.For the rockfalls/rock avalanche, the following 6 parameters were considered; uniaxial compressive strength of rock, rockmass modulus of rock, Poisson's ratio of rock, 'mi' parameter, stiffness ratio, and field stress coefficient.The 'mi' is a Generalized Hoek-Brown (GHB) parameter that is equivalent to the angle of friction of Mohr-coulomb (M-C) criteria.

3 Geomorphic indices
Considering the possibility of landslide dam formation in case of slope failure, the following geomorphic indices were also used;

4 Rainfall and Earthquake regime
Precipitation in the study area is related primarily to the Indian Summer Monsoon (ISM) and Western Disturbance (WD) and varies spatially-temporally due to various local and regional factors (Gadgil et al. 2007;Hunt et al. 2018).Therefore, we have taken the TRMM_3B42 (Huffman et al. 2016) daily rainfall data of years 2000-2019 at four different locations; Moorang, Kalpa, Nachar, and Rampur (Locations mentioned in Fig. 1).The dataset of earthquake events (2<M<8) in and around study area during the years 1940-2019 was retrieved from the ISC catalogue (http://www.isc.ac.uk/iscbulletin/search/catalogue/, retrieved on 02 March 2020) to determine the spatio-temporal pattern.

5 Run-out modelling
Since the study area has witnessed many disastrous landslides, mostly rainfall triggered, and flash floods in past (Gupta and Sah 2008;Ruiz-Villanueva et al. 2016), run-out analysis was performed to understand the post-failure scenario.Such run-out predictions will also be helpful to ascertain the possibility of damming because various studies have noted river damming by the debris flows (Li et al. 2011;Braun et al. 2018;Fan et al. 2020).The landslides that have potential to form dams based on the indices (sec.3.3) are evaluated for such run-out analysis.
The RAMMS for debris flow uses the Voellmy friction law and divides the frictional resistance into a dry-Coulomb type friction (μ) and viscous-turbulent friction (ξ).The frictional resistance S (Pa) is : where = hgcos() is the normal stress on the running surface, ρ is density, g is gravitational acceleration,  is slope angle, h is flow height and u= (ux, uy), consisting of the flow velocity in the x-and y-directions.In this study, a range of friction (µ) and turbulence (ξ) values, apart from other input parameters, are used to evaluate the uncertainty in output (Table 4).Generally, the values for μ and ξ are determined using the reconstruction of real events through the simulation and subsequent comparison between the dimensional characteristics of real and simulated event.However, the landslides in the study area merge with the river floor and/or are in close proximity and hence there is no failed material left from the previous events to reconstruct.Therefore, the µ and ξ values were taken from a range in view of topography of landslide slope and run-out path, landslide material, similar landslide events/material, and results from previous studies/models (H¨urlimann et al. 2008;Rickenmann and Scheidl 2013;RAMMS v.1.7.0).Since these landslides are relatively deep in nature and during the slope failure, irrespective of type of trigger, entire loose material might not slide down, the depth of landslide is taken as only ¼ (thickness) in the run-out calculation.Further, a release area concept (for unchanneled flow or block release) was used for the run-out simulation.During the field visits, no specific flow channels (or gullies) were found on the landslide slopes except a few centimeters deep seasonal flow channels for S. N. 5 and S.N. 15 landslides (Table 1).
However, the data pertaining to the spatial-temporal pattern of discharge at these two landslides was not available.Therefore, the release area concept was chosen because it has been more appropriate when the flow path (e.g.gully) and its possible discharge on the slope is uncertain

Slope instability regime and parametric output
Out of the 44 landslides studied here, 31 are in meta-stable state (1 ≤FS≤ 2) and 13 in unstable state (FS <1) (Fig. 4).Most of the unstable landslides are debris slides, whereas the majority of the meta-stable landslides are rock fall/rock avalanche.Debris slides constitute ~ 90 % and ~99 % of the total area and volume, respectively, of the unstable landslides.Notably, about ~70 % of the total human population along the study area resides in the vicinity (~500 m) of these unstable debris slides (Fig. 4).Rock falls/Rock avalanches constitute ~84 % and ~78 % of the area and volume, respectively, of the meta-stable landslides.Out of total 20 debris slides, 12 debris slides are found to be in unstable stage, whereas 8 in the meta-stable condition (Fig. 4).These 20 debris slides occupy ~1.9 ±0.02 x 10 6 m 2 area and ~ 26 ±6 x 10 6 m 3 volume.When comparing the Factor of Safety (FS) with the Total Displacement (TD) and Shear Strain (SS), nonlinear poor correlation is achieved (Fig. 5).Since the TD and SS are a relatively good correlation (Fig. 5), only the TD is used further along with the FS.The TD ranges from 7.4± 8.9 cm to 95.5± 10 cm for the unstable debris slides and ~18.8 cm for meta-stable landslides (Fig. 4).Out of 13 rockfalls, 1 belongs to the unstable state and 12 to the meta-stable state (Fig. 4).The TD varies from 0.4 to 80 cm with the maximum for Bara Kamba rockfall (S.N. 31).
Out of 11 rock avalanches, 1 belongs to the unstable state and 10 to the meta-stable state (Fig. 4).The TD varies from 6.0 to 132.0 cm with the maximum for the Kandar rock avalanche (S.N.

25
). Relatively higher TD is obtained by the rock fall and rock avalanche of the Lesser Himalaya Crystalline region (Fig. 4).The landslides of the Higher Himalaya Crystalline (HHC), Kinnaur Kailash Granite (KKG) and Tethyan Sequence (TS), despite being only 17 out of the total 44 landslides, constituted ~ 67 % and ~ 82 % of the total area and total volume of the landslides.
The Factor of Safety (FS) of debris slides is found to be relatively less sensitive to the change in the value of input parameters than the Total Displacement (TD) (Fig. 6).In case of Akpa (Fig. 6a) and Pangi landslide (Fig. 6b), soil friction and field stress have more influence on the FS.However, for TD, field stress, elastic modulus and Poisson's ratio of the soil are relatively more controlling parameters.The FS and TD of the Barauni Gad landslide (Fig. 6c) are relatively more sensitive to soil cohesion and 'mi' parameter.Therefore, it can be inferred that the FS of debris slides is more sensitive to soil friction and field stress, whereas TD is mostly controlled by the field stress and deformation parameters i.e, elastic modulus and Poisson's ratio.Similar to the debris slides, the FS of rock falls and rock avalanche are found to be relatively less sensitive than TD to the change in the value of input parameters (Fig. 7).Tirung Khad rock fall (Fig. 7a) and Baren Dogri rock avalanche (Fig. 7b) show dominance of 'mi' parameter and field stress in the FS as well as in TD.In case of Chagaon rock fall (Fig. 7c), Poisson's ratio and UCS have relatively more influence on FS and TD.Thus, it can be inferred that the rock fall/rock avalanche are more sensitive to 'mi' parameter and field stress.

Potential landslide damming
Based on the MOI, out of total 44 landslides, 5 (S.N. 5, 7, 14, 15, 19) are observed to be in the formation domain, 15 in uncertain domain, and 24 in non-formation domain (Fig. 8a).Thefive landslides that have potential to dam the river in case of slope failure comprise ~26.3 ± 6.7 x 10 6 m 3 volume (Fig. 9 a-e).In terms of temporal stability (or durability), out of these five landslides, only one landslide (S.N. 5) is noted to attain the 'uncertain' domain, whereas the remaining four show 'instability' (Fig. 8b,d).The lacustrine deposit in the upstream of Akpa landslide (S.N. 5) in Fig. 9a shows signs of landslide damming in the past (Fig. 10).The 'uncertain' temporal stability indicates that the landslide dam may be stable or unstable depending upon the stream power and landslide volume, which in turn are dynamic factors and may change owing to the changing climate and/or tectonic event.The landslides that have been observed to form the landslide dam but are noted to be in temporally unstable category (S.N. 7, 14, 15, 19) are still considerable owing to the associated risks of lake-impoundment and generation of secondary landslides.Urni landslide (S.N. 19) (Fig. 9e) that damaged the part of National Highway road (NH)-05 has already partially dammed the river since year 2016 and holds potential for the further damming (Kumar et al. 2019a).Apart from the S.N. 5 and S.N.
19 landslides, remaining landslides (S.N. 7, 14, 15) belong to the Higher Himalaya Crystalline (HHC) region that has been observed to accommodate many landslide dams and subsequent flash floods events in the geological past (Sharma et al. 2017).

Rainfall and Earthquake regime
In order to explain the spatio-temporal variation in rainfall, the topographic profile of the study area is also plotted along with the rainfall variation (Fig. 11a).The temporal distribution of rainfall is presented at annual, monsoonal i.e., Indian Summer Monsoon (ISM): June-September and non-monsoonal i.e., Western Disturbance (WD): Oct-May (Fig. 11b-d) level.
Rainfall data of the years 2000-2019 revealed a relative increase in the annual rainfall since the year 2010 (Fig. 11b).The Kalpa region (orographic barrier) received relatively more annual rainfall than the Rampur, Nachar and Moorang regions throughout the time period, except the year 2017.The rainfall dominance at Kalpa is more visible in the non-monsoonal season (Fig. 11d).This difference may be due to the orographic influence on the saturated winds of the WD (Dimri et al. 2015).Further, the rainfall during the monsoon season that was dominant at the Rampur region till year 2012 gained dominance at Kalpa region since the year 2013 (Fig. 11c).However, the contribution of the ISM and WD associated rainfall was variable in those years (Fig. 11).Such frequent but inconsistent rainfall events that possess varied (temporally) dominance of the ISM and WD are noted to owe their occurrence to the El-Nino Southern Oscillation (ENSO), Equatorial Indian Ocean Circulation (EIOC), and planetary warming (Gadgil et al. 2007;Hunt et al. 2018).The orographic setting is noted to act as a main local factor as evident from the relatively more rainfall (total precipitation=1748±594 mm/yr.) at Kalpa region (orographic barrier) in the non-monsoon and monsoon season from the year 2010 onwards (Fig. 11).Prediction of the potential landslide damming sites in the region revealed that four (S.N. 7, 14, 15, 19) out of five landslides that are predicted to be able to form dams belong to this orographic barrier region.Therefore, in view of the prevailing rainfall trend since the year 2010, regional factors, discussed above, and orographic setting, precipitation triggered slope failure events may be expected in the future.Such slope failure events, if they occur, at the predicted landslide damming sites may certainly dam the river.

Extreme rainfall events of
The seismic pattern revealed that the region has been hit by 1662 events during the years 1940-2019 with the epicenters located in and around the study area (Fig. 12a).However, ~99.5 % of these earthquake events had a magnitude of less than 6.0 and only 8 events are recorded in the range of 6.0 to 6.8 Ms (International Seismological Centre 2019).Out of these 8 events, only one event i.e., at 6.8 Ms (19 th Jan. 1975), has been noted to induce widespread slope failures in the study area (Khattri et al. 1978).The majority of the earthquake events in the study area occurred in the vicinity of the N-S oriented trans-tensional Kaurik -Chango Fault (KCF) that accommodated the epicenter of 19 th Jan. 1975 earthquake (Hazarika et al. 2017; http://www.isc.ac.uk/iscbulletin/search/catalogue/, retrieved on 02 March 2020).About 95% of the total 1662 events had their focal depth within 40 km (Fig. 12b).Such a relatively low magnitude -shallow seismicity in the region has been related to the Main Himalayan Thrust (MHT) decollement as a response to the relatively low convergence (~14±2 mm/yr) of India and Eurasia plates in the region (Bilham 2019) (Fig. 12c).Further, the arc (Himalaya)perpendicular Delhi-Haridwar ridge that is under thrusting the Eurasian plate in this region has been observed to be responsible for the spatially varied low seismicity in the region (Hazarika et al. 2017).Thus, though the study area has been subjected to frequent earthquakes, chances of earthquake-triggered landslides have been relatively low in comparison to rainfall-triggered landslides and associated landslide damming.For this reason and the lack of reliable dynamic load of major earthquake event, we have performed the static modelling in the present study.
However, we intend to perform the dynamic modelling in the near future if the reliable dynamic load data will be available.

Run-out analysis
All five landslides (S.N. 5, 7, 14, 15, 19 in Fig. 9) that are predicted to form potential landslide dams in case of slope failure were also used for the run-out analysis to evaluate expected runout distances in the event of reactivation and failure in the future.Results are as follows; 4. 4.1 Akpa landslide (S.N. 5) Though it is difficult to ascertain the extent to which the predicted debris flow might contribute in the river blockage, it will certainly block the river in view of ~38 m high debris material with ~50 m wide run-out across the channel in this narrow part of river valley (Fig. 9a) even at maximum value of coefficient of friction (i.e., µ =0.3) (Fig. 13a).Notably, not only the run-out extent but flow height also decreases on increasing the friction value (Fig. 13a.1-13.a.3).The maximum friction takes into account the shear resistance by slope material and the bed-load on the river channel.However, apart from the frictional characteristics of run-out path, turbulence of a debris flow also controls its dimension and hence consequences like potential damming.

Baren dogri landslide (S.N. 7)
At the maximum friction value (µ =0.4), the Baren dogri landslide would attain a peak value of flow height i.e., ~30 m at the base of central part of landslide (Fig. 13b).Similar to the valley configuration around the Akpa landslide (sec 4.4.1), the river valley attains a narrow/deep gorge setting here also (Fig. 9b).The maximum value of debris flow height obtained using the different µ and ξ values is 25.6 ± 2.1m (Fig. 14b).Flow material is also noted to attain more run-out in upstream direction of river (~1100 m) than in the downstream direction (~800 m).
This spatial variability in the run-out length might exist due to the river channel configuration as river channel in upstream direction is relatively narrower than the downstream direction.

Pawari landslide (S.N. 14)
The Pawari landslide attains maximum flow height of ~20 m at the maximum friction of runout path (µ=0.4) (Fig. 13c).The resultant debris flow that is achieved using the different values of µ and ξ parameters attains a peak value of 24.8 ± 2.7 m and decreases gradually with a runout of ~1500 m in upstream and downstream direction (Fig. 14c).This landslide resulted in the relatively long run-out of ~1500 in the upstream and downstream direction.Apart from the landslide volume that affects the run-out extent, valley morphology also controls it as evident from the previous landslides.The river channel in upstream and downstream direction from the landslide location is observed to be narrow (Fig. 9c).

Telangi landslide (S.N. 15)
The Telangi landslide would result in peak debris flow height of ~24 m at the maximum friction (µ=0.4) (Fig. 13d).On increasing the friction of run-out path, flow run-out decreased along the river channel but increased across the river channel resulting into possible damming.The debris flow after taking into account different values of µ and ξ parameters attains a peak value of 25.0± 4.0 m (Fig. 14d).Similar to Baren dogri landslide (S.N. 7), material attained more runout in the upstream direction of river (~1800 m) than in the downstream direction (~600 m) ; this difference can be attributed to a narrower river channel in upstream than the downstream direction.The downstream side attains wider river channel due to the Main Central Thrust (MCT) fault in the proximity (Fig. 1).Since the Pawari and Telangi landslides (S.N 14 &15) are situated ~500 m from each other, their respective flow run-outs might mix in the river channel resulting into disastrous cumulative effect.

Urni landslide (S.N. 19)
The Urni landslide is predicted to attain a peak value of ~44 m of debris flow height at the maximum friction value (µ=0.4) (Fig. 13e).After considering different values of µ and ξ parameters, the debris flow would attain a height of 26.3± 1.8 m (Fig. 14e).The relatively wider river channel in the downstream direction (Fig. 9e) results in longer run-out in downstream direction than in the upstream.

DISCUSSION
This study aimed to determine the potential landslide damming sites in the Satluj River valley, NW Himalaya.In order to achieve this objective, 44 active landslides were considered.First, slope stability evaluation of all the slopes, at these landslides sites was performed alongwith the parametric evaluation.Then the geomorphic indices, i.e., Morphological Obstruction Index (MOI) and Hydro-morphological Dam Stability Index (HDSI), were used to predict the formation of potential landslide dams and their subsequent stability.Rainfall and earthquake regime were also explored in the study area.Finally, run-out analysis was performed for those landslides that have been observed to form the potential landslide dam.
The MOI revealed that out of 44 active landslides in the Satluj valley, five (S.N. 5, 7, 14, 15, 19) have the potential to form the landslide dam (Fig. 8,9).On evaluating the stability of such potential dam sites using the HDSI, one landslide (S.N. 5) is predicted to attain an 'uncertain' domain (5.74<HDSI<7.44) in terms of dam stability.The uncertain term implies that the resultant dam may be stable or unstable depending upon the landslide/dam volume, upstream catchment area (or water discharge) and slope gradient (sec 3.3).Since this landslide (S.N.5) presents clear signs of having already formed a dam in the past, as indicated by the alternating fine-coarse layered sediment deposit (or lake deposit) in the upstream region (Fig. 10), recurrence is expected in the future.Further, run-out analysis of this landslide has predicted a 39.8± 4.0m high debris flow in the event of failure that will block the river completely (Fig. 13a,14a).However, the durability of the blocking cannot be ascertained as it will depend on the volume of landslide that will be retained in the channel and river discharge.
The remaining four landslides (S.N. 7, 14, 15, 19), though showing instability i.e., HDSI <5.74 at present, may form dams in the near future as the region accommodating these landslides has been affected by such damming and subsequent flash floods in the past (Sharma et al. 2017).
The last one of these i.e., S.N. 19 (Urni landslide) has already dammed the river partially and holds potential to completely block the river in near future (Kumar et al. 2019a).Run-out analysis of these landslides (S.N. 7, 14, 15, 19) has predicted 25.6 ± 2.1m, 24.8 ± 2.7m, 25.0± 4.0m and 26.3± 1.8m flow height, respectively that will result in temporary blocking of the river (Fig. 13,14).These findings of run-out indicate the blocking of river in the event of slope failure, irrespective of durability, despite the conservative depth as input because only ¼ of landslide thickness is used in the run-out analysis (sec.3.5).
Stability evaluation of these five landslide slopes (S.N. 5, 7, 14, 15, 19) that have potential to form landslide dams revealed that one landslide (S.N.7) is meta-stable, while the other four belong to the unstable category (Fig. 4).Further, these four unstable landslide slopes are debris slides in nature.It is noteworthy to discuss the implications of FS<1.The Factor of Safety (FS) in the Shear Strength Reduction (SSR) approach is a factor by which the existing shear strength of material is divided to determine the critical shear strength at which failure occurs (Zienkiewicz et al. 1975;Duncan 1996).Since the landslide represents a failed slope i.e., critical shear strength > existing shear strength, FS<1 is justifiable.Further, the failure state of a slope in the FEM can be defined by different criteria; the FS of the same slope may vary a little depending upon the usage of failure criteria and the convergence threshold (Abramson et al. 1996;Griffiths and Lane 1999).
The possible causes of instability (FS<1) may be steep slope gradient, rockmass having low strength, and joints.Three (S.N. 7, 14, 15) out of the five landslides that have potential to form dams belong to the tectonically active Higher Himalaya Crystalline (HHC).The notion of steep slope gradient cannot be generalized because the HHC accommodates voluminous (~10 5 -10 7 m 3 ) landslides (Fig. 4).These deep seated landslides must require smooth slope gradient to accommodate the voluminous overburden.Further, the HHC comprises gneiss having high compressive strength and Geological Strength Index (Supplementary Table 2, Kumar et al. 2021), therefore the notion of low strength rockmass also may not be appropriate.However, the jointed rock mass that owes its origin to numerous small-scale folds, shearing, and faults associated with the active orogeny process can be considered as the main factor for relatively more instability of debris slide type landslides.Since, the study area is subjected to the varied stress regime caused by the tectonic structures (Vannay et al. 2004), thermal variations (Singh et al. 2015), and anthropogenic cause (Lata et al. 2015), joints may continue to develop and destabilize the slopes.Apart from this inherent factor like joints, external factors like rainfall and exhumation rate may also contribute to instability of these landslides.This region receives relatively more annual rainfall owing to orographic barrier setting (Fig. 11) and is subjected to relatively high exhumation rate of 2.0-4.5 mm/yr (Thiede et al. 2009).
Two landslides (S.N. 5, 19) that are also capable of forming potential landslide dams (Fig. 8,9a;e) and are also unstable (FS<1) in nature (Fig. 4) do not belong to the HHC.The first landslide (S.N. 5) exists at the lithological contact of schist of the Tethyan Sequence and Kinnaur Kailash Granite rockmass.A regional normal fault, the Sangla Detachment (SD), passes through this contact.Some prior studies suggest that the SD is an outcome of reactivation of a former thrust fault that has resulted in intense rockmass shearing (Vannay et al. 2004;Kumar et al. 2019b).Owing to its location in the orographic interior region, hillslopes receive very low annual rainfall (Fig. 11) and thus have the least vegetation on the hillslopes in this region.The lack of vegetation on hillslopes has been observed to result in low shear strength of material and hence in the instability (Kokutse et al. 2016).Thus, lithological contrast, rockmass shearing, and lack of vegetation are thought to be the main reasons of instability of S.N. 5 landslide.The second landslide (S.N. 19) belongs to the inter-layered schist/gneiss rockmass of the Lesser Himalaya Crystalline (LHC) and is situated at the orographic front where rainfall increases suddenly (Fig. 11).Further, this region is also subjected to the high exhumation rate of 2.0-4.5 mm/yr (Thiede et al. 2009).Therefore, lithological contrast, high rainfall and high exhumation rate are considered as the main reasons of instability of this landslide slope.
The landslides that are inferred not to result in the river damming are mostly in the LHC and Lesser Himalaya Sequence (LHS) region.These regions consist of a majority of the rock fall and rock avalanches that are generally of meta-stable category (Fig. 4).Despite the narrow valley setting, landslides in these regions may not form the potential landslide dam, at present, owing to the relatively small landslide volume.The possible causes of their meta-stability may be high compressive strength and geological strength index of gneiss (Kumar et al. 2021), dense vegetation on the hillslopes (Chawla et al. 2012), relatively less sheared rock mass in comparison to the HHC region, and relatively less decrease in land use/landcover (Lata et al. 2015).Maximum Total Displacement (TD) is also associated with the rock fall and rock avalanche of this region (Fig. 4).
In the parametric study, soil friction and in-situ stress are noted to affect the FS most in case of the debris slide, whereas the FS of rock fall and rock avalanche are mainly controlled by the 'mi' and the in-situ stress.The 'mi' is a GHB criteria parameter that is equivalent to the friction in the M-C criteria.For the TD of the debris slides, field stress, elastic modulus and Poisson's ratio, whereas for rock falls and rock avalanches, the 'mi' parameter and in-situ stress played the dominant role (Fig. 6,7).Friction has been a controlling factor for the shear strength, and its decrease has been observed to result in the shear failure of slope material (Matsui and San 1992).Since rainfall plays an important role in decreasing the friction of slope material by changing the pore water pressure regime (Rahardjo et al. 2005), frequent extreme rainfall events in the study area since the year 2013 (Kumar et al. 2019a) amplify the risk of hillslope instability.Furthermore, the in-situ field stress that has been compressional and/or extensional owing to the orogenic setting in the region may also enhance the hillslope instability (Eberhardt et al. 2004;Vannay et al. 2004).Deformation parameters, e.g.elastic modulus and Poisson's ratio, are also observed to affect the displacement in slope models of the debris slides.Similar studies in other regions have also noted the sensitivity of the elastic modulus and Poisson's ratio on the slope stability (Zhang and Chen 2006).
The study area has been subjected to extreme rainfalls since the year 2010 and received widespread slope failures and flash-floods (Fig. 11b).Three (S.N. 7,14,15 in Fig. 9) out of five potential landslide dams belong to the Higher Himalaya Crystalline (HHC) that receives relatively more rainfall (Fig. 11).Contrary to the along 'Himalayan' arc distribution of earthquakes, the study area has received most of the earthquakes around the N-S oriented Kaurik-Chango Fault (Fig. 12a).However, the only major earthquake event has been Mw 6.8 earthquake on 19 th Jan. 1975 that resulted in the widespread landslides (Khattri et al. 1978).
The low-magnitude recent seismicity in the region has been attributed to the northward extension of the Delhi-Haridwar ridge (Hazarika et al. 2019), whereas the shallow nature is attributed to the MHT ramp structure in the region that allows strain accumulation at shallow depth (Bilham 2019).Thus, earthquakes have not been a major landslide triggering process in the region in recent times.Finally, the word "active landslide" refers to the hillslope that is still subjected to the slope failures caused by the various factors.The word "landslide" can be perceived in the following three ways; pre-failure deformations, failure itself, and post-failure displacement (Terzaghi 1950;Cruden & Varnes, 1996;Hungr et al., 2014).Landslide slopes in this study pertain to the post-failure state that are categorized into "unstable" and "metastable" stages based on their existing FS.If an active landslide is not categorized as "unstable", it means that the existing slope geometry provides it a "meta-stable" stage that might transform into an unstable stage with time due to the stability controlling parameters (Sec.4.1).Though the field visits were performed in different seasons to cover all the landslides along the study area, there might be a possibility of vegetation growth on the failed slopes, particularly in the LHC and LHS.However, the landslides in the LHC and LHS are mostly rockfall/rock avalanche type because of the deep gorge setting, whereas the vegetation growth generally requires the debris laden hillslopes.Nonetheless, such aspect will be explored in the future prospects.The HHC and the TS region are subjected to the semi-humid to semi-arid climate and hence the vegetation type is mostly scattered trees/shrubs.Therefore, such possibility might not exist in these regions.
A supplementary table involving all the details like landslides dimensions, factor of safety, and geomorphic indices output of each landslide is provided in the data repository (Kumar et al. 2021).
In view of the possible uncertainties in the predictive nature of the study, the following assumptions and simplifications were made; • To account the effect the spatial variability in the slope geometry, 3D models have been in use for the last decade (Griffiths and Marquez 2007).However, the pre-requisite for the 3D models involves the detailed understanding of slope geometry and material variability in the subsurface that was not possible in the study area considering steep and inaccessible slopes.Therefore, multiple 2D sections were chosen, wherever possible.To account the effect of sampling bias and material variability, a range of values of input parameters was used (sec.4.1).
• Determination of the debris thickness has been a major problem in the landslide volume measurement particularly in the steep, narrow river valleys of the NW Himalaya.
Therefore, the thickness was approximated by considering the relative altitude of the ground on either side of the deposit, as also performed by Innes (1983).It was assumed that the ground beneath the deposit is regular.
• The resultant dam volume could be different from the landslide volume due to the entrainment, rockmass fragmentation, pore water pressure, size of debris particles, and washout of landslide material by the river (Hungr and Evans 2004;Dong et al. 2011;Yu et al. 2014).Therefore, dam volume is presumed to be equal to landslide volume for the worst-case scenario (sec.3.3).Stream power is manifested by the upstream catchment area and local slope gradient in the geomorphic indices.It may also vary at temporal scale owing to the temporally varying water influx from glaciers and precipitation systems i.e., ISM and WD (Gadgil et al. 2007;Hunt et al. 2018).Though our study is confined to the spatial scale at present, the findings remain subjected to the change at temporal scale.
• The RAMMS model (Voellmy 1955;Salm 1993;Christen et al. 2010) requires the calibrated friction and turbulence values for the run-out analysis.Though the previous debris flow events have not left evidence in the study area owing to the convergence of the landslide toes with the river channel, a range of µ and ξ values were used in the study in view of the material type and run-out path characteristics.
Despite these uncertainties, studies such as this one are required to minimize the risk and avert the possible disasters in the terrain where human population lives in the proximity of unstable landslides.

CONCLUSION
Out of forty-four landslides that are studied in the Satluj valley in the NW Himalaya, five landslides are noted to form the potential landslide dam, if failure occurs.Though the blocking duration is difficult to predict, upstream and downstream consequences of these damming  Fig. 3 The FEM configuration of some of the slope models.S.N. refers to the serial no. of landslides in Table 1.The joint distribution in all the slopes was parallel-statistical with the normal distribution of joint spacing.Barauni Gad_I_S (S.N. 38).S. N. refers to the serial no. of landslides in Table 1.          3 The FEM configuration of some of the slope models.S.N. refers to the serial no. of landslides in Table 1.The joint distribution in all the slopes was parallel-statistical with the normal distribution of joint spacing.

Fig. 4
The FEM analysis of all forty-four landslides.Grey bar in the background highlights the Higher Himalaya Crystalline (HHC) region that comprises relatively more unstable landslides, landslide volume and human population..TS, KKG, HHC, LHC and LHS are Tethyan Sequence, Kinnaur Kailash Granite, Higher Himalaya Crystalline, Lesser Himalaya Crystalline and Lesser Himalaya Sequence, respectively.Barauni Gad_I_S (S.N. 38).S. N. refers to the serial no. of landslides in Table 1.
and ultrasonic tests (CATS Ultrasonic (1.95) of Geotechnical Consulting & Testing Systems).The ultrasonic test was conducted to determine the density, elastic modulus, and Poisson's ratio of rock samples.The soil samples were tested for grain size (IS: 2720-Part 4-1985), UCS test (IS: 2720-Part 10-1991), and direct shear test (IS: 2720-Part 13 (dam volume)= Vl (landslide volume), m 3 ; Ab is upstream catchment area (km 2 ); Wv is width of the valley (m) and S is local slope gradient of river channel (m/m).Though the resultant dam volume could be higher or lower than the landslide volume owing to slope entrainment, rockmass fragmentation, retaining of material at the slope, and washout by the river(Hungr and Evans 2004;Dong et al. 2011), dam volume is assumed to be equal to landslide volume for the worst case.By utilizing the comprehensive dataset of ~300 landslide dams of Italy,Stefanelli et al. (2016) have classified the MOI into (i) non-formation domain: MOI <3.00, (ii) uncertain evolution domain: 3.00 <MOI >4.60, and (iii) formation domain: MOI >4.60.By utilizing the same dataset,Stefanelli et al. (2016) defined the HDSI into following categories (i) instability domain: HDSI <5.74, (ii) uncertain determination domain: 5.74<HDSI >7.44, and (iii) Stability domain: HDSI>7.44.
June 2013 that resulted in the widespread slope failure in the NW Himalaya also caused landslide damming at places (National Disaster Management Authority, Govt. of India, 2013; Kumar et al. 2019a).Similar to the year 2013, the years 2007, 2010 and 2019 also witnessed enhanced annual rainfall and associated flash floods and/or landslides in the region (hpenvis.nic.in,retrieved on March 1, 2020; sandrp.in,retrieved on March 1, 2020).

Fig. 4
Fig. 4The FEM analysis of all forty-four landslides.Grey bar in the background highlights the Higher Himalaya Crystalline (HHC) region that comprises relatively more unstable landslides, landslide volume and human population..TS, KKG, HHC, LHC and LHS are Tethyan Sequence, Kinnaur Kailash Granite, Higher Himalaya Crystalline, Lesser Himalaya Crystalline and Lesser Himalaya Sequence, respectively.

Fig. 5 Fig. 6
Fig. 5 Relationship of Factor of Safety (FS), Total Displacement (TD) and Shear Strain (SS).DS, RF, and RA refer to Debris slide, rock fall and rock avalanche, respectively.

Fig. 10
Fig. 10 Field signatures of the landslide damming near Akpa_III landslide.(a) Upstream view of Akpa landslide with lacustrine deposit at the left bank; (b) enlarged view of the lacustrine deposit with an arrow indicating the lacustrine sequence; (c) alternating finecoarse sediments.F and C refer to fine (covered by yellow dashed lines) and coarse (covered by green dashed lines) sediments, respectively.

Fig. 12
Fig. 12 Earthquake distribution.(a) Spatial variation of earthquakes.The transparent circle represents the region within 100 km radius from the Satluj River (blue line).The black dashed line represents the seismic dominance around the Kaurik-Chango fault;(b) earthquake magnitude vs. focal depth.The red dashed region highlights the concentration of earthquakes within 40 km depth; (c) Cross section view (Based on Hazarika et al. 2017; Bilham 2019).Red dashed circle represents the zone of strain accumulation caused by the Indian and Eurasian plate collision (Bilham 2019).ISC: International Seismological Centre.HFT: Himalayan Frontal Thrust.

Fig. 1
Fig. 1 Geological setting.WGC: Wangtu Gneissic Complex.The red dashed circle in the inset represents the region within 100 km radius from the Satluj River (marked as blue line) thatwas used to determine the earthquake distribution in the area.Yellow lines represent the regional faults in the region.KCF in inset refers to Kaurik-Chango Fault.The numbers 1-44 refer to serial number of landslides in Table1.

Fig.
Fig.3The FEM configuration of some of the slope models.S.N. refers to the serial no. of

Fig. 5
Fig. 5 Relationship of Factor of Safety (FS), Total Displacement (TD) and Shear Strain (SS).DS, RF, and RA refer to Debris slide, rock fall and rock avalanche, respectively.

Fig. 10
Fig. 10 Field signatures of the landslide damming near Akpa_III landslide.(a) Upstream view of Akpa landslide with lacustrine deposit at the left bank; (b) enlarged view of the lacustrine deposit with an arrow indicating the lacustrine sequence; (c) alternating fine-coarse sediments.

F
and C refer to fine (covered by yellow dashed lines) and coarse (covered by green dashed lines) sediments, respectively.

Fig. 12
Fig. 12 Earthquake distribution.(a) Spatial variation of earthquakes.The transparent circle represents the region within 100 km radius from the Satluj River (blue line).The black dashed line represents the seismic dominance around the Kaurik-Chango fault;(b) earthquake magnitude vs. focal depth.The red dashed region highlights the concentration of earthquakes within 40 km depth; (c) Cross section view (Based onHazarika et al. 2017;Bilham, 2019).Red dashed circle represents the zone of strain accumulation caused by the Indian and Eurasian plate collision(Bilham, 2019).ISC: International Seismological Centre.HFT: Himalayan Frontal Thrust.

Fig. 13
Fig. 13 Results of the run-out analysis.µ refers to coefficient of friction.

Fig. 14
Fig. 14 Results of run-out analysis at different values of µ and ξ. µ and ξ refer to coefficient of friction and turbulence, respectively.
Wulf, H.,Bookhagen, B. and Scherler, D. 2012."Climatic and geologic controls on suspended sediment flux in the Sutlej River Valley, western Himalaya".Hydrology and Earth Geological setting.WGC: Wangtu Gneissic Complex.The red dashed circle in the inset represents the region within 100 km radius from the Satluj River (marked as blue line) that was used to determine the earthquake distribution in the area.Yellow lines represent the regional faults in the region.KCF in inset refers to Kaurik-Chango Fault.The numbers 1-44 refer to serial number of landslides in Table1.
Fig. 14 Results of run-out analysis at different values of µ and ξ. µ and ξ refer to coefficient of Details of the landslides used in the study.Details of the satellite imagery.Criteria used in the Finite Element Method (FEM) analysis.Details of input parameters used in the run-out analysis.Details of landslides used in the study.