remote sensing

Article
Reconstruction and Characterisation of Past and the Most
Recent Slope Failure Events at the 2021 Rock-Ice Avalanche Site
in Chamoli, Indian Himalaya
Anshuman Bhardwaj * and Lydia Sam

School of Geosciences, University of Aberdeen, Meston Building, King’s College, Aberdeen AB24 3UE, UK;
lydia.sam@abdn.ac.uk
* Correspondence: anshuman.bhardwaj@abdn.ac.uk

Abstract: Frequent ice avalanche events are being reported across the globe in recent years. On
the 7 February 2021, a flash flood triggered by a rock-ice avalanche with an unusually long runout
distance, caused significant damage of life and property in the Tapovan region of the Indian Himalaya.
Using multi-temporal satellite datasets, digital terrain models (DTMs) and simulations, here we
report the pre-event and during-event flow characteristics of two large-scale avalanches within a
5-year interval at the slope failure site. Prior to both the events, we observed short-term and long-
term changes in surface velocity (SV) with maximum SVs increasing up to over 5 times the normal
values. We further simulated the events to understand their mechanical characteristics leading to
long runouts. In addition to its massive volume, the extraordinary magnitude of the 2021 event can
partly be attributed to the possible remobilisation and entrainment of the colluvial deposits from
previous ice and snow avalanches. The anomalous SVs should be explored further for their suitability
as a possible remotely observable precursor of ice avalanches from hanging glaciers. This sequence
of events highlights that there is a need to take into account the antecedent conditions, while making





 a holistic assessment of the hazard.
Citation: Bhardwaj, A.; Sam, L.
Reconstruction and Characterisation Keywords: ice avalanche; Himalaya; Chamoli; debris flow; rapid mass movement simulation (RAMMS);
of Past and the Most Recent Slope co-registration of optically sensed images and correlation (COSI-Corr); glacier velocity; high mountain;
Failure Events at the 2021 Rock-Ice disaster; hazard sequence; mass movement mechanics
Avalanche Site in Chamoli, Indian
Himalaya. Remote Sens. 2022, 14, 949.
https://doi.org/10.3390/rs14040949

Academic Editor: Ulrich Kamp 1. Introduction

Received: 2 January 2022

Slope failure events in high-mountains globally hamper infrastructure and economic
Accepted: 10 February 2022
development, and cause hundreds of fatalities each year [1–3]. Contemporary research [4–7]
Published: 16 February 2022
linked the increasing mass movement hazards in high mountains to the increase in temper-
atures and erratic precipitation conditions. However, many of these studies also attribute
Publisher’s Note: MDPI stays neutral
the cascading effects of rapidly changing periglacial landscape and local topography for
with regard to jurisdictional claims in
these mountain disasters [3–6]. Therefore, it is relevant to have a basic understanding of
published maps and institutional affil-
various slope failure processes in high mountains [8]. Mass movements in high moun-
iations.
tains can be characterised based on three main parameters [9]: (1) constituents involved
(regolith, rock, glacial ice, and snow), (2) propagation velocity (slow, intermediate, and
fast), and (3) characteristics of the movement (debris cloud, slurry, and coherent mass). An
Copyright: © 2022 by the authors.
understanding of these characterising parameters further allows for a distinct classification
Licensee MDPI, Basel, Switzerland. of high-mountain slope failure events within categories, such as mudflow, debris flow,
This article is an open access article avalanche, landslide, rockfall, and debris fall. Falls are abrupt movements of rocks or
distributed under the terms and debris after detaching from steep slopes along structural discontinuities, such as joints,
conditions of the Creative Commons faults, bedding planes, and other fractures, and movement can occur in form of free-fall,
Attribution (CC BY) license (https:// bouncing, and rolling. A mudflow is a very rapid movement of saturated, fine-grained
creativecommons.org/licenses/by/ deposits with water contents equal to or higher than the liquid limit, changing the phase of
4.0/). the deposits from solid to plastic phase [10]. Debris flows are also very rapid movements

Remote Sens. 2022, 14, 949. https://doi.org/10.3390/rs14040949 https://www.mdpi.com/journal/remotesensing

Remote Sens. 2022, 14, 949 2 of 32

of rock debris, in which the constituents can be in the form of tens to even hundreds of
meters-wide boulders [11]. A debris flow usually forms deposits at the mouth of a canyon
in form of an alluvial or debris cone consisting of a poorly sorted mixture of particles, and
these deposits are steeper than the usual alluvial fans [12]. Avalanches are characterised by
turbulent clouds of debris mixed with air that rush down steep slopes [3]. An avalanche
can be named as an ice, snow, and rock avalanche based on its major constituent. While wet
avalanches rich in liquid water behave like a viscous slurry and move relatively slowly, dry
avalanches consisting of cold powdery snow move rapidly (~60 m/s) above the ground
surface on a layer of pressurised air [13]. On the other hand, debris avalanches, consisting
primarily of rocks and ice slabs, can move with velocities reaching even in excess of 80 m/s,
depending on the slope angle [14]. Landslides are defined by rockslides and debris slides.
While rockslides move along a plane parallel to the surface as a coherent mass, debris slides
display rotational sliding typically in the upper part. Debris flows differ from landslides
as debris flows propagate within numerous layers of the medium, whereas a landslide, as
the name implies, occurs only along one or several interfaces or beds [15]. While disaster
chains are initiated by landslides and avalanches, debris flow occur along a valley bottom
typically after the debris deposition from an avalanche or a landslide. The deposits from
landslides and avalanches are known to form barrier lakes in high mountains, developing
the risk of future outburst floods.
The Hindu Kush Himalaya (HKH) is a hotspot for land degradation and mass move-
ment events [7], such as snow avalanches, ice avalanches, and debris flows, due to the high
seismicity along the Indian and Eurasian plate boundary and the extreme precipitation
conditions [16]. The HKH countries are some of the most densely populated ones globally
with limited land resources, and this makes the commercial activities (e.g., [17]) in the
high mountains prone to mass movement hazards, both unavoidable and unsustainable
at the same time. Snow and ice avalanches and landslides are frequent occurrences in
high Himalayan mountains [18,19]. Nevertheless, the hazardous nature is defined by
their magnitude and proximity with the local population [3]. Ice avalanches are typical
of hanging glaciers. However, in the past couple of decades, a remarkable frequency of
sudden large-volume detachments even from low-angle mountain glaciers has also been
reported [20]. Ice avalanches have three main modes of release, i.e., frontal block failure,
ice slab failure, and ice-bedrock failure [3]. Even minor ice avalanches can lead to a se-
rious secondary hazard risk [3] and the Himalayan regions are mainly affected by these
secondary hazard risks through the formation of dangerous ice-dammed lakes or through
glacial lake outburst floods (GLOFs). Another secondary hazard risk is the entrainment of
the snowpack in winter generating a snow avalanche; even a minor ice avalanche could
trigger a huge snow avalanche. However, the studied event, i.e., the flash flood in Dhauli-
ganga valley in Chamoli District of Uttarakhand State of India, on the morning of the
7 February 2021 was caused by a rock-ice avalanche with an extraordinarily long runout of
over ~12 km, before it turned into an extremely fluidised debris flow in Rishiganga River.
Within minutes, the Rishiganga River was gushing with water and debris damaging the
13.2 MW Rishiganga hydropower plant, while the debris flow further continued in the
Dhauliganga River and damaged the under-construction 520 MW Tapovan–Vishnugad
hydropower project, causing fatalities and injuries to over 200 people. The scale of the
destruction could have magnified, had it not happened in winter, when the water-level
in the river was low. A study [21] provided a detailed overview of the 7 February 2021
avalanche hazard cascade and modelled this event to be consisting of 80% rock and 20% ice
volume, generating a total avalanche volume of ~27 × 106 m3 . Two more studies [22,23]
also adopted similar approaches to put a volumetric constraint on the 7 February 2021 event
and described the possible chain of events that led to the extensive destruction. Another
recent study [24] used satellite datasets to investigate and comprehend the possible causes
for progressive destabilisation at the ice-rock interface.
In this paper, we intend to supplement these published studies (e.g., [21–24]) by
focusing on the pre-event and during-event cascading flow characteristics using multi-
Remote Sens. 2022, 14, 949 3 of 32

temporal satellite datasets and robust rock-ice avalanche modelling approaches. We focus
on quantifying and discussing the mechanical characteristics (momentum, pressure, and
shear stress) of this mass movement leading to its translation into an enormous debris flow.
We selected two massive avalanche events, temporally spaced by a ~4.5-year gap, which
ran through the same valley stretch, but were constituently different. The first event of
September 2016, hereafter referred to as “the first slope failure (SF1)”, was a frontal ice-block
failure of a hanging glacier. The second event of the 7 February 2021, hereafter referred
to as “the second slope failure (SF2)”, was a rock-ice avalanche and was over ~2.5 times
more voluminous than the 2016 ice avalanche. The successive occurrence of these two
major events in the same topographical setting provides us an excellent opportunity to
understand the difference in their mass movement behaviour owing to varying degree
of ice content within them. The manifestation of these two massive slope failure events,
SF1 and SF2, originating from nearly the same release height and hitting the same valley
within a period of 5 years, is unique and offers us an unprecedented natural testbed to
understand how frequent mass movements in the coming years can vary in terms of their
runout and destructive nature. A thermomechanical simulation of SF1 offers us the ability
to modify the valley topography to include deposits from SF1 in analysing their possible
influence on the debris flow runout of SF2 [25]. In the following sections, we provide a
brief overview of the avalanche sites, and separately present our findings on the estimated
surface velocity (SV) and the modelled accounts of the two slope failures. Considering that
SF1 was a frontal ice-block failure, while SF2 was a landslide causing ice-bedrock failure, in
this paper, we used SV as an all-inclusive term to account for all the movements (sliding
and deformational) occurring at the avalanche sites. Thus, this paper investigates three
main research questions: (1) are there any quantifiable pre-event changes in SVs? (2) how
do SV changes vary between two constituently different slope failures?, and (3) what are
the differences in mechanical characteristics of two large slope failures originating from the
same release height and hitting the same valley?

2. Study Area
The study area is situated in Central Himalaya within a transitional zone between the
wetter Southern Himalayan ranges dominated by Indian Summer Monsoons (ISMs) and
the cold arid Tibetan Highlands in the north [26]. Although the predominant precipitation
source is ISM during the warmer months, the region also receives considerable precipitation
from mid-latitude westerlies in the form of winter snowfall [27] between December and
March. The average annual precipitation is ~1100 mm with sub-zero winter temperatures
and summer temperature oscillating between 5 ◦ C and 29 ◦ C [28]. These weather systems
display a significant control over regional cryosphere, ecosystems, and landscape evolution,
ultimately influencing the socioeconomics of the region [26]. The rugged topography
of this region lies between ~2000 and ~7800 m asl [29]. Topographically, the region is a
hotbed for high-velocity mass movement events as it is characterised by high peaks, narrow
river valleys, and deep gorges [26]. The study area partly falls within the Alaknanda
River catchment and is abundant in several mountain streams and tributaries fed by
meltwaters [26].
Geologically, the study area is located between the Vaikrita Thrust to the south and
the Malari Fault to the North [26], and consists of the Vaikrita Group of high-grade meta-
morphic rocks with both acidic and basic intrusions, which are well-exposed along the
Rishiganga River valley [30]. These metamorphic rocks are dominated by schist, amphi-
bolite, and gneisses [31], and such a stratified geology is considered prone to episodes
of glaciation, deglaciation, and neotectonics [32,33]. Thus, the geology of the study area
is structurally fragile and exposed to the contrasting regional climatic conditions [26].
The elevation range of the detached SF1 ice slab was ~5200–5400 m asl while that of the
avalanched SF2 rock-ice body was ~4900–5600 m asl. Interestingly, the equilibrium line
altitude (ELA), i.e., the dividing altitude between the accumulation and ablation zones of a
glacier, reported for glaciers in this region is at ~5400 m asl [28]. At ELA, the annual snow
 Remote Sens. 2022, 14, x FOR PEER REVIEW 4 of 33

Remote Sens. 2022, 14, 949 lanched SF2 rock-ice body was ~4900–5600 m asl. Interestingly, the equilibrium line 4 ofalti-
32
tude (ELA), i.e., the dividing altitude between the accumulation and ablation zones of a
glacier, reported for glaciers in this region is at ~5400 m asl [28]. At ELA, the annual snow
accumulationequals
accumulation equalsthe theannual
annualsnow
snowablation
ablationover
overaaperiod
periodofofseveral
severalyears,
years,making
makingthis this
altitude very sensitive to climate fluctuations [28]. The fact that both these
altitude very sensitive to climate fluctuations [28]. The fact that both these events occurredevents occurred
atatthe
theslopes
slopesclose
closetotothe
thereported
reportedELA,ELA,indicates
indicatestheir
theirparticular
particularsensitivity
sensitivitytowards
towardsbothboth
winteraccumulation
winter accumulationand andsummer
summerablation.
ablation.Further
Furtherdetails
detailson
onthe
thegeomorphic
geomorphicand anddimen-
dimen-
sionalcharacteristics
sional characteristicsof ofthese
theseslope
slopefailures
failuresare
areprovided
providedininthe theResults
Resultsand
andDiscussion
Discussion
sectionfor
section forbetter
bettercontinuity
continuityand andunderstanding
understandingofofthe theevents.
events.TheTheavalanched
avalancheddebris
debrishit hit
the valley bottom of Ronti Gad stream, which merges in the Rishiganga
the valley bottom of Ronti Gad stream, which merges in the Rishiganga River before Raini River before Raini
Village.Rishiganga
Village. Rishigangacontributes
contributestotothe theDhauliganga
DhauligangaRiverRiverafter
aftercrossing
crossingRaini
RainiVillage.
Village.SF2SF2
followedthis
followed thisroute
routeforfor~32
~32km kmasasaaflash
flashflood
flood[29],
[29],causing
causingaawidespread
widespreaddestruction
destructionofof
lifeand
life andproperty
propertyandandvisibly
visiblychanging
changingthe thegeomorphology
geomorphologyof ofthe
thevalley.
valley.AAshort
shortYouTube
YouTube
video[34]
video [34]captured
capturedthe theenormity
enormityand andthethedestructive
destructivenature
natureofofthe
theflash
flashflood
floodtriggered
triggered
bythe
by theSF2
SF2event
eventupuptoto~8 ~8km kmdownstream
downstreamfrom fromthethedeposition
depositionsite.
site.Figure
Figure11provides
providesaa
detailedoverview
detailed overviewofofthe thestudy
studyarea.
area.

Figure1.1.Location
Figure Locationmap:
map: (a)
(a)satellite
satelliteimage
imagehighlighting
highlighting the
thestudy
studyarea
areaand
and(b)
(b)an
anoverview
overviewmapmap
showingthe
showing thecontextual
contextuallocation
locationofofthe
thestudy
studysite
site(red
(redrectangle)
rectangle)ininnorth
northIndia
India(the
(themain
mainriver
riverbasins
basins
ofofthe
theIndian
IndianSubcontinent
Subcontinentare
arealso
alsonamed
namedininblue
bluetext,
text,while
whilethe
thecountries
countriesand
andtheir
theircapitals
capitals(green
(green
stars) are mentioned in black and red texts, respectively); (c) zoomed-in 3D perspective views (north
is downward) of the two slope failure sites (the red arrow shows the snow in the upper reaches of the
slopes joining SF1 and SF2 zones and the dashed red curve shows the indentation from where repetitive
Remote Sens. 2022, 14, 949 5 of 32

slope failure events take place at SF1 site); (d) 3D perspective views of the starting zone (yellow
rectangle in the top panel) of SF1 before and after the event, and bottom panel showing the de-
position zone (orange rectangle) of SF1 in the Ronti Gad stream (the 9 October 2016 Sentinel-2 ID:
S2A_OPER_MSI_L1C_TL_SGS__20161009T053631_A006780_T44RLU); and (e) 3D perspective views
of the starting zone (yellow rectangle) of SF2 on 7 February 2021 (the 13 September 2020 Sentinel-
2 ID: S2B_MSIL1C_20200913T051649_N0209_R062_T44RLU, and the 10 February 2021 Sentinel-2
ID: S2B_MSIL1C_20210210T051939_N0209_R062_T44RLU). Sentinel-2 images used in (c–e) are ac-
quired within the Copernicus Programme of the European Space Agency (ESA); the source for (b) is
world grey scale layer of Environmental Systems Research Institute (ESRI) online base map; the
source of the HKH river basin boundary used in (b) is the International Centre for Integrated Moun-
tain Development [35]; and the bottom panel of (d) comprises of the Google Earth images (Credit:
CNES/Airbus).

3. Materials and Methods

In addition to the high-resolution images obtained from Google Earth for visual
analysis, we mainly used two remote sensing datasets for performing quantitative analyses:
(1) multi-temporal satellite images for SV estimation and (2) high-resolution multi-temporal
digital terrain models (DTMs) for Rapid Mass Movement Simulation (RAMMS) modelling.
The details of these datasets are provided in Table 1. The image IDs of the satellite images
used for visual analysis are provided in the respective figure captions. The overall workflow,
designed to investigate the pre-event and during-event flow characteristics, is depicted in
Figure 2.

Table 1. Data used for SV estimations and RAMMS modelling.

Satellite Images
Acquisition Date Scene ID/Band Satellite/Sensor Spatial Resolution (m/pixel)
1 September 2015 LC81450392015244LGN01/Panchromatic Landsat 8/OLI 15
17 September 2015 LC81450392015260LGN01/Panchromatic Landsat 8/OLI 15
3 September 2016 LC81450392016247LGN01/Panchromatic Landsat 8/OLI 15
19 September 2016 LC81450392016263LGN01/Panchromatic Landsat 8/OLI 15
9 October 2016 S2A_OPER_MSI_L1C_TL_SGS__20161009T053631_A006780_T44RLU/B2 Sentinel-2/MSI 10
14 October 2017 S2A_OPER_MSI_L1C_TL_SGS__20171014T104205_A012071_T44RLU/B2 Sentinel-2/MSI 10
24 October 2018 S2B_OPER_MSI_L1C_TL_MTI__20181024T090955_A008525_T44RLU/B2 Sentinel-2/MSI 10
24 October 2019 S2A_OPER_MSI_L1C_TL_EPAE_20191024T082325_A022653_T44RLU/B2 Sentinel-2/MSI 10
18 October 2020 S2A_OPER_MSI_L1C_TL_VGS1_20201018T072817_A027801_T44RLU/B2 Sentinel-2/MSI 10
27 December 2020 S2A_OPER_MSI_L1C_TL_EPAE_20201227T062601_A028802_T44RLU/B2 Sentinel-2/MSI 10
16 January 2021 S2A_OPER_MSI_L1C_TL_EPAE_20210116T063033_A029088_T44RLU/B2 Sentinel-2/MSI 10
26 January 2021 S2A_OPER_MSI_L2A_TL_EPAE_20210126T074559_A029231_T44RLU/B2 Sentinel-2/MSI 10
5 February 2021 S2A_OPER_MSI_L2A_TL_VGS1_20210205T082130_A029374_T44RLU/B2 Sentinel-2/MSI 10
DTM Details
Date Range for Source
Satellite Spatial Resolution (m/pixel)
Stereopairs Link/Reference
1 September 2015–
WorldView-1 and WorldView-2 [21,36–38] 2
5 October 2015
10–11 February 2021 WorldView-2, WorldView-3 and GeoEye-1 [21,37,39,40] 2

3.1. Surface Velocity (SV) Change Estimation

We used the Co-Registration of Optically Sensed Images and Correlation (COSI-Corr)
tool [41], which works as a plug-in module to the Environment for Visualizing Images
(ENVI) software for estimating horizontal displacement fields. COSI-Corr implements a
phase correlation method based on the Fourier Shift Theorem and the retrieved relative
displacement between an image pair is a resultant of an iterative, unbiased process, esti-
mating the phase difference of their Fourier transform [41]. The standard outputs are in the
form of East–West (EW) displacement, North–South (NS) displacement, and signal-to-noise
ratio (SNR), and the absolute displacement is calculated as the square root of the sum of
squares of the EW and NS displacement components. This displacement can be divided by
various units of time, i.e., day in the present study, to obtain the SVs. The typical correlation
accuracy achieved at sub-pixel resolution by COSI-Corr is as high as 1/20–1/10 of the
spatial resolution of an image pixel [42].
 ysis, we mainly used two remote sensing datasets for performing quantitative analyses:
(1) multi-temporal satellite images for SV estimation and (2) high-resolution multi-tem-
poral digital terrain models (DTMs) for Rapid Mass Movement Simulation (RAMMS)
modelling. The details of these datasets are provided in Table 1. The image IDs of the
satellite images used for visual analysis are provided in the respective figure captions. The
Remote Sens. 2022, 14, 949 6 of 32
overall workflow, designed to investigate the pre-event and during-event flow character-
istics, is depicted in Figure 2.

Figure2.2.AAworkflow
Figure workflowfor
forinvestigating
investigatingthe
thepre-event
pre-eventand
andduring-event
during-eventflow
flowcharacteristics.
characteristics.

The panchromatic band (Band 8) of Landsat 8 Operational Land Imager (OLI) and
blue band (B2) of Sentinel-2 MultiSpectral Instrument (MSI) are the preferred bands for
obtaining glacier and mass movement SVs [43,44]. Landsat 8 images were downloaded
from the United States Geological Survey’s (USGS’s) EarthExplorer portal [45] and Sentinel-
2 images were downloaded from the Copernicus Open Access Hub [46]. A 15 m/pixel
resolution Landsat 8 Panchromatic band was selected for the SV estimation of SF1, as
appropriate Sentinel-2 data was not available during September 2015. We found cloud-free
Landsat 8 images covering 1–17 September 2015 and 3–19 of September 2016, with identical
surface conditions for comparing the same duration SVs in the two consecutive years prior
to SF1. We did not go back farther than one year before SF1 as it was primarily a frontal
block failure and the critical ice load conditions for such events vary significantly from one
year to another in this region [47]. The selected image pairs clearly highlighted the effect of
this aspect on the estimated SVs. In fact, SF1 is a recurring frontal ice-block failure and we
identified a previous similar event in early 2000 (Figure 3). In contrast, SF2 being primarily a
landslide causing an ice-bedrock failure, we opted to perform yearly SV observations for it,
starting from 2016 with good temporal availability of Sentinel-2 images. We again ensured
cloud-free and similar surface conditions by selecting all the images from October, and this
approach significantly reduced uncertainties. For SV observations ~40 days prior to SF2,
we again opted for Band 2 of Sentinel-2. The presence of seasonal snow was unavoidable in
this case. In order to effectively deal with mountain shadow, which was covering the entire
SF2 hanging glacier uniformly during these ~40 days, we first clipped and then contrast-
stretched the SF2 site in all the scenes by comparative histogram matching. Afterwards,
we performed COSI-Corr processing only for this contrast-stretched part and the obtained
results were considerably good with acceptable uncertainties (Table 2).
 vations ~40 days prior to SF2, we again opted for Band 2 of Sentinel-2. The presence of
seasonal snow was unavoidable in this case. In order to effectively deal with mountain
shadow, which was covering the entire SF2 hanging glacier uniformly during these ~40
days, we first clipped and then contrast-stretched the SF2 site in all the scenes by compar-
Remote Sens. 2022, 14, 949
ative histogram matching. Afterwards, we performed COSI-Corr processing only for7this of 32
contrast-stretched part and the obtained results were considerably good with acceptable
uncertainties (Table 2).

Figure 3. Similar temporal events at SF1 site: (a) 3D perspective views generated using the October
Figure 3. Similar temporal events at SF1 site: (a) 3D perspective views generated using the Oc-
2000 (post-event) and October 1999 (pre-event) Landsat 7 Panchromatic Band (15 m/pixel) (the event
tober 2000 (post-event) and October 1999 (pre-event) Landsat 7 Panchromatic Band (15 m/pixel)
occurred between January and April of 2000); (b) 3D perspective views generated using the October
(the event
2016 occurred and
(post-event) between January2016
September and (pre-event)
April of 2000); (b) 3D perspective
Sentinel-2 true colour views generated
composite (TCC)using
(10
the October
m/pixel). The2016
event(post-event) and September
occurred between the 19 and2016 (pre-event)
24 September Sentinel-2
2016. true
The yellow colour composite
rectangles highlight
the starting
(TCC) zone while
(10 m/pixel). Thethe orange
event rectangles
occurred show the
between the 19
deposition zone, before2016.
and 24 September and after the events.
The yellow rect-
Landsat
angles 7 image the
highlight IDs:starting
LE71450392000275SGS00 and LE71450391999288EDC00;
zone while the orange Sentinel-2zone,
rectangles show the deposition imagebefore
IDs:
S2A_OPER_MSI_L1C_TL_SGS__20160919T104034_A006494_T44RKU
and after the events. Landsat 7 image IDs: LE71450392000275SGS00 and LE71450391999288EDC00; and
S2A_OPER_MSI_L1C_TL_SGS__20161009T053631_A006780_T44RLU.
Sentinel-2 image IDs: S2A_OPER_MSI_L1C_TL_SGS__20160919T104034_A006494_T44RKU and
S2A_OPER_MSI_L1C_TL_SGS__20161009T053631_A006780_T44RLU.
Table 2. Estimated uncertainties in SVs. Please note that Cosi-Corr provides the displacement val-
ues, which are divided by the number of days between the image pairs to obtain velocities, ex-
Table 2. Estimated uncertainties in SVs. Please note that Cosi-Corr provides the displacement values,
plaining the corresponding, seemingly low, uncertainty values.
which are divided by the number of days between the image pairs to obtain velocities, explaining the
corresponding, seemingly low, uncertainty values. Cumulative Uncertainty
Correlated Pair (Represented by Acquisition Dates)
(±cm/day)
Correlated Pair (Represented by Acquisition Dates) Cumulative Uncertainty (±cm/Day)
1–17 September 2015 2.13
1–17 3–19
September 2015 2016
September 2.13 1.25
3–19 September 2016 1.25
17 September 2015–19 September 2016 1.48
17 September 2015–19 September 2016 1.48
9 October
9 October 2016–14
2016–14 October
October 2017 2017 1.19 1.19
14 October
14 October 2017–24
2017–24 October
October 2018 2018 1.24 1.24
24 October 2018–24
24 October October
2018–24 2019 2019
October 1.56 1.56
24 October 2019–18 October 2020 1.23
27 December 2020–16 January 2021 2.53
16–26 January 2021 3.21
26 January 2021–5 February 2021 31.70

Lacroix et al. [43] suggested to keep the window size for COSI-Corr processing as a
compromise between the robustness of large windows and the ability of smaller windows
to detect motions of rather small objects. Thus, we tried several window sizes and the most
effective window size for our study site with the least uncertainties for both Landsat 8 and
Sentinel-2 was found to be 32 × 32 window. We used this window for further processing
all the data presented in this study. In the post-processing stage, we applied two of the
filters effectively used by previous studies (e.g., [48]) to make the results more robust.
Remote Sens. 2022, 14, 949 8 of 32

Using the SNR filter, pixels with low SNR values (SNR < 0.90) were masked as poorly
correlated pixels. Using the magnitude filter, we masked the abrupt values, assuming that
the movement may not change abruptly but gradually [42], thus obtaining spatially smooth
SV maps. However, we intentionally did not apply the magnitude filter for observing yearly
movements as we wanted to identify the localised zones showing anomalous movements
on a yearly scale. This approach perfectly served the purpose, as we could see the gradual
spread of the stress field along the indentation for SF2 hanging glacier with each passing
year until 2019.

3.2. Surface Velocity Uncertainty Estimation

The uncertainty introduced due to possible orthorectification errors may result in some
minor sub-pixel-scale horizontal shifts in the Landsat 8 and Sentinel-2 images. However, in
case of imageries for the same path and row used within this study, the orthorectification
uncertainties affecting SV estimations are negligible [49]. Eventually, for estimating the
cumulative uncertainties in our SV estimates, we adopted the approach suggested by previ-
ous studies (e.g., [48,50]). We sampled all the pixels representing the stable terrain devoid
of ice and vegetation. For these stable pixels, the overall mean and standard deviation of
displacements were summed up and divided by the temporal separation of the correlated
images to obtain the overall uncertainty per correlated scene-pair. The estimated uncertain-
ties are provided in Table 2 and were of the same order of magnitude as reported by the
previous studies (e.g., [48,50]), except for the last pair of 26 January 2021–5 February 2021,
owing to extreme snow cover in 5 February 2021 image.

3.3. Rapid Mass Movement Simulation (RAMMS)

The RAMMS tool for mass movement modelling, developed at the Swiss Institute for
Snow and Avalanche Research, was successfully implemented for high mountain regions,
such as the Alps (e.g., [51–53]), Andes (e.g., [54]), Himalaya (e.g., [55]), Southern Alps
(e.g., [56]), Tatras (e.g., [57]), and others for realistic thermomechanical simulations. The
RAMMS tool [13,25,58] has a user-friendly interface that allows for the integration of
state-of-the-art numerical solution methods with input parameters derived from DTMs
and satellite images or maps. An efficient second-order numerical solution of the depth-
averaged avalanche dynamics equations forms the core of this programme [25]. More
specifically, the model relies on the two-parameter Voellmy model, which is actually field
calibrated using data from a real scale avalanche test site in Vallée de la Sionne, canton
Valais, Switzerland [25]. The usual input parameters consist of topography, release area,
release height, forest extent, and calibrated frictional coefficients, while the modelled output
includes all types of mechanical parameters, such as avalanche velocity, momentum, flux,
pressure, runout, and depositional height. Generating input data for this programme
is straightforward as all the input parameters can easily be derived from the DTM and
satellite images in the form of geographical information system (GIS) compatible shape
files. The software manual provides sufficient details on how the tool works and how
various frictional coefficients and input parameters can be manipulated to obtain the most
realistic simulations. Detailed information on RAMMS and its input and output parameters
can be downloaded from the software manual webpage [59]. More importantly, not long
ago, a rock/ice avalanche model capable of modelling both extreme cases, i.e., avalanches
containing mainly ice (e.g., SF1 in the present study) as well as avalanches that mainly
consist of rock debris (e.g., SF2 in the present study), was embedded in the RAMMS
software and was already successfully used by several case studies worldwide [60].
The hazard potential of many of the recent rock/ice avalanches, including SF2 in
Chamoli, has been complex to comprehend, due to the extreme flow transitions among
the constituents, i.e., ice, snow, and debris, and the temperature of the material. The
proportion of the constituents (ice, snow, and debris) usually decides whether the dry
fluidised (powder avalanches) flow regime would prevail, or a fluid-lubricated (debris
flows) flow regime. The thermomechanical properties of the material entrained by the
Remote Sens. 2022, 14, 949 9 of 32

avalanche often determine its flow regime and modelling rock/ice avalanche speed and
runout, and therefore requires a thermomechanical model accounting for the melting and
lubricating properties of the entrained material, which is usually a mixture of snow, ice,
water, and debris [13,58]. RAMMS is capable of modelling the thermomechanical regime of
various types of avalanches, and therefore we used it as a tool to model during-event flow
characteristics of SF1 and SF2.
High-resolution DTMs are extremely useful for geomorphology and mountain hazard
research [61,62] and are most preferred for RAMMS simulations, as they inform the model
about all the obstructions and natural dams to employ the most appropriate frictional pa-
rameters for a particular event. Fortunately, we have the availability of both pre-event and
post-event high-resolution (2 m/pixel) DTMs covering SF1 and SF2 (Table 1). These DTMs
were used in a study by Shugar et al. [21] and were derived by Bhushan and Shean [36]
and Shean et al. [40] from the stereopairs of WorldView-1, WorldView-2, WorldView-3, and
GeoEye-1 satellite images. Table 1 provides all the related references detailing the genera-
tion and characteristics of these DTMs. These DTMs allowed us to precisely demarcate and
estimate the release area and volumes of SF1 and SF2 through DTM differencing as input to
RAMMS. These areal and volumetric estimates are detailed in the Results and Discussion
Sections. Another model input has to be the avalanche return period, which is pre-defined
in the model as 10 years, 30 years, 100 years, and 300 years, representing various known
or hypothetical recurrences. For SF1, we examined historical data and found it to be a
recurring event, once the hanging glacier reaches its critical geometry. We could confirm
that nearly 16 years ago, an event of same areal scale occurred at the SF1 site (Figure 3).
Therefore, an average return period of 10 years was selected for SF1. This also translated in
the simulation results in which the simulated SF1 avalanche followed the exact depositional
boundaries that were observed in the after-event satellite images. We further ran the model
using 30 years as the return period and the results were not equally comparable to the
after-event satellite images that we derived using 10 years as the return period. On the
other hand, SF2 predominantly consisted of rock debris [21] and was more suited for
simulation within debris flow module of RAMMS, without a need to define the return
period. Considering that both SF1 and SF2 were over 10 × 106 m3 , the avalanche volume
category within RAMMS was set to “Large”. The model also takes into account the forest
cover present along the avalanche path as an important input parameter, since vegetation
provides resistance to the avalanche front. We digitised the forest mask using satellite
images acquired prior to both the events. SF1 was mostly uninterrupted by any forest
cover, while SF2, being a larger event, encountered dense vegetation along its path in lower
reaches of the valley. Considering that SF1 was majorly an ice avalanche, for the avalanche
core at rest, we used a density of 850 kg m−3 as recommended by Margreth et al. [53]. For
SF2, considering the suggested rock:ice proportion [21] of 80:20 and the predominance of
metamorphic rocks at the avalanche site [23], we inputted a density of ~2050 kg m−3 [63].
As the final vital input to RAMMS, the Coulomb friction (µ) and the turbulent friction
(ε) are defined. The friction values µ and ε strongly depend on the return period and
avalanche volume, and, in our study, owing to high-resolution DTMs and temporal image
analysis, we were certain of these input parameters. The friction values, µ and ε, depend
on the configurational energy content of the avalanche core [13], and in the avalanche
module used for SF1 simulation in RAMMS, if a calculation is obtained with constant
friction values, no terrain undulations and forest areas are considered. Therefore, for SF1,
we opted for the variable friction values, which are automatically estimated by RAMMS,
based on DTM-derived parameters (slope angle, altitude, and curvature), and we further
provided forest information and the global parameters of the return period and avalanche
volume. In case of SF2, the debris flow module of RAMMS needed to be initialised through
a carefully selected combination of µ and ε, calibrated using the documented [21–24] time-
and flow-related details of the event. This calibration was achieved through the observed
runout extent on the post-disaster satellite imagery and the reported debris flow zones, flow
heights, and velocities [21–24]. Thus, the model was precisely calibrated through extensive
Remote Sens. 2022, 14, 949 10 of 32

iterations to define the optimal friction values of µ and ε as 0.13 and 200 ms−2 , respectively.
These values correspond well with the suggested range of values in the RAMMS manual.
A general topic of previous investigations for SF2 has been to explain the possible
source of entrainment material, which led to such an enormous debris flow. A recent
study [29], based on field observations, suggests that the vast deposits of the 2016 event
might have been a possible source of material, which significantly enlarged the volume and
extent of SF2. To investigate this aspect, we further used the erosion module in RAMMS
while simulating SF2. We provided variable depths of sediment erosions between 2–8 m,
based on the simulations by Jiang et al. [22], in the upper reaches of Ronti Gad, in order
to model the increase in volume and velocity of SF2 as it travelled along the channel.
This erosion corresponds to the material that aids the entrainment of the avalanche front.
Supplementary Materials Files (Video S1: Simulated SF1 event, Video S2: Simulated SF2
event) of this paper provide a 3D visualisation of the RAMMS modelling results.

4. Results and Discussion

We present and discuss the main results of this research under the following subhead-
ings, focussing on the pre-event and during-event flow characteristics of SF1, the long-term
and short-term pre-event flow characteristics of SF2, and the simulated during-event runout
characteristics of SF2.

4.1. SF1: Key Observations and Pre-Event and During-Event Flow Characteristics
SF1 site (Figures 1 and 3) appears to be a hotspot for ice avalanching from the same
initiation zone. The maximum cross-sectional lengths and area of the broken part of the
hanging glacier during SF1 were ~0.57 km (across the flow), ~0.92 km (along the flow), and
~0.28 km2 , respectively. SF1 was a frontal block failure and the multi-temporal DTMs of 2021
and 2015 provided by Shugar et al. [21] helped us to estimate a volume of ~10 × 106 m3 ,
corresponding predominantly to the glacier ice and possibly some seasonal snow. Using
multi-temporal satellite images, we observed at least two major and comparable slope
failure events at the SF1 site during the past two decades. These events were spaced
~16.5 years apart and they produced exactly the same indentation (Figure 3), marked by the
dashed red curve in Figure 1c. For the earlier event, the 3D perspective views in Figure 3a
depict the avalanched zone, and by using all the available Landsat images, we could narrow-
down the timing of this event between January and April 2000. Similarly, the 3D perspective
views in Figure 3b highlight the repetition of another avalanche (called SF1 in this paper) at
the same indentation site. Using freely accessible Sentinel-2 and PlanetScope [64] images,
we could narrow-down the timing of SF1 between 19 and 24 September 2016.
To understand the enormity of SF2, there is a need to consider these events from 2000
and 2016 at the SF1 site. These events were large-volume (~10 × 106 m3 ) ice avalanche
events leading to large deposits of ice and debris (orange rectangles in Figure 3) in the
avalanche deposition zone in a narrow mountain stream named Ronti Gad. This deposition,
following the 2000 and 2016 events, covered a ~3.5 km length of the stream (Figures 3 and 4).
The thick ice-debris deposits from these events tend to persist for multiple years after the
event. For example, some of the possible remnants of the deposits from the 2000 event
could be observed even in 2005, more than 5 years after the event (Figure 5b).
Another observation made on multiple years images (Figure 5) is the significantly
seasonally active nature of the avalanche slope (Figure 5h). In addition to the ice avalanches,
these slopes are always susceptible to snow avalanches during the winter months (Figure 5).
A ~500 m-long stretch of the valley, directly beneath the avalanche slopes, is mostly covered
by compacted snow, ice, and debris (Figure 5). These deposits also often hinder the full
flow of the narrow stream, as is evident in Figure 4b, building up the temporary zones
of increased hydraulic gradient and developing the dammed water pockets along the
stream. We identified one such seasonally active avalanche slope below SF1 and SF2
sites (Figures 5 and 6). While in the 3 November 2013 image of Figure 5 we observe a
clear stream devoid of any avalanche deposit from the previous event of 2000 or from
Remote Sens. 2022, 14, 949 11 of 32

subsequent snow avalanches, in the 3 June 2014 and the 7 May 2015 images, we see thick
fresh avalanche deposits covering the stream valley. Figure 5h highlights the seasonally
active avalanche slope below the SF1 and SF2 sites, which accumulates significant snow
mass during winter, often leading to slope failures relatively smaller than SF1 and SF2,
but still large enough to form a deposition zone of ~500 m to 1 km in length. This large
volume of snow, ice, and debris can accentuate any subsequent slope failure events, such
as SF2, which occurred in the winter of 2021. We observed the available images of these
Remote Sens. 2022, 14, x FOR PEER REVIEW 11 of 33
slopes, days before SF2. Figure 6 highlights the significant snow accumulation on the
seasonally active slope on 2 February 2021, while on 5 February 2021 Sentinel-2 image,
it could be seen that a fresh cumulative snowfall of ~112 mm [22] significantly added to
events
the snowleading to large
cover. Figure deposits of shows
6 additionally ice andthe debris (orange
developed rectanglesfrom
indentation in Figure
where 3)thein the
avalanche
hanging deposition
glacier broke on zone in a narrow
7 February 2021, mountain
causing SF2. stream named
Terrain Ronti
profiles forGad.
SF1 andThisSF2
deposi-
tion, following
slopes (Figure 7),the 2000 and 2016
representing events,before
the terrain covered a ~3.5
either kmreported
of the length ofevents
the stream (Figures
occurred,
3 and 4).
estimate The
the thick ice-debris
maximum deposits
slope reaching fromas
as high ~37◦events
these tend to average
. The reported persist for multiple
slope for theyears
entire ◦
after SF2 mass [24]
the event. Forwas even higher
example, some ofat ~40 . These slope
the possible values
remnants ofare
theentirely
depositswithin
from the
the 2000
most
eventavalanche-prone slopeeven
could be observed range
inreported
2005, morefromthan
this5region
years [18].
after the event (Figure 5b).

Figure4. 4.
Figure Large
Large deposition
deposition of and
of ice ice and debris
debris during
during previous
previous slope failures
slope failures in 2016
in 2016 and 2000 and
from2000
SF1 from
site: (a) ice-debris deposition reaching ~3.5 km down-valley after the 2016 event; (b) an ice-debris daman ice-
SF1 site: (a) ice-debris deposition reaching ~3.5 km down-valley after the 2016 event; (b)
debris
(red arrowdam (red image),
in 2003 arrow in~22003 image), ~2 km
km downstream downstream
from the startingfrom
pointthe starting
of the pointzone,
deposition of thewhich
deposition
zone, which eventually disappeared in the 2005 image. Source: Google Earth images (Credit:
eventually disappeared in the 2005 image. Source: Google Earth images (Credit: CNES/Airbus).
CNES/Airbus).
Remote
RemoteSens. 2022,14,
Sens.2022, 14,949
x FOR PEER REVIEW 12 12
ofof3332

Figure 5.5. Multi-temporal

Figure Multi-temporal (a–g) status
status of
of the
the avalanche
avalanche deposition
depositionzone
zone(dashed
(dashedyellow
yellowpolygon),
polygon),
right beneath the active slopes. The image in (e) shows a fresh avalanche and the imageinin(h)
right beneath the active slopes. The image in (e) shows a fresh avalanche and the image (h)shows
shows
the seasonally active avalanche slope (red rectangle) below SF1 and SF2, with a long runout marked
the seasonally active avalanche slope (red rectangle) below SF1 and SF2, with a long runout marked
by the red arrow. Source: Google Earth images (Credit: CNES/Airbus).
by the red arrow. Source: Google Earth images (Credit: CNES/Airbus).

Another observation
COSI-Corr made ontool
[41] is an efficient multiple years surface
to estimate images deformations
(Figure 5) is the andsignificantly
mass move-
seasonally active nature of the avalanche slope (Figure
ments through the correlations between temporal images, even in the absence 5h). In addition to theofice ava-
ground-
lanches, these slopes are always susceptible to snow avalanches
control points (GCPs). COSI-Corr was conclusively established in estimating the motion induring the winter months
(Figure
mass 5). A ~500[43,65–68],
movements m-long stretch of the
glaciers [44],valley,
and sanddirectly
dunes beneath the avalanche
[69], with the relativeslopes,
accuracy is
mostly covered by compacted snow, ice, and debris (Figure 5).
of co-registration and correlation up to ~1/20th of a pixel [41]. Thus, there were two These deposits also often
hinderobjectives
prime the full flow of the
for us narrow
to assess thestream, as is evident inofFigure
flow characteristics 4b, building
SF1. First, up the tem-
to use COSI-Corr for
porary zones of increased hydraulic gradient and developing
quantifying the SV changes leading to SF1, and, second, to perform the simulation and the dammed water pockets
alonganalyses
flow the stream. Webroken
of the identifiedice one such during
volume seasonallythe active
SF1 eventavalanche
using slope
RAMMS below SF1 and
model. We
used the cloudless and snow-free Landsat 8 satellite Panchromatic band (15 m/pixel) ato
SF2 sites (Figures 5 and 6). While in the 3 November 2013 image of Figure 5 we observe
clear stream
estimate the SVdevoid
during of1–17
any September
avalanche 2015 deposit frombefore
(a year the previous
the 2016 event
event)ofand2000
for or from
a similar
subsequent snow avalanches, in the 3 June 2014 and the 7 May
period during 3–19 September 2016 (Figure 8a,b), immediately before SF1 and at the end 2015 images, we see thick
fresh
of the avalanche
yearly ablationdeposits covering
season. We alsotheestimated
stream valley. FigureSVs
the yearly 5h during
highlights the seasonally
17 September 2015
active avalanche slope below the SF1 and SF2 sites, which accumulates
to 19 September 2016, highlighting the higher SVs in the avalanched zone of the hanging significant snow
mass during
glacier (Figurewinter, often leading
8a,b). Compared to 1–17to slope failures
September relatively
2015, the maximumsmallerSV than SF1avalanched
in the and SF2,
but still large enough to form a deposition zone of ~500 m to
part of the hanging glacier is ~5 times higher during 3–19 September 2016 (Figure 8b) 1 km in length. This large
and
volume of snow, ice, and debris can accentuate any subsequent slope
the formed stress is strikingly evident in the region that leads to SF1 (middle and rightmost failure events, such
as SF2, of
images which
Figureoccurred
8a). Theinaverage
the winter of 2021.
velocity Weyearly
at the observed scalethe available
(~12.3 ± 1.5images
cm/day) of for
thesethe
avalanched part of the hanging glacier is also ~2 times higher than the average short-term
 sonally active slope on 2 February 2021, while on 5 February 2021 Sentinel-2 image, it
could be seen that a fresh cumulative snowfall of ~112 mm [22] significantly added to the
snow cover. Figure 6 additionally shows the developed indentation from where the hang-
ing glacier broke on 7 February 2021, causing SF2. Terrain profiles for SF1 and SF2 slopes
Remote Sens. 2022, 14, 949 13 of 32
(Figure 7), representing the terrain before either of the reported events occurred, estimate
the maximum slope reaching as high as ~37°. The reported average slope for the entire
SF2 mass [24] was even higher at ~40°. These slope values are entirely within the most
velocities (~6.8 ± 1.7 cm/day) (Figure 8b). The maximum SVs show a distinctly increasing
avalanche-prone
trend before theslope
event,range
while reported
the overallfrom this
average region
SVs do not[18].
follow a pattern.

Figure 6. 3D perspective views of the slope failure sites and the seasonal avalanche slope (red rec-
Figure 6. 3D perspective views of the slope failure sites and the seasonal avalanche slope (red
tangle). The avalanche
rectangle). slope represents
The avalanche the same
slope represents terrain
the same as the
terrain as red rectangle
the red in Figure
rectangle 5: 5:
in Figure (a) 2 Feb-
ruary(a)2021 PlanetScope image [64]; (b) 5 February 2021 Sentinel-2 image. The yellow
2 February 2021 PlanetScope image [64]; (b) 5 February 2021 Sentinel-2 image. The yellow rect- rectangle
showsangle shows the developed indentation from where the major part of the hanging glacier broke on 7 Feb-
the developed indentation from where the major part of the hanging glacier broke on
ruary7 2021, causing
February SF2. The
2021, causing images
SF2. are contrast-stretched
The images totocounter
are contrast-stretched counter the shadowand
the shadow and sensor satu-
sensor
ration due to the
saturation due snow. PlanetScope
to the snow. image
PlanetScope ID:ID:20210202_052653_11_227c_3B;
image Sentinel-2
20210202_052653_11_227c_3B; Sentinel-2 image ID:
image
S2A_OPER_MSI_L2A_TL_VGS1_20210205T082130_A029374_T44RLU.
ID: S2A_OPER_MSI_L2A_TL_VGS1_20210205T082130_A029374_T44RLU.

The visualisation
COSI-Corr [41] is an of efficient
the surface deformations
tool to estimate caused by the
surface load effects, through
deformations and mass an move-
estimation of the surface displacements, has proven effective in observing the development
ments through the correlations between temporal images, even in the absence of ground-
and trajectory of mass movements [67,70]. As early as 1973, SVs were measured for the first
control
timepoints (GCPs).cold-based
on an unstable COSI-Corr was conclusively
hanging glacier in order established incollapse
to predict its estimating[71]. the
The motion
in mass
observed velocities were shown to increase as a power-law function of time, up to infinity, at accu-
movements [43,65–68], glaciers [44], and sand dunes [69], with the relative
racythe
of theoretical
co-registration
time of and correlation
failure [71]. This isup
alsotothe
~1/20th of a pixel
characteristic [41].ofThus,
signature there
the stress were two
build-
up as a critical phenomenon in the case of various other naturally occurring
prime objectives for us to assess the flow characteristics of SF1. First, to use COSI-Corr forruptures, such
as earthquakes
quantifying the SV [72], landslides
changes [73], and
leading to rockfalls
SF1, and, [74]. Velocitytomeasurements
second, perform theperformed
simulation and
on cold-based Weisshorn and Mönch hanging glaciers in Switzerland further confirmed
flow analyses of the broken ice volume during the SF1 event using RAMMS model. We
this behaviour [75]. The validity of an SV monitoring approach was demonstrated in 2014
usedbythe
thecloudless
successful and snow-free
prediction Landsat
of a hanging 8 satellite
glacier Panchromatic
“break-off” bandface
from the south (15ofm/pixel)
the to
 Remote Sens. 2022, 14, x FOR PEER REVIEW 14 of 33

Remote Sens. 2022, 14, 949 estimate the SV during 1–17 September 2015 (a year before the 2016 event) and for a14sim- of 32

ilar period during 3–19 September 2016 (Figure 8a,b), immediately before SF1 and at the
end of the yearly ablation season. We also estimated the yearly SVs during 17 September
2015 to 19
Grandes September
Jorasses (Monte2016, highlighting
Bianco, Italy), 10 the
dayshigher SVs
prior to theinevent
the avalanched
[76]. However,zoneowing
of theto
hangingissues,
logistic glaciersuch
(Figure 8a,b).
efforts wereCompared to 1–17 September
understandably limited in 2015,
their the maximum SVdomain,
spatiotemporal in the
avalanched
and all reliedpart of the
on field hanging
data glacier
collection for is
SV~5estimations.
times higher Ourduring 3–19 September
observations 2016of
in the cases
(Figure
SF1 and 8b)
SF2,and the formed
highlight thatstress
remote is strikingly evidentSVs
sensing-derived in the region
have that
large leadscoverage
spatial to SF1 (mid-
and
dle and rightmost
sufficient continuityimages
and canof represent
Figure 8a).a The average velocity
cost-effective and saferat the
wayyearly scale (~12.3
to monitor hanging ±
1.5 cm/day)
glaciers, andfor the avalanched
possibly part of
even predict the hanging
large glacier isice-rock
and dangerous also ~2 times higherin
avalanches than the
certain
average
cases. short-term
However, anyvelocities (~6.8 ±predictability
possible future 1.7 cm/day) (Figure
certainly8b).requires
The maximum SVs show
more research a
in this
distinctlyto
direction increasing trend
identify the before
types the event,
of slope while
failures the overall
which average SVssignificant
show statistically do not follow a
trends
pattern.
in SVs.

Figure 7. Elevation profile along the transects covering the avalanche slopes. The elevation is above
Figure 7. Elevation profile along the transects covering the avalanche slopes. The elevation is above
the mean sea level and the source is the Shuttle Radar Topography Mission (SRTM) DEM from 2000,
the mean sea level and the source is the Shuttle Radar Topography Mission (SRTM) DEM from 2000,
representing the terrain before any of the reported events. Source: Google Earth image (Credit:
representing
CNES/Airbus).the terrain before any of the reported events. Source: Google Earth image (Credit:
CNES/Airbus).
The visualisation of the surface deformations caused by the load effects, through an
The RAMMS
estimation simulation
of the surface further revealed
displacements, the flow
has proven characteristic
effective of this
in observing the entire SF1
develop-
event (Figures 8 and 9), which lasted for ~3 min as per the modelled results.
ment and trajectory of mass movements [67,70]. As early as 1973, SVs were measured for Based on the
mechanical characteristics,
the first time on an unstablethe cold-based
event can behanging
dividedglacier
into three prominent
in order phases
to predict its(Figure
collapse9):
(1) break-off and avalanching in transect AB, (2) impact and shearing in transect
[71]. The observed velocities were shown to increase as a power-law function of time, BB’, and
up
(3) deposition in transect B’A’. Overall, the maximum velocity, maximum flow
to infinity, at the theoretical time of failure [71]. This is also the characteristic signature height,
of
and maximum
the stress pressure
build-up reachedphenomenon
as a critical were ~90 m/s, ~100
in the m,of
case and ~6000other
various kPa, respectively.
naturally occur-The
avalanching following the frontal break-off continued for ~2300 m before it struck the
valley. The flow velocity, momentum, and pressure reached their peaks during this phase.
Immediately after the impact, the maximum shear stress reached its peak (~1200 kPa),
and this zone of ~200 m signifies the shearing and fragmentation of the ice core (Figure 9).
Finally, a depositional analysis in RAMMS simulates the deposition zone within a ~3000 m
Remote Sens. 2022, 14, 949 15 of 32
Remote Sens. 2022, 14, x FOR PEER REVIEW 15 of 33

length along the Ronti Gad stream (Figure 8c). The maximum deposition height is as much
as ~48.10
ring m insuch
ruptures, the middle regions[72],
as earthquakes of the deposition
landslides [73], zone (Figure [74].
and rockfalls 8c), highlighting
Velocity meas- the
magnitude of the SF1 on
urements performed event. This also
cold-based points towards
Weisshorn and Mönchthe possibility that SF1ininSwitzer-
hanging glaciers 2016 and
the
landprevious event in 2000
further confirmed this might have[75].
behaviour preconditioned
The validity theof an valley significantly
SV monitoring in terms
approach
of
wassedimentation
demonstratedavailability
in 2014 by and topographic
the successful smoothness,
prediction leadingglacier
of a hanging to the“break-off”
enormity of
debris
from theflow during
south faceSF2of in
the2021. We further
Grandes generated
Jorasses the deposition
(Monte Bianco, Italy), 10outline for SF1
days prior to and
the it
remarkably matched owing
event [76]. However, the ice-debris
to logisticextent observed
issues, in 9were
such efforts October 2016 post-SF1
understandably Sentinel-
limited in
2their spatiotemporal
image. domain, context,
In the Himalayan and all relied on field data
the modelled flowcollection for SV estimations.
characteristics Our
and mechanical
observations
parameters incomparable
are the cases of to SF1those
and SF2, highlight
reported that
for the remote
Gyari sensing-derived
glacier avalanche inSVs thehave
Shyok
large spatial
Basin coverage
of the Indus andwhich
River, sufficient
killedcontinuity
139 peopleand[55].
can represent
The volume a cost-effective
of the Gyariand safer
avalanche
way4.31
was × 106 m
to monitor 3 , as opposed
hanging glaciers,toand × 106 m
~10possibly 3 ofpredict
even SF1, but large
theand dangerous
maximum ice-rock
height of the
avalanches
deposits in certain
reached cases.asHowever,
as high 40 m [55].any possible
The future
difference predictability
between certainlymaximum
the simulated requires
more research
pressure in this
for Gyari direction
(~2500 kPa)toandidentify the types
SF1 (~6000 of slope
kPa) failures
was also which show
proportional to statisti-
the mass
cally significant trends in SVs.
difference between the two events [55].

Figure 8. Surface velocity (SV) and RAMMS simulation for SF1: (a) 3D perspective views with over-
Figure 8. Surface velocity (SV) and RAMMS simulation for SF1: (a) 3D perspective views with
laid short-term and long-term SVs. The Sentinel-2 image of 9 October 2016 (after SF1) is in the back-
overlaid
ground. short-term andramp
The SV colour long-term
for allSVs. The Sentinel-2
the time image ofto
periods is adjusted 9 October
the same2016range (after SF1)m/day)
(0–1.27 is in the
background. The SV colour
for visual comparisons; (b) ramp for plots
detailed all theoftime
SVsperiods is adjusted
for the entire to the
hanging same(i.e.,
glacier range (0–1.27
HG), m/day)
including
for
its visual
intact comparisons;
zone (IZ) and(b) detailed plots
avalanched zone of SVs For
(AZ). for the
the entire
IZ andhanging
the AZ,glacier
the SVs(i.e.,
areHG),
alsoincluding
providedits
separately
intact withand
zone (IZ) associated uncertainties;
avalanched zone (AZ).(c) ForRapid
the IZMass Movement
and the Simulation
AZ, the SVs are also (RAMMS) model-
provided separately
ling of SF1. The panel left-to-right shows the maximum velocity field, flow height
with associated uncertainties; (c) Rapid Mass Movement Simulation (RAMMS) modelling of SF1. The and avalanche
deposition
panel zone, and
left-to-right the transect
shows and respective
the maximum velocityavalanche
field, flowdeposit
heightprofile.
and avalanche deposition zone,
and the transect and respective avalanche deposit profile.
The RAMMS simulation further revealed the flow characteristic of this entire SF1
event
4.2. (Figures
SF2: 8 andand
Long-Term 9), Short-Term
which lasted for ~3Velocity
Surface min as Evolution
per the modelled results. Based on the
mechanical characteristics, the event can be divided into three prominent phases (Figure
We further monitored the temporal evolution in the SVs of SF2 hanging glacier. For the
9): (1) break-off and avalanching in transect AB, (2) impact and shearing in transect BB’,
long-term yearly SV estimation, we used available best quality multi-temporal Sentinel-2
and (3) deposition in transect B’A’. Overall, the maximum velocity, maximum flow height,
images from the October months of all the years, to ensure similar shadow and surface
and maximum pressure reached were ~90 m/s, ~100 m, and ~6000 kPa, respectively. The
conditions. Cloud-free images from the September months were not available for some
avalanching following the frontal break-off continued for ~2300 m before it struck the val-
of the years. 10 m/pixel bands of Sentinel-2 images have been reported to be extremely
ley. The flow velocity, momentum, and pressure reached their peaks during this phase.
 Immediately after the impact, the maximum shear stress reached its peak (~1200 kPa), and
this zone of ~200 m signifies the shearing and fragmentation of the ice core (Figure 9).
Finally, a depositional analysis in RAMMS simulates the deposition zone within a ~3000
Remote Sens. 2022, 14, 949 m length along the Ronti Gad stream (Figure 8c). The maximum deposition height 16 of 32 is as
much as ~48.10 m in the middle regions of the deposition zone (Figure 8c), highlighting
the magnitude of the SF1 event. This also points towards the possibility that SF1 in 2016
and the previous
useful event
in detecting in 2000 motions
precursory might have preconditioned
in the mass movements theusing
valley significantly
COSI-Corr [43]. in
Weterms
of generated
sedimentation availability
a time-series and topographic
of the yearly smoothness,
SVs (Figure 10a) to observeleading to thechanges
the long-term enormityandof de-
bris
theflow
zones during
of the SF2 in 2021.
hanging Weshowing
glacier furtheranomalous
generatedSVs. the During
deposition outline
October for SF1 and it
2016–October
2017, we can
remarkably observethe
matched anomalous ice extent
ice-debris movement fields developing
observed in 9 Octoberin SF2
2016break-off
post-SF1 zone of
Sentinel-2
the hanging glacier, below its bergschrund (a large crevasse often
image. In the Himalayan context, the modelled flow characteristics and mechanical marking the head of a pa-
mountain glacier), with the maximum SV reaching as high as ~39.2 ±
rameters are comparable to those reported for the Gyari glacier avalanche in the Shyok 1.2 cm/day, while
the average SV for the entire hanging glacier is ~2.9 ± 1.2 cm/day (Figure 10a). The aver-
Basin of the Indus River, which killed 139 people [55]. The volume of the Gyari avalanche
age SV kept6on 3fluctuating through 2017–2019; ~5.9 ± 1.2 cm/day during 2017–2018 and
was 4.31 × 10 m , as opposed to ~10 × 106 m3 of SF1, but the maximum height of the de-
~3.1 ± 1.6 cm/day during 2018–2019 (Figure 10a). However, the most interesting observa-
posits reached
tion is as high
the distinct stressas 40 m
field [55]. The
growing difference
across the entire between the simulated
break-off zone, maximum pres-
below the bergschrund
sure for hanging
of the Gyari (~2500
glacier,kPa)
withand SF1 (~6000SVs
the maximum kPa) was also
reaching proportional
up to to the mass
~51.3 ± 1.2 cm/day during differ-
ence between
2017–2018 the two
(Figure 10a).events [55].

Figure 9. Mechanical characteristics of SF1 during the event: the left panel shows the transect of SF1
Figure 9. Mechanical characteristics of SF1 during the event: the left panel shows the transect of SF1
along which various mechanical characteristics are plotted in the right panel. The knickpoints where
along which various mechanical characteristics are plotted in the right panel. The knickpoints where
the avalanche transformed its mechanical properties are marked by dashed white lines perpendicular
to the runout. (MFH = Maximum Flow Height, MFV = Maximum Flow Velocity, MFM = Maximum
Flow Momentum, MFP = Maximum Flow Pressure, and MSS = Maximum Shear Stress).
 ing the head of a mountain glacier), with the maximum SV reaching as high as ~39.2 ± 1.2
cm/day, while the average SV for the entire hanging glacier is ~2.9 ± 1.2 cm/day (Figure
10a). The average SV kept on fluctuating through 2017–2019; ~5.9 ± 1.2 cm/day during
2017–2018 and ~3.1 ± 1.6 cm/day during 2018–2019 (Figure 10a). However, the most inter-
esting observation is the distinct stress field growing across the entire break-off zone, be-
Remote Sens. 2022, 14, 949 17 of 32
low the bergschrund of the hanging glacier, with the maximum SVs reaching up to ~51.3
± 1.2 cm/day during 2017–2018 (Figure 10a).

Figure 10. SV evolution for SF2 hanging glacier: (a) 3D perspective views with overlaid yearly SVs
Figure 10. SV evolution for SF2 hanging glacier: (a) 3D perspective views with overlaid yearly SVs in
in the top panel and the bottom panel shows the graphical representation of evolving SVs. The SV
the top ramp
colour panel for
andall
thethe
bottom panel shows
time periods the graphical
is adjusted representation
to the same of evolving
range (0–0.6 m/day) SVs. The SV
for visual colour
compar-
ramp for all the time periods is adjusted to the same range (0–0.6 m/day) for visual comparisons;
(b) 3D perspective views with overlaid short-term SVs and the bottom panel shows the graphical
representation of evolving SVs. The SV colour ramp for all the time periods is adjusted to the same
range (0–6 m/day) for visual comparisons.

These estimations are in accordance with the conclusions provided by Shugar et al. [21],
where they report the destabilisation in SF2 mass starting from 2016, i.e., after SF1, and
displaying maximum displacements during the summer months of 2017 and 2018. These
SV values indicate that SF2 mass was already unstable post-SF1, and experienced significant
displacement by 2020. The indentation that was almost invisible in 2016, grew to be as large
as ~60–80 m by October 2020 (Figure 11), highlighting the major displacements during this
period. The heavy precipitation [26,29] in the early weeks of February 2021 further added
to this large, destabilised volume (~27 × 106 m3 ) of rock, ice, and snow, leading to SF2. In
fact, as per the Supplementary Data provided by Jiang et al. [22], the cumulative snowfall
only on 5 February 2021 was ~112 mm, considerably higher than any of the days of that
winter season.
SF1 in itself was a fairly large event that exerted a huge pressure of ~6000 kPa at
its impact point, and the localised tremor could have initiated the destabilisation in the
bedrock of SF2 hanging glacier below its bergschrund. A valid question here might be
that why a similar intensity 2000 event at SF1 site did not trigger the destabilisation in the
bedrock of SF2 hanging glacier. The possible answer lies in the timing of both these events
at SF1 site. The 2000 event occurred in the winter months, while the 2016 SF1 occurred
in September, i.e., at the end of summer, when the bedrock is more prone to detachment.
The powerful impact of SF1 might have been sufficient to cause such detachment in the
Remote Sens. 2022, 14, 949 18 of 32

bedrock of SF2 hanging glacier, which continues to grow during the following years
(Figure
Remote Sens. 2022, 14, x FOR PEER REVIEW 11). However, it will need a separate detailed investigation to conclusively prove
18 of 33
if impacts such as SF1 can trigger further glacier detachments. Here, we also report the
short-term SV changes during the ~40 days prior to SF2 (Figure 10b). While the average
SV during
isons; (b) 3D27 December
perspective views2020–16 January
with overlaid 2021 was
short-term ~11.9
SVs ± 2.5
and the cm/day,
bottom panel it nearly
shows thedoubles
graph-
to ~23.2 ± 3.2 cm/day during 16–26 January 2021, and reaches over
ical representation of evolving SVs. The SV colour ramp for all the time periods is adjusted eight times thistofigure
the
to
same range ±
~105.2 31.7
(0–6 cm/day
m/day) during
for visual 26 January–5 February 2021 (Figure 10b). However,
comparisons.
an even more dramatic increase is recorded in the maximum SVs, showing a ~10-fold
These estimations
rise during these ~40 days,are in with
accordance with the SVs
the maximum conclusions ~481.8 ±by
reachingprovided Shugar
31.7 cm/day et al.by
5[21], where they
February 2021,report
2 days theprior
destabilisation
to SF2 (Figurein SF210b).
mass However,
starting from 2016,
these i.e., aftershould
estimates SF1,
andinterpreted
be displaying with maximum
care due displacements during the summer
to the proportionally months of 2017
high uncertainties and 2018.
associated with
These
the SVSV values indicate
measurements forthat
the SF2 mass was February
26 January–5 already unstable
2021 periodpost-SF1, and10b),
(Figure experienced
owing to
significant displacement
significantly high snow cover by 2020.
and The
poorindentation
illumination that was almost
conditions in 5invisible
Februaryin20212016,Sentinel-
grew
to be as large as ~60–80 m by October 2020 (Figure 11), highlighting
2 image (Figures 6b and 11f). Such conditions are reported to significantly affect the quality the major displace-
ments
and during this
accuracy of SVperiod. The heavy
estimations [44],precipitation
and therefore,[26,29]
weincannot
the early weeks of February
conclusively say if the
2021 further
derived added
SVs for the to
26this large, destabilised
January–5 February 2021 volume
period(~27 × 10to6 m
refer 3) of rock,
actual ice, and snow,
displacements or are
leading atoresult
majorly SF2. Inof fact,
image asartifact.
per the Nevertheless,
supplementary datashort-term
these provided SVs by Jiang
provideet al. [22], the
a sequential
cumulative
picture of thesnowfall only
pre-event on 5characteristics,
flow February 2021 during
was ~112 mm, considerably
a period of approximatelyhigherone thanmonth,
any
of the days
leading of that winter season.
to SF2.

Figure 11. Multi-temporal (a–f) true colour composites (TCCs) of Sentinel-2 images showing the
Figure 11. Multi-temporal (a–f) true colour composites (TCCs) of Sentinel-2 images showing the
evolving break-off zone (red ellipses) for SF2.
evolving break-off zone (red ellipses) for SF2.
SF1 in itself was a fairly large event that exerted a huge pressure of ~6000 kPa at its
4.3. SF2: Key Observations and During-Event Flow Characteristics
impact point, and the localised tremor could have initiated the destabilisation in the bed-
2 with
rock The total
of SF2 area ofglacier
hanging the avalanched part of SF2 hanging
below its bergschrund. glacier was
A valid question here~0.27
mightkm be that
the
whymaximum length of2000
a similar intensity ~850event
m. The destabilisation
at SF1 site did notintrigger
the bedrock started afterinSF1
the destabilisation thein
2016 at an elevation of ~5600 m asl (Figure 11), and on 7 February 2021, this sliding
bedrock of SF2 hanging glacier. The possible answer lies in the timing of both these events mass
of rocksite.
at SF1 andThe ice2000
crashed
eventinto the Ronti
occurred in theGad valley
winter at anwhile
months, elevation of ~3800
the 2016 m asl. in
SF1 occurred We
simulated this SF2 event using high-resolution DTMs and RAMMS model.
September, i.e., at the end of summer, when the bedrock is more prone to detachment. We contrast-
The powerful impact of SF1 might have been sufficient to cause such detachment in the
bedrock of SF2 hanging glacier, which continues to grow during the following years (Fig-
ure 11). However, it will need a separate detailed investigation to conclusively prove if
Remote Sens. 2022, 14, 949 19 of 32

stretched a high-resolution image of 11 February 2021 (Figure 12e), and the scar certainly
suggested a significant removal of bedrock, as also modelled by Shugar et al. [21]. This event
started as a landslide and turned into a rock-ice avalanche. Real-time images provided by
PlanetScope [64] suggest that before Raini village, SF2 turned into a more fluidised debris
flow, causing a flash flood in Rishiganga River. This flash flood took nearly 27 min to cover
the ~6 km of the distance between Raini village and Tapovan town (Figure 13), clocking
a disastrous velocity of ~4 m/s. This velocity is comparable to those of several reported
flash floods in mountain catchments and is capable of carrying a significant amount of
debris and sediment to a large distance, destroying any obstructions in the path [77,78].
Remote Sens. 2022, 14, x FOR PEER REVIEW 20 of 33
The large debris and sedimentation spread along the river channel after SF2 is visible on
the 10 February 2021 Sentinel-2 image of Figure 14d.

Figure 12.
Figure 12. Pre-SF2
Pre-SF2 (a,c)
(a,c) and
and post-SF2
post-SF2 (b,d)
(b,d) images
images of
of the
the avalanche
avalanche scar
scar (zoomed-in view in
(zoomed-in view in (e))
(e)) and
and
valley. The
valley. The extent
extent of
of coverage
coverage of
of SF1
SF1 deposits
deposits can
can be
be visualised
visualised by
by comparing
comparing (c,d).
(c,d). The
The red
red ellipse
ellipse
shows the spur that was obliterated by SF2. Source: Google Earth images (Credit: CNES/Airbus).
shows the spur that was obliterated by SF2. Source: Google Earth images (Credit: CNES/Airbus).

The RAMMS simulation provided us an informative account of the mechanical char-

acteristics of SF2 (Figure 15), which further helped us in deciphering the various zones of
this entire event (Figure 16). The reported accounts [21–24,29] based on field validation, air
surveys, and the visual interpretation of satellite images suggest that the total distance be-
tween SF2 site and the end of the rock-ice avalanche deposition zone was ~13 km (~2 km of
mountain slope and ~11 km of river valley). From this point, the deposits started becoming
more fluidised and turned into a debris flood [21]. The same published sources [21–24,29],
using seismic modelling and videographic evidence provided by the local population,
also put a constraint on the timeline of this event. The initial slope failure occurred at
10:21 Indian Standard Time (IST) (4:51 Greenwich Mean Time (GMT)) and it took slightly
less than a minute to hit the valley [21]. PlanetScope image taken at 5:01:51 GMT (Figure 13),
shows the SF2 front at the exact bend before Raini village where the rock-ice avalanche
started turning into a debris flow. This indicates that SF2 rock-ice avalanche front took
10 min to travel to the transition point. However, this transition bend was also one of the
sharpest and narrowest along the valley, stopping the runout and depositing the majority of
Remote Sens. 2022, 14, 949 20 of 32

the avalanched and entrained mass. This further consolidation of rock debris and ice along
Rishiganga valley could have taken a few more minutes, as a long trail of avalanche can
still be seen in Figure 13a. Thus, we calibrated our RAMMS model for the aforementioned
distance and time to reach the optimal friction values of µ and ε. The final model outcome
presented us with a time of ~13.5 min for the entire SF2 avalanche event, from slope failure
until complete deposition. This is in accordance with the estimated timeline above and the
simulated timeline provided by Jiang et al. [22]. This consistency in the spatiotemporal
premise for the final outcomes gave us the needed confidence in the simulated mechanical
parameters for SF2 (Figures 15 and 16). The improving satellite remote sensing capabilities,
such as the temporal resolutions, help to study natural hazards more efficiently, as in this
case too, PlanetScope was able to capture the first image within ~10 min of the
Remote Sens. 2022, 14, x FOR PEER REVIEW
event, giving
21 of 33
us a reliable timestamp for calibrating the model.

Figure 13. PlanetScope images showing the flood and debris-cloud in real-time (~10 min after the
Figure 13. PlanetScope imagesavalanche
landslide): (a) rock-ice showingfrontthe flood
beyond and debris-cloud
the confluence in Rishiganga
of Ronti Gad and real-time(Image(~10 min after
the landslide):ID:(a)
20210207_050151_1039_3B
rock-ice avalanche andfront
20210207_050152_1039_3B);
beyond the confluence(b) Debris of
flowRonti
leadingGad
to flash floodRishiganga
and
beyond Raini village (Image ID: 20210207_052833_01_2413_3B and 20210207_050152_1039_3B).
(Image ID: 20210207_050151_1039_3B
Courtesy: Planet Team [64]. and 20210207_050152_1039_3B); (b) Debris flow leading to flash
flood beyond Raini village (Image ID: 20210207_052833_01_2413_3B and 20210207_050152_1039_3B).
Courtesy: Planet Team [64].
Remote Sens. 2022, 14, 949 21 of 32
Remote Sens. 2022, 14, x FOR PEER REVIEW 22 of 33

Figure 14. Temporal changes (a–d) in the SF1 deposition zone, as seen on Sentinel-2 images. Yellow
Figure 14. Temporal changes (a–d) in the SF1 deposition zone, as seen on Sentinel-2 images. Yellow
outlines in the 18 October 2020 image highlight the regions along the channel, which significantly
outlines in the 18 October 2020 image highlight the regions along the channel, which significantly
changed after SF1 and accumulated significant ice-debris deposits. Red arrows mark the avalanche
changed after SF1image
slope. Sentinel-2 and accumulated significant ice-debris deposits. Red arrows mark the avalanche
IDs: S2A_OPER_MSI_L1C_TL_SGS__20160919T104034_A006494_T44RKU,
slope. Sentinel-2 image IDs: S2A_OPER_MSI_L1C_TL_SGS__20160919T104034_A006494_T44RKU,
S2A_OPER_MSI_L1C_TL_SGS__20161009T053631_A006780_T44RLU,
S2A_OPER_MSI_L1C_TL_SGS__20161009T053631_A006780_T44RLU,
S2A_OPER_MSI_L1C_TL_VGS1_20201018T072817_A027801_T44RLU, S2A_OPER_MSI_L1C_TL_ and
S2B_MSIL1C_20210210T051939_N0209_R062_T44RLU.
VGS1_20201018T072817_A027801_T44RLU, and S2B_MSIL1C_20210210T051939_N0209_R062_T44RLU.

The maximum
RAMMS simulation velocity reached
providedbyusSF2 was ~65 m/s,
an informative whileof
account most of the flow height
the mechanical char-
remains
acteristicsbelow
of SF2120 m (Figure
(Figure 16). These
15), which furthervalues
helped also correspond
us in deciphering well thewith other
various studies;
zones of
Shugar et al. [21] reported flow height of up to 120 m and the maximum
this entire event (Figure 16). The reported accounts [21–24,29] based on field validation, velocity of 60 m/s
for this event,
air surveys, andwhile Jianginterpretation
the visual et al. [22] modelled the
of satellite maximum
images suggest flow
thatheight of distance
the total 100 with
the maximum
between SF2 site velocity
and the of end
65 m/s.
of theThe maximum
rock-ice pressure
avalanche reacheszone
deposition ~6000 kPa,
was ~13while
km (~2 the
momentum reaches ~4000 −1 in the GG’ transect (Figure 16). Based on the modelled
km of mountain slope andkgms
~11 km of river valley). From this point, the deposits started
mechanical
becoming more characteristics,
fluidised and SF2turned
as a rock-ice avalanche
into a debris floodcan
[21].beThe
divided
same into five prominent
published sources
phases before
[21–24,29], usingit turns
seismic into a more fluidised
modelling debris flow
and videographic in Dhauliganga
evidence provided by River
the (Figure
local pop- 16):
(1) slopealso
ulation, failure
put aand slidingon
constraint asthe
slope runout
timeline (initial
of this ~1800
event. m of transect
The initial GG’,occurred
slope failure i.e., GA);
(2) impactIndian
at 10:21 and shearing
Standard(~1800–~2400
Time (IST) (4:51 m ofGreenwich
transect GG’, Mean i.e.,Time
AB); (3) entrainment
(GMT)) and it tookfrom
valley
slightlysediments
less than a(~2400–~5000
minute to hit the m of transect
valley GG’, i.e., BC);
[21]. PlanetScope (4) liquefaction
image taken at 5:01:51 andGMT rock-
ice avalanche
(Figure 13), shows frontthepropagation
SF2 front atasthe valley
exactrunout (~5000–~8900
bend before Raini villagem of transect
where GG’, i.e.,
the rock-ice
CD); and (5)
avalanche deposition,
started turningdetrainment,
into a debris and flow.transition into athat
This indicates more SF2fluidised
rock-icedebris flood
avalanche
(~8900–~12,500
front took 10 min m toof travel
transect GG’,
to the i.e., DG’).point. However, this transition bend was also
transition
one of the sharpest and narrowest along the valley, stopping the runout and depositing
 sistency in the spatiotemporal premise for the final outcomes gave us the needed confi-
dence in the simulated mechanical parameters for SF2 (Figures 15 and 16). The improving
satellite remote sensing capabilities, such as the temporal resolutions, help to study natu-
ral hazards more efficiently, as in this case too, PlanetScope was able to capture the first
Remote Sens. 2022, 14, 949 22 of 32
image within ~10 min of the event, giving us a reliable timestamp for calibrating the
model.

Figure 15. 3D perspective views of the mechanical parameters defining SF2: (a) height, (b) velocity,
Figure 15. 3D perspective views of the mechanical parameters defining SF2: (a) height, (b) velocity,
(c) pressure, and (d) momentum.
(c) pressure, and (d) momentum.
 phases before it turns into a more fluidised debris flow in Dhauliganga River (Figure 16):
(1) slope failure and sliding as slope runout (initial ~1800 m of transect GG’, i.e., GA); (2)
impact and shearing (~1800–~2400 m of transect GG’, i.e., AB); (3) entrainment from valley
sediments (~2400–~5000 m of transect GG’, i.e., BC); (4) liquefaction and rock-ice ava-
lanche front propagation as valley runout (~5000–~8900 m of transect GG’, i.e., CD); and
Remote Sens. 2022, 14, 949 (5) deposition, detrainment, and transition into a more fluidised debris flood (~8900– 23 of 32
~12,500 m of transect GG’, i.e., DG’).

Figure 16. Mechanical characteristics of SF2 during the event, revealing its distinct phases: the left
panel (a) shows
showsthethetransect
transectofofSF2
SF2along which
along various
which variousmechanical characteristics
mechanical areare
characteristics plotted in the
plotted in
right panel.
the right TheThe
panel. knickpoints where
knickpoints the the
where mass movement
mass movementtransformed its mechanical
transformed properties
its mechanical are
properties
marked by dashed
are marked white
by dashed lineslines
white perpendicular to thetorunout.
perpendicular The middle
the runout. panel (b)
The middle shows
panel the differ-
(b) shows the
ence DTM highlighting the obliterated spur (within yellow polygon) during SF2 propagation, as
difference DTM highlighting the obliterated spur (within yellow polygon) during SF2 propagation,
shown by the red ellipse in Figure 12. The red dashed line in the right panel (c) shows the dip in all
as shown by the red ellipse in Figure 12. The red dashed line in the right panel (c) shows the dip in all
the observed (i.e., difference DTM) and simulated profiles, corresponding to the breakage of this spur.
The transects in the deposition zone, i.e., DG’, provide the spatial context for Figure 17. Numbers
in the top panel of (c) represent various mechanical zones of SF2, as explained in the text. Please
note that MFH and MFV, and MFM and MFP are of the same order of magnitude, and are therefore
plotted on the same axis using the same scale in their respective plots. (MFH = Maximum Flow
Height, MFV = Maximum Flow Velocity, MFM = Maximum Flow Momentum, MFP = Maximum
Flow Pressure, and MSS = Maximum Shear Stress).
 Video S2) confirms that while the avalanche was airborne briefly, it soon became
grounded and that could explain the obliteration of this spur. We perceive that the pres-
ence of this spur offered a topographic barrier to the lowermost part of the eastern glacier,
causing a slight pause in SF2 propagation (hence the dip in the mechanical parameters),
containing it within a channel, and becoming obliterated in turn. This topographic barrier
Remote Sens. 2022, 14, 949 24 of 32
and friction also caused some deformation and possible breakage of the SF2 mass, as is
evident from the spikes in the shear stress profile in Figure 16.

Figure17.17.Depositional
Figure Depositionalprofiles
profilescompared
comparedtotoDTM
DTMdifference
differencebetween
between2021
2021and
and 2015
2015 DTMs.The
DTMs. The
yellow
yellow curve
curve representing
representing the
the modelled
modelled deposition
deposition height
height along
along the
the channel
channel shows
shows the
the terrainright
terrain right
afterSF2
after SF2onon 7 February,
7 February, while
while thethe DTM
DTM difference
difference curve
curve represents
represents thethe terrain
terrain onon
11 11 February.
February.

In the first phase of the landslide, the detached mass travelled ~1800 m before it struck
the valley. Our model results highlight an interesting and previously unreported aspect of
this phase. At nearly 750 m below the terminus of the SF2 hanging glacier, we observed a
dip in all the mechanical parameters, and, remarkably, we could also observe a dip in the
difference DTM, highlighting the removal of significant rock volume by the SF2 rock-ice
mass before it struck the valley (Figure 16). SF2 majorly consisted of rock volume and
started as a landslide, and, therefore, when the large rock body hit the spur, it eroded
it. The red ellipse in Figure 12 highlights this spur in the path of SF2 and the middle
panel of Figure 16 shows the difference DTM highlighting the obliterated spur during SF2
propagation. This spur had a height range of ~4–74 m with an average height of ~25 m;
The total area covered by it was ~31,078 m2 and the total volume of this spur removed by
SF2 was ~795,810 m3 . This observation signified two aspects of this phase: (1) the model
worked significantly well in predicting a hindrance in the flow, leading to the possible
removal of this spur owing to the mechanical characteristics (high momentum and power)
of the landslide, and (2) a significant rock volume of ~0.8 × 106 m3 was further added to
the avalanched volume of ~27 × 106 m3 even before it reached the valley. Shugar et al. [21]
postulated that the reason behind the intact lowermost part of the neighbouring larger
eastern glacier post-SF2, might be the possibility that the avalanche became airborne for
a while. A 3D visualization of our modelling result (Supplementary Materials Video S2)
confirms that while the avalanche was airborne briefly, it soon became grounded and that
could explain the obliteration of this spur. We perceive that the presence of this spur offered
a topographic barrier to the lowermost part of the eastern glacier, causing a slight pause
in SF2 propagation (hence the dip in the mechanical parameters), containing it within
a channel, and becoming obliterated in turn. This topographic barrier and friction also
caused some deformation and possible breakage of the SF2 mass, as is evident from the
spikes in the shear stress profile in Figure 16.
In the impact and shearing phase, the SF2 mass hit the valley floor and all the mechan-
ical parameters reached closer to their peak values during this phase. Immediately after
the impact, the maximum shear stress reaches its peak (~1200 kPa) and this zone of ~600 m
signifies the shearing and fragmentation of the rock and ice core (Figure 16). While the
ice-debris from SF1 in 2016 mostly remains within ~3 km of the deposition zone, covering
the main channel and surroundings, the SF2 event, starting from slightly steeper slopes,
releases over ~2.5 times more of the rock-ice volume on the valley floor with previously
deposited avalanche debris, significant seasonal snow (Figure 6b), and possibly dammed
water ponds (Figures 4 and 5). This impact could have led to considerable ice and snow
Remote Sens. 2022, 14, 949 25 of 32

melt through the substantial change in momentum (Figure 16), leading to the transmission
of energy by the falling and sliding rock and ice along the stretch. This immense rate of
change in momentum and power, once the avalanche hit the valley, also caused the large
burst whose bang could reportedly be heard from as far as ~15 km from the avalanche
site, and was later identified in the seismic readings from a station located over 100 km
away [21]. For a similar event at Kolka in the North Ossetia–Alania region of the Caucasus
Mountains in 2002, where a sliding hanging glacier of ~20 × 106 m3 failure volume hit and
eroded a small glacier, a broad estimate, considering that the substantial amount of energy
was dissipated during the fall and impact, suggested that ~0.44 × 106 m3 of ice could be
instantly melted by a falling volume of 8 × 106 m3 [79]. Considering that SF2 had a total
volume of over ~27 × 106 m3 and it travelled an almost similar drop height as Kolka, it
could have melted at least 3 times of ~0.44 × 106 m3 of ice instantly, i.e., ~1.32 × 106 m3 of
ice, corresponding to the ~1/4th of the total SF2 ice volume, due to the direct impact and
frictional sliding. This can explain why in contrast to most previously documented rock
avalanches, very little debris was found at the base of the failed slope [21]. The Kolka event
too, similar to SF2, did not allow the fallen ice and debris to deposit in the vicinity of the
fall, but led to an unprecedented runout distance of ~20 km, followed by flash mudflows.
In such voluminous avalanches, granular friction can effectively be reduced by ice melt
and associated fluidisation effects, leading to much larger travel distances than usual [79].
The next phase of SF2 consisted of entrainment from the valley sediments, mostly
within the ~2400–~5000 m of transect GG’, i.e., in the BC region (Figure 16). Point C in
Figure 16 marks the sharp turn in the channel where previous ice avalanches of 2000 and
2016 from SF1 site had stopped. While the majority of the ice-debris deposit from SF1
appears to have melted away before the SF2 event, we cannot ignore the enormity of SF1
deposits (Figures 4, 5 and 8), when the maximum deposition height reached as high as
~48 m. The SF1 event had not just altered the valley profile, but had also enriched it with
significant volume of sediments and percolated meltwater in subsequent years. A mystery
in explaining the long runout of SF2 has been that how only 20% of ice volume could
displace nearly the entire 80% of rock volume, so far down the valley [29]. We propose
that the long runout of SF2 was facilitated by the considerable entrainment from the valley
sediments of past ice and snow avalanches (Figures 4 and 5) during this phase of SF2 in the
BC transect (Figure 16). To evaluate this hypothesis, we modelled two scenarios: (1) with
the inclusion of differential erosion areas and erosion depth based on the modelled results
of Jian et al. [22], and (2) without specifying the erosion areas and depth. By default, the
model includes a critical shear stress that permits erosion only when the shear stress in
any given cell exceeds the critical shear stress value for the onset of erosion. However,
with the available information, i.e., in this case, the extent of SF1 deposits and modelled
erosion by Jiang et al. [22], the inclusion of defined erosion zones can improve the model
accuracy. RAMMS also provides the option to input erosion area parameters, and these
were appropriately adjusted as per the instructions in RAMMS manual. Scenario 2 of
not specifying the erosion areas and depth also meant that we had to consider a smaller
turbulent friction (ε) parameter. While, for Scenario 1, after the model calibration, we opted
for a value of ε = 200 m/s2 , representing a transition from a solid-dominated granular
flow to a muddier flow owing to sediment entrainment, for Scenario 2, we chose a value
of ε = 100 m/s2 , representing a relatively drier solid-dominated flow. As we expected,
compared to Scenario 1, the rock-ice avalanche in Scenario 2 took nearly twice the time
with an unrealistically low avalanche front momentum (only the ~1/3rd of Scenario 1), for
reaching point G’ in Figure 16. This supports our assumption that the role of past avalanche
deposits in entraining SF2 was considerably significant.
The subsequent phase of SF2 between ~5000 m and ~8900 m of transect GG’ (i.e., CD)
does not only show a relatively constricted valley profile, but is also marked by dense
vegetation along the channel, which offered considerable resistance to SF2 propagation. The
rapidly changing shear stress, flow pressure, and maximum flow height values in Figure 16
highlight this aspect. This friction led to further liquefaction of the avalanche core and the
Remote Sens. 2022, 14, 949 26 of 32

front took the shape of a well-defined rock-ice avalanche between this ~5000–~8900 m of
transect GG’, i.e., in section CD. Point D represents another sharp turn in the valley and
allowed for the temporary accumulation of material, leading to rise in the flow heights
(Figure 16).
The last phase of SF2 is marked by the deposition of transported material between
~8900 and ~12,500 m of transect GG’, i.e., section DG’. Transect HH’ in Figure 16 repre-
sents the place where Ronti Gad meets Rishiganga River and, from this point onwards, a
vast volume of sediments and rocks were deposited, temporarily damming the river and
forming a ~700 m long lake [21]. A deposition analysis in RAMMS resulted in the total
deposition volume of ~15 × 106 m3 . The maximum flow height reached up to ~120 m.
While Shugar et al. [21] also reported the maximum deposition height in the same range,
based on DTM differencing, they reported the total deposition volume of ~8 × 106 m3 ,
between the Ronti Gad–Rishiganga River confluence and point G’ in Figure 16. There
are two reasons why our deposition volume is larger: (1) we observed and quantified the
deposition on a larger stretch of transect GG’, i.e., DG’, and (2) our deposition estimates
are in real-time, representing the terrain just after the rock-ice avalanche stopped at point
G’. Figure 17 clarifies this point further. The 2021 post-SF2 DTM was generated using
stereopairs acquired on 10 and 11 February and, by then, a significant volume of deposits
fluidised and escaped as the debris flood. After 24 h of the event, the sediment plume could
be seen on satellite images ~150 km downstream of the source, and ~80 times of the normal
turbidity could be observed ~500 km downstream within 8 days of the event [21]. This
means that the deposited avalanche material continuously kept on fluidising, contributing
to a constant, albeit gradually diminishing, debris flow. To evaluate our modelled depo-
sition results, we hypothesised that while DTM difference elevations would show lower
deposition heights, considering that they were generated using 10/11 February stereopairs,
the profile pattern should match our modelled deposition heights. When we performed
the profile analysis along four transects in the deposition region (Figure 16), we found that
there was a remarkable similarity between the profile patterns between the modelled curve
representing the deposition height immediately after the SF2 event on 7 February and the
DTM difference curve representing the terrain on 11 February (Figure 17). This further
proved the authenticity of our simulated results.
The SF2 event could be a result of the combined contribution of both the gradually
destabilising slope and a deformational stress building up due to the winter snow accu-
mulation load. A recent paper [47] suggested a possible polythermal regime for glaciers
north of the study area, promoting both sliding and deformational forms of motion. Both
pressure melting and creep-rate enhancement driven by stress concentration act as prime
factors in producing noticeable surface motion under polythermal regime [80]. An ap-
proximate range for the existence of cold, polythermal, and temperate ice based on the
relation between the mean annual air temperature (MAAT) and critical slope for failure
on ramp-type glaciers [81] suggests that, for the average slope of SF2 site and assuming
the polythermal nature of ice, a slope failure is possible within an MAAT range of −2.5 ◦ C
to −5 ◦ C. The MAAT for SF1 and SF2 corresponds to ~−5 ◦ C, based on the published
data for Chamoli district [82], extrapolated to the SF2 site using a constant environmental
lapse rate of −6.5 ◦ C [83]. This means that such slope failures can become increasingly
possible, irrespective of the seasons, with changing climatic conditions. In their review
on avalanching glacier instabilities, Faillettaz et al. [84] identified three different types of
instabilities depending on the thermal properties of the ice/bedrock interface. If cold, the
glacier experiences a critical acceleration following the power-law up to the break-off, and
a prediction of the final break-off is possible. The type of motion observed in case of SF1 is
closer to a possible cold regime, indicating a possibility of predicting its future break-offs.
In case of a temperate regime, water plays a key role in the development of instability, and
critical conditions promoting the final instability can still be identified [84]. However, in the
case of a partly temperate or polythermal regime, predicting the actual break-off becomes
difficult. While SF1 and SF2 are different in their mechanical characteristics, in both these
Remote Sens. 2022, 14, 949 27 of 32

cases, we observe certain abrupt but characteristically different changes in SVs, which
raises another research question concerning whether their thermal regimes are different
despite being in proximity, or if they experience a transition/evolution in their thermal
regimes. The temporal evolution observed in the case of SF2 might represent the instability
in a glacier experiencing a rapid transition from a cold to a temperate glacier bed due
to percolation through the detached headwall. This possibility makes our SV estimates
relevant to possibly open-up a new research frontier on understanding the transition in
thermal regimes of hanging glaciers in these mountains, and how this transition might
make them more prone to slope failures.

5. Conclusions
This research started with targeting three key research questions (Figure 2) mentioned
in the Introduction Section, and we successfully derived several useful conclusions for them.
The change in satellite-derived SVs should further be explored as an important precursor
of ice-rock avalanching, and its temporal observation in high mountains can prove to
be extremely helpful for predicting several SF1-type future slope failure events. Our SV
estimates prove that continuously improving temporal and spatial resolution of satellite
datasets can be effectively employed in performing such observations at unprecedented
spatiotemporal scales. We also observed how SV changes vary between two constituently
different slope failures, i.e., SF1 and SF2. In the case of the landslide causing ice-bedrock
failure, i.e., SF2, SV anomalies developed and disappeared gradually and could be observed
even at yearly scales. However, that makes predicting the actual time of failure difficult.
The observed increase in SVs days-to-weeks prior to SF2 should be interpreted with caution,
as relatively higher uncertainties are associated with them, due to high snow cover and
poor illumination in satellite images. More of similar case studies are needed to infer if
other SF2-type events also experience an immediate surge in SVs prior to their occurrence.
On the other hand, for the frontal block failure, i.e., SF1, the SV anomalies are more localised
and pronounced closer to the time of slope failure, as the glacier attains its critical geometry
(Figure 18). We further performed calibrated RAMMS simulations to infer the differences in
the mechanical characteristics of the two large slope failures, i.e., SF1 and SF2, originating
from the same release height and hitting the same valley floor. The simulated mechanical
parameters were able to describe the during-event flow characteristics of both the events,
and the changes in these parameters further helped in identifying the various zones of
these mass movements.
Although the past events at SF1 site had not inflicted any direct damage to life and
property, it is difficult to predict their future impacts in combination with other glacial
hazards, such as SF2. While SF2 is an ice-bedrock failure, SF1 shows the characteristics of a
frontal block failure as it tends to periodically replicate itself. More importantly, we can
observe new ice accumulation at the avalanched slope of SF1 site since its last failure in
2016 (Figure 18). The hanging glacier is recovering its previous critical geometry before
SF1, and ~0.052 km2 (~20% of the total lost area in SF1) was recovered by October 2020
(Figure 18). Based on the available high-resolution images on Google Earth, we can observe
that the full recovery usually takes 12–15 years, as happened after the 2000 event at SF1
site. Our results provide the first-hand account of noticeable SV changes days before SF1,
and indicate that temporal SV monitoring of hanging glaciers might help to predict similar
future events However, it is more important to first develop statistically robust “anomaly”
characterisation routines for a variety of ice avalanche events to identify the cases with
the most potential to be predicted through remote sensing. This chain of events clearly
highlights the complexity of high mountain hazards, where considering the hazardous
nature of an event in isolation is no longer sufficient. Instead, there is an increasing need
to take into account the antecedent conditions, while making a holistic assessment of any
high mountain hazard.
 similar future events However, it is more important to first develop statistically robust
“anomaly” characterisation routines for a variety of ice avalanche events to identify the
cases with the most potential to be predicted through remote sensing. This chain of events
clearly highlights the complexity of high mountain hazards, where considering the haz-
Remote Sens. 2022, 14, 949
ardous nature of an event in isolation is no longer sufficient. Instead, there is an increasing
28 of 32
need to take into account the antecedent conditions, while making a holistic assessment
of any high mountain hazard.

Figure 18.
Figure 18. Landsat
Landsat 77 and
and Sentinel-2
Sentinel-2 images
images showing
showing ice
ice loss
loss (red
(red polygon)
polygon) and
and subsequent
subsequent gain
gain
(shaded blue polygon) for 2000 and 2016 events at SF1 site. Landsat 7 image IDs:
(shaded blue polygon) for 2000 and 2016 events at SF1 site. Landsat 7 image IDs: LE71450392000275-
LE71450392000275SGS00 and LE71450391999288EDC00; Sentinel-2 image IDs:
SGS00 and LE71450391999288EDC00; Sentinel-2 image IDs: S2A_OPER_MSI_L1C_TL_SGS__201609-
S2A_OPER_MSI_L1C_TL_SGS__20160919T104034_A006494_T44RKU,
19T104034_A006494_T44RKU, S2A_OPER_MSI_L1C_TL_SGS__20161009T053631_A006780_T44RLU,
S2A_OPER_MSI_L1C_TL_SGS__20161009T053631_A006780_T44RLU, and
and S2A_OPER_MSI_L1C_TL_VGS1_20201018T072817_A027801_T44RLU.
S2A_OPER_MSI_L1C_TL_VGS1_20201018T072817_A027801_T44RLU.
Supplementary Materials: The simulated Graphics Interchange Format (.gif) files for SF1 and SF2
can be downloaded at: https://www.mdpi.com/article/10.3390/rs14040949/s1, Video S1: Simulated
SF1 event, Video S2: Simulated SF2 event.
Author Contributions: Conceptualisation, A.B. and L.S.; methodology, A.B. and L.S.; software,
A.B. and L.S.; validation, A.B. and L.S.; formal analysis, A.B. and L.S.; investigation, A.B. and L.S.;
resources, A.B. and L.S.; writing—original draft preparation, A.B.; writing—review and editing, L.S.;
visualisation, A.B. and L.S.; project administration, A.B. and L.S.; funding acquisition, A.B. and L.S.
All authors have read and agreed to the published version of the manuscript.
Remote Sens. 2022, 14, 949 29 of 32

Funding: This research was funded from the Interdisciplinary Pump Priming Fund (grants nos.:
SF10237-19 and SF10206-67) and Global Challenge Research Internal Funds (grant no.: SF10206-78)
granted by the University of Aberdeen, U.K., and the Scottish Funding Council.
Institutional Review Board Statement: Not Applicable.
Informed Consent Statement: Not Applicable.
Data Availability Statement: All the data used in this study is available free-of-cost in open ac-
cess and the data sources have been cited. The simulation graphics have been provided as the
Supplementary Materials.
Acknowledgments: We acknowledge NASA, USGS, ESA, Planet Labs, and previous studies [36,40]
for providing free-of-cost medium-to-high resolution satellite images and DTMs. We acknowledge
the support provided by RAMMS Team at the Swiss Institute for Snow and Avalanche Research,
Davos, Graubünden, in offering the RAMMS tool. We thank the reviewers for their constructive
suggestions, which improved the quality of this paper.
Conflicts of Interest: The authors declare no conflict of interest.

References
1. Froude, M.J.; Petley, D.N. Global Fatal Landslide Occurrence from 2004 to 2016. Nat. Hazards Earth Syst. Sci. 2018, 18, 2161–2181.
[CrossRef]
2. Petley, D.N.; Hearn, G.J.; Hart, A.; Rosser, N.J.; Dunning, S.A.; Oven, K.; Mitchell, W.A. Trends in Landslide Occurrence in Nepal.
Nat. Hazards 2007, 43, 23–44. [CrossRef]
3. Richardson, S.D.; Reynolds, J.M. An Overview of Glacial Hazards in the Himalayas. Quat. Int. 2000, 65, 31–47. [CrossRef]
4. Haeberli, W.; Schaub, Y.; Huggel, C. Increasing Risks Related to Landslides from Degrading Permafrost into New Lakes in
De-Glaciating Mountain Ranges. Geomorphology 2017, 293, 405–417. [CrossRef]
5. Huggel, C. Recent Extreme Slope Failures in Glacial Environments: Effects of Thermal Perturbation. Quat. Sci. Rev. 2009, 28,
1119–1130. [CrossRef]
6. Huggel, C.; Clague, J.J.; Korup, O. Is Climate Change Responsible for Changing Landslide Activity in High Mountains? Earth
Surf. Processes Landf. 2012, 37, 77–91. [CrossRef]
7. Kirschbaum, D.; Kapnick, S.B.; Stanley, T.; Pascale, S. Changes in Extreme Precipitation and Landslides over High Mountain Asia.
Geophys. Res. Lett. 2020, 47, e2019GL085347. [CrossRef]
8. Fischer, L.; Huggel, C.; Kääb, A.; Haeberli, W. Slope Failures and Erosion Rates on a Glacierized High-Mountain Face under
Climatic Changes. Earth Surf. Processes Landf. 2013, 38, 836–846. [CrossRef]
9. Martinsen, O. Mass Movements. In The Geological Deformation of Sediments; Maltman, A., Ed.; Springer: Dordrecht,
The Netherlands, 1994; pp. 127–165; ISBN 978-94-011-0731-0.
10. Lee, S.H.-H.; Widjaja, B. Phase Concept for Mudflow Based on the Influence of Viscosity. Soils Found. 2013, 53, 77–90. [CrossRef]
11. Abe, K.; Konagai, K. Numerical Simulation for Runout Process of Debris Flow Using Depth-Averaged Material Point Method.
Soils Found. 2016, 56, 869–888. [CrossRef]
12. Takahashi, T. Debris Flow. Annu. Rev. Fluid Mech. 1981, 13, 57–77. [CrossRef]
13. Bartelt, P.; Buser, O.; Vera Valero, C.; Bühler, Y. Configurational Energy and the Formation of Mixed Flowing/Powder Snow and
Ice Avalanches. Ann. Glaciol. 2016, 57, 179–188. [CrossRef]
14. Schaub, Y.; Huggel, C.; Cochachin, A. Ice-Avalanche Scenario Elaboration and Uncertainty Propagation in Numerical Simulation
of Rock-/Ice-Avalanche-Induced Impact Waves at Mount Hualcán and Lake 513, Peru. Landslides 2016, 13, 1445–1459. [CrossRef]
15. Wang, Z.-Y.; Lee, J.H.; Melching, C.S. Debris Flows and Landslides. In River Dynamics and Integrated River Management; Springer:
Berlin/Heidelberg, Germany, 2015; pp. 193–264. [CrossRef]
16. Kirschbaum, D.; Stanley, T.; Zhou, Y. Spatial and Temporal Analysis of a Global Landslide Catalog. Geomorphology 2015, 249, 4–15.
[CrossRef]
17. Sam, L.; Bhardwaj, A.; Sinha, V.S.P.; Joshi, P.K.; Kumar, R. Use of Geospatial Tools to Prioritize Zones of Hydro-Energy Potential
in Glaciated Himalayan Terrain. J. Indian Soc. Remote Sens. 2016, 44, 409–420. [CrossRef]
18. Snehmani; Bhardwaj, A.; Pandit, A.; Ganju, A. Demarcation of Potential Avalanche Sites Using Remote Sensing and Ground
Observations: A Case Study of Gangotri Glacier. Geocarto Int. 2014, 29, 520–535. [CrossRef]
19. Snehmani; Bhardwaj, A.; Singh, M.K.; Gupta, R.D.; Joshi, P.K.; Ganju, A. Modelling the Hypsometric Seasonal Snow Cover Using
Meteorological Parameters. J. Spat. Sci. 2015, 60, 51–64. [CrossRef]
20. Kääb, A.; Jacquemart, M.; Gilbert, A.; Leinss, S.; Girod, L.; Huggel, C.; Falaschi, D.; Ugalde, F.; Petrakov, D.; Chernomorets,
S.; et al. Sudden Large-Volume Detachments of Low-Angle Mountain Glaciers—More Frequent than Thought? Cryosphere 2021,
15, 1751–1785. [CrossRef]
Remote Sens. 2022, 14, 949 30 of 32

21. Shugar, D.H.; Jacquemart, M.; Shean, D.; Bhushan, S.; Upadhyay, K.; Sattar, A.; Schwanghart, W.; McBride, S.; de Vries, M.V.W.;
Mergili, M. A Massive Rock and Ice Avalanche Caused the 2021 Disaster at Chamoli, Indian Himalaya. Science 2021, 373, 300–306.
[CrossRef]
22. Jiang, R.; Zhang, L.; Peng, D.; He, X.; He, J. The Landslide Hazard Chain in the Tapovan of the Himalayas on 7 February 2021.
Geophys. Res. Lett. 2021, 48, e2021GL093723. [CrossRef]
23. Martha, T.R.; Roy, P.; Jain, N.; Kumar, K.V.; Reddy, P.S.; Nalini, J.; Sharma, S.; Shukla, A.K.; Rao, K.D.; Narender, B. Rock Avalanche
Induced Flash Flood on 07 February 2021 in Uttarakhand, India—A Photogeological Reconstruction of the Event. Landslides 2021,
18, 2881–2893. [CrossRef]
24. Mao, W.; Wu, L.; Singh, R.P.; Qi, Y.; Xie, B.; Liu, Y.; Ding, Y.; Zhou, Z.; Li, J. Progressive Destabilization and Triggering Mechanism
Analysis Using Multiple Data for Chamoli Rockslide in 7 February 2021. Geomat. Nat. Hazards Risk 2022, 13, 35–53. [CrossRef]
25. Christen, M.; Bühler, Y.; Bartelt, P.; Leine, R.; Glover, J.; Schweizer, A.; Graf, C.; McArdell, B.W.; Gerber, W.; Deubelbeiss, Y.; et al.
Integral hazard management using a unified software environment: Numerical simulation tool “RAMMS” for gravitational
natural hazards. In Proceedings of the 12th Congress Interpraevent, Grenoble, France, 23–26 April 2012; Koboltschnig, G., Hübl,
J., Braun, J., Eds.; 2012; pp. 77–86.
26. Pandey, P.; Chauhan, P.; Bhatt, C.M.; Thakur, P.K.; Kannaujia, S.; Dhote, P.R.; Roy, A.; Kumar, S.; Chopra, S.; Bhardwaj, A. Cause
and Process Mechanism of Rockslide Triggered Flood Event in Rishiganga and Dhauliganga River Valleys, Chamoli, Uttarakhand,
India Using Satellite Remote Sensing and in Situ Observations. J. Indian Soc. Remote Sens. 2021, 49, 1011–1024. [CrossRef]
27. Shekhar, M.; Bhardwaj, A.; Singh, S.; Ranhotra, P.S.; Bhattacharyya, A.; Pal, A.K.; Roy, I.; Martín-Torres, F.J.; Zorzano, M.-P.
Himalayan Glaciers Experienced Significant Mass Loss during Later Phases of Little Ice Age. Sci. Rep. 2017, 7, 10305. [CrossRef]
[PubMed]
28. Kumar, V.; Mehta, M.; Mishra, A.; Trivedi, A. Temporal Fluctuations and Frontal Area Change of Bangni and Dunagiri Glaciers
from 1962 to 2013, Dhauliganga Basin, Central Himalaya, India. Geomorphology 2017, 284, 88–98. [CrossRef]
29. Mehta, M.; Kumar, V.; Sain, K.; Tiwari, S.K.; Kumar, A.; Verma, A. Causes and Consequences of Rishiganga Flash Flood, Nanda
Devi Biosphere Reserve, Central Himalaya, India. Curr. Sci. 2021, 121, 1483–1487.
30. Mauro, Y. Geology and Metamorphism of the Nanda Devi Region, Kumaun Higher Himalaya. Himal. Geol. 1979, 9, 3–17.
31. Valdiya, K.S.; Goel, O.P. Lithological Subdivision and Petrology of the Great Himalayan Vaikrita Group in Kumaun, India. Proc.
Indian Acad. Sci.-Earth Planet. Sci. 1983, 92, 141–163. [CrossRef]
32. Sahoo, P.K.; Kumar, S.; Singh, R.P. Neotectonic Study of Ganga and Yamuna Tear Faults, NW Himalaya, Using Remote Sensing
and GIS. Int. J. Remote Sens. 2000, 21, 499–518. [CrossRef]
33. McColl, S.T. Paraglacial Rock-Slope Stability. Geomorphology 2012, 153, 1–16. [CrossRef]
34. Massive Flood as Glacier Breaks Off at Uttarakhand’s Joshimath, 150 Labourers Missing. Available online: https://www.youtube.
com/watch?v=DoWivEFpbsE (accessed on 23 January 2022).
35. ICIMOD. Major River Basins in the Hindu Kush Himalaya (HKH) Region [Data Set]; ICIMOD: Lalitpur, Nepal, 2021. [CrossRef]
36. Bhushan, S.; Shean, D. Chamoli Disaster Pre-Event 2-m DEM Composite: September 2015 (1.0) [Data Set]; Zenodo: Brussels,
Belgium, 2021. [CrossRef]
37. Shean, D.E.; Alexandrov, O.; Moratto, Z.M.; Smith, B.E.; Joughin, I.R.; Porter, C.; Morin, P. An Automated, Open-Source Pipeline
for Mass Production of Digital Elevation Models (DEMs) from Very-High-Resolution Commercial Stereo Satellite Imagery. ISPRS
J. Photogramm. Remote Sens. 2016, 116, 101–117. [CrossRef]
38. Shean, D.E.; Bhushan, S.; Montesano, P.; Rounce, D.R.; Arendt, A.; Osmanoglu, B. A Systematic, Regional Assessment of High
Mountain Asia Glacier Mass Balance. Front. Earth Sci. 2020, 7, 363. [CrossRef]
39. Deschamps-Berger, C.; Gascoin, S.; Berthier, E.; Deems, J.; Gutmann, E.; Dehecq, A.; Shean, D.; Dumont, M. Snow Depth Mapping
from Stereo Satellite Imagery in Mountainous Terrain: Evaluation Using Airborne Laser-Scanning Data. Cryosphere 2020, 14,
2925–2940. [CrossRef]
40. Shean, D.; Bhushan, S.; Berthier, E.; Deschamps-Berger, C.; Gascoin, S.; Knuth, F. Chamoli Disaster Post-Event 2-m DEM Composite
(10–11 February 2021) and Difference Map (1.0) [Data Set]; Zenodo: Brussels, Belgium, 2021. [CrossRef]
41. Leprince, S.; Barbot, S.; Ayoub, F.; Avouac, J.-P. Automatic and Precise Orthorectification, Coregistration, and Subpixel Correlation
of Satellite Images, Application to Ground Deformation Measurements. IEEE Trans. Geosci. Remote Sens. 2007, 45, 1529–1558.
[CrossRef]
42. Scherler, D.; Leprince, S.; Strecker, M.R. Glacier-Surface Velocities in Alpine Terrain from Optical Satellite Imagery—Accuracy
Improvement and Quality Assessment. Remote Sens. Environ. 2008, 112, 3806–3819. [CrossRef]
43. Lacroix, P.; Bièvre, G.; Pathier, E.; Kniess, U.; Jongmans, D. Use of Sentinel-2 Images for the Detection of Precursory Motions
before Landslide Failures. Remote Sens. Environ. 2018, 215, 507–516. [CrossRef]
44. Sam, L.; Bhardwaj, A.; Singh, S.; Kumar, R. Remote Sensing Flow Velocity of Debris-Covered Glaciers Using Landsat 8 Data. Prog.
Phys. Geogr. Earth Environ. 2016, 40, 305–321. [CrossRef]
45. USGS EarthExplorer. Available online: https://earthexplorer.usgs.gov/ (accessed on 10 March 2021).
46. Copernicus Open Access Hub. Available online: https://scihub.copernicus.eu/dhus/#/home (accessed on 10 March 2021).
47. Sam, L.; Bhardwaj, A.; Kumar, R.; Buchroithner, M.F.; Martín-Torres, F.J. Heterogeneity in Topographic Control on Velocities of
Western Himalayan Glaciers. Sci. Rep. 2018, 8, 12843. [CrossRef]
Remote Sens. 2022, 14, 949 31 of 32

48. Garg, P.K.; Shukla, A.; Jasrotia, A.S. On the Strongly Imbalanced State of Glaciers in the Sikkim, Eastern Himalaya, India. Sci.
Total Environ. 2019, 691, 16–35. [CrossRef]
49. Heid, T.; Kääb, A. Repeat Optical Satellite Images Reveal Widespread and Long Term Decrease in Land-Terminating Glacier
Speeds. Cryosphere 2012, 6, 467–478. [CrossRef]
50. Scherler, D.; Strecker, M.R. Large Surface Velocity Fluctuations of Biafo Glacier, Central Karakoram, at High Spatial and Temporal
Resolution from Optical Satellite Images. J. Glaciol. 2012, 58, 569–580. [CrossRef]
51. Dreier, L.; Bühler, Y.; Ginzler, C.; Bartelt, P. Comparison of Simulated Powder Snow Avalanches with Photogrammetric Measure-
ments. Ann. Glaciol. 2016, 57, 371–381. [CrossRef]
52. Frank, F.; McArdell, B.W.; Huggel, C.; Vieli, A. The Importance of Entrainment and Bulking on Debris Flow Runout Modeling:
Examples from the Swiss Alps. Nat. Hazards Earth Syst. Sci. 2015, 15, 2569–2583. [CrossRef]
53. Margreth, S.; Funk, M.; Tobler, D.; Dalban, P.; Meier, L.; Lauper, J. Analysis of the Hazard Caused by Ice Avalanches from the
Hanging Glacier on the Eiger West Face. Cold Reg. Sci. Technol. 2017, 144, 63–72. [CrossRef]
54. Rodríguez-Morata, C.; Villacorta, S.; Stoffel, M.; Ballesteros-Cánovas, J.A. Assessing Strategies to Mitigate Debris-Flow Risk in
Abancay Province, South-Central Peruvian Andes. Geomorphology 2019, 342, 127–139. [CrossRef]
55. Gilany, N.; Iqbal, J. Simulation of Glacial Avalanche Hazards in Shyok Basin of Upper Indus. Sci. Rep. 2019, 9, 20077. [CrossRef]
56. Allen, S.K.; Schneider, D.; Owens, I.F. First Approaches towards Modelling Glacial Hazards in the Mount Cook Region of New
Zealand’s Southern Alps. Nat. Hazards Earth Syst. Sci. 2009, 9, 481–499. [CrossRef]
57. Chrustek, P.; Świerk, M.; Biskupić, M. Snow Avalanche Hazard Mapping for Different Frequency Scenarios, the Case of the Tatra
Mts., Western Carpathians. In Proceedings of the 2013 International Snow Science Workshop, Grenoble–Chamonix Mont-Blanc,
France, 7–11 October 2013; pp. 745–749.
58. Bartelt, P.; Valero, C.V.; Feistl, T.; Christen, M.; Bühler, Y.; Buser, O. Modelling Cohesion in Snow Avalanche Flow. J. Glaciol. 2015,
61, 837–850. [CrossRef]
59. RAMMS Downloads. Available online: https://ramms.slf.ch/ramms/index.php?option=com_content&view=article&id=53
&Itemid=70 (accessed on 30 July 2021).
60. New Thermomechanical Model for Rock/Ice Avalanches. Available online: https://www.slf.ch/en/news/2018/04/new-
thermomechanical-model-for-rockice-avalanches.html#tabelement1-tab2 (accessed on 30 July 2021).
61. Bhardwaj, A.; Sam, L.; Martín-Torres, F.J.; Zorzano, M.-P. Are Slope Streaks Indicative of Global-Scale Aqueous Processes on
Contemporary Mars? Rev. Geophys. 2019, 57, 48–77. [CrossRef]
62. Singh, M.K.; Snehmani; Gupta, R.D.; Bhardwaj, A.; Joshi, P.K.; Ganju, A. High Resolution DEM Generation for Complex Snow
Covered Indian Himalayan Region Using ADS80 Aerial Push-Broom Camera: A First Time Attempt. Arab. J. Geosci. 2015, 8,
1403–1414. [CrossRef]
63. Smithson, S.B. Densities of Metamorphic Rocks. Geophysics 1971, 36, 690–694. [CrossRef]
64. Planet Team. Planet Application Program Interface: In Space for Life on Earth; Planet Team: San Francisco, CA, USA, 2017.
65. Fernandez, P.; Whitworth, M. A New Technique for the Detection of Large Scale Landslides in Glacio-Lacustrine Deposits Using
Image Correlation Based upon Aerial Imagery: A Case Study from the French Alps. Int. J. Appl. Earth Obs. Geoinf. 2016, 52, 1–11.
[CrossRef]
66. Peppa, M.V.; Mills, J.P.; Moore, P.; Miller, P.E.; Chambers, J.E. Brief Communication: Landslide Motion from Cross Correlation of
UAV-Derived Morphological Attributes. Nat. Hazards Earth Syst. Sci. 2017, 17, 2143–2150. [CrossRef]
67. Türk, T. Determination of Mass Movements in Slow-Motion Landslides by the Cosi-Corr Method. Geomat. Nat. Hazards Risk 2018,
9, 325–336. [CrossRef]
68. Turner, D.; Lucieer, A.; De Jong, S.M. Time Series Analysis of Landslide Dynamics Using an Unmanned Aerial Vehicle (UAV).
Remote Sens. 2015, 7, 1736–1757. [CrossRef]
69. Sam, L.; Gahlot, N.; Prusty, B.G. Estimation of Dune Celerity and Sand Flux in Part of West Rajasthan, Gadra Area of the Thar
Desert Using Temporal Remote Sensing Data. Arab. J. Geosci. 2015, 8, 295–306. [CrossRef]
70. Gili, J.A.; Corominas, J.; Rius, J. Using Global Positioning System Techniques in Landslide Monitoring. Eng. Geol. 2000, 55,
167–192. [CrossRef]
71. Flotron, A. Movement Studies on a Hanging Glacier in Relation with an Ice Avalanche. J. Glaciol. 1977, 19, 671–672. [CrossRef]
72. Bowman, D.D.; Ouillon, G.; Sammis, C.G.; Sornette, A.; Sornette, D. An Observational Test of the Critical Earthquake Concept. J.
Geophys. Res. Solid Earth 1998, 103, 24359–24372. [CrossRef]
73. Sornette, D.; Helmstetter, A.; Andersen, J.V.; Gluzman, S.; Grasso, J.-R.; Pisarenko, V. Towards Landslide Predictions: Two Case
Studies. Phys. A Stat. Mech. Its Appl. 2004, 338, 605–632. [CrossRef]
74. Amitrano, D.; Grasso, J.R.; Senfaute, G. Seismic Precursory Patterns before a Cliff Collapse and Critical Point Phenomena. Geophys.
Res. Lett. 2005, 32. [CrossRef]
75. Faillettaz, J.; Pralong, A.; Funk, M.; Deichmann, N. Evidence of Log-Periodic Oscillations and Increasing Icequake Activity during
the Breaking-off of Large Ice Masses. J. Glaciol. 2008, 54, 725–737. [CrossRef]
76. Faillettaz, J.; Funk, M.; Vagliasindi, M. Time Forecast of a Break-off Event from a Hanging Glacier. Cryosphere 2016, 10, 1191–1200.
[CrossRef]
77. Asano, Y.; Uchida, T. Detailed Documentation of Dynamic Changes in Flow Depth and Surface Velocity during a Large Flood in a
Steep Mountain Stream. J. Hydrol. 2016, 541, 127–135. [CrossRef]
Remote Sens. 2022, 14, 949 32 of 32

78. Zhang, G.; Cui, P.; Yin, Y.; Liu, D.; Jin, W.; Wang, H.; Yan, Y.; Ahmed, B.N.; Wang, J. Real-Time Monitoring and Estimation of the
Discharge of Flash Floods in a Steep Mountain Catchment. Hydrol. Processes 2019, 33, 3195–3212. [CrossRef]
79. Huggel, C.; Zgraggen-Oswald, S.; Haeberli, W.; Kääb, A.; Polkvoj, A.; Galushkin, I.; Evans, S.G. The 2002 Rock/Ice Avalanche
at Kolka/Karmadon, Russian Caucasus: Assessment of Extraordinary Avalanche Formation and Mobility, and Application of
QuickBird Satellite Imagery. Nat. Hazards Earth Syst. Sci. 2005, 5, 173–187. [CrossRef]
80. Weertman, J. On the Sliding of Glaciers. J. Glaciol. 1957, 3, 33–38. [CrossRef]
81. Huggel, C.; Haeberli, W.; Kääb, A.; Bieri, D.; Richardson, S. An Assessment Procedure for Glacial Hazards in the Swiss Alps. Can.
Geotech. J. 2004, 41, 1068–1083. [CrossRef]
82. Yadav, R.; Tripathi, S.K.; Pranuthi, G.; Dubey, S.K. Trend Analysis by Mann-Kendall Test for Precipitation and Temperature for
Thirteen Districts of Uttarakhand. J. Agrometeorol. 2014, 16, 164.
83. Barry, R.G. Mountain Weather and Climate; Routledge: London, UK, 2008.
84. Faillettaz, J.; Funk, M.; Vincent, C. Avalanching Glacier Instabilities: Review on Processes and Early Warning Perspectives. Rev.
Geophys. 2015, 53, 203–224. [CrossRef]