Reconstructing historic Glacial Lake Outburst Floods through numerical modelling and geomorphological assessment: Extreme events in the Himalaya

Recession of high‐mountain glaciers in response to climatic change frequently results in the development of moraine‐dammed glacial lakes. Moraine dam failure is often accompanied by the release of large volumes of water and sediment, termed a Glacial Lake Outburst Flood (GLOF). Chukhung Glacier is a small (~3 km2) receding valley glacier in Mt. Everest (Sagarmatha) National Park, Nepal. Unlike many Himalayan glaciers, which possess a thick mantle of supraglacial debris, its surface is relatively clean. The glacier terminus has receded 1.3 km from its maximum Holocene position, and in doing so provided the space for an ice‐contact moraine‐dammed lake to develop. The lake had a maximum volume of 5.5 × 105 m3 and drained as a result of breaching of the terminal moraine. An estimated 1.3 × 105 m3 of material was removed from the terminal moraine during breach development. Numerical dam‐breach modelling, implemented within a Generalised Likelihood Uncertainty Estimation (GLUE) framework, was used to investigate a range of moraine‐dam failure scenarios. Reconstructed outflow peak discharges, including failure via overtopping and piping mechanisms, are in the range 146–2200 m3 s‐1. Results from two‐dimensional hydrodynamic GLOF modelling indicate that maximum local flow depths may have exceeded 9 m, with maximum flow velocities exceeding 20 m s‐1 within 700 m of the breach. The floodwaters mobilised a significant amount of material, sourced mostly from the expanding breach, forming a 300 m long and 100 m wide debris fan originating at the breach exit. moraine‐dam. These results also suggest that inundation of the entire floodplain may have been achieved within ten minutes of initial breach development, suggesting that debris fan development was rapid. We discuss the key glaciological and geomorphological factors that have determined the evolution of a hazardous moraine‐dammed lake complex and the subsequent generation of a GLOF and its geomorphological impact. © 2014 The Authors. Earth Surface Processes and Landforms published by John Wiley & Sons Ltd.

The aim of this paper is to quantify GLOF-generating scenarios and the downstream passage of the outburst flood waves resulting from different dam breach mechanisms using a combination of remote sensing, digital elevation modelling and dam-breach and hydrodynamic modelling. Specifically, this approach involved: (i) reconstruction of the post-GLOF characteristics (topography and the bathymetry) of a drained moraine-dammed glacial lake using terrestrial photogrammetry; (ii) numerical dam-breach and hydrodynamic modelling, executed within a probabilistic framework, to quantify a spectrum of flood generation scenarios and downstream flood pathways associated with this event; and (iii) to use the results of (i) and (ii) in combination with field-and remote sensing-based investigation to interpret the glaciological, hydrological and glaciological evolution of a drained moraine-dammed glacial lake complex. A further aim of the article is to place the development of the GLOF hazard at Chukhung Glacier in its glaciological and geomorphological context at both the regional and global scale.
Studies using both numerical modelling and geomorphological and sedimentological characterisation to reconstruct sudden-onset outburst floods from glacial sources are rare Wohl, 2001, 2003;Carrivick, 2006Carrivick, , 2007aAlho and Aaltonen, 2008;Worni et al., 2012;Dunning et al., 2013). The research presented in this article is among the first to present a truly holistic interpretation of the key factors that have determined the evolution of a hazardous moraine-dammed lake complex, as well as the subsequent generation and geomorphological impact of a GLOF. The glacial geomorphology of neighbouring glaciers has already been described in detail (Benn and Owen, 2002;Hambrey et al., 2009) but Chukhung Glacier has received relatively little attention to date.

Study Site
Glaciological and hydrological context The Himalayan Range comprises the most extensively glacierised area outside of the polar regions and contains an estimated 60 000 km 2 of glacier or perennial surface ice  and widespread mountain permafrost (Jakob, 1992;Ishikawa et al., 2001;Fukui et al., 2007;Regmi, 2008). Most glaciers, particularly in the eastern and central Himalaya, are of the summer-accumulation type and gain mass predominantly through Asian summer-monsoon snowfall . This seasonal distribution of precipitation results in the coincidence of the maximum accumulation and ablation periods (Ageta and Higuchi, 1984;Benn and Owen, 1998), with 80% of the total annual precipitation typically occurring between June and September (Ageta, 1976).
Many Himalayan glaciers possess heavily debris-covered tongues ( Figure 1); a consequence of oversteepened slopes and active tectonics which release considerable quantities of frost-shattered rock through frequent ice-and snow-avalanche activity (Benn et al., 2003;Hambrey et al., 2009;Quincey et al., 2009;Bolch et al., 2012). Supraglacial debris affects glacier response to climate change by altering surface (or subdebris) ablation rates (Østrem, 1959;Benn, 2006, 2012;Mihalcea et al., 2008;Brock et al., 2010;Reid and Brock, 2010;Benn et al., 2012). The maintenance of a glacier surface gradient of less than 2°, combined with a stagnating flow regime has been shown previously to be conducive to the development of a supraglacial pond network (Reynolds, 2000;Quincey et al., 2007). In turn, this may lead to the formation of a fully-formed moraine-dammed glacial lake. Alternatively, a moraine-dammed lake may develop in the proglacial area between a terminal moraine and a receding glacier terminus as glacial meltwater ponds in the newly deglaciated basin. However, this mode of lake formation requires that the moraine is stable enough to resist seepage-or piping-induced failure, and to maintain a spillway on the dam crest (Costa and Schuster, 1988;Clague and Evans, 2000;Korup and Tweed, 2007).
Most Himalayan glaciers have been receding since the mid-19th century, resulting in widespread mass loss and a reduction in total glacierised area (Bolch et al., 2008. Kääb et al. (2012) estimated a Himalayan-wide thinning of 0.26 ± 0.06 m (w.e.) yr -1 (equivalent to a sea-level rise of 0.035 ± 0.01 mm (w.e.) yr -1 ). Glaciers in eastern Nepal and Bhutan are estimated to have thinned by 0.3 ± 0.09 m yr -1 . The total glacierised area in the Khumbu Himal in 2005 was estimated as 87.39 km 2 , of which 58.1% and 41.9% was 'clean' and debris-covered, respectively.
In recent decades, a number of moraine-dammed lakes have developed in the Khumbu Himal (Watanabe et al., 1994(Watanabe et al., , 1995Thompson et al., 2012). Nepal's largest moraine-dammed glacial lake, Tsho Rolpa, is located in the neighbouring Rolwaling Himal, and contained an estimated~10 7 m 3 of water at the turn of the century, when it was still growing rapidly (Reynolds, 1998;Rana et al., 2000;Sakai et al., 2000). A rapidly expanding moraine-dammed lake is developing on the snout of Ngozumpa Glacier (Benn et al., 2001;Thompson et al., 2012). If expansion continues unchecked, this lake has the potential to attain a volume in excess of~10 8 m 3 within the next two to three decades Thompson et al., 2012).

Chukhung Glacier and moraine complex
Chukhung Glacier is located in the north-east of Sagarmatha (Mt. Everest) National Park (Figure 1) in the Khumbu Himal, Nepal. The park has an area of~1450 km 2 and ranges in elevation from 2800 m to 8848 m at the summit of Mt. Everest. The glacier is approximately 2.5 km long, has an altitudinal range of approximately 5100 to 6000 m, and is situated in a wide, north-facing amphitheatre. It is connected to Ama Dablam Glacier in its upper south-western reaches (Hambrey et al., 2009). A breached latero-terminal moraine-dam complex is located approximately 1 km in front of the contemporary glacier terminus ( Figure 2). The moraine is topographically complex and contains at least two drained lake beds, additional breaches, and bedrock outcrops (Figure 3). The moraine dam was breached prior to the acquisition of the first satellite imagery of the region in 1962 by the Corona mission, in which the breach and outburst fan are visible. An alluvial debris fan is located at the exit of the main breach (Br_1), and extends for a distance of~200 m onto the floodplain. A brief description of the fan is provided by Hambrey et al. (2009). In summary, fan lithofacies are characterised by a high proportion of subangular and subrounded bouldergravel and sandy boulder-gravel, with the latter prevailing towards the lower extent of the fan. Incision of the fan by an outflow channel reveals the fan to be moderately well-sorted.
Beyond~200 m, superficial sediment coverage becomes increasingly sparse towards the confluence with the Imja Khola, at a distance of~700 m from the breach. The floodplain is 130 m wide along its length, and is bounded to the east by a mountainside spur, and to the west by lateral moraines associated with Ama Dablam Glacier. Outburst flood deposits including large boulder-jams are confined to the floodplain, and deposits further downstream are indistinguishable from other in-channel and overbank material along the Imja Khola.

Structure-from-Motion photogrammetry
Structure-from-Motion (SfM) photogrammetry (Snavely et al., 2008) was used to generate 3-D digital terrain models (DTM) of the moraine and floodplain. This innovative photogrammetric method, which combines multi-view stereo and computer vision algorithms is capable of producing 3-D point clouds from multiple overlapping images with spatial resolutions, precisions and accuracies commensurate with airborne LiDAR (James and Robson, 2012;Westoby et al., 2012;Fonstad et al., 2013;Javernick et al., 2014). In this study we utilise SfM-derived DTMs to: (a) reconstruct the pre-GLOF bathymetry of Chukhung Tsho; (b) aid geometric characterisation of the moraine-dam, both required for input to higher-order numerical dam-breach modeling; (c) quantify the volume of material removed from the moraine during breach development, and d) represent (post-GLOF) floodplain topography for two-dimensional hydrodynamic simulation.
Digital terrain model generation A total of 2056 photographs of the proglacial area were taken using a consumer-grade 12-Megapixel digital camera (Panasonic DMC-G10). Automatic focusing and exposure settings were enabled. Photographs were divided into seven batches to improve  (Astre, 2010;Westoby et al., 2012) was used to produce sparse and dense point-cloud datasets. Forty-nine 1 m 2 numbered yellow targets with a centred red cross were distributed in a quasiuniform pattern across the site (34 across the moraine, 15 across the floodplain) and were used as ground control points (GCP). GCP locations were surveyed using a Leica GPS1200 and postprocessed using data from a static Trimble R7 GNSS base station situated at the northern-most extent of the site. Three-dimensional GPS point qualities averaged 0.478 m and 0.002 m in xyz for the moraine and floodplain GCPs, respectively (Table I).
Following the identification of ground-control points (GCP) in the dense 3-D point-cloud, these data were georeferenced to an absolute co-ordinate system (UTM Zone 45N) using a rigid body transformation (Horn, 1987). The Topographic Point-Cloud Analysis Toolkit (ToPCAT; Brasington et al., 2012;Rychkov et al., 2012) was used to extract detrended minimum grid-cell elevations at 1 m 2 spatial resolution, in order to improve the computational efficiency of surface interpolation in ArcGIS®, while also preserving sub-grid statistics.
Moraine and floodplain topography The final DTM extends from the ice-distal southern lateral moraine of Lhotse Glacier and the adjacent Imja Khola, to the breached terminal moraine of Chukhung Glacier (Figure 4). Topographic features resolved in the DTM include the breach  Figure 4. A two-dimensional perspective of the final hill-shaded, Structure-from-Motion-derived digital terrain model of the Chukhung moraine dam complex and floodplain. These data were used to aid geometric characterisation of the moraine dam, to reconstruct the bathymetry of Chukhung Tsho (inset), and to serve as the topographic domain for hydrodynamic GLOF simulation. Data are georegistered to UTM Zone 45N geographic coordinate system. This figure is available in colour online at wileyonlinelibrary.com/journal/espl through the terminal moraine, multiple, mostly abandoned, anabranching channels, and moraine ridges that are likely to have formed during a Lateglacial advance of Ama Dablam Glacier. Reconstruction and georegistration results are shown in Table I. For the moraine, the dense reconstruction generated >16 million points, while 44 million points were recovered to quantify the floodplain topography. Six of the 15 GCPs that were surveyed across the floodplain, and 10 of the 34 GCPs that were surveyed across the moraine prior to photograph acquisition were used to generate independent SfM-DTM quality metrics using standard point-in-grid estimates, returning residual vertical errors of 0.085 m and 0.814 m for the floodplain and moraine surfaces respectively. The average Euclidean point spacing of the final merged SfM-DTM was 0.054 m (σ = 0.068 m).
The increased reconstruction density and accuracy of the floodplain is explained by the acquisition of input photographs from elevated positions along the margins of the floodplain. Elevated perspectives, coupled with relatively short camerato-feature and camera-to-camera baselines, minimised line-ofsight losses and facilitated the successful 3-D reconstruction of floodplain topography, which is supported by typically higher point densities (>100 points per m 2 ) and improved accuracies when compared to the moraine model. Following decimation, the final number of points used as input to SfM-DTM surface interpolation for the moraine was 300,000, and for the floodplain numbered~460 000. Point densities were highest on ice-proximal faces of the terminal moraine ( Figure 5; >500 per m 2 ), a fact explained by the increased photographic density in these areas. Increased point densities on the floodplain are attributed to an abundance of texturally complex terrain (i.e. cobbles and large boulders) and dense photographic coverage. Point densities of 7, 38, 144, 313 and 443 points per m 2 correspond with the 10th, 25th, 50th, 75th, and 90th percentiles of the point density distribution ( Figure 5). Locally, densities commonly exceeded 200 points per m 2 on the valley floor, while densities across the lowest sections of the valley flanks exceeded 500 points per m 2 . The textural homogeneity of the distal latero-terminal moraine faces and a general lack of elevated and oblique photograph perspectives led to comparatively sparse (<30 points per m 2 ) recovery of 3-D point data in these areas; a methodological shortcoming which also applied to the topographic reconstruction of large areas of the deglaciated moraine basin floor ( Figure 4).

Moraine-dammed lake reconstruction
To calculate the pre-GLOF elevation of the dam crest and reconstruct pre-GLOF moraine geometry, the breach was artificially filled in ArcGIS. All 3-D points located within the breach boundaries were manually removed, and surface point-interpolation was used to 'fill' the edited SfM-DTM. The reconstructed spillway, or crest elevation for the moraine was 4906 m. Direct raster-based differencing of pre-and post-GLOF DTMs permitted the calculation of the volume of material removed from the moraine dam, which is estimated as~125 600 m 3 . Given the inevitable post-GLOF slumping and reworking of the breach side-walls, this figure represents a conservative estimate of the breach volume immediately after its formation.
Methods for quantifying glacial lake bathymetry include the interpolation of spot depths calculated from manual 'plumb line' surveys (Yamada and Sharma, 1993;Fujita et al., 2009)  or GPS-integrated SONAR (Lamsal et al., 2011;Sawagaki et al., 2012). Such methods are appropriate for the survey of extant glacial lakes, but quantification of the bathymetry of drained lakes necessitated an alternative approach. The bathymetry of Chukhung Tsho was reconstructed through the interpolation of SfM-DTM contours spaced at 1 m elevation intervals ( Figure 4). Assuming zero freeboard prior to dam failure, the maximum lake surface area was 65.2 × 10 4 m 2 , and is associated with a maximum lake volume of 5.5 × 10 5 m 3 . This volumetric estimate is based on the post-GLOF SfM-DTM, which includes the breach that connects L3 and L2, and that we infer was formed as water flowed from sub-basin L2, through the expanding Br_2, before exiting the main breach (Br_1) via L3 ( Figure 3). Although we have not calculated the volume of Br_2, it is assumed that it contributes relatively little to our estimate of the volume of water that comprised Chukhung Tsho. Accordingly, the reconstructed volume represents an upper estimate.
The maximum water-depth was calculated as 24 m and volume-stage relationships ( Figure 6) were derived for inputs to hydraulic simulations using HR BREACH (see 'Twodimensional hydrodynamic modeling' below). The reconstructed extent of Chukhung Tsho corresponds well with the trace of a former lake bed, which is observable on Corona satellite imagery from 1962 (Figure 7), and with well-defined palaeo-shorelines that were observed during field investigation (Figure 2(e); Figure 3).
Geomorphological mapping and sedimentological investigation Fine-resolution (1 pixel = 0.5 m), georeferenced GeoEye satellite imagery acquired on 3 November 2009 was used as a base layer for local-scale geomorphological mapping (Figure 3). Mapping was augmented with SfM-derived DTMs and field photographs. The sedimentological characteristics of the moraine and basin-floor sediments were sampled at three locations; the side walls of the breach through the terminal moraine (Br_1), a second breach located~75 m south of Br_1 (Br_2), and a third breach (Br_3) which links the upper (L3) and central (L2) lake basins (Figure 2(f); Figure 3).

Dam-breach modelling
Dam-breach model The current generation of numerical dam-breach models are fully physically based, and employ a combination of geomechanical, hydraulic and erosion or sediment-transport equations to simulate breach development and calculate the breach hydrograph (Mohamed et al., 2002;Hanson et al., 2005;Westoby et al., 2014). These models require extensive parameterisation to define both the initial conditions and material characteristics of the dam, many of which are quantifiable only through rigorous field investigation. Studies linking advanced numerical dam-breach models to modes of morainedam failure are rare (Worni et al., 2012).  Here, we employ the numerical dam-breach model, HR BREACH (Mohamed et al., 2002) within a Generalised Likelihood Uncertainty (GLUE) simulation framework (Beven and Binley, 1992;Kuczera and Parent, 1998;Romanowicz and Beven, 1998;Beven and Freer, 2001;Aronica et al., 2002;Blazkova and Beven, 2004;Beven, 2005;Hunter et al., 2005) to explore range of failure mechanisms for the Chukhung moraine-dam. A novel feature of HR BREACH is an in-built stochastic parameter sampling module, which enables Monte Carlo analysis of the parameter space of specified input dammaterial parameters. These parameters include material density, cohesion, porosity, median diameter (D 50 ), roughness, and the internal angle of friction (Table II). A weir coefficient and an erosion width-to-depth ratio also require specification. We adopted this probabilistic approach because of the inherent uncertainty that surrounds the dam-material properties due to the heterogeneous nature of moraine composition, and the lack of detailed observations of the breaching process.
Initial dam geometry and lake bathymetry were extracted from the SfM-DTM (Figure 4). Reconstruction of pre-GLOF dam geometry facilitated the estimation of dam crest width (20 m) and proximal and distal face slope ratios (1:1.9 m/m and 1:2.8 m/m, respectively), which represent mean values for five elevation profiles taken across the moraine dam and were used for geometric characterisation of the moraine dam during the setup of the numerical dam-breach model. A dam length of 150 m was assumed, and a downstream valley slope of 1:0.09 m/m was extracted from the SfM-DTM.
Input parameter ranges (Table II) were derived through a combination of field investigation and published values for moraine dams (Xin et al., 2008;Worni et al., 2012). A uniform probability distribution for each parameter was assumed in the absence of any a priori knowledge of the dam-material composition. HR BREACH employs a one-dimensional hydrodynamic model to simulate breach water-profiles for an overtopping failure (Mohamed et al., 2002;Morris et al., 2008) and includes a broad-crested weir discharge equation to represent flow through the breach and a simplified version of the Saint-Venant equations to simulate breach outflow as it descends the distal face of the dam. The model includes various options for sediment transport and erosion equations to model the erosion of the dam material during breach development. In this study, Chen and Anderson's (1986) erosion equation for the calculation of continuous breach erosion through non-cohesive sediments was used to model sub-aerial breach development, such as that associated with failure that has been initiated by wave overtopping. The equation takes the form: where ε r is the detachment rate per unit area (e.g. m 3 s -1 /m -2 ), τ e is the flow shear stress (kN) at the breach boundary, τ c is the critical shear stress required to initiate particle detachment (kN), and K d and a are coefficients based on the sediment properties. HR BREACH accounts for bending (or 'tensional') and 'shear' type failures (Mohamed et al., 2002). The former are represented by a moment of force, M o : where W e is the weight of a dry block of soil (kN), W s is the weight of a saturated block of soil (kN), W u is the weight of a submerged block of soil (kN), H 1 is the hydrostatic pressure force in the breach channel (kN), H 2 is the hydrostatic pressure force inside the dam structure (kN), e, e s and e u are weight force eccentricities (m), e 1 and e 2 are hydrostatic pressure force eccentricities (m) (Hassan and Morris, 2012). For non-cohesive material, there is a higher likelihood of shear failure of the breach walls occurring during breach development. HR BREACH calculates shear failure through the analysis of Factor of Safety (FoS) values using the following equation: where c is soil cohesion (kN/m 2 ), L is the length of the failure plane (m), and Φ is the soil angle of friction (°). HR BREACH can also simulate piping-style breach development. Piping is initiated when water seeping through a dam structure removes the finest sediment fractions, eventually forming a pipe, through which water drains. Flow through the resulting orifice is calculated as: with: where Q b is discharge (m 3 s -1 ), g is acceleration due to gravity (m/s 2 ), A is the cross-sectional area of the pipe (m 2 ), H is the water level in the dam (m), H p is the pipe centre-line elevation (m), h l represents losses due to friction and contraction, f is a friction coefficient determined as a function of the material D 50 (mm), L is the pipe length (m), and D is the pipe diameter (m) (Mohamed et al., 2002). The initial pipe may then evolve by erosion which is assumed to occur uniformly along the pipe. The volume of sediment eroded (V s ) per time step is taken to be: where Q s is the sediment transport rate (m 3 s -1 ), and Δt the model time step (s). This represents a highly simplified approximation of pipe expansion mechanisms (Fread, 1988), and further progress is necessary to improve the representation of the processes driving pipe-expansion processes (Morris, 2009;Westoby et al., 2014). Material from the downstream dam face that falls into the expanding pipe is removed by the flow. Transition from breach enlargement through pipe-flow to sub-aerial, overtopping-type flow occurs when the force of the water in the lake overcomes the shear strength of the pipe overburden, or this material collapses under its own weight (Mohamed et al., 2002).
Experimental design Four breaching scenarios were parameterised to capture a range of dynamical breach conditions. In the first, CK_control, breaching is initiated as water overtops the moraine dam through a predefined notch (1 m wide, 1 m deep). In the remaining scenarios (CK_ 4890, CK_4895 and CK _4900), breaching is initiated via piping. The numerical value in each of the scenario IDs denotes the elevation at which the initial pipe is located. The reconstructed elevation of the dam crest was 4906 m, while the elevation of the dam toe was set to 4854 m. An initial pipe diameter of 0.3 m was used; this was the minimum diameter which was found to produce a numerically stable solution. One thousand Monte Carlo simulations were executed per scenario. This number of runs was deemed an appropriate compromise, which satisfied the demands of a comprehensive search of the model parameter space and practical limitations of the computational time (each 1000 runs requiring 8 h of simulation time with a 10 s model timestep.) In the absence of observed hydrodynamic data describing the actual breach, the performance of each simulation was evaluated using two reconstructed morphometric descriptors and classified using simple geometrical likelihood functions. First, simulations that reproduced the final upstream breach depths of 20 m ±1 m were deemed acceptable, or behavioural, and were retained. This error margin was intended to account for observed errors in the GCP positional accuracies and transformation residuals of the SfM data (Table I). A triangular likelihood function was used to quantify simulation performance on a linear scale of 0-1, centred on an observed final breach depth of 20 m.
Secondly, the residual sum of squared errors (RSS) of the modelled versus the observed final elevation profile of the breach centreline was calculated for each behavioural simulation: where Z obs is the observed elevation profile and Z pred. is the modelled elevation profile. Here, a linear function was used to classify the likelihood of each simulation (between 0 and 1), scaled to reflect the range between a perfect simulation of the breach profile and the worst case behavioural simulation. These two likelihood values were combined using Bayesian updating (after Lamb et al., 1998) to produce a multi-criteria, unified likelihood value. It should be noted that this method assumes that the elevation profile of the breach floor has undergone no subsequent alteration after breaching had ceased. The method also assumes that Chukhung Tsho was drained in a single, largely uninterrupted dam-breaching event, and that modelled peak discharges are representative of such a mode of drainage. Behavioural breach-hydrograph data were used to construct cumulative probability distribution functions (CDFs) for each scenario. For each time step, discharges were sorted by magnitude, and weighted by their final likelihood values, to support the derivation of time step-specific CDF curves quantifying the potential range of discharge values at each time step.

Two-dimensional hydrodynamic modelling
Hydrodynamic model Two-dimensional models based on the depth-averaged shallow water equations (Chanson, 2004;Hervouet, 2007) provide a proven framework for predictions of distributed flow hydraulics.
The advantages of using 2-D numerical models for hydrodynamic simulation of sudden-onset floods include their ability to simulate multi-directional and multi-channel flow, superelevation of flow around channel bends, hydraulic jumps, transcritical flow regimes, and flow recirculation (see for example, Bohorquez and Darby, 2008). A number of two-dimensional models include sediment transport functionality, making them ideally suited for quantifying the geomorphological impact of such flood flows.
For GLOF simulation, we employed a MacCormack-Total Variation Diminishing (TVD) scheme in ISIS 2D (Halcrow, 2012). This uses predictor and corrector steps to compute flow depth and discharge across a regular grid. A TVD term is added at the corrector step to suppress numerical oscillations near sharp gradients, therefore making it suitable for the simulation of rapidly evolving, transcritical and supercritical flows (Liang et al., 2006Liao et al., 2007).

GLOF reconstruction
Breach-outflow hydrograph data were used as the upstream hydraulic boundary condition for two-dimensional hydrodynamic modelling. A single breach hydrograph (likelihood value = 0.95) was used to provide a representative assessment of the near-field hydrodynamic characteristics of a GLOF from Chukhung Tsho. Floodplain topography was represented by a SfM-derived DTM with a spatial resolution of 4 m 2 , which encompassed the width (~300 m) and length (~700 m) of the floodplain. A model time step of 0.04 s was specified to satisfy the standard Courant-Friedrich-Levy stability criterion. At the downstream model boundary, flow was able to exit the model domain freely, through linkage with an ASTER GDEM-derived DTM, thereby eliminating artificial backwater effects. Sediment transport functionality was not available, and so the GLOF was simulated as a clearwater-type flow.

Dam-breach modelling
Simulation retention statistics revealed that 8.8%, 6.3%, 6.1% and 1.7% of each set of 1000 Monte Carlo simulations were retained as behavioural for CK_control, CK_4890, CK_4895 and CK_4900, respectively. These results indicate that the simulation of breach initiation through piping, and decreasing pipe elevation, leads to a decrease in the number of retained parameter ensembles. Maximum non-behavioural values of Q p ranged from 2000-2200 m 3 s -1 , with a corresponding T p range of 12-16 min. The lowest non-behavioural values of Q p were in the range 146-164 m 3 s -1 , within a T p range of 50-65 min.
Model evaluation resulted in a significant refinement of Q p and T p across all scenarios. Analysis of the likelihood response surfaces for each input parameter revealed a high degree of scatter, and therefore no obvious structural relationships in the parameter space. The exception is Manning's n, which displayed an underlying linear relationship with increasing likelihood. The behavioural range of Manning's n for all scenarios was 0.027-0.049, and is representative of the pipingstyle failure scenarios. Behavioural Manning's n values for CK_control were 0.042-0.049.
A decrease in maximum and minimum behavioural values of material roughness is observed following the introduction of piping, and decreasing pipe-elevation. The lowest values were attributed to CK_4890 (0.021-0.038). This observation is consistent as lower roughness coefficients serve to offset rapid pipe-development and breach enlargement, which results from an increasing pressure-head in the pipe with decreasing pipe elevation. This feedback was confirmed by a relatively homogeneous spread of behavioural Manning's n values along the parameter axis for CK_4895. Final simulation-specific likelihood values ranged from 0.01-0.95 for CK_control, CK_4890 and CK_4900. The likelihood range for CK_4895 was 0.81-0.94.
Behavioural hydrograph envelopes of the breach and percentile hydrographs are displayed in Figure 8. Behavioural hydrographs for CK_4890 and CK_4900 show a tendency towards a dual-peak form, with a dominant early peak. By comparison, the form of the behavioural hydrographs for CK_control and CK_4895 are broadly similar, though the time to concentration is longer for the control scenario. A rapid rise-to-peak and increased maximum peak discharges are associated with piping-style failures (specifically CK_4890 and CK_4895). These characteristics may be the result of a large initial pipe diameter (0.3 m), and a rapid transition from confined to unconfined breach-development. To improve the physical robustness of the model, stability should be ensured when specifying extremely small initial pipe diameters, as our results provide little indication as to the precise breach initiation time associated with piping-style failures.

GLOF modelling
Inundation of the entire floodplain occured within~10 minutes ( Figure 9). Maximum inundation corresponds with the outburst hydrograph peak (1915 m 3 s -1 ). The relatively short length (~0.7 km) of the floodplain was conducive to the rapid downstream translation of the breach hydrograph. Maximum flowdepths (10-12 m) occur at~15 min, and are largely confined to a natural depression located approximately 550 m from the breach. Upstream, increased flow-depths are concentrated in a single channel that traces the eastern side of the floodplain. We reiterate that the underlying topography is a post-GLOF reproduction of the floodplain, and so the observed flow-concentration is a reflection of the contemporary channel, as opposed to a true reflection of the characteristics of the undated GLOF. This caveat applies to interpretations of other flood characteristics (e.g. flow velocities, clast mobility), and is unavoidable in the absence of pre-GLOF topographic data. Flow depths gradually decrease with a decrease in breach-hydrograph discharge. However, the entire floodplain remains submerged following the cessation of breach outflow as a consequence of the restricted exit of the reach (Figure 9). In addition, The highest flow-velocities (16-20 m s -1 ) are generally confined to the eastern edge of the floodplain, but are more equally distributed than the associated flow-depths (Figure 9), and occur between~9 and 12 min. Increases in flow-velocity cease immediately following the breach-hydrograph peak, and for the remainder of the simulation remain comparatively low (<4 m s -1 ). Increased flow velocities (8-12 m s -1 ) are observed at the downstream boundary of the model domain, and are associated with the significant constriction in flow width as the valley floor narrows between the lateral moraine of Lhotse Glacier, and the eastern lateral moraine of Ama Dablam Glacier (Figure 9).
The simulated GLOF represents a clear water modelling scenario because of the numerical limitations of the hydrodynamic model used. Simulating GLOFs as clearwater flows is a questionable practice, as they are widely documented to entrain significant volumes of morainic material during breach development and downstream flow propagation, which induces fundamental changes in their flow rheology and hydraulic behaviour (O'Connor et al., 2001;Kershaw et al., 2005;Carrivick et al., 2011;Westoby et al., 2014). The geomorphology of the valley indicates that the majority of sediment removed from the breach was deposited within a few hundred metres of the breach (Figure 2(h)). Comparatively limited reworking of the distal floodplain environment indicates that the GLOF was unlikely to be heavily sediment laden at this point. Despite the modelling limitations described above, the calculation of stream power, in combination with the application of an indicative scale to illustrate the maximum sediment fraction that would be mobilised as bedload for each grid-cell (after Carrivick et al., 2010), permitted a first-pass assessment of the potential morphodynamic impact of a GLOF from Chukhung Tsho (Figure 9). Our results indicate that the inundation of the floodplain may be associated with the mobilisation of gravel-and cobble-sized material. However, the rapid progression of the breach-hydrograph would result in the mobilisation of boulders in excess of 2.56 m diameter within a distance of~500 m from the breach, and across a zone measuring up to 100 m in width. Hydrograph recession is associated with a decrease in stream power and an accompanying decrease in sediment mobility. At 18 minutes, the GLOF remains capable of mobilising clasts ≤ 64 mm at its downstream extent (Figure 9).

Palaeocompetence-based GLOF reconstruction
The results of a palaeocompetence-based reconstruction of ata-point GLOF discharges along the length of the floodplain are displayed in Figure 10. The length of the intermediate axis of the 10 largest boulders encountered along transects spaced at 50 m intervals in a down-valley direction across the floodplain from the breach to 150 m upstream of the confluence with the Imja Khola were recorded, and the arithmetic mean of these data using for input to palaeoGLOF reconstruction. Our reconstruction is based on equations presented by Costa (1983), which enable the estimation of flow velocity, palaeodepth and discharge from particle size data: where v is flow velocity (m s -1 ), d i is intermediate (b-axis) boulder diameter, ω is unit stream power (N m s -1 ), D is flow depth (m), γ f is the specific weight of water (9800 N m 3 ) and S is channel slope (m m -1 ). Wetted channel perimeter and crosssectional area were calculated from the SfM-DTM, and Manning's n of 0.05 was assumed. Reconstructed discharges (Figure 10) largely mirror the modelled pattern of a down-valley decrease in the sediment transport capacity of the escaping floodwaters (Figure 9). These results reflect the dissipation of flow energy with lateral flow expansion and travel distance and are characteristic of the sedimentology of outburst flood deposits observed both in the Himalaya (Cenderelli and Wohl, 2003) and further afield (Kershaw et al., 2005;Carrivick, 2007aCarrivick, , 2007bBreien et al., 2008).
The GLOF from Chukhung Tsho was capable of transporting significant volumes of sediment, and it is assumed that the entire volume of material removed from the breach was deposited within~0.5 km of the terminal moraine from which it was mobilised. However, the absence of a pre-GLOF topographic dataset of the floodplain, and a hydrodynamic model with morphodynamic capability precludes an in-depth assessment of the temporal and spatial characteristics of landscape change in the proglacial area (Breien et al., 2008;Carrivick et al., 2010Worni et al., 2012;Dunning et al., 2013).

Discussion and Interpretation
Evolution of the GLOF hazard at Chukhung Glacier Moraine and basin construction The lack of glaciotectonic structure in the exposed terminal moraine indicates a single period of moraine-building, most likely attributable to glacier re-advance during the Little Ice Age (LIA). There are no older Holocene or Late Glacial Maxima moraines associated with Chukhung Glacier in the immediate vicinity. It is conceivable that Chukhung Glacier overrode older, smaller moraines during its LIA advance. The tapering planform of the Chukhung moraine loop is explained as the result of lateral confinement to the west by the eastern lateral moraine and bulk of Ama Dablam Glacier, and to the east by the western flanks of a ridge that extends northwards from Amphu peak (5663 m; Figure 3). The deposition of a continuous moraine loop, followed by late-Holocene glacier recession allowed a significant volume of glacial meltwater to pond in the deglaciated moraine basin.
Character of palaeo-glacial lakes and GLOF trigger mechanisms The recession of Chukhung Glacier from its maximum LIA position may have been rapid, since short and steep glaciers are far more responsive to climatic change than longer glaciers with shallower surface gradients (Benn and Evans, 1998). In addition, the expansion of an ice-contact proglacial lake is likely to have been associated with the development of lake-terminating ice cliffs. Such cliffs will have accelerated the rate of terminus retreat as discrete subaerial and subaqueous ice-calving and progressive backwasting became the dominant modes of terminus retreat (Robertson et al., 2012;Thompson et al., 2012).
Three-dimensional reconstruction of the bathymetry of Chukhung Tsho indicates that the moraine-dammed lake would have extended as far south as the head of the central sub-basin L2 (Figure 4). A bedrock outcrop that spans the central section of the basin would have limited further expansion. Volumetric and areal estimates assume zero dam freeboard prior to failure, and also that the glacier terminus had receded beyond the head of the reconstructed lake. These estimates therefore represent an upper limit.
The southern-most breach (Br_3) and the debris fan that extends from its exit into L2 are absent in the 1962 image  Figure 7). As Chukhung Tsho had drained prior to 1962, the formation of Br_3 is unrelated to the undated GLOF, and was most likely formed as a result of overspill from the uppermost lake (L3), and in turn possibly generated a 'mini' GLOF. Possible trigger mechanisms might include the generation of a displacement wave from collapsing ice séracs perched on the bedrock outcrop above L3, or a sudden rise in lake level as the result of the rapid input of a meltwater pulse from Chukhung Glacier. Formation of this breach by prolonged down-cutting through the till drape is discounted because of the presence of the debris fan, which comprises predominantly gravel-and cobble-sized material that could only have been mobilised by significant flows. The cascading of water from L1 into Chukhung Tsho, with an associated rise in lake level, and the overtopping of the terminal moraine is therefore discounted as a trigger mechanism for the initiation of terminal moraine breaching, and generation of the undated GLOF.
Undated moraine-dam failure Numerical modelling Perhaps the most striking finding of the dam-breach modelling results is the short duration of the piping-failure hydrographs (Figure 8), of which the longest duration of was 13 min (CK_4890). In comparison, behavioural overtopping hydrographs lasted between 12.5 and 18 min. This difference is explained as the result of a rapid initial transition from confined to unconfined breach-enlargement during piping simulations. During this transition the dimensions of the initial pipe increased exponentially from the beginning of each simulation, with pipes quickly collapsing (typically within~1 min) under the weight of their overburden. Whether or not this is an accurate representation of the real-world characteristics of pipe expansion is questionable, but is explained here as the result of the specification of a relatively large initial pipe diameter, and large pressure-heads in the pipe at the onset of simulation. Combined, these produced extremely high rates of sediment evacuation and pipe expansion. The rapid transition to an overtopping mode of breach development also accounts for the increased peak discharges which are associated with the piping-simulation results, as the initial dimensions of the newly formed unconfined breach were conducive to the evacuation of elevated volumes of water.
The dominance of the roughness coefficient, Manning's n, in determining peak discharge and time-to-peak is unsurprising. Increasing the boundary resistance to flow will simultaneously increase local flow depths, and therefore rates of subaqueous erosion and side-wall undercutting, thereby accounting for the positive relationship that was observed between this parameter and model output. Given that the input range for Manning's n encompassed values characteristic of channels composed of low-lying vegetation and weeds (0.02), to cobbles and large boulders (0.05; Chow, 1959), the behavioural parameter range for overtopping failure (CK_control) appears to reflect the robustness of HR BREACH for the simulation of this particular mode of breach initiation. The extended range of behavioural roughness values observed for the piping scenarios (0.027-0.049) is unusual, and might suggest that functional relationships between the remaining input parameters exert an overriding influence over breach development for this pipe-initiated breaching.
These hydrodynamic results are broadly consistent with the simulation of similar high-magnitude outburst floods, whereby the highest near-field flow-depths, velocities, and shear stresses are largely confined to a single channel (Alho and Aaltonen, 2008;Carrivick et al., 2010). Beyond~0.5 km from the breach, the reducing definition of a single channel, combined with a reduction in the gradient of the valley floor and increasing lateral flow expansion as a result of increased upstream discharges, produced a widening zone of increased flow depths (Figure 9).
Comparison with other documented moraine-dam failures Geomorphological reworking of valley-floors downstream of moraine-dammed lakes by GLOFs has been widely documented (Clague and Evans, 2000;Richardson and Reynolds, 2000;Cenderelli and Wohl, 2003;Korup and Tweed, 2007;Osti and Egashira, 2009;Worni et al., 2012). In many instances, distinct alluvial debris-fans similar to that observed beneath the main breach at Chukhung are located immediately downstream of moraine (Vuichard and Zimmerman, 1987;Evans and Clague, 1994;Kershaw et al., 2005) and landslide dambreaches (Dunning et al., 2006). In both settings, these fans comprise material predominantly sourced from the dam structures themselves, as rapid energy losses immediately downstream of a dam-breach result in the deposition of the coarsest sediment fraction. Downstream, reaches variously undergo erosion and deposition both during and following GLOF passage. The morphodynamic regime is primarily dependent on the longitudinal distribution of channel expansion and contraction (Clague and Evans, 2000;Cenderelli and Wohl, 2003;Kershaw et al., 2005).
The volume of water released by the reconstructed Chukhung GLOF (5-6.5 × 10 6 m 3 ) is smaller than estimated for the Dig Tsho; Vuichard and Zimmerman, 1987), Sabai Tsho (18 × 10 6 m 3 ; Osti and Egashira, 2009), or Luggye Tsho (18 × 10 6 m 3 ; Watanabe and Rothacher, 1996) GLOFs, and is more comparable with estimates for the Nare GLOF (0.5-5 × 10 6 m 3 ; Buchroithner et al., 1982;Fushimi et al., 1985), and for several historical GLOFs in British Columbia (Clague and Evans, 2000;Kershaw et al., 2005). The volume of impounded water is a direct product of moraine-basin morphology and mode of lake development, and the pre-failure water surface elevation. In contrast, lakes that originate on the tongues of long, low-gradient, debris-covered glaciers have the potential to impound water volumes in excess of 10 7 m 3 . The latter is determined by various factors including dam crest spillway elevation, or the elevation of intra-moraine seepage pathways and pipes

Regional and global context
The Chukhung moraine complex is regionally unique, in that it represents an example of a deglaciated, breached and fully drained moraine-dammed lake. It provides an indication of how other glacierised basins might behave as glaciers recede. Breached moraine-dammed lake complexes are uncommon in the Khumbu Himal, the most obvious exception being Dig Tsho, whose moraine-dam was breached on 4 August 1986 (Vuichard and Zimmerman, 1987). A GLOF also originated from the vicinity of the Nare Glacier on 3 September 1977 (Buchroithner et al., 1982;Wohl, 2001, 2003), and has been attributed to the failure of a series of ice-cored moraine dams. In the neighbouring Hinku valley, a GLOF was produced by the failure of the Tam Pokhari (Sabai Tsho) moraine-dammed glacial lake (Osti and Egashira, 2009;Osti et al., 2011). In comparison, extant moraine-dammed glacial lakes in the Khumbu Himal currently include those found in front of Imja Glacier (Watanabe et al., 1994(Watanabe et al., , 2009Reynolds, 2006;Hambrey et al., 2009), Ngozumpa Glacier Thompson et al., 2012) and Lumding Glacier (Bajracharya, 2008).
Geomorphologically, the Chukhung Glacier-moraine complex resembles those in British Columbia which are typically 1688 M. J. WESTOBY ET AL.
situated in cirques or the upper reaches of steep-walled valleys with low width-to-height dam ratios and steep proximal and distal moraine slopes (Blown and Church, 1985;Evans, 1986;Clague and Evans, 2000;McKillop and Clague, 2007). These characteristics are also typical of many glacier-lake systems in the South American Cordillera Blanca (Lliboutry et al., 1977;Hubbard et al., 2005). The initiation of moraine-dam failure has been described as the result of wave-overtopping in both regions (Lliboutry et al., 1977;Reynolds, 1992;Clague and Evans, 1994;Richardson and Reynolds, 2000), and the development of intra-morainal springs and seepage, which may serve as a precursor to pipe-initiated failure, has also been documented (Clague and Evans, 2000;Lliboutry et al., 1977), although direct observations of piping style moraine-dam failure are scarce. The Chukhung moraine-dam complex represents an endmember in a landform continuum proposed by Clague and Evans (2000) which begins with the development of supraglacial ponds and subsequent development of a 'proto' moraine-dammed lake (cf. Benn et al., 2012), and ends with the failure of a moraine dam and the generation of a GLOF (Benn and Owen, 2002). By comparison, many debris-covered glacier tongues in the region (Hambrey et al., 2009;Benn et al., 2012), and in the wider Himalaya (Sakai and Fujita, 2010) currently host a network of supraglacial ponds, or are undergoing a transition to fully-formed moraine-dammed lake development. The progression of the Chukhung Glacier-moraine-dam complex to the end of this continuum is explained by the nature of its formation; that it proceeded along an alternative and accelerated pathway associated with comparatively steep and largely debris-free parent glaciers ( Figure 11). Glaciers of this type typically respond to increased air temperatures through rapid recession of the glacier terminus. Morainedammed lake development in the proglacial zone of such glaciers is comparatively straightforward and rapid when compared with their heavily debris-mantled counterparts (Clague and Evans, 2000). When combined, rapid lake development and impoundment by a moraine-dam with a high width-toheight ratio, and proximity to slopes which generate frequent ice-and rock-avalanches are conducive to GLOF generation.

Conclusions
The Chukhung Glacier moraine-dam complex in the Khumbu Himal, Nepal, comprises multiple breaches and interconnected lake basins. A numerical dam-breach model, executed within a probabilistic framework, was used to simulate moraine breaching. The results suggest that breaching scenarios associated with both piping and overtopping styles of breach initiation were capable of reproducing observed contemporary breach morphology. The maximum lake volume was 5.5 × 10 5 m 3 with modelled peak discharges in the range 146 to 2200 m 3 s -1 . Two-dimensional hydrodynamic simulations indicate that inundation of the entire floodplain was achieved within 10 min of breach development. The temporal and spatial variability in the sediment transport capacity of the floodwaters supports the deposition of an alluvial fan, which is also observed in the field. Although numerous examples of breached moraine-dam complexes from other regions have been described in the literature, Chukhung Glacier represents a rare example of a Khumbu Himal landsystem that has progressed to the end of a landform continuum which begins with glacier recession and moraine-dammed lake development, and ends with moraine-dam failure and the generation of a GLOF. In contrast, the majority of glaciers in the Khumbu Himal, and many similar adjacent catchments, and those further afield in the Himalaya, currently host fully formed or 'proto' moraine supraglacial or proglacial lakes or a network of supraglacial melt ponds. The uniqueness of Chukhung Glacier is attributed to topographic controls on glacier surface gradient, and the lack of a thick mantle of supraglacial debris that characterises the majority of glaciers in the region. In combination, this means that the glacier is more sensitive to the effects of climatic change than the majority of other glaciers in the region, which are largely debris-covered and possess low surface gradients along most of their length.