Rapid fluvial remobilization of sediments deposited by the 2021 Chamoli disaster, Indian Himalaya Open Access

Mass movement–driven sediment cascades are hazardous to life, property, and infrastructure (e.g., Dowling and Santi, 2014) and are a key mechanism of mountain landscape evolution (e.g., Stock and Dietrich, 2003). The nature of individual mass-flow elements is well studied (e.g., Mergili et al., 2017; Geertsema et al., 2022), but the manner in which flow deposits are remobilized and the sediment flux response remain poorly understood. This knowledge gap needs addressing due to the relationship between mass-flow activity and climate change (e.g., Chiarle et al., 2021; Zheng et al., 2021), the uncertain role of mass flows in the geological record (Bovis and Jakob, 1999; Toby et al., 2022), and the impacts that elevated post-event sediment loads can have for hydropower projects (HPPs; Li et al., 2021) and ecosystems (e.g., Danehy et al., 2011).

On 7 February 2021, ~27 Mm³ of glacier ice and bedrock detached from Ronti Peak (6063 m), Chamoli district, Uttarakhand, India (Fig. 1), and transitioned into a debris flow then to a debris-laden flood along the Rishiganga-Dhauliganga-Alaknanda rivers (Shugar et al., 2021; Tiwari et al. 2022). The flow caused >200 fatalities and severely damaged key roads, bridges, and HPP installations―one under-construction (520 MW) and one operational (13.2 MW) HPP, with estimated damages to HPP infrastructure totaling ~$200 million USD (Dutta, 2021). The flow caused extensive geomorphic work (Fig. 2) including widespread channel infill as much as 100 m thick (Shugar et al., 2021), scour of hillslopes 200 m above the valley floor, temporary damming of tributary valleys, and exacerbated slope instability at Raini village, prompting rehabilitation (UDRP-AF, 2021). Rapid field reconnaissance identified dewatering debris drapes and proto-channels incising into new deposits (Fig. 2; Fig. S2 in the Supplemental Material¹). Qualitative analysis of sediment deposited by the flow reveals clasts ≥64 mm supported by a fine-grained matrix that may constitute as much as 50% by volume (Fig. S3), in contrast to pre-event channel-floor sediments, which likely comprised coarse alluvium without a supporting matrix (Fig. S4). We quantify the short-term geomorphological response of the 7 February 2021 mass flow by conducting a fluvial morphodynamic assessment for ~1 yr post-event (10 February 2021–27 January 2022). By reconstructing post-event geomorphological adjustment and sediment yields (Figs. 3 and 4), we deliver new quantitative insights into fluvial recovery following an intense landscape-scale perturbation.

A comprehensive description of our methods and their primary limitations is found in the Supplemental Material. Volumetric change on the valley floor was established from a 2 m satellite digital elevation model (DEM) of difference (DoD; 10 February 2021–27 January 2022, 2σ error = ±1.1 m). We generated high-temporal-resolution (daily) detail of fluvial erosion and deposition using the hydrodynamic-landscape evolution model CAESAR-Lisflood (Coulthard et al., 2013; Fig. 3), permitting novel insights into fast-acting erosive processes and their relationship to river discharge; CAESAR-Lisflood can simulate a wide range of flows, erosion, and deposition for as many as nine grain sizes, suspended and bedload sediment transport, and failure of unstable slopes. We specified boundary conditions and parameters including DEMs representing bedrock and surface topography, an input discharge time series for all rivers entering the model, and grain-size information. Due to a lack of local flow measurements, we drove CAESAR-Lisflood using reanalysis-derived streamflow predictions (https://geoglows.ecmwf.int). Flow was routed across the 10 February 2021 DEM, resampled to a 10 m grid to improve computational efficiency while preserving the ability to simulate channel adjustment at a scale appropriate to the aim of the study. We used a satellite DoD for 10 February 2021–25 December 2021 (2σ error = ±1.1 m) to calibrate the model, and for analysis, we subdivided the study area into a series of geomorphological “segments.” In the absence of field observations of sediment export, the timing and relative magnitude of modeled suspended-sediment export was validated with in-channel turbidity estimates retrieved from short-temporal-repeat satellite imagery for locations 25 km (Batula) and 85 km (Srinagar) downstream of the study area (Fig. 1).

Overall, 67.2% of sediment deposited along the primary flow pathway was removed in the first 12 months; 7.0 ± 1.5 Mm³ was removed from a 10.4 ± 1.6 Mm³ deposit. Of the removed volume, 66% was sourced from a thick (mean 36.4 m) deposit located immediately downstream of the Ronti Gad–Rishiganga confluence (Fig. 2A). The median DoD-derived erosion rate for this period was 2.3 ± 1.1 m a^–1. The Dhauliganga River, between Raini and Tapovan (Figs. 2B and 2C), experienced a median erosion rate of 2.2 ± 1.1 m a^–1. High erosion rates (mean 2.8 m a^–1) on the valley floor in the vicinity of the damaged Tapovan Vishnugad HPP were driven by fluvial action, while other elevation change was artificial (berm reconstruction). Downstream of Tapovan, the mean incision rate was 3.2 m a^–1 and as much as 10 m a^–1 locally at Joshimath (Fig. 1).

Morphodynamic modeling with CAESAR-Lisflood shows 6.4 Mm³ net volume loss, consistent with the satellite DoD analysis (–8.7% difference). However, at the segment scale, there were differences between the model and DoD results (Fig. 4B). From CAESAR-Lisflood modeling, we quantified four sediment discharge (Q_s) phases: (1) immediately post-event (early February–late March), (2) pre-monsoon (late March–mid-June), (3) summer monsoon (mid-June–mid-September), and (4) post-monsoon (mid-September onward) (Fig. 3A).

Immediately post-event, the mean erosion rate was 2.4 m a^–1, and 4% of the new sediment was evacuated. At this time, the system remained inefficient, with a sediment delivery ratio (SDR; a measure of sediment transport efficiency) of ~0.3 as channels re-established (Fig. 3B). The pre-monsoon phase was characterized by short-lived, high-magnitude sediment export, predominantly of fines (Fig. 3A); sediment export accounted for 62% (4.3 × 10⁶ m³) of annual Q_s, and the mean erosion rate was 5.3 m a^–1. This phase also saw a switch to an efficient, albeit short-lived, hydro-geomorphic system (SDR >1). Peak Q_s (13.2 × 10⁴ m³ day^–1) occurred in mid-May, preceding peak water discharge (Q_w) by one month.

During the summer monsoon, exhaustion of fine-grained sediment, perhaps compounded by bed armoring and the development of perched deposits due to rapid incision (e.g., Fig. 2C, Fig. S4B), markedly reduced Q_s while Q_w remained high (Fig. 3A). The falling limb of the Q_s curve occurred primarily in this period, during which total Q_s was 1.8 × 10⁶ m³ (mean erosion rate 2.0 m a^–1) or 26% of the cumulative annual load. The remaining 4% of annual Q_s was exported in the post-monsoon phase (mean erosion rate 0.7 m a^–1), when SDR reduced to <1. During 17–19 October, Uttarakhand experienced unseasonably heavy rainfall (Rawat et al., 2022). The consequent rapid rise in sediment transport was resolved by CAESAR-Lisflood (Fig. 3A), triggering a short-lived increase in sediment transport efficiency (Fig. 3B; SDR = 1.7) that is evident in our independent turbidity analyses (Figs. 3C and 3D).

The CAESAR-Lisflood morphodynamic modeling enables us to reconstruct complex and transient spatio-temporal patterns of post-event erosion and deposition (Fig. 4C). Some valley segments have long periods when they are primarily erosional or depositional, and we clearly observed the passage of at least three sediment waves with a velocity of 0.1–0.3 km day^–1. The waves initiated in segments 1, 10, and 20 and were most active from mid-March (week 6; pre-monsoon) to late August (week 25; monsoon).

The initial event produced a sudden and extreme increase in turbidity (Figs. 3C and 3D). We detected a maximum of 12,500 FNU (formazin nephelometric units) at Srinagar, the site of a 330 MW HPP, which exceeded pre-monsoon norms (mean 24 FNU equivalent; Sharma and Sharma, 2016) by three orders of magnitude. Turbidity reduced to an order of magnitude above pre-disturbance background values for the remainder of the pre-monsoon phase and sustained this difference through the monsoon. In April (at Batula) and June (at Srinagar), turbidity approached or exceeded that observed during the initial flood. The velocity of the post-event suspended-sediment pulse from Joshimath to Srinagar was ~5.7 km day^–1. Modeled suspended sediment yields at Joshimath were three and four orders of magnitude higher than observed pre-disturbance values for the monsoon (2014) and non-monsoon (2012, 2014) periods, respectively (Panwar et al., 2017).

Our results indicate a fluvial system with a short relaxation time (time taken to recover toward a pre-disturbance state) and a relatively high degree of resilience (ability to recover toward a pre-disturbance state) (e.g., Pitlick, 1993; Phillips, 2009; Rathburn et al., 2018). There was an exponential decay in sediment yield in the first 12 months, indicating a high system efficiency in evacuating initially abundant sediment along with rapid incision isolating sediment deposits (cf. Adams, 1980; Major et al., 2000; Nelson and Dubé, 2016) (Fig. S4). Our observations imply a sediment exhaustion half-life of ~1 yr in the short-term, which is comparable with those of catchments choked with sediment from the 1980 Mount St. Helens eruption (Washington, USA; Major et al., 2000) or the 1991 eruption of Mount Pinatubo (Philippines; Gran and Montgomery, 2005). The half-life tentatively extends to ~2–10 yr with a longer-term outlook and is sensitive to future hydrological regime (Fig. S5).

Importantly, our method of using numerical modeling to fill gaps between satellite DEM surveys reveals spatio-temporal nuances of sediment transport, including the role and speed of large-scale sediment waves. As channels affected by the Chamoli ice-rock avalanche–debris flow continue to relax, ongoing geomorphological work will be progressively restricted to being achieved by episodic flows that exceed seasonal norms (Wolman and Miller, 1960; Major et al., 2021), such as during mid-October 2021 (Fig. 3; Fig. S6).

Understanding rates and processes by which mountain fluvial systems respond to large-scale perturbations is essential for effective resource management and hazard mitigation. Extreme increases in fluvial sediment yields are detrimental for HPPs because sediment pulses can overload desilting chambers and increase the infill rate of storage reservoirs, ultimately decreasing HPP efficiency and lifetime (Li et al., 2022). Elevated sediment yields over the longer term may also need to be considered by operators of new HPPs (e.g., the 444 MW Vishnugad Pipalkoti HPP, 8.5 km downstream of Joshimath). Our modeling implies that ~25% of the sediment yield in the year post-event was suspended load, comprising grain fractions that are effective at damaging turbines (Singh et al., 2013). The bedload component of the modeled annual sediment yield (~75%) is high relative to the upper bound of bedload estimates from relatively undisturbed mountain rivers (e.g., Turowski et al., 2010) but is characteristic of a post-disturbance system with initially abundant and accessible sediment and with ample capacity for transport by seasonal high flows.

This field- and modeling-based study has characterized a high-mountain fluvial cascade that appears to be geomorphologically resilient to an extreme perturbation and that has a quantifiably short relaxation time. This system response must be understood for present and future HPPs, which are increasingly sited within high-altitude river reaches (Li et al., 2022) and are exposed to geomorphological impacts of extreme mass flows. The spatio-temporal sequencing of erosion and deposition in the form of potential “memory effects” (de Haas et al., 2020), whereby erosion is enhanced in zones of deposition in the post-flood system (and vice versa), may need to be considered in the context of future flood risk, while impacts on aquatic ecosystems caused by the reorganization of channel substrate and modifications to nutrient load could manifest for years to decades (e.g., Danehy et al., 2011).

Data involved in this research are in the process of being archived by the UKRI Environmental Information Data Centre and in the interim are available from the corresponding author on request. Westoby, Coulthard, Dunning, and Carrivick acknowledge support via the UK Natural Environment Research Council grant NE/W002930/1. Westoby and Berthier acknowledge DINAMIS/Centre National d’Études Spatiales (France) for facilitating the acquisition of Pléiades imagery. Shugar acknowledges Natural Sciences and Engineering Research Council of Canada (DG2020-0407) and Alberta Innovates. Shean and Haritashya respectively acknowledge support via NASA High-Mountain Asia grants 80NSSC20K1595 and 80NSSC19K0653, and Haritashya also acknowledges NASA Interdisciplinary Research in Earth Science grant 80NSSC18K0432 and the University of Dayton’s Mann Endowed Chair in the Sciences. Sain acknowledges support from the J.C. Bose National Fellowship program (SERB-DST, New Delhi). We thank Planet Labs (San Francisco, USA) for provision of satellite data, and Shashank Bhushan for assistance with DEM generation. The authors thank the three anonymous reviewers whose comments helped to improve the manuscript.