Himalayan hazard cascades – modern and medieval outburst floods in Pokhara, Nepal

In May 2012, a sediment‐laden flood along the Seti Khola (= river) caused 72 fatalities and widespread devastation for > 40 km in Pokhara, Nepal's second largest city. The flood was the terminal phase of a hazard cascade that likely began with a major rock‐slope collapse in the Annapurna Massif upstream, followed by intermittent ponding of meltwater and subsequent outburst flooding. Similar hazard cascades have been reported in other mountain belts, but peak discharges for these events have rarely been quantified. We use two hydrodynamic models to simulate the extent and geomorphic impacts of the 2012 flood and attempt to reconstruct the likely water discharge linked to even larger medieval sediment pulses. The latter are reported to have deposited several cubic kilometres of sediment in the Pokhara Valley. The process behind these sediment pulses is debated. We traced evidence of aggradation along the Seti Khola during field surveys and from RapidEye satellite images. We use two steady‐state flood models, HEC‐RAS and ANUGA, and high‐resolution topographic data, to constrain the initial flood discharge with the lowest mismatch between observed and predicted flood extents. We explore the physically plausible range of simplified flood scenarios, from meteorological (1000 m3 s−1) to cataclysmic outburst floods (600,000 m3 s−1). We find that the 2012 flood most likely had a peak discharge of 3700 m3 s−1 in the upper Seti Khola and attenuated to 500 m3 s−1 when arriving in Pokhara city. Simulations of larger outburst floods produce extensive backwater effects in tributary valleys that match with the locations of upstream‐dipping medieval‐age slackwater sediments in several tributaries of the Seti Khola. Our findings are consistent with the notion that the medieval sediment pulses were linked to outburst floods with peak discharges of >50,000 m3 s−1, though discharge may have been an order of magnitude higher.

of kilometres beyond the source. Several recent destructive landslides have received extensive attention. For example, the 2017 Xinmo rock avalanche in China has been the focus of nearly 80 scientific papers (e.g., Fan et al., 2017;Huang et al., 2019). Other examples include the 2000 Yigong rock avalanche, south-eastern Tibetan Plateau (Shang et al., 2003); the 1987 Parraguirre rock avalanche, Chile (Hauser, 2002); the 1970 Nevados Huascaran rock avalanche, Peru (Evans et al., 2009); and the 2021 Chamoli rock-ice avalanche, India . Sedimentary evidence is often used to reconstruct the magnitude and behaviour of past events (Ely & Baker, 1985;Toonen et al., 2020;Wilhelm et al., 2018). However, mass flow deposits are reworked swiftly, and limited exposure make it challenging to reconstruct past flows. Here we utilize fresh sedimentary deposits and hydrodynamic modelling to place constraints on the terminal flow phase of a recent hazard cascade in and below the Annapurna Massif in Nepal. We then reconstruct the discharge associated with older and potentially flood-derived deposits from the medieval period by linking hydrodynamic model results with preserved sedimentary evidence.
On 5 May 2012 a hyperconcentrated flow impacted the Seti Khola (= river) in the Pokhara Basin, Nepal's second most densely populated area. The flood claimed at least 72 lives and incurred a loss of some 50 million Nepalese Rupee (370,000 Euro) (Gurung et al., 2015; see Gurung et al., 2021 for a list of impacts). International media initially attributed the flood to a glacial lake outburst, whereas subsequent field and remote sensing investigations identified a series of rock-slope failures in the Annapurna Massif > 20 km upstream as the cause (Kargel et al., 2013). Field reconnaissance (e.g., Gurung et al., 2015;Kargel et al., 2013;Oi et al., 2012;SANDRP, 2014) and scientific studies (e.g., Dwivedi & Neupane, 2013;Hanisch et al., 2013) were conducted in the immediate aftermath to infer the mechanism(s) that triggered the sediment-laden flood. Debris generated by the landslides proposedly blocked the narrow and steep headwater gorge of the Seti Khola, impounding meltwater from the partly glacier-covered Sabche Cirque (Gurung et al., 2015;Kargel et al., 2013). A subsequent rock-and ice-avalanche of 33 million m 3 (Oi et al., 2014) fell from the south-western ridge of Annapurna IV (7525 m above sea level [a.s.l.]), and caused the rockslide dam to burst, although the sequence of events remains debated (Hanisch et al., 2013). The sediment-laden flood released by the proposed dam break rushed through a steep bedrock gorge in the headwaters of the Seti Khola, and destroyed Kharapani village (1100 m a.s.l.) some 18 km downstream. The flood had more than 20 reported pulses involving some estimated 7.5 million m 3 of water and sediment and a peak discharge Q p of 8400 m 3 s À1 (Gurung et al., 2015;Kargel et al., 2013). There are no stream gauge data to infer discharge, hence previous studies estimated Q p from video footage and empirical flow equations at a limited number of river cross-sections (Table 1). Sedimentation associated with the flood has also been evaluated at a limited number of locations; based on data from five cross-sections, Gurung et al. (2021) suggest there was net deposition during the May 2012 flood. We expand on these local assessments by combining hydrodynamic modelling and sedimentary evidence gathered along the length of the Seti Khola to quantify discharge and compare our values against previous estimates.
Our analysis utilizes fresh flood deposits to infer flood stage, and hence discharge, an approach that has been used to reconstruct preinstrumental and prehistoric floods (Wohl, 1995). The spatial extent of the 2012 hazard cascade is much smaller than that of at least three medieval (12th to 14th century CE) catastrophic sediment pulses that originated in the Annapurna Massif and are recorded in the valley fill of the Pokhara Basin (Fort, 2010;Schwanghart et al., 2016). These sediment pulses were likely derived from a combination of earthquake-triggered rock-ice avalanches, dam-break floods, debris flows, and intermittent fluvial reworking (Schwanghart et al., 2016;Stolle et al., 2017). However, the exact nature of the underlying depositional processes remains unclear. Although sedimentary evidence suggests an important role of fluvial processes, flow dynamics have not been evaluated quantitatively. Although a number of large prehistoric outburst events or hazard cascades are known from the Himalayas (Coxon et al., 1996;Srivastava et al., 2017;Turzewski et al., 2019), their flow characteristics are challenging to reconstruct.
The Late Pleistocene (26.9 to 43.4 ka BP) outburst of former glacial lake Batal in the Lahul Himalaya, India, impacted the Chandra Valley with an estimated Q p of 27,000 m 3 s À1 and a total flood volume of 1.5 km 3 (Coxon et al., 1996;Richardson & Reynolds, 2000). Based on observations of large, exotic boulders in the Trishuli and Sunkoshi Rivers of the Central Himalaya, Huber et al. (2020) suggested that mid-Holocence outbursts along these rivers had peak discharges between 10 3 and 10 5 m 3 s À1 . Even larger Quaternary glacier-and landslidedammed outburst floods (Q p > 10 5 m 3 s À1 ) occurred along the Yigong River, China (Korup & Montgomery, 2008;Turzewski et al., 2019).
Here we use numerical hydrodynamic models to reconstruct flow dynamics along the Seti Khola during the 2012 flood and the proposed medieval outburst floods. Hence, we focus solely on the downstream flow phases of these hazard cascades. We pursue two objectives. First, we use a HEC-RAS model calibrated using mapped sedimentary evidence of the 2012 flood to assess whether and how well our simulations match estimates of peak discharge largely derived from eyewitness accounts and amateur video footage. Second, we T A B L E 1 Reported flood estimates of the May 2012 flood along the Seti Khola (for a comparison with this study's results also see Figure 10 Figure 1).
The most recent sediment fill in the Pokhara Basin is a large alluvial fan with an estimated volume of > 1 km 3 (Blöthe & Korup, 2013) and at least three depositional units: the Tallakot, Ghachok, and Pokhara Formations (Fort, 2010;Schwanghart et al., 2016). The conglomeratic deposits of the stratigraphically oldest Tallakot and the overlying Ghachok Formations are inferred to have been deposited by Pleistocene to post-LGM (last glacial maximum) mass flows from the Sabche Cirque (Fort, 2010). The stratigraphically youngest Pokhara Formation is 60-100 m thick (Fort, 2010), with decimetre-to metrethick cobble to boulder beds of HHC provenance (Fort, 2010;Schwanghart et al., 2016). Radiocarbon dates of the Pokhara Formation are medieval and coincide with the timing of at least three documented M > 8 earthquakes (1100( , 1255( , 1344, such that the deposits may have arisen from seismically triggered, long-runout mass flows (Schwanghart et al., 2016). The processes responsible for these sediment pulses along the Seti Khola remain unclear and both lake outburst floods and long-runout ice-rock avalanches have been discussed as likely candidates (Schwanghart et al., 2016;Stolle et al., 2017). Due to this ambiguous nature (potential terminology: mass flows, debris flows, or [outburst] floods, etc.), we use the more neutral term '(medieval) sediment pulses', following the terminology of Schwanghart et al. (2016). However, as we test the plausibility of whether these sediment pulses were outburst floods representing the terminal phase of a larger hazard cascade with our models, they will also be referred to as '(medieval) outburst floods' in the context of our simulations.
Several gorges less than 1 km long but up to 90 m deep occur in the indurated, calcareous Ghachok Formation and the LHS bedrock (Fort, 2010;Stolle et al., 2019). Several lakes (Phewa, Rupa, and Begnas) formed when the medieval and older mass-flow deposits of the Seti Khola dammed several tributary mouths (Fort, 2010;Stolle et al., 2019).
The climate is seasonal with heavy summer monsoon rainfall on the southern flank of the Annapurna Massif from May to October, when > 80% of the mean annual precipitation of 4000 mm occurs (Ross & Gilbert, 1999). The central Pokhara Basin has a humid subtropical to humid temperate climate with mean monthly temperatures of 12.8 to 25.8 C (Ross & Gilbert, 1999), whereas the Annapurna Massif has temperate to alpine climate.
Pokhara is Nepal's second largest city with an estimated population of 523,000 in 2020 that tripled since the 1990s (Rimal et al., 2015). In past decades, the city and the surrounding basin have seen rapid growth in population and tourism. Urban areas increased  (Rimal, 2012;Rimal et al., 2015).
Urbanization also led to an increase in construction of informal settlements along the active channel of the Seti Khola (Rimal et al., 2015(Rimal et al., , 2018  We used our field data to manually correct the DEM-derived crosssections along the Seti Khola's deep and narrow gorges as the resolution of the AW3D DEM was not able to accurately capture the topography there. We estimated Manning's n, that is a hydraulic loss coefficient, for the channel at 61 field sites following the method of Arcement & Schneider (1984) and Chow (1959). For the floodplains, we linked Manning's n to land-cover types that we manually mapped from a 1-m resolution Maxar satellite image from 2020 that is available as ESRI army.mil/software/hec-ras/), which applies a step method for simulating 1D steady, that is gradually varied but constant, channel flow F I G U R E 2 Manning's n values estimated in the field (or from field-photographs) with the Arcement and Schneider (1984) method (grey circles) and as remotely-sensed land cover class mapping (coloured polygons) from ESRI basemap Maxar satellite imagery acquired in 2020 (ESRI and Maxar Technologies, 2022). [Color figure can be viewed at wileyonlinelibrary.com] (Brunner, 2020a;Klimeš et al., 2014). We followed the recommendations by Brunner (2020b) and choose HEC-RAS 1D modelling over the computationally more expensive two-dimensional (2D) modelling as the May 2012 flood along the steep Seti Khola was highly gravity driven and flow generally followed the river's path. In 1D form, HEC-

RAS computes water surface profiles by iteratively solving
Equation (1), in which flow energy losses are due to friction, contraction, and expansion in the natural channel geometry (Brunner, 2020a): where Z 1 and Z 2 are the elevations of the main channel bed at crosssections 1 and 2, Y 1 and Y 2 are the corresponding flow depths, V 1 and V 2 are the mean flow velocities, and a 1 and a 2 are weighting coefficients; g is the gravitational acceleration and h e is the energy head loss (Equation 2): where L is the distance weighted reach length, S f is the energy gradient, and C is the expansion (or contraction) loss coefficient. The discharge Q for each cross-section is calculated using Manning's equation (Equation 3): where A is the cross-sectional area of the flow and R is the hydraulic radius.
HEC-RAS allows for modelling steady flow in subcritical, supercritical, or mixed flow regimes (Brunner, 2020a). For modelling supercritical flow in HEC-RAS, the necessary critical depth for each crosssection is iteratively solved via the total energy head H (Equation 4): where WS is the water surface elevation and aV 2 2g is the velocity head. We computed water-surface-profiles in HEC-RAS for a total of We first simulated previous estimates of peak discharge (Table 1) with HEC-RAS. We simulate a Q p of 8400 m 3 s À1 in the in the upper Seti Khola reach, which is based on estimates for Kharapani from Gurung et al. (2015) and Oi et al. (2014). Further downstream, we used the detailed model at Seti dam to test the flood reconstructions from SANDRP (2014), including the measured flood stage of 2.5 m and Q p of 935 m 3 s À1 .
We then iteratively simulated Q p values between 1000 and 10,000 m 3 s À1 . We approximated local Q p during the May 2012 flood by determining which of the 29 modelled discharges produced flood inundation and stage that best matched the mapped flood inundation and field evidence of high water.
For all simulations, we distinguished between areas where predictions underestimated or overestimated the observed extent of sediment deposition. We assumed a mixed flow regime with critical flow depth as an upper boundary condition and normal depth with an approximated energy slope of 0.0065 as a lower boundary condition, given the high channel gradient and frequent alterations in channel geometry (Klimeš et al., 2014). In the detailed model at Seti dam and Kaseri, however, we replaced critical flow depth by the recorded water level as an upper boundary condition, while maintaining all other model specifications. Discharge data are unavailable, thus we assumed an approximated steady base flow of 100 m 3 s À1 (i.e., < 10% of our Q p scenarios) in the tributaries for all scenarios.

| Numerical flood routing with ANUGA
For simulating the larger medieval outburst floods we used ANUGA, which has been used to simulate outburst floods from breaches of man-made (Mungkasi et al., 2013) and natural dams (David et al., 2022;Larsen & Lamb, 2016;Lehnigk & Larsen, 2022) and meteorological floods (Chen et al., 2021). We chose ANUGA over HEC-RAS for computational efficiency when simulating higher magnitude flows. The purpose of these simulations is to determine the smallest discharge that inundates slackwater deposits, rather than to simulate a flood hydrograph. ANUGA is a 2D, finite-volume hydrodynamic modelling software that solves the depth-averaged shallow water equations (Roberts et al., 2015). For our simulations in ANUGA, we

| Discharge constraints from modelling and field evidence
Several aspects of our approach and the inherent assumptions might put constraints on the flood reconstructions. For example, the results of both models are sensitive to the choice of Manning's n (Westoby et al., 2014;Wohl, 1998). To assess the impact of Manning's n on our HEC-RAS simulations, we compared inundation areas and depths between models with both spatially varying and uniform Manning's n = 0.1 for a mixed flow for the scenario with the highest Q p = 10,000 m 3 s À1 . This roughness value has been used in steep gravel-bed rivers (Cenderelli & Wohl, 2001). We also estimated the sensitivity of 1D flow simulations in HEC-RAS with respect to contributions from tributaries (Mardi, Kali, and Phusre Khola) compared to a model run for the Seti Khola without any tributaries.
We used the spatial extent and elevation of the Pokhara Formation in tributary valleys to estimate the Q p of the medieval outburst floods. The field evidence for flooding included slackwater deposits mapped by Stolle et al. (2017) and two outcrops of mainly breccious conglomerates we examined near the Mardi and Phusre confluences.
We computed the horizontal distance between mapped slackwater deposits and simulated flood extents for Q p > 50,000 m 3 s À1 . We assume that the discharge that emplaced the deposits was associated with the simulated flood that resulted in the minimum mean horizontal distance. This analysis provides conservative estimates of flood size, as the elevation of slackwater deposits in the tributary valleys place minimum constraints on flood stage due to potential erosion of the fan deposits after the Pokhara Formation had been emplaced.
The topography of the Pokhara Basin prior to the catastrophic infill is unknown, but was dominated by an older, potentially less dissected alluvial fan of Ghachok Formation material (Fort, 1987). In order to test the sensitivity of our simulations of possible medieval outburst floods with respect to the currently dissected topography, we performed ANUGA modelling with the specifications described in At the Seti dam, our detailed model predicts supercritical flow of 500 m 3 s À1 and flow velocities of 0.2-6.7 m s À1 (Figure 6b). Magdi Khola for Q p > 50,000 m 3 s À1 and match sedimentary evidence best again for the scenario of Q p = 500,000 m 3 s À1 .

| Model applicability
The resolution and quality of our topographic data affect to first order the accuracy of our hydraulic models. The 5-m ALOS DEM appears suitable for HEC-RAS modelling of outburst floods (Westoby et al., 2014;Zhang & Liu, 2015), but only partly resolves the steep gorges along the Seti Khola fan. Earlier attempts at field surveys with terrestrial laser scanning were hampered by haze and water vapour, while airborne mapping using unmanned aerial vehicles (UAVs) is legally restricted in Nepal. Hence, we found that a handheld laser rangefinder proved most useful and economic during field work. While our topographic data are more detailed than in most previous lake- outburst studies (e.g., Mergili et al., 2011;Somos-Valenzuela et al., 2014;Wang et al., 2018), the DEM does not reflect the topography of the Pokhara Basin at the time of the medieval sediment pulses.
Hence, volumetric estimates are limited, which is a common issue in reconstructing outburst floods (Westoby et al., 2014).
There are several simplifying assumptions inherent to our modelling approach that may influence the interpretation of the results. The sedimentary record of both the May 2012 flood and the medieval sediment pulses indicates that they were likely sediment-laden and transient flows with high erosive potential. However, steady-flow models in HEC-RAS (version 5.0.7) and ANUGA both simulate clearwater flow in stable channels and, hence, ignore potential geomorphic effects during large floods such as changing topography by scouring and deposition (David et al., 2022;Larsen & Lamb, 2016;Lehnigk & Larsen, 2022) or an increase in fluid viscosity due to entrainment of sediments (Mungkasi et al., 2013). The latter might also alter flow mobility and runout (Westoby et al., 2014). The entrainment of high volumes of sediment by outburst floods increases flow volumes (Frank et al., 2015;Westoby et al., 2014). Thus, when using paleostage indicators for event reconstructions, clear-water models might overestimate both peak discharge and the volume of flood water. Further, hazard assessments relying on clear-water models might significantly underestimate the downstream impacts for outbursts of smaller bodies of water, whose momentum, peak flow duration, and flow depth might significantly increase downstream, depending on the amount of sediments available for entrainment (Carrivick, 2010;Frank et al., 2015). However, direct observations and data on the studied events are few, and based on proxies, particularly with regard to the medieval sediment pulses. To find a compromise between uncertainty and simplicity for our preliminary reconstructions we mostly site-specific (Dussaillant et al., 2010;Walder & Costa, 1996) and need additional estimates such as the breach rate and depths or potential water volumes impounded by dams in the headwaters of the Seti Khola, which is unconstrained for the 2012 flood. Given these limitations, our reconstructions of the peak discharge of the May 2012 flood remain first-order estimates, which are likely closer to the upper limit of the potential discharge range. However, our numerical simulations expand on previous empirical assessments of the flood (Kargel et al., 2013;Oi et al., 2014;SANDRP, 2014), which consistently described higher peak discharges (Table 1). Moreover, regression-based estimates of peak discharge from a given dam break may vary by more than an order of magnitude (Walder & Costa, 1996). Our sensitivity analysis for Manning's n in HEC-RAS showed that the model with uniform roughness had a slightly larger flood extent on average ($4%). These differences are more pronounced in areas with for our study area low surface roughness (Manning's n < 0.08) such as grassland, bare floodplains or gravel mining sites, consistent with findings from previous work (e.g., Jha & Khare, 2016;Wang et al., 2018;Wohl, 1998). Although these differences seem minute and a varying Manning's n might be less important for larger Q p scenarios, we argue that roughness may capture impor-
Based on sedimentary evidence, Stolle et al. (2017) (Oi et al., 2012;SANDRP, 2014)  The key limitations to peak discharge of lake outbursts are the dam height and the geomorphology of the upstream catchment, which constrain the maximum water volume to be stored in a paleo-lake (Costa & Schuster, 1988). Based on the height of three relict landslide dams in the uppermost Seti Khola gorge and the distinctly steep-walled and bowl-shaped morphology of the Sabche Cirque, Schwanghart et al. (2016) estimated that former lakes could have stored up to 1 km 3 of meltwater in the headwaters.
The reconstructed discharge of the May 2012 flood is of similar magnitude to reported historic flood discharges, including the Nepalese Dig Tsho GLOF of 1985 (1.6 Â 10 3 m 3 s À1 Q p ; Vuichard & Zimmermann, 1987) and the Tam Phokari GLOF of 1998 (1 Â 10 4 m 3 s À1 Q p ; Osti & Egashira, 2009)  where most of the fatalities and damage occurred, was 3700 m 3 s À1 .
Some 15 km downstream, at the Seti dam in Pokhara city, the estimated, strongly attenuated, discharge is 500 m 3 s À1 . We predict peak discharges that are about a factor of two lower than several (empirical) flood estimates reports (Gurung et al., 2015;Oi et al., 2012;SANDRP, 2014). F I G U R E 1 1 Estimated peak discharges of outburst floods in the Pokhara Valley compared to known or reconstructed discharges of glacial lake outburst floods (GLOFs) in the greater Himalayan region derived from Veh et al. (2022). [Color figure can be viewed at wileyonlinelibrary.com]