Geomorphological and hydrological controls on sediment export in earthquake-affected catchments in the Nepal Himalaya

Large earthquakes can contribute to mountain growth by building topography but also contribute to mass removal from mountain ranges through widespread mass wasting. On annual to decadal or centennial timescales, large earthquakes also have the potential to significantly alter fluvial sediment dynamics if a significant volume of the sediment generated reaches the fluvial network. In this contribution, we focus on the Melamchi–Indrawati and Bhote Koshi rivers in central Nepal, which have both experienced widespread landsliding associated with the 2015 Gorkha (Nepal) earthquake. Using a time series of high-resolution satellite imagery, we have mapped exposed sediment along the rivers from 2012–2021 to identify zones of active channel deposition and document changes over time. Counter to expectations, we show negligible increases in coarse-sediment accumulation along both river corridors since the Gorkha earthquake. However, an extremely high-concentration flow event on 15 June 2021 caused an approximately 4-fold increase in exposed sediment along a 30 km reach of the channel with up to 12 m of channel aggradation in the Melamchi–Indrawati rivers; this event was localised and did not impact the neighbouring Bhote Koshi catchment. Based on published reports, new helicopter-based photography, and satellite data, we demonstrate that this event was sourced from a localised rainfall event between 4500 and 4800 m and that a significant fraction of the sediment was supplied from sources that were unrelated to the landslides generated by the Gorkha earthquake.

1 Introduction

By building topography, large earthquakes contribute to the growth of mountain ranges (Avouac, 2007), but they also play a role in their erosion by producing sediment through coseismic landsliding, which will eventually be evacuated from the mountain range (e.g. Keefer, 1994; Dadson et al., 2004; Larsen et al., 2010; Hovius et al., 2011; Parker et al., 2011; Egholm et al., 2013). In some cases, the volume remobilised through mass wasting can equal or even outweigh that gained through surface uplift (e.g. Hovius et al., 2011; Parker et al., 2011; Marc et al., 2016, 2019; Francis et al., 2020; Li et al., 2014). Malamud et al. (2004) suggest that earthquakes above a threshold moment magnitude of M_w=4.3 trigger landslides through ground shaking, with landslides usually concentrated in areas of highest peak ground acceleration (PGA) (Khazai and Sitar, 2004) and near the top of hillslopes (Densmore and Hovius, 2000; Yin et al., 2009). In the short term, the sediment generated by coseismic landsliding either remains on the hillslopes or is delivered to fluvial systems, where it can cause aggradation and decrease channel capacity, thereby increasing flood risk downstream and causing rapid changes in channel morphology (Korup et al., 2004; Sims and Rutherfurd, 2017). For example, the 1999 Chi-Chi (Taiwan) and 2008 Wenchuan (China) earthquakes were associated with up to 18 and 10 m of channel bed aggradation, respectively (Chen and Petley, 2005; Whadcoat, 2011; Yanites et al., 2010). Coseismic landslide sediment on hillslopes often remains as loose sediment until hillslopes re-vegetate or re-stabilise, which may take years. Consequently, large earthquakes are frequently followed by a period of increased landsliding as (i) the ground shaking destabilises hillslope material, leaving loose sediment readily mobilised during subsequent rainfall events and (ii) coseismic landslides are remobilised (Dadson et al., 2004; Chen and Petley, 2005; Huang and Fan, 2013; Marc et al., 2015). In the 2 years following the Kashmir earthquake, Saba et al. (2010) documented an increase in landsliding, succeeded by a period of relative stability, while Li et al. (2022) show a growth in landslide area for 5 years after the 2008 Wenchuan earthquake. Similarly, Kincey et al. (2021) noted elevated rates of landsliding for at least 3 years in the wake of the 2015 Gorkha earthquake.

Coseismic landslides are thought to affect the fluvial network for decades to centuries (Croissant et al., 2017; Hovius et al., 2011; Wang et al., 2015; Yanites et al., 2010; Wang et al., 2017). Wang et al. (2015) estimate that the residence time of suspended sediment supplied by coseismic Wenchuan landslides ranges from a year to more than a century, which ties in with previous estimates of sub-annual to centennial timescales, as suggested by a range of studies from across the globe (Pain and Bowler, 1973; Pearce and Watson, 1986; Koi et al., 2008; Howarth et al., 2012). The residence time of coseismic bedload is more challenging to constrain due to the difficulty of monitoring bedload transport (Li et al., 2014; Croissant et al., 2017). However, based on a case study of a large landslide in the Bhote Koshi valley, Croissant et al. (2017) estimate that half the coarse fraction of a given sediment volume can be removed within 5–25 years by channel narrowing to increase sediment transport capacity.

While sediment can be generated through landsliding, sediment transfer downstream will depend on the ability of the river network to transport the newly generated sediment. In areas that experience widespread landsliding, river systems tend to become transport limited (e.g. Pratt-Sitaula et al., 2004; Yanites et al., 2010), whereby the quantity and calibre of sediment transported by rivers is limited by the magnitude of subsequent floods. Extreme flood events such as those generated by glacial lake outburst floods (GLOFs) and unusually intense rainfall have the capability to mobilise very large amounts of sediment. Channel gradients and valley morphologies exert a control on where sediment is mobilised and where it is deposited during extreme flood events (Belmont et al., 2011; Devrani et al., 2015; Gran and Czuba, 2017). While widespread landsliding has the potential to turn supply-limited systems into transport-limited systems, we hypothesise that the exported sediment flux and the loci of sediment entrainment will be dominated by both the magnitude of subsequent hydrological events and the inherited nature of the river valley in terms of morphology and sediment fill.

Our study focuses on the Melamchi–Indrawati and Bhote Koshi valleys in central Nepal that were both affected by the 2015 Gorkha earthquake and a series of geomorphically significant hydro-meteorological events, making them ideal locations to detect the relative roles of inherited valley form, remobilisation of seismically induced landslides, and the origins of the hydrological events. The earthquake had a similar moment magnitude to the 2008 Wenchuan earthquake (M_w 7.8 compared to M_w 7.9) and caused widespread landsliding (e.g. Collins and Jibson, 2015; Gnyawali and Adhikari, 2017; Roback et al., 2018). Here, we use optical satellite imagery to identify zones of sediment input and aggradation along the Melamchi–Indrawati and Bhote Koshi rivers and assess how they change over time, in particular in response to notable geomorphic events. Our analysis focuses on the period 2012–2021, which includes the 2015 Gorkha earthquake, as well as a series of major hydrological events such as the 2016 Bhote Koshi GLOF and the 2021 Melamchi flood (these are presented in detail below). The aim of this work is to test for the respective roles played by the remobilisation of hillslopes by the earthquake, inherited valley morphology, and the nature of major hydrological events on the flux of coarse bedload out of these high Himalayan catchments. In the following section, we present the study area, the Gorkha earthquake and associated landslides, and the major hydrological events that have affected the study area. We then present our methods and results, before discussing the implications of this work in terms of understanding sediment generation and transport processes in mountainous areas and their associated hazards.

2 Study area

Our study focuses on the Melamchi–Indrawati and Bhote Koshi catchments. Both are located in Nepal's Sindhupalchok district, which was one of the worst-affected areas in the aftermath of the Gorkha earthquake in terms of landslide density (e.g. Collins and Jibson, 2015; Kargel et al., 2016; Roback et al., 2018). Of the 87 km² of landslides triggered by the Gorkha earthquake and mapped by Roback et al. (2018), 12.3 km² are located in the Melamchi–Indrawati catchment and 24.5 km² in the Bhote Koshi (Roback et al., 2018).

Figure 1 (a) Catchment geology adapted from Yin (2006) and Dhital (2015). Gorkha landslides mapped by Roback et al. (2018) are shown as black polygons. The two studied catchments are delineated by the white outline: Melamchi–Indrawati to the west and Bhote Koshi to the east. (b) Melamchi–Indrawati catchment with main rivers and localities mentioned in text labelled. (c) Southern portion of the Bhote Koshi catchment with main rivers and localities mentioned in text labelled. Studied river reaches are shown in orange. The background of all three panels is a shaded relief map from a 30 m Copernicus DEM. B is Bahrabise, BT is Bremthang, CB is Chanaute Bazaar, D is Dolalghat, K is Kodari, MB is Melamchi Bazaar (location of Fig. 2), MG is Melamchigaon, N is Nakote, T is Timbu, MCT is the Main Central Thrust, and STD is the South Tibetan Detachment.

The Melamchi–Indrawati catchment is located ca. 20 km east of Kathmandu (Fig. 1), and the Indrawati River is fed by three main tributaries, the Melamchi, Yangri, and Larke rivers (Fig. 1b), originating in the Ganja La range. We focus on the Melamchi River and the Indrawati River downstream of its confluence with the Melamchi River, herein referred to as the Melamchi–Indrawati rivers. Elevations in the Indrawati catchment range from around 600 m at Dolalghat (D in Fig. 1b) to >5000 m, and the climate spans temperate to polar tundra environments (Pandey et al., 2021). Above Nakote (N in Fig. 1b), the Melamchi River is steep and narrow, confined into a “classic” V-shaped valley; it widens downstream of Nakote to include fertile agricultural land on its banks (Baskota et al., 2021) (Fig. 1b). Downstream of Melamchi Bazaar (MB in Fig. 1b), the Indrawati River widens into a braided system with a floodplain reaching up to a kilometre in width. The catchment also hosts the Melamchi Water Supply Project, currently under construction, aiming to supply water to households as far away as Kathmandu (e.g. Rest, 2019).

The Bhote Koshi, called Poiqu in China, originates in China's Nyalam County (Liu et al., 2020). At 3400 km², its catchment is almost 3 times the size of the Melamchi–Indrawati and sources roughly half of its area from the Tibetan Plateau (Fig. 1c). As in the Melamchi–Indrawati, the climate in the Bhote Koshi catchment ranges from temperate to polar (Karki et al., 2016), and the lower, temperate part of the catchment features distinct wet and dry seasons, with ∼80 % of the annual precipitation occurring during the Indian Summer Monsoon (Bookhagen and Burbank, 2010). Annual rainfall is high, ranging between 1300 and 4100 mm (Tanoli et al., 2017), and the Bhote Koshi is prone to glacial lake outburst floods (GLOFs) – Khanal et al. (2015) reported six GLOFs in the Bhote Koshi since 1935. The Araniko Highway, an important trade and transport route and one of the two routes linking Nepal with China (Liu et al., 2020; Whitworth et al., 2020), runs alongside the Bhote Koshi. In their study of landslide susceptibility along the Araniko Highway from Dolalghat to Kodari (see Fig. 1c for locations), Nepal et al. (2019) identified the area around Kodari as particularly prone to landsliding.

2.1 Geologic and geomorphic setting

Both the Melamchi–Indrawati and Bhote Koshi catchments lie within the gneisses, schists, quartzites, and Miocene granitic intrusions of the Higher Himalayan Crystallines (HHC) and the lower-grade metamorphic rocks (phyllites, quartzites, and slates) of the Lesser Himalayan Sequence (LHS), which are separated by the Main Central Thrust (MCT) (Gansser, 1964; Upreti, 1999). The Bhote Koshi catchment additionally extends into the carbonate-dominated sedimentary rocks of the Tethyan Himalayan Sequence, separated from the HHC by the South Tibetan Detachment (STD) (Yin, 2006), as shown in Fig. 1a. Our studied reach along the Melamchi–Indrawati rivers (Fig. 1b) extends through the HHC, crosses the MCT, and then flows predominantly through LHS units, including the quartzites of the Syangja Formation, as well as the Seti Formation's phyllites and quartzites. The Bhote Koshi study reach (Fig. 1c) lies mostly within the LHS, predominantly the phyllites and quartzites of the Seti Formation, and also flows through LHS shales, slates, and limestones in the MCT zone near the Nepal–China border.

2.2 The 25 April 2015 M_w 7.8 Gorkha (Nepal) earthquake and associated landslides

The April 2015 Gorkha earthquake was the result of a 140 km long rupture along the Main Himalayan Thrust (Avouac et al., 2015; Hayes et al., 2015; Elliott et al., 2016). A sequence of aftershocks followed the main event, the largest of which occurred on 12 May 2015 in the Dolakha district of Nepal (Fig. 1), reaching a moment magnitude of 7.3 (Collins and Jibson, 2015). A total of 9000 deaths and billions of US dollars in economic losses resulted from the mainshock alone (Government of Nepal National Planning Commission, 2015). While initial surveys (e.g. Kargel et al., 2015) reported a significantly lower number of landslides than expected for an earthquake of moment magnitude 7.8, later inventories identified up to 25 000 coseismic landslides covering a total area of ca. 87 km² (Martha et al., 2017; Roback et al., 2018), which is more consistent with expectations (Malamud et al., 2004; Marc et al., 2016). Most landslides observed in the months following the earthquake had been initiated in weathered and fractured surface material (Collins and Jibson, 2015). In terms of regional landslide distribution, Roback et al. (2018) report no regional lithological control – instead, landslides are concentrated where high peak ground acceleration (PGA), steep slopes, and high mean annual precipitation coincide. Regarding landslide threat to human populations, Collins and Jibson (2015) documented over 69 partial, full, or temporary landslide dams in the aftermath of the earthquake, although formation of a landslide lake was observed for only half the dams. At the time of writing, a wide range of Gorkha landslide inventories have been published, mapping coseismic landslides as points (Kargel et al., 2016; Tiwari et al., 2017); a mixture of points, polygons, and polylines (Williams et al., 2018); or polygons (Regmi et al., 2016; Gnyawali and Adhikari, 2017; Martha et al., 2017; Roback et al., 2018; Xu, 2018; Kincey et al., 2021; Valagussa et al., 2021). Of these inventories, Roback et al. (2018) appears the most complete, including 24 915 landslides covering a total area of 87 km² and separating the full landslide area from the landslide source.

2.3 Significant hydro-geomorphic events

This section briefly presents significant hydro-geomorphic events in both catchments that occurred in the time frame covered by our study (2012–2021).

Figure 2 Impact of the 2021 Melamchi flood in Melamchi Bazaar (MB in Fig. 1b). (a) Picture from a foot bridge towards Melamchi Bazaar, looking NW, taken on 28 October 2019. The colourful building in the centre of the photo is Hotel Roj, used as a reference in the other photos (position: 27.82917^∘ N, 85.57714^∘ E). Note gabions along the gravel bar to limit bank erosion during monsoon floods. (b) Picture from a road bridge looking NE, upstream the Indrawati River, taken on 28 October 2019. Hotel Roj is to the left. Note gabions on the gravel bar. (c) Google Maps image of Melamchi Bazaar before the 2021 Melamchi flood (© 2021 Maxar Technologies). Road and foot bridges are indicated, as well as road junction and Hotel Roj, for reference. (d) Aerial image of Melamchi Bazaar after the main flood event of June 2021, looking NW (source: Geovation Nepal). Both bridges have been destroyed, and only the top two floors of Hotel Roj emerge from the sediment, indicating aggradation in excess of 10 m.

In August 2014, collapse of the hillside at Jure village along the Bhote Koshi produced 6 million cubic metres of debris, killing 156 people and creating a deposit 100 m thick that also reached the opposite bank of the river (van der Geest, 2018). The landslide created a 55 m high dam blocking the river and forming a 3 km long lake which filled and overflowed in 12 h (Acharya et al., 2016). The landslide dam eventually breached on 7 September, mostly due to overtopping (Lamichhane et al., 2021). According to van der Geest (2018), the landslide was triggered by heavy rainfall during the preceding 2 d; the slope had been noted as being unstable in the years before the event (Shrestha and Nakagawa, 2016; Ao et al., 2020).

Originating from Gongbatongshacuo Lake in Tibet Autonomous Region, China, a moraine-dammed lake that broke its dam as a result of a rainfall-induced debris flow, a GLOF swept along the Bhote Koshi on 5 July 2016 (Cook et al., 2018), destroying or damaging parts of the Araniko highway, buildings in Kodari and Tatopani, and the intake dam of a hydropower project.

A total of 4 years later, on 9 July 2020, heavy monsoon rainfall caused a debris flow along the Bhote Koshi which killed two people, destroyed or damaged dozens of houses, swept away a 700 m stretch of the Araniko Highway, and damaged the Middle Bhotekoshi Hydropower Project (Lamichhane et al., 2021).

Along the Melamchi River, a flash flood occurred on 15 June 2021, lasting for 10 h and depositing up to 12 m of mostly sand and silt in the river channel (Fig. 2) (Pandey et al., 2021). A total of 5 people were killed, with 20 still missing as of September 2021, and 525 families were displaced by the event. The flood damaged 337 houses, severely damaged the headworks of the Melamchi Water Supply Project, delaying delivery of water to households in Kathmandu, and disrupted access to several villages by blocking roads and damaging or destroying more than a dozen bridges (Baskota et al., 2021; Maharjan et al., 2021; Pandey et al., 2021). Maharjan et al. (2021) report that more than 1.75 km² of agricultural land, crucial for subsistence farming, was lost in the wake of the event within the Melamchi Municipality alone. The event is thought to have been caused by a combination of glacial moraine collapse, landslide dam outburst, and heavy precipitation (Baskota et al., 2021; Maharjan et al., 2021). On 31 July 2021, a second debris flow scoured the area previously devastated by the prior event (Baskota et al., 2021).

The above events are well constrained in space and time, facilitating taking into account their impact during interpretation of our data.

Figure 3 Examples of mapped sections, with the mapped sediment area polygons shown in yellow. Panels (a) and (b) show a section around Timbu overlain on imagery obtained through Google Earth Pro (© 2022 Maxar Technologies). Panels (c)  and (d) show a section around Melamchi Bazaar overlain on Planetscope imagery (Planet Team, 2017). See Fig. 1b for location of the sections within Melamchi–Indrawati catchment. (a) Sediment area on 3 May 2015. (b) Sediment area on 1 February 2016. Note the influx of sediment from the tributary to the NE. (c) Sediment area on 2 June 2021. (d) Sediment area on 29 July 2021. Note the dramatic increase in sediment area extent following the 15 June 2021 Melamchi disaster.

3 Methods

In order to quantify changes in sediment deposition patterns, which could then be linked to events supplying or evacuating sediment and valley morphology, we map exposed sediment along the Melamchi–Indrawati and Bhote Koshi rivers over a time period spanning 10 years (2012–2021); these data are used to constrain zones of sediment input and quantify planform changes in channel morphology. We then measure channel steepness through long river profile analysis and valley width within the studied catchments. Lastly, in light of the significant impact of the 2021 Melamchi event on the morphology of the Melamchi–Indrawati rivers, we further use optical satellite data and photographic images to identify the sources of sediment for the Melamchi flash floods and assess the likely contribution of the remobilisation of landslide sediment generated during the Gorkha earthquake.

Figure 4 Steps for processing the data shown using an example section near Melamchi Bazaar (see Figs. 1b and 3c, d for the location). (a) Original mapped channel and sediment areas. (b) The polygon representing the sediment area is converted to a raster, and the line representing the channel is split into points. Note that the rasterisation process creates apparent gaps where sediment areas are narrower than the resolution of the raster res_r; this does not affect the trends in area along the channel (Figs. A1a and A2a). (c) A window of width d_o and length d_p is created around each channel point, and the raster area within each window is calculated. (d) The raster area calculated within each window along the river is smoothed using a rolling window and plotted as a function of distance along the channel.

3.1 Determining sediment area along the river channel

3.1.1 Mapping and data processing

We manually mapped exposed sediment along the Melamchi–Indrawati and Bhote Koshi rivers every year from 2012–2021 to document changes after each monsoon. We used sub-metre resolution imagery from Maxar Technologies and CNES/Airbus, available in Google Earth Pro, for the time period between 2012 and 2019 inclusive, and 3 m resolution Planetscope imagery, obtained through Planet's Education and Research Program, for 2020 onward (Fig. 3) (see Appendix A for the full breakdown of imagery dates and sources). From 2012–2019 inclusive, each river was mapped once per year so that eight post-monsoon mapping epochs were available. For 2020 and 2021, when sufficient imagery was not available in Google Earth Pro, we switched to Planetscope imagery. In 2021, we took advantage of Planetscope's higher temporal resolution to map both rivers monthly from June–December inclusive. In the Melamchi–Indrawati catchment, mapping was carried out from the confluence with the Bhote Koshi at Dolalghat up to 3.5 km upstream of Nakote (Fig. 1b). Along the Bhote Koshi, we mapped from the confluence with the Indrawati at Dolalghat to 10 km upstream of the Nepal–China border crossing at Kodari (Fig. 1c). As it was not possible to pick images from the same month or season every year, we mapped exposed sediment and water combined to avoid changes in exposed sediment area due to seasonal water level variations (Fig. 3).

Figure 5 (a) Five versions of both the Google Earth (G1-5) and Planetscope (P1-5) mapped area for Melamchi–Indrawati, 22 April 2018. (b) Histogram of the standard deviations of the 10 different measured areas for each data point in (a). (c) Five versions of both the Google Earth (G1-5) and Planetscope (P1-5) mapped area for Bhote Koshi, 7 and 9 December and 10 December 2017, respectively. Visibility in the Planetscope imagery was too poor to map the section from 0–18 km. (d) Histogram of the standard deviations of the 10 different measured areas for each data point in (c), with the 75th percentile indicated by the dashed black line. G is Google Earth, while P is Planetscope.

Total area from each mapping epoch is delineated as a single polygon shapefile. In order to calculate and visualise changes in sediment area along the river profile and thereby pinpoint the location of any observed changes, the data are processed as follows (Fig. 4).

1. Each polygon shapefile of sediment area is converted into a raster of a resolution res_r=40 m.

2. For each river, a polyline shapefile representing the river channel within the studied reach is manually delineated. This polyline is not intended to represent the actual centreline of the river as it is only used to segment the river longitudinally. For comparability, the same polyline is used for all mapping epochs. This polyline shapefile is then split into points with a spacing d_s=50 m.

3. A window of dimensions d_p×d_o is created centred on each of the points, with d_p representing the window distance along the channel and d_o the window distance orthogonal to the channel. We set d_p=500 m and d_o=1000 m in the following. The sensitivity of the results to changing the size of these windows is explored in Appendix A2.

4. The area of the raster representing the mapped combined channel area (i.e. water and visible sediment) is calculated within each window.

5. This “windowed area” is then plotted against distance along the channel, additionally smoothed using a rolling mean. In the following, we chose a window size k=7, which achieves the best balance in terms of providing reach-scale information while limiting noise (see sensitivity analysis in Appendix A2).

3.1.2 Constraining uncertainty in mapping and limitations

All the mapping was carried out by one researcher. The rigorous uncertainty calculations following Fan et al. (2018), who had five independent mappers and assessed the degree of matching between each of the five data sets and a reference mapping data set, could therefore not be applied here. Instead, we selected dates on which imagery was available from both Planetscope and Google Earth Pro and mapped the full study reach of each image 5 times, i.e. 10 times in total for each river. For the Melamchi–Indrawati rivers, 22 April 2018 was selected, and for the Bhote Koshi, we chose 7/9 and 10 December 2017 for Google Earth Pro and Planetscope, respectively. Each of the repeatedly mapped areas was then processed as outlined in Sect. 3.1.1, and the median and standard deviation of the 10 mapped area values were calculated along each point along the river (Fig. 5). The 75th percentile of the collection of standard deviations was subsequently added as an uncertainty envelope in all plots showing area along profile and also carried forwards to plots showing expansion or removal of sediment along the profile (Figs. 6 and 7). The uncertainty values are $\mathrm{6.2}\times {\mathrm{10}}^{-\mathrm{3}}$ km² for the Melamchi–Indrawati rivers and $\mathrm{5.3}\times {\mathrm{10}}^{-\mathrm{3}}$ km² for the Bhote Koshi.

Figure 6 Changes in sediment area along the Melamchi–Indrawati rivers: (a) long profile of the Melamchi–Indrawati rivers within the study area, (b) mapped sediment area per 0.5 km² window along profile, (c) difference in mapped area along profile between each time period, and (d) valley width extracted with LSDTopoTools plotted alongside mapped expansion in sediment between June and July 2021 (i.e. difference in mapped sediment area before and after the Melamchi event). Areas until November 2019 or April 2020 inclusive are mapped in Google Earth Pro; areas from June 2021 onward are mapped in QGIS from Planetscope imagery. Locations mentioned in the text are labelled in (a) and shown with vertical grey lines. The dashed grey line indicates the approximate location of the Main Central Thrust (MCT), separating the Higher Himalayan Crystallines (HHC) from the Lesser Himalayan Sequence (LHS).

Figure 7 Changes in sediment area along the Bhote Koshi: (a) long profile of Bhote Koshi within the study area, (b) mapped sediment area per 0.5 km² window along profile, (c) difference in mapped area along profile between each time period, and (d) valley width extracted with LSDTopoTools. Areas from October 2016 and October 2018 onward are mapped in QGIS from Planetscope imagery; all other areas are mapped in Google Earth Pro. Locations mentioned in the text are labelled in (a) and shown with vertical grey lines. The dashed grey line indicates the approximate location of the Main Central Thrust (MCT), separating the Higher Himalayan Crystallines (HHC) from the Lesser Himalayan Sequence (LHS).

A limitation of our approach is that it does not account for the vertical dimension; all changes are documented in plan view. For a given amount of aggradation, channel reaches with steep valley sides will exhibit less discernible widening of the visible water and sediment extent. The method will not be able to detect changes in a gorge with vertical walls. However, even in the most confined reaches, valley walls are rarely completely vertical, so we believe that aggradation, while not as clear as in wide valley sections, would still be visible. Furthermore, the amount of sediment that would fill up a wide valley by a couple of metres would cause much thicker aggradation if deposited in a narrow reach, increasing the likelihood of detection. We will use documented local aggradation locations to assess this bias. Finally, it is important to note that where local constrictions occur, aggradation does tend to occur upstream of the constriction rather than in the gorge where shear stress will be the greatest (e.g. Turzewski et al., 2019), except if the gorge itself were dammed by a landslide, which would likely lead to detectable changes at the reach scale.

3.2 Valley morphology and channel steepness

We use the topographic analysis package LSDTopoTools (Mudd et al., 2021) to extract river profiles, normalised channel steepness, and valley width within the Melamchi–Indrawati and Bhote Koshi catchments. The aim of this analysis is to further identify potential zones of sediment storage. For valley width, we use the automated method developed by Clubb et al. (2017), which first identifies floodplains based on channel relief and local gradient before calculating the width of the extracted floodplain at each point in the channel. We apply the floodplain and valley width extraction to a 30 m resolution Copernicus Digital Elevation Model (DEM). In order to limit the likelihood of overestimating channel width, which typically occurs when channels become narrower than a couple of pixels, all tributary channels below a Strahler stream order of 3 are ignored. We calculate normalised channel steepness k_sn from the equation $S=k_{\mathrm{sn}}{A}^{-{\mathit{\theta }}_{\mathrm{ref}}}$, where S is the channel gradient (dimensionless), A is the drainage area (m²), and θ_ref is the reference concavity index (dimensionless), using the statistical approach described in Mudd et al. (2014). The value of θ_ref can significantly affect the values of k_sn (Gailleton et al., 2021), and we therefore constrain θ_ref using the method outlined by Gailleton et al. (2021). This yielded a best-fit θ_ref value of 0.15 for the Melamchi–Indrawati and Bhote Koshi catchments combined. This value is relatively low compared to the global data set shown in Gailleton et al. (2021), which may reflect the importance of debris flows in controlling the long-term shape of the long profiles in these catchments (e.g. Gasparini et al., 2004; Whipple et al., 2013). Additionally, Gailleton et al. (2021) suggest that the presence of glaciers and spatial variations in tectonic processes can obscure the optimal concavity value. However, as we are using k_sn as a metric to describe relative spatial variations in channel steepness within the catchments, the low value of the reference concavity does not present a major issue.

4 Results

4.1 Determining sediment area along the river channel

The sediment area mapped along the Melamchi–Indrawati rivers and the valley width extracted from digital topographic data highlight the transition from a confined river system to a wide and alluvial channel ca. 5 km downstream of Melamchi Bazaar (Fig. 6b and d), with valley width reaching 1000 m at about 45 km. Sediment area shows temporal noise across all mapping epochs; this is most pronounced in the reach extending from 40 to 60 km (Fig. 6b and c). Following the Gorkha earthquake (i.e. from the May 2015 to January/February 2016 mapping epochs), the Melamchi–Indrawati rivers experienced localised increases in sediment area. Although the mainshock of the Gorkha earthquake series pre-dates the May 2015 imagery, the time interval between the earthquake and the imagery is short enough that we can assume the majority of coseismic landslide sediment had not yet been delivered to the river network or had been transported a significant distance in May 2015. The most prominent increase in sediment area following the earthquake is seen at Timbu between May 2015 and January/February 2016, at a distance of around 16 km downstream (Figs. 1 and 6a–c), where the Melamchi River received an influx of sediment from a debris flow in a tributary channel, as shown in Fig. 3a and b; this peak persists until early 2020. We confirmed that the deposit leading to this increase in sediment area is from a debris flow during a field visit to Timbu in October 2019. As Timbu is located in one of the narrowest steep-sided sections of the Melamchi–Indrawati valley (Fig. 6a and d), the fact we are able to detect this change gives us confidence that our mapping method works even in steep-sided valleys, despite the limitations discussed in Sect. 3.1.2.

No other significant change is detected along the Melamchi–Indrawati rivers until June to July 2021, when there is a marked, consistent increase in sediment area from just upstream of Timbu to ca. 5 km downstream of Melamchi Bazaar (Fig. 6). The maximum increase in sediment area observed within a 0.5 km² window is 0.18 km² at ca. 37 km (i.e. downstream of Melamchi Bazaar), which coincides with the furthest downstream extent of the continuous increase in sediment area between the two mapping epochs and a marked increase in valley width (Fig. 6b–d). Although the spatially continuous increase in sediment area between mapping epochs stops at the transition to a wider alluvial valley, localised increases of up to 0.08 km² are observed along the remaining study reach downstream until Dolalghat. The location and timing of this signal is consistent with field reports documenting the aftermath of the 15 June 2021 Melamchi disaster (Baskota et al., 2021; Pandey et al., 2021) (Sect. 2.3). Immediately upstream of Melamchigaon, a local increase in sediment area of ca. 0.1 km² is observed between June and July 2021, and the resulting peak persists until the last mapping epoch in December 2021. In the same time period, Baskota et al. (2021) and Pandey et al. (2021) document a reactivation of a landslide at Melamchigaon, which briefly dammed the Melamchi River, thereby contributing to the scale of the June 2021 Melamchi disaster. Between July and December 2021, we observe additional increases in sediment area between Timbu and Melamchi Bazaar, and the July 2021 peak translates several kilometres more downstream to the approximate location of the MCT at 40 km along our study reach (Fig. 6b and c).

Valley width is much less variable along the Bhote Koshi, with values not exceeding 200 m for the full length of the studied reach (Fig. 7b and d). Sediment area along the Bhote Koshi is lower and also less variable spatially than along the Melamchi–Indrawati rivers. At 37 km downstream along the Bhote Koshi, an increase in sediment area after February 2013 is visible (Fig. 7 and c). This peak persists over the subsequent mapping epochs, attenuating slightly beginning in October 2020. The timing and location of this signal is consistent with the 2014 Jure landslide, which dammed the river and led to subsequent aggradation upstream (Acharya et al., 2016; van der Geest, 2018). Following the Gorkha earthquake (i.e. from the May 2015 to January/February 2016 mapping epochs), localised increases in sediment area are observed as along the Melamchi–Indrawati rivers, but they are smaller in magnitude than the Jure landslide. No signal can be clearly linked in space or time to the hydro-geomorphic events outlined in Sect. 2.3, including the July 2016 GLOF and the July 2020 debris flow. While it is known that the 2016 GLOF caused valley widening of up to 35 m in places (Cook et al., 2018), the spatial and temporal resolution of our analysis does not allow us to unambiguously detect these changes.

In summary, no significant continuous signals in either the spatial or temporal dimension are observed along the Bhote Koshi or the Melamchi–Indrawati rivers in the aftermath of the 2015 Gorkha earthquake. The first and only such signal is observed along the Melamchi–Indrawati rivers in July 2021, following the June 2021 Melamchi disaster, and persists until the end of our mapping period. While some changes in sediment area along both rivers could be attributed to the initial impact of additional sediment yield in the 6 years following the Gorkha earthquake, these are localised and small in comparison to the impact of the 2021 Melamchi event.

4.2 Valley width and channel steepness

As shown in Figs. 6d and 7d and further highlighted in Fig. 8a and b, the Melamchi–Indrawati and Bhote Koshi rivers show clear differences in valley width along the profile. While the Melamchi–Indrawati transitions from a confined river system to a wide, alluvial channel at ca. 37 km, with valley width reaching 1000 m at 45 km, valley width remains comparatively constant along the studied reach of the Bhote Koshi and does not exceed 200 m. The spatial distribution of normalised channel steepness in the Melamchi–Indrawati catchment follows the expectations established by the pattern of total valley width: the headwaters are steepest, and the wide alluvial section of the Indrawati River downstream of Melamchi Bazaar has the lowest gradient. In the Bhote Koshi catchment, the pattern of normalised channel steepness is slightly more complex. The lowest gradients are present in the headwaters on the Tibetan Plateau as well as the lower reaches near the confluence with the Indrawati River. In the section around the Nepal–China border, the tributary channels and trunk channel show the steepest gradients. The spatial distribution of normalised channel steepness in the Balephi River (adjacent to the Melamchi–Indrawati catchment) mirrors that seen in the Melamchi–Indrawati catchment, with the highest gradients in the headwaters and smaller tributaries and the gentlest channel slopes in the lower reaches. In the headwaters of the Melamchi River, Fig. 8a and b show a wide (>300 m) section with low gradient compared to the adjacent reaches. Downstream of this section, the valley narrows to <50 m again, and normalised channel steepness abruptly increases. There is a similar configuration (i.e. a wide, flat section abruptly followed by a narrow and steep channel) in the neighbouring Yangri valley, also indicated in Fig. 8a and b. Our preliminary analysis shows no such configurations in the Bhote Koshi catchment or elsewhere in the Melamchi–Indrawati catchment.

Figure 8 (a) Total valley width in the Melamchi–Indrawati catchment. (b) Total valley width in the Bhote Koshi catchment. (c) Normalised channel steepness k_sn in the Melamchi–Indrawati and Bhote Koshi catchments, using a best-fit value of θ_ref=0.15 for the combined catchments. Boxes in the headwaters of the Melamchi and Yangri rivers in panels (a) and (c) outline a similar anomalous configuration: a wide, flat reach followed by a narrow and steep reach (see Sect. 6 and Fig. 10). To extract the channel network for this figure, we selected a drainage area threshold of 15 000 pixels (or 13.5 km², as pixels are 30×30 m). The minimum and maximum values of the colour bars in all panels are set to the 10th and 90th percentile, respectively, of the data shown. The background of all three panels is a shaded relief map derived from a 30 m Copernicus DEM.

5 Exploration of the Melamchi event

Our mapped sediment areas along the Melamchi–Indrawati and Bhote Koshi rivers show temporal noise throughout the mapping period from 2012–2021, indicating that the two river systems are not static and gravel bars shift over time. Although we note several localised increases in sediment area following the Gorkha earthquake, potentially linked to input of coseismic landslide sediment, these changes are small compared to the significant increase in sediment area along the central reaches of the Melamchi–Indrawati rivers in the months following the June 2021 Melamchi disaster or in the Bhote Koshi following the Jure landslide (Figs. 6 and 7). In the following analysis, we assess whether the Gorkha earthquake was a necessary pre-condition for the 2021 Melamchi disaster. We do this by considering the causes of the Melamchi disaster and the source of its sediment.

5.1 Reporting of Melamchi disaster

At the time of writing, three detailed reports documenting the impact of the Melamchi disaster and exploring its causes have been published (Baskota et al., 2021; Maharjan et al., 2021; Pandey et al., 2021). We summarise these reports and support them with photographs taken during a subsequent helicopter flight to the region in December 2021. The initial report from Pandey et al. (2021) describes the event as an outburst flood from a landslide dam triggered by intense localised rainfall. The Department of Hydrology and Meteorology's (DHM) Sermathang station, 5 km ESE of Timbu, did record higher than normal cumulative rainfall in the 6 d preceding the event (Maharjan et al., 2021) and Global Precipitation Measurement (GPM) data shows high precipitation around Melamchi Bazaar on 14 June (Pandey et al., 2021); however, neither of these measurements are considered extreme enough to explain the scale of the event.

Maharjan et al. (2021) and Baskota et al. (2021), with the benefit of additional information from satellite imagery and further field excursions, rule out a single trigger, instead listing multiple factors as a cause of the event. They propose a chain of events culminating in the Melamchi disaster: first, collapse of glacial moraine in the Pemdang Khola, which is a tributary of the Melamchi River, where heavy rainfall or high snowmelt released a wave of glacial lake water and glacial sediments. The resulting debris accumulated in Bremthang, a wide sandy plain downstream of the confluence between the Pemdang Khola and Melamchi Khola (Baskota et al., 2021; Maharjan et al., 2021), which itself has been interpreted as a valley fill behind an old landslide dam (Maharjan et al., 2021). The floodwaters and debris released in the Pemdang Khola incised into the Bremthang old dam, mobilising additional sediment and causing a debris flow that contributed to bank erosion and hillslope toe cutting further downstream. At Melamchigaon, an old landslide which was reactivated first by the Gorkha earthquake, and subsequently during this event, is also thought to have dammed the river; the nearby gauging station at Nakote recorded a decrease in water level followed by a rapid increase before the station was destroyed during the event and stopped recording (Baskota et al., 2021; Pandey et al., 2021). Baskota et al. (2021) point out that the Melamchigaon landslide dam could have retained 32 400–43 200 m³ of water, taking into account the decrease in water level recorded at Nakote station and assuming a discharge of 9–12 m³ s⁻¹. This amount is small compared to the 70 million cubic metres of debris deposited around Melamchi, highlighting that while the Melamchigaon landslide dam outburst likely added to the disastrous scale of the event, it cannot be the sole cause of the disaster.

While Maharjan et al. (2021) cannot link the Melamchi event to the 2015 Gorkha earthquake with certainty, they suggest that the earthquake destabilised the hillslopes in the area affected by coseismic and post-seismic landsliding, as corroborated by Kincey et al. (2021), who show increased landslide activity since the earthquake.

To summarise, the previous studies have demonstrated the compound nature of the Melamchi disaster but as yet have not tackled the question of whether these processes were conditional on sediment released from hillslopes by the Gorkha earthquake.

5.2 Sediment source derived from satellite data

Planet Explorer 3.7 m optical satellite imagery was used to explore the likely sources of sediment during the Melamchi disaster. The satellite data confirm the break-out of the Pemdang Lake and the accumulation of fresh sediment observed in the Bremthang valley. However, the data also show widespread sediment excavation in a number of neighbouring valleys in the catchment, suggesting increased discharge over a broad area rather than a point source (Fig. 9). For example, many of the channels in the upper Yangri Valley, immediately east of Melamchi, also show evidence of remobilisation and expansion of sediment accumulation, as does a smaller valley to the west. The location at which the headwaters of these channels exhibited sediment remobilisation falls between an elevation range of approximately 4500–4800 m elevation and extends along the range for at least 10 km between the valleys in an WSW–ENE orientation. Critically, the source locations for the initiation of sediment mobilisation tend to be upstream of the areas that experienced intense landsliding in response to the Gorkha earthquake. While Gorkha-related landslides tend to concentrate immediately downstream of the Bremthang, sediment mobilisation is seen to occur upstream (e.g. the Pemdang Lake outburst and nearby gulleys) but also in valleys immediately east and west that did not experience any landsliding in response to the earthquake (Fig. 9).

Figure 9 Planet Explorer 3.7 m optical satellite imagery (Planet Team, 2017) of the upper reaches of the Melamchi and Yangri valleys from October 2020 (a) and October 2021 (b). The images are overlain by the coseismic landslide data (orange) by Roback et al. (2018). The lighter grey portions along the channels in (b) represent the remobilisation of sediment during the Melamchi disaster (see Sect. 5.1). Panel (c) is an interpretation of (b) highlighting the channels where sediment has been remobilised during the event (in red), where extensive accumulation of sediment occurred in the Bremthang valley (yellow) and where gorge incision has taken place downstream of the Bremthang valley (blue shading). The pink dots approximate the position of the headwaters that record the first evidence of sediment remobilisation. We hypothesise that an extreme storm event initiated the disaster: the pink dots delineate the extent of the area affected by this hypothetical storm. The star labelled P shows location of Pemdang Lake where a glacial outburst flood has originated.

These observations support the interpretations from previous reports that suggest a single dam failure or outburst event is insufficient to explain the volumes of discharge and sediment recorded in the Melamchi valley. Instead, the evidence of sediment mobilisation in multiple catchments along the range implies an intense precipitation event that was too localised to be picked up by any of the meteorological stations. Overtopping and release of the Pemdang Lake accompanied by rapid incision of the old landslide that dammed the Bremthang valley were a part of the response to this rainfall event that added to its impact. A similar “cloudburst” event caused devastation in the Ladakh Himalaya in 2010, where both the meteorological stations and radar satellite data missed the rainfall as it was localised over a 6 km wide zone that runs between multiple small catchments along strike at an elevation of approximately 4500 m (Hobley et al., 2012). Jones et al. (2021) identified similar extreme cloudburst events resulting in widespread sediment mobilisation in southern Nepal in 1993 and 2002. Similarly, in 2018, a localised band of extreme precipitation in the Californian Sierra Nevada has been shown to initiate hundreds of landslides, increasing fluvial sediment yield to up to 16-fold relative to the local annual sediment yield (Collins et al., 2020). An additional consequence of these observations is that the Melamchi disaster appears to have been independent of the Gorkha earthquake, as there is no evidence that the main sources of sediment were obviously pre-conditioned by the Gorkha earthquake. Most sediment appears to have been mobilised upstream of the locus of coseismic landslides, with a large part of the sediment sourced from valley excavation in the Bremthang area (Fig. 9), although whether surrounding hillslopes were weakened facilitating sediment mobility during intense rainfall is not evaluated here. However, incorporation of landslide-derived sediment downstream of the Bremthang valley and the river damming and dam outburst at Melamchigaon likely amplified the impact of the flood.

6 Discussion

6.1 Landslide connectivity following the 2015 Gorkha earthquake

It is possible that large-scale channel aggradation such as that observed following the 1999 Chi-Chi (Taiwan) and 2008 Wenchuan (China) earthquakes (Chen and Petley, 2005; Whadcoat, 2011; Yanites et al., 2010) did not happen in the aftermath of the Gorkha earthquake because the landslides were not connected to active river channels to the same extent (Collins and Jibson, 2015; Cook et al., 2016). A landslide connected to the river network will deliver sediment to the channel and therefore contribute to an increase in sediment load downstream, whereas a landslide isolated on a hillslope and not connected to the fluvial domain will not immediately contribute to the fluvial sediment budget. Roback et al. (2018) found that connectivity values compiled across a broad area affected by the Gorkha earthquake were similar to those found for landslides triggered by the Wenchuan earthquake, and it is therefore unlikely that low connectivity is the reason for the lack of coseismic landslide sediment in the river channels. In the following, we explore potential causes for the lack of large-scale aggradation in our study area following the Gorkha earthquake.

Whether a landslide ends up contributing to the sediment budget will also depend on how quickly it becomes re-stabilised, mostly through re-vegetation, although the landslide mass will likely remain more easily re-mobilised than an intact hillslope. Re-vegetation of hillslopes can occur rapidly in tropical and/or wet climates. For example, Lin et al. (2007) reported 66 % of vegetation recovering in 4 years following the 1999 Chi-Chi earthquake, and Saba et al. (2010) demonstrated that hillslopes had re-vegetated and stabilised within 2 years of the 2005 Kashmir earthquake – much earlier than expected – facilitating a rapid return of landslide activity to pre-earthquake conditions. Fan et al. (2018) and Gan et al. (2019) showed that hillslopes returned to pre-earthquake conditions between 7 and 9 years after the 2008 Wenchuan earthquake, with Fan et al. (2018) citing re-vegetation, grain coarsening and densification of the mass wasting deposits as the main processes promoting rapid hillslope stabilisation. Based on these observations, Fan et al. (2018) suggest that long-term hillslope processes contribute significantly more to mass removal following large earthquakes than short-term mass wasting processes.

In the area affected by coseismic landsliding during the 2015 Gorkha earthquake, Kincey et al. (2021) compiled a multi-temporal landslide inventory over 11 mapping epochs, from pre-monsoon 2014 to post-monsoon 2018. For each 1 km² grid cell in their mapping area, landslide “birth” was assigned to the first epoch in which the cell was first impacted by a landslide, and landslide “death” indicated the epoch in which the cell was last impacted by a landslide. Although Kincey et al. (2021) did not assess hillslope stability or re-vegetation patterns, they showed that in ca. 75% of the cells registering a coseismic landslide birth, the landslides persist until the final mapping epoch, after the 2018 monsoon. This tentatively suggests that hillslopes had not yet extensively re-vegetated in the 3 years following the 2015 Gorkha earthquake. These findings therefore suggest that the lack of large-scale aggradation following the Gorkha earthquake cannot be attributed to anomalously low connectivity or rapid re-vegetation in the studied area.

While the landslide connectivity values following the Gorkha earthquake are similar to those of the Wenchuan earthquake regionally, spatial differences may explain different geomorphic responses locally. We chose the study catchments for their high density of coseismic landslides. A visual inspection of the mapped coseismic landslides (Fig. 1a) shows clusters of landslides along all of the rivers investigated in this study, with no obvious difference in density between the Melamchi–Indrawati and the Bhote Koshi. All the landslides are located in the upper part of the study reaches, suggesting that aggradation would have occurred should the landslide material have been mobilised following the earthquake. The only clear difference between the Melamchi–Indrawati and the Bhote Koshi landslides is that they are mostly located in the upper half of the Melamchi–Indrawati catchment, while they are focused south of the Nepal–China border in the Bhote Koshi catchment, which extends far upstream beyond the border. Despite clear differences in the locus of coseismic landslides within their catchments, both river systems show little evidence for aggradation following the Gorkha earthquake, suggesting that the specific location of clusters along the studied river systems is unlikely to be responsible for the lack of aggradation as a result of coseismic landsliding. However, this specific location of clusters may have played a role in the occurrence of the Melamchi disaster.

The coseismic landslides themselves have not led to significant downstream aggradation. We will therefore explore the importance of alternative controls on sediment transfer, including major hydrological events and valley morphology, in the following sections.

6.2 Importance of hydro-meteorological events for sediment export

Following the 1999 Chi-Chi earthquake in Taiwan, Typhoon Toraji (2001) resulted in a 3-fold increase in landslide area compared to pre-earthquake Typhoon Herb even though the latter brought more rainfall (Lin et al., 2004, 2006). Previous studies following large earthquakes note elevated suspended sediment loads for several years following the event itself, suggesting that further triggering events beyond the earthquake itself are needed to flush out coseismic landslide sediment (Dadson et al., 2004; Yanites et al., 2010; Hovius et al., 2011). In their study focusing on the Wenchuan earthquake, Wang et al. (2015) additionally found that the residence time of coseismic landslide sediment in catchments is reduced by the prevalence of intense rainfall.

Figure 10 Morphological anomalies and potential sediment stores in the headwaters of the Melamchi and Yangri rivers. (a, b) Google Map images of areas characterised by wide, low-gradient valleys with braided channels (L) abruptly transitioning downstream into steep narrow channels (S) (see boxes in Fig. 8) (© 2021 Maxar Technologies). The stars indicate the transition. (a) Bremthang in the Melamchi Valley, interpreted to be a valley fill behind a landslide dam (represented by fanning arrows); the breach of the landslide dam and re-incision of the valley fill is believed to have contributed a significant amount of sediment during the Melamchi disaster. (b) A similar configuration in the Yangri valley immediately east of the Melamchi valley. (c) River profiles showing the location and extent of the low-gradient reaches shown in (a) and (b) and the abrupt transition to a steep reach downstream (indicated by the stars). (d) Photo looking north (upstream) taken from the approximate original location of the knickpoint at Bremthang (star in a), showing rapid incision into valley fill and remobilisation of large amounts of sediment during the Melamchi disaster. The approximate location of the star in (a) is 28.0930^∘ N, 85.5461^∘ E. The approximate location of the star in (b) is 28.0658^∘ N, 85.5963^∘ E.

The triggers of the Melamchi event proposed by Baskota et al. (2021) and Maharjan et al. (2021) indicate that the majority of the sediment involved in the June 2021 event was sourced from an area in the catchment that was located higher than the area that had been worst affected by Gorkha earthquake-triggered landslides, suggesting that hillslope pre-conditioning by the earthquake likely played only a minor role in facilitating the Melamchi disaster. Analysis of the Planet Explorer images (Fig. 9) confirms this origin and highlight that sediment mobilisation likely occurred as a result of an extreme hydro-meteorological event focused on an area across the 4500–4800 m elevation range. The associated glacial lake outburst in the Pemdang Khola, valley damming by the Melamchigaon landslide, and subsequent dam outburst likely contributed to the severity of the event.

6.3 Importance of valley morphology for sediment export

Our study focused on the Melamchi–Indrawati and Bhote Koshi catchments in central Nepal, both of which experienced severe coseismic landsliding during the 2015 Gorkha earthquake sequence. Additionally, the studied reaches in both catchments flow predominantly across the rocks of the Higher Himalayan Crystallines (HHC) and Lesser Himalayan Sequence (LHS) and are crossed by the Main Central Thrust (MCT). However, the evolution of sediment area during our study period from 2012–2021 is dramatically different between the two catchments, with the Melamchi–Indrawati rivers experiencing a dramatic increase in sediment area associated with a catastrophic high-concentration flow event in June 2021, which persists until the end of our mapping period. In this section, we explore the potential influence of valley morphology and geomorphic history on the occurrence of this event in the Melamchi–Indrawati but not the Bhote Koshi catchment. Baskota et al. (2021), Maharjan et al. (2021), and Pandey et al. (2021) all agree that Bremthang, a wide, sandy plain in the headwaters of the Melamchi River, was a key component of the Melamchi disaster. As Baskota et al. (2021) point out, the low gradient of Bremthang allowed water and debris to accumulate. We therefore hypothesise that the configuration of Bremthang facilitated the Melamchi disaster, with a flat and wide reach upstream acting as a sediment trap connected to a steep, confined channel allowing rapid evacuation of sediment once the old landslide dam was breached, adding to the catastrophic impact of the event further downstream. As shown in Figs. 8 and 10, Bremthang is visible in maps of channel steepness and channel width as a wide and flat section abruptly followed by a transition to a steep and narrow channel downstream. We also note that the headwaters of the Yangri River, immediately to the east of the Melamchi River, show a similar configuration, albeit less pronounced than Bremthang along the Melamchi River (Fig. 10). Landslide dams such as Bremthang are made of heterogeneous material with a wide range of grain sizes, up to very large boulders. They can therefore persist for decades to centuries. By contrast, valley fill accumulating behind a landslide dam is composed of finer-grained alluvium and lake sediment (Weidinger, 2006), meaning it can be mobilised relatively easily once the landslide dam is incised. As reported by Baskota et al. (2021), the Melamchi River is currently incising into the downstream end of the Bremthang valley fill, forming a steep scarp whose collapse led to a second debris flow event on 31 July 2021 (Fig. 10d). This further highlights the potential of landslide-dammed valley fill as sediment stores that can be released during extreme hydrological events and ties in well with the Devrani et al. (2015) observation from the upper Ganga basin that low-gradient reaches aggrade and high-gradient reaches incise during extreme events. Our insights from assessing the patterns of valley width and normalised channel steepness in the Melamchi–Indrawati and Bhote Koshi catchments indicate that topographic analysis can help identify sediment stores that may pose a similar risk to the valley fill in Bremthang, which significantly contributed to the scale of the 2021 Melamchi disaster.

7 Conclusions

By producing sediment through widespread mass wasting, large earthquakes are expected to influence fluvial and sediment dynamics for years to centuries. In this study, we have mapped sediment area in the Melamchi–Indrawati and Bhote Koshi catchments in central Nepal from 2012 to late 2021 to assess the impact of the 2015 Gorkha earthquake, which triggered over 25 000 landslides and was expected to result in channel aggradation in a similar manner to that following the 2008 Wenchuan (China) earthquake, which had a similar moment magnitude. However, our data show only localised and small increases in sediment area that can be directly linked to landslides triggered by the Gorkha earthquake. The only clear reach-scale signal in either catchment is a sharp increase in sediment area from June to July 2021 along a 30 km stretch of the Melamchi–Indrawati rivers caused by an extreme high-concentration flow event on 15 June 2021. We explore the reasons for the absence of a clear signal attributable to the Gorkha earthquake and assess the influence of the Gorkha earthquake on the 2021 Melamchi disaster and the potential control of valley morphology on the occurrence of similar events. From a review of field reports following the Melamchi disaster, combined with new analysis of optical satellite data, we find that the majority of sediment transported during the event was likely sourced from a rainfall event localised in elevation but stretching between neighbouring catchments; the sediment supplied by this event was further upstream than the coseismic landslides. This indicates that a hydro-meteorological trigger played the dominant role in the Melamchi disaster rather than the hillslope pre-conditioning that caused the 2015 Gorkha earthquake. We suggest that extreme climatic events and the existence of sediment stores within valleys, rather than the instability of hillslopes, limit the export of sediment from this part of the Himalaya. We have demonstrated that maps of normalised channel steepness and valley width can help to identify such sediment stores, and we have shown that in addition to coseismic landslides, such sediment stores can pose a risk to populations in a region increasingly exposed to extreme hydro-meteorological events.

Appendix A: Mapping

A1 Imagery dates and sources

Table A1 Imagery dates and sources for the Melamchi–Indrawati rivers. Extent refers to distance along the polyline representing the channel and is consistent with the distances shown in Fig. 6. The asterisk indicates imagery used to constrain mapping uncertainty.

Table A2 Imagery dates and sources for the Bhote Koshi. Extent refers to distance along the polyline representing the channel and is consistent with the distances shown in Fig. 7. The asterisk indicates imagery used to constrain mapping uncertainty.

A2 Parameter selection and sensitivity analysis

We tested the influence of each of the five parameters involved in processing the mapping data on the resulting area values. The five parameters and their selected ranges are (i) the resolution of the raster res_r, varying from 10–70 m at an interval of 10 m; (ii) the distance between the points along the channel line d_s, varying from 10–100 m at an interval of 10 m; (iii) the window distance along the channel d_p, varying from 200 to 1100 m at an interval of 100 m; (iv) the window distance orthogonal to the channel d_o, varying from 200 to 2000 m at an interval of 200 m; and (v) the size k of the rolling window to smooth the area for plotting, for which we used values of 2, 5, and 7 and then ranging from 10–70 at an interval of 10.

We varied each parameter in turn by the ranges given above, while keeping all other parameters constant. We then processed all areas mapped from Planetscope imagery along the Melamchi–Indrawati rivers, and for each data point along the river we plotted the corresponding area calculated for the lowest chosen parameter value against that calculated using the highest parameter value (Fig. A1). Assuming that a given parameter does not influence the results, the data plotted should lie along the 1:1 line. Additionally, we show the variation in sediment area along the Melamchi–Indrawati river in one mapping epoch, October 2020, for the full range of each parameter (Fig. A2).

Figure A1 Sensitivity analysis results for (a) raster resolution res_r, (b) distance between the points along the channel line d_s, (c) window distance along the channel d_p, (d) window distance orthogonal to the channel d_o, and (e) rolling window size k. For each of the five parameters we varied in the sensitivity analysis, the area calculated at each channel point using the highest value in the parameter range is plotted against the area calculated using the lowest value in the parameter range. The choice of parameter value is considered to not influence the results if the data points lie along the dashed grey 1:1 line diagonally crossing the plot, as is the case with raster resolution in panel (a). Areas in these plots are mapped from Planetscope imagery along the Melamchi–Indrawati rivers. See the text for a detailed discussion of the results shown here, including the special cases of d_o and d_p.

Figure A2 Area along the profile for one mapping epoch showing the variation between parameter values. Areas are mapped along the Melamchi–Indrawati rivers from Planetscope imagery of 30 October 2020.

Figure A3 Total sediment area calculated along the Melamchi–Indrawati rivers as a function of d_o. The roll-off on the graph helps to identify the value of d_o above which we include the full width of the sediment areas. The more that value is exceeded, the more likely it is to get overlap between windows at bends in the river. Areas in these plots are mapped from Planetscope imagery of 30 October 2020 along the Melamchi–Indrawati rivers.

Raster resolution res_r has no discernible influence on the results, with the data points in Fig. A1 lying along the 1:1 line, and we therefore chose 40 m as the final value as this allows for low computational cost while also faithfully reproducing the shape of the mapped areas. The step distance for resampling the channel line to points, d_s, follows a roughly linear trend, albeit with some scatter. It is interesting that, as opposed to raster resolution, the epochs do not all follow the same trend. However, Fig. A2b shows that changes in the sediment area along the river channel as a function of d_s are negligible, and we could freely choose the value to be used in the processing. We decided on 50 m as it breaks the channel into enough points to accurately sample the river course while also not taking up too much computing time. Window distance along the river, d_p, is not expected to follow a 1:1 trend since the longer the window, the longer the stretch of river that is being sampled in any given window, and therefore the larger the area within this window. However, the data points here closely follow a linear trend that remains consistent across mapping epochs, and Fig. A2c shows that the overall trend does not change. Window distance orthogonal to the river, d_o, is a special case, as evidenced by the cut-off seen on the plot, indicating that for d_o=200, the area calculated within each window never exceeds 0.15 km². This is because the widest parts of the mapped channel will exceed 200 m in width, and the window will therefore not be able to capture the full extent of the combined channel area (Fig. A2d). To choose the best value for d_o, we therefore needed to make sure that it is at least the width of the widest part of the channel in any time period in order to not lose information. However, setting d_o too high causes more overlap between sections of the river that are not adjacent, e.g. at bends in the river course, and would therefore likely overestimate the area within a given window. Figure A3 shows the change in total calculated area along the river as a function of d_o. The “roll-off” apparent in this figure indicates the value of d_o beyond which no information is lost due to missing the edges of the active channel, which correlates with the manually derived value of the widest channel section, 800 m. In order to have an extra safety margin, we chose a value of d_o=1000 m. The size of the rolling window for smoothing the data, k, also shows a complex impact on the results, as evidenced by the different trends in each mapping epoch (Fig. A1e). For this case, we again refer to Fig. A2e, which shows that while high values of k remove local signals and translate peaks in sediment area downstream, this does not happen for k≤10. The final value chosen (k=7) was based on visual estimation of the smallest value which would sufficiently reduce noise in the data while not removing local signals.

A3 Handling incomplete imagery

For some dates, imagery either did not cover the entire study area or visibility was not good enough to reliably map the rivers (see Tables A1 and A2). The latter was particularly a problem for Planetscope images in the upper reaches of the study area, and especially in the Bhote Koshi, since the NE–SW orientation and steep flanks of the valley create shadows. Where imagery was not available or not sufficiently clear, these areas were not included in the data set. Since the figures in the Sect. 4 show averaged mapped areas along the profile, it was necessary to ensure that the shorter sections of mapped area did not produce erroneous spikes. To that end, we removed the mapped area within a single unit of d_p (i.e. 500 m) at either end of the mapped reach prior to smoothing the data with a rolling window.