Numerical modeling of glacial lake outburst floods using physically based dam-breach models

. The instability of moraine-dammed proglacial lakes creates the potential for catastrophic glacial lake outburst ﬂoods (GLOFs) in high-mountain regions. In this research, we use a unique combination of numerical dam-breach and two-dimensional hydrodynamic modelling, employed within a generalised likelihood uncertainty estimation (GLUE) framework, to quantify predictive uncertainty in model outputs associated with a reconstruction of the Dig Tsho failure in Nepal. Monte Carlo analysis was used to sample the model parameter space, and morphological descriptors of the moraine breach were used to evaluate model performance. Multiple breach scenarios were produced by differing parameter ensembles associated with a range of breach initiation mechanisms, including overtopping waves and mechanical failure of the dam face. The material roughness coefﬁcient was found to exert a dominant inﬂuence over model performance. The downstream routing of scenario-speciﬁc breach hydrographs revealed signiﬁcant differences in the timing and extent of inundation. A GLUE-based methodology for constructing probabilistic maps of inundation extent, ﬂow depth, and hazard is presented and provides a useful tool for communicating uncertainty in GLOF hazard assessment.


Moraine-dammed glacial lakes
Glacier recession is occurring globally as a result of recent climatic change (Oerlemans, 1994;Kaser et al., 2006;Zemp et al., 2009;Bolch et al., 2011. The exposure of terminal and lateral moraine complexes is becoming increasingly commonplace as a result of glacier recession, particularly in high-mountain regions (Hambrey et al., 2008). Moraines reflect the historical maximum extent of a given glacier, and are typically composed of poorly consolidated glacial material. A latero-terminal moraine can present a physical barrier to drainage of glacial meltwater and, in such cases, result in the formation of a moraine-dammed lake (Costa and Schus-ter, 1988) and create the potential for a glacial lake outburst flood (GLOF) hazard (e.g. Clague and Evans, 2000;Westoby et al., 2014a;Worni et al., 2014).
Published by Copernicus Publications on behalf of the European Geosciences Union.
Breaching of a moraine dam can result in the generation of a GLOF (Lliboutry 1977;Vuichard and Zimmerman, 1987;Clague and Evans, 2000;Kershaw et al., 2005;Harrison et al., 2006;Osti and Egashira, 2009;Worni et al., 2012Worni et al., , 2014Westoby et al., 2014a, b). These sudden-onset floods represent high-magnitude, low-frequency catastrophic phenomena that have enormous potential for geomorphological reworking of channel and floodplain environments (Cenderelli and Wohl, 2003;Worni et al., 2012;Westoby et al., 2014b). These floods can pose significant hazards, resulting in the destruction of in-channel and riparian assets, including hydroelectric power facilities, trekking routes, and impact on settlements with ensuing loss of life (Vuichard and Zimmerman, 1987;Lliboutry et al., 1977;Watanabe and Rothacher, 1996). Of the various glacial hazards, GLOFs have far-reaching, distal impacts, with destruction often reported tens or hundreds of kilometres downstream of their source (Vuichard and Zimmerman, 1987;Clague and Evans, 2000;Richardson and Reynolds, 2000).

Modelling glacial lake outburst floods
A number of approaches to model GLOFs have been presented; these are summarised by Westoby et al. (2014a). At the outset, it is worth highlighting that it is currently not expedient to employ our modelling workflow for GLOF analysis at catchment or basin scales but rather as a logical next step based on the results from a broader hazard assessment screening exercise. Such exercises should, in the first instance, identify and quantify the risk posed by individual glacial lakes in a region through the use of multi-criteria hazard analysis techniques (e.g. RGSL, 2003;Reynolds, 2014) and may also involve the application of DEM-based floodrouting models to provide rapid, first-pass assessments of likely patterns of inundation (e.g. Huggel et al., 2004) that form the basis for subsequent detailed GLOF analysis.
Most GLOF modelling approaches typically employ an empirical or numerical dam-breach model to derive a breach outflow hydrograph, and couple this to a numerical hydrodynamic model to simulate the downstream propagation of the flood wave. Such loose model coupling, or "process chain" simulation, is now a well-established practice and has been applied for reconstructive (e.g. Osti and Egashira, 2009;Klimeš et al., 2014;Westoby et al., 2014b;Worni et al., 2014) and predictive GLOF modelling studies (e.g. Xin et al., 2008;Worni et al., 2013).
Dam-breach modelling typically relies on an empirical or analytical formulation (e.g. Walder and O'Connor, 1997). Empirical models return a single breach hydrograph descriptor, usually peak discharge or flood time to peak, to which a hydrograph is fitted. Such models are derived from historical case study data and use simple geometric descriptors of the moraine or lake, such as dam height and lake surface area, as the sole inputs. Their robustness may be questionable because of their explicit reliance on the case-study data and their representativeness. Analytical models are similar to their empirical counterparts, but require that the user know the final breach dimensions and formation time, to which an analytical solution for the rate of breach growth is fitted. This approach suffers equally where limited data are available to describe the breach and its rate of formation (i.e. in particular where dams have yet to fail and the breach form and dimensions must be estimated a priori).
Physically based numerical models represent the current state of the art in dam-breach modelling, and combine numerical erosion and sediment transport models in combination with a 1-D or 2-D flow hydraulics solver to simulate breach expansion and channel flow (Worni et al., 2012;Westoby et al., 2014a, b). Despite their advantages, these models are not yet in widespread use by the glacial hazard communities, with recent studies still preferring to adopt established empirical approaches (e.g. Byers et al., 2013). In large part, this slow rate of adoption may reflect the high data requirements of physically based models. For example, these physics-based simulation require knowledge of the mean particle size (or a full particle-size distribution), internal angle of friction, cohesion, and porosity, amongst other physical attributes of the material. Such data are scale-dependent and at best require detailed field or laboratory investigation, which may be logistically challenging in remote, high-altitude environments and may be impossible to obtain for some reconstruction studies.

Motivation for this study and modelling approach
Numerous sources of predictive uncertainty exist in the parameterisation of contemporary numerical dam-breach and hydrodynamic models when used in either reconstructive or predictive GLOF simulation (Westoby et al., 2014a). Reconstructions of historic events where no field-based or published data exist to describe the geometric and material characteristics of the dam represent the most poorly constrained case. However, the characteristic compositional heterogeneity of moraines implies that even data-rich scenarios are likely to undersample the actual system complexity. Moreover, reconciling the often disparate spatial scales of processes, field observations, and effective model parameters presents a persistent ambiguity even under optimal conditions. For such ill-conditioned problems, an increasingly popular approach is to accept uncertainty in the a priori model parameters, and rather than seek the optimal parameter set through calibration, numerical methods are used to quantify the associated predictive uncertainty and present simulation output probabilistically (Beven, 2005).
To date, few studies have sought to explore the predictive uncertainty of numerical dam-breach models and link this to the downstream consequences in terms of flood wave propagation and inundation (e.g. Wang et al., 2012;Worni et al., 2012;Westoby et al., 2014b). In this paper we address this research gap through the development and demon-Earth Surf. Dynam., 3,2015 www.earth-surf-dynam.net/3/171/2015/ stration of a unifying framework for cascading uncertainty in a GLOF model chain, illustrated for a reconstruction of a major historical outburst flood in the Nepalese Himalaya. We demonstrate how numerical models representing components of the GLOF process can be coupled to provide probabilistic predictions of the breaching process and downstream flood propagation using generalised likelihood uncertainty estimation (GLUE; Beven and Binley, 1992).

Uncertainty and equifinality in numerical modelling
2.1 Types, sources, and the potential significance of uncertainty Pervasive uncertainty surrounds almost all aspects of environmental modelling, and may be broadly classified as aleatory or epistemic . Aleatory uncertainties relate to fundamental random variability in the modelling process. Epistemic uncertainties are associated with a lack of knowledge of system characteristics or architecture. Specific, epistemic uncertainties relate to boundary conditions, initial conditions, parameter values, and model structures (Beven et al., 2014) -in essence, elements that are common to almost all environmental modelling problems. Similar sources of uncertainty exist in the dam-breach and flood-routing components of the GLOF model chain (Westoby et al., 2014a). In the case of the former, appreciable, and predominantly epistemic, sources of uncertainty surround the establishment of initial conditions (e.g. dam geometry, reservoir hypsometry), dam material parameterisation (e.g. grainsize distribution curves, material porosity, density, cohesion, roughness coefficients, internal angle of friction) and the establishment of computational constraints (e.g. model time step and grid discretisation).
The construction of high-resolution digital terrain models (DTMs) using novel geomatics technologies such as terrestrial laser scanning and low-cost, structure-from-motion photogrammetry have the potential to help constrain these uncertainties (e.g. Westoby et al., 2012Westoby et al., , 2014b. For example, these data can be used to extract accurate models of the crosssectional geometry of a moraine dam and the bathymetry of a (drained) lake basin. However, the accurate quantification of the material characteristics of the dam structure is more difficult to sample effectively and requires logistically challenging fieldwork (e.g. Hanson and Cook, 2004;Osti et al., 2011;Worni et al., 2012). Furthermore, the heterogeneity of moraine sediments is so high that field observations inevitably undersample this variability, so that even the most detailed field measurements must result in significant spatial averaging The predictive performance of linked hydrodynamic models is strongly influenced by model dimensionality (e.g. Alho and Aaltonen, 2008;Bohorquez and Darby, 2008), grid discretisation strategies (e.g. Sanders, 2007;Huggel et al., 2004), DEM quality or the frequency of cross-section data (e.g. Castellarin et al., 2009), and the parameterisation of inchannel and floodplain roughness coefficients (Wohl, 1998;Hall et al., 2005;Pappenberger et al., 2005Pappenberger et al., , 2006. Recent studies have undertaken various forms of sensitivity analyses and uncertainty estimation to quantify the effects of uncertainty on numerical dam-breach model output (e.g. IM-PACT, 2004;Xin et al., 2008;Dewals et al., 2011;Zhong et al., 2011;Worni et al., 2012). The uncertainties surrounding model parameterisation are likely to translate into varying characteristics of the breachoutflow hydrograph, specifically peak discharge, time to peak, and the hydrograph form. When these data are then used subsequently as the upstream boundary condition for hydrodynamic modelling, their effects can be manifest in significantly different simulations of downstream flood propagation. This can give rise to strong variations in the celerity of the flood wave, time-varying flow depths and velocities, and variations in inundation extent. Importantly, it is these data that form the basis for flood hazard assessment. Consequently, it is essential not just to better quantify the predictive uncertainty associated with the component models but also to extend this to consider how uncertainty propagates through the model chain and ultimately might influence the development of effective flood mitigation and management strategies (Beven et al., 2014).

Equifinality in numerical modelling
One consequence of the parametric and structural complexity of many environmental models is the possibility to obtain similar or even identical model outputs through different combinations of input parameters, initial conditions, and model structures (Beven and Binley, 1992;Beven and Freer, 2001;Beven, 2005). Beven and Binley (1992) termed this behaviour model "equifinality". The concept is linked to the more generic form of landscape equifinality used in geomorphology, in which interactions between different processes can give rise to similar landscapes or landform assemblages through differing but equally plausible genetic mechanisms (e.g. Nicholas and Quine, 2010;Stokes et al., 2011). Both forms of equifinality have their origins in systems theory (von Bertalanffy, 1968) and has been identified and quantified in a range of geoscience settings (e.g. Beven and Binley, 1992;Kuczera and Parent, 1998;Romanowicz and Beven, 1998;Beven and Freer, 2001;Blazkova and Beven, 2004;Hunter et al., 2005;Hassan et al., 2008;Vasquez et al., 2009;Vrugt et al., 2009;Franz and Hogue, 2011).
For event-specific reconstructions, calibration of input parameters is often undertaken to identify optimal model parameter sets (e.g. Kidson et al., 2006;Cao and Carling, 2002). Such calibration typically involves a search of the model parameter space to identify optimal model performance with respect to an observed set of observations (e.g. Refsgaard, 1997;Beven and Freer, 2001;Hunter et al., www.earth-surf-dynam.net/3/171/2015/ Earth Surf. Dynam., 3, 171-199, 20152005Westerberg et al., 2011). However, a problem emerges when the observed calibration data contain insufficient information to uniquely identify a single optimal parameter set, and may highlight the presence of multiple optima within the parameter space that may correspond to different mechanistic solutions to achieve the same end results obtained in the calibration data. This "equifinality" is particularly problematic in the case of complex, spatially distributed models where the number of model parameters vastly outweighs the information content of calibration data, which are often spatially averaged model responses, such as flood wave celerity of peak discharge. Multiple strategies can be adopted to minimise these effects, and range from reducing model complexity to increasing the variety of calibration criteria (for example involving multiple data streams). However, where an imbalance in between model parameterisation and calibration is persistent, Beven and Binley (1992) argued that any combination of model input that reproduces the observed outcome, within acceptable limits, must be considered equally likely as a simulator of the system under investigation (Beven and Binley, 1992). This principle has been used to advocate a quantitative method to assess model performance, and present results in a probabilistic framework, based on a method termed GLUE (Beven and Binley, 1992;Beven, 2005). In this article, we adopt this numerical approach to quantify the predictive uncertainty associated with a physically based dambreach model for the reconstruction of a historical outburst flood in the Nepalese Himalaya. In light of the considerable uncertainty that surrounds dam-breach model parameterisation, the application of the GLUE method enables the quantification of the influence which this uncertainty exerts over derivation dam-breach outflow hydrographs. Our approach to probabilistic dam-breach modelling is described below.

Numerical dam-breach modelling
Simulations of incipient breach development were based on HR BREACH, a physically based, numerical dam-breach model. This simulation tool predicts the progressive growth of a dam breach initiated by either the overtopping or piping of non-cohesive and cohesive embankment materials (Morris et al., 2008;Westoby, 2013;Westoby et al., 2014b). The model employs physically based hydraulic, sediment-erosion and discrete embankment stability modules to calculate the evolving breach geometry and associated drainage outflow hydrograph (Morris et al., 2008). This approach offers a significant advance over simplified, empirically derived breaching models or analytical methods that fit the pattern of breach expansion to user-defined final breach dimensions and formation time.
Breach enlargement is simulated through the interaction of two mechanisms: (i) continuous hydraulic erosion based on either equilibrium sediment-transport equations or erosiondepth equations and (ii) discrete mass failures as a consequence of side-slope instability (Mohamed et al., 2002). In this study, erosional processes were modelled using the erosion-depth equation for non-cohesive embankments, after Chen and Anderson (1986): 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 dimensionless coefficients based on the sediment properties. Bending (tensional) failure is represented by a moment of force, M o : where W 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); and e 1 and e 2 are hydrostatic pressure-force eccentricities (m) (Hassan and Morris, 2012). The occurrence of shear-type failure is evaluated 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 ( • ). Subaerial flows across and through the evolving breach are represented using a steady-state one-dimensional flow model (Mohamed et al., 2002;Morris et al., 2008). This model combines a weir discharge equation and a simplified version of the Saint-Venant equations (Chanson, 2011) to simulate breach flow and flow descending the distal face of the dam, respectively. A variable weir discharge coefficient (C d ) is used to reflect the changing profile of the dam crest and is incorporated into calculations of breach discharge (Q b ) as follows (Hassan and Morris, 2012): where L and H are the length and height of the dam crest, respectively. While a significant simplification, this approach accounts for some of the complexity that arises during breach development. Primary model output is time-stepspecific breach discharge, which reflects water-only breach flow.
Earth Surf. Dynam., 3, 171-199, 2015 www.earth-surf-dynam.net/3/171/2015/ Application of the model requires definition of the following key boundary and initial conditions: (i) geometric descriptors of the dam, specifically proximal and distal dam face slope, dam length, crest width, spillway dimensions (width and height), and downstream valley slope; (ii) lake or reservoir hypsometry based on either stage-area or stagevolume curves. In addition, parameter range estimates and their likelihood distributions (either uniform, linear, or triangular), and computational settings (model time step, total run time, and minimum wetted depth threshold) must be defined. The version of the model used in this study incorporates a Monte Carlo parameter sampling routine to automatically search the model parameter set and generate multiple realisations of the potential drainage hydrographs and predicted breach geometry at specific time steps. The input scenarios modelled are described in Sect. 5.3.

Study site
The Dig Tsho moraine-dammed lake complex is located at the head of the Langmoche Valley in the western sector of Sagarmatha (Mt. Everest) National Park in the Khumbu Himal, Nepal ( Fig. 1; 27 • 52 23.89 N, 86 • 35 14.91 E). A relict moraine-dammed lake, oriented west-east, is impounded by a latero-terminal moraine complex to the north and east. The basin is bounded by a near-vertical bedrock face to the south. The parent Langmoche Glacier has receded approximately 2 km since its maximum Neoglacial terminal position and is now confined to the upper north-east face of Tangi Ragi Tau (6940 m).
The moraine possesses steep (25-30 • ) ice-distal and iceproximal faces along the northern lateral moraine. In contrast, ice-proximal faces of the terminal moraine are generally shallower, but they are more topographically complex than their distal counterparts. The moraine is composed of sandy boulder gravel, and includes rare large (> 10 m diameter) boulders. A 200 m wide, 400 m long, 60 m deep breach dissects the northern sector of the terminal moraine ( Fig. 1; Vuichard and Zimmerman, 1987;Westoby et al., 2012). Breach dimensions conform to those for other documented Nepalese moraine-dam failures (Mool et al., 2001).
On 4 August 1985 an ice avalanche from the receding Langmoche Glacier traversed the steep (∼ 30 • ) avalanche and debris cone into the proglacial lake and triggered a major displacement wave. This initial wave is believed to have overtopped the moraine dam and initiated its failure by scour of the moraine crest. Although no directly documented observations are available, anecdotal accounts suggest the breach formation occurred within 4-6 h (Vuichard and Zimmerman, 1987). The geomorphological effects of the Dig Tsho GLOF have been well documented elsewhere (Vuichard andZimmerman, 1986, 1987;Cenderelli andWohl, 2001, 2003). Downstream, the flood created a major hazard, destroying a newly completed hydroelectric power installation at Thamo (∼ 11.5 km from source) and all trekking bridges for a distance of > 30 km. Up to five people were reported to have been killed (Galay, 1985).
Existing observational data for the site include first-hand descriptions of the condition of the glacial lake prior to, and immediately following, moraine dam failure in 1985, as well as a description of the pre-GLOF nature of the valley-floor sedimentology immediately downstream (Vuichard and Zimmerman, 1987). From these observations, an estimate of the total mass of material removed from the breach was calculated as 9 × 10 5 m 3 , and that of the drained volume of the moraine basin as 5 × 10 6 m 3 . The vast majority of the dam material eroded during breach development was deposited within 2 km of the moraine. More recent reconstructions of the moraine and lake basin involving high-precision 3-D photogrammetric reconstructions by Westoby et al. (2012) have challenged these initial estimates and suggest the volume of the contemporary breach (surveyed in 2010) to be approximately 5.8 × 10 5 m 3 , and the volume of water released by the 1985 GLOF to be ∼ 7.3 × 10 6 m 3 . Sedimentological (D 50 ) sampling of matrix material exposed in the breach sidewalls was undertaken during field investigation, and was used to aid dam-breach model parameterisation.

Method
In order to quantify the predictive uncertainty associated with the coupled modelling system used here, a sequential workflow was developed that involved the following steps: i. Topographic data were derived from digital terrain modelling using terrestrial structure-from-motion (SfM) photogrammetry to reconstruct pre-existing moraine and floodplain topography and extract glacial lake bathymetry and dam geometry data.
ii. A priori parameter values for the numerical dam-breach model were estimated using material properties obtained from field investigation and from the published literature.
iii. A series of potential breach initiation scenarios were hypothesised to account for uncertainty in the principal driving factors behind failure (overtopping and mechanical failure).
iv. For each scenario, the multiple model outputs were weighted using a multi-criterion likelihood statistic based on a comparison of the modelled and observed final breach geometries. This provided a means to weight the ensemble of predicted discharges at each time step and derive probabilistic hydrograph forecasts for each scenario.
v. Two-dimensional hydrodynamic modelling was used to simulate the downstream propagation of the probabilistic outflow hydrographs, using upper and lower 95 % confidence intervals to describe the range of potential flows. These simulations provide a range of estimates for the potential inundation extent, flow velocity, and flow depth data for use as input to probabilistic GLOF hazard mapping.
The following sections describe each method in detail.

Topographic data acquisition
Three topographic models are required for this analysis: (i) a reconstruction of the pre-failure moraine dam and lake basin; (ii) detailed models of the breach topography to determine the position, width, length, and slope of the breach and reconstruct the drainage volume; and (iii) the downstream floodplain topography. These models were derived using terrestrial photogrammetry, based on novel structure-from-motion and multi-view stereo (SfM-MVS). A full description of the image processing is beyond the scope of this paper and the reader is referred to more comprehensive accounts elsewhere area-stage relationships, extracted from a fine-resolution, structure-from-motion-derived digital terrain model of the Dig Tsho moraine-dammed lake basin (Westoby et al., 2012). (James and Robson, 2012;Westoby et al., 2012;Fonstad et al., 2013;Javernick et al., 2014). Briefly, SfM-MVS involves the reconstruction of 3-D point cloud data from highly redundant photographic data sets that have a high degree of image overlap. Unlike traditional photogrammetry, SfM uses feature tracking between images and bundle adjustment methods to reconstruct the camera alignment and reconstruct the 3-D scene geometry (Lowe, 2004;Snavely, 2008;Snavely et al., 2008).
In this study, we use published high-resolution reconstructions of the pre-and post-flood moraine topography derived by Westoby et al. (2012). These terrain models were used to quantify the moraine geometry and lake basin hypsometry ( Fig. 2) for the establishment of initial conditions for numerical dam-breach modelling.
A topographic model of the downstream floodplain was obtained using similar methods. The analysis was based on two photosets, incorporating 226 and 303 photographs obtained from south-and north-facing transects along the valley flanks. Photographs were acquired with a Panasonic DMC-G10 (12 MP) digital camera. Automatic focusing and exposure settings were enabled during photograph acquisition. The freely available software bundle SFMToolkit3 (Astre, 2010) was used to process the input images using feature extraction (Lowe, 2004) and sparse and dense SfM-MVS reconstruction algorithms. Following manual outlier removal and editing, final dense point clouds numbered 4.8 × 10 6 and 1.3 × 10 6 points for north-and south-facing photosets, respectively. Data transformation was achieved through the identification of clearly identifiable boulders from kite aerial photography and an existing SfM-DTM of the easternmost sector of the floodplain domain, as well as boulders visible in the georeferenced DTM of the moraine dam. Averaged georegistration residual errors for both dense point clouds were 1.37, 0.30,and 0.06 m for x, y, and z dimensions, respectively.
Georeferenced point cloud data were merged and decimated to improve data handling whilst preserving sub-grid statistics using the C++/Python-based Topographic Point Cloud Analysis Toolkit point cloud decimation toolkit Rychkov et al., 2012). A bare-earth DTM Earth Surf. Dynam., 3, 2015 www.earth-surf-dynam.net/3/171/2015/ at 8 m spatial resolution was extracted for hydrodynamic GLOF routing. Sensitivity analyses have previously revealed this particular grid discretisation to produce time-varying inundation extents and wetting-front travel times of comparable (relative) accuracies to those produced using finer grid resolutions, and of a higher accuracy than coarser grids (Westoby et al., 2014a).

Model parameterisation
The first step in the GLUE workflow is to establish the parameters for inclusion and their respective ranges. In the absence of any prior information regarding parameter and range choice, all available input parameters and their entire simulation range should be included. In practice, complete uncertainty regarding parameter and range choice is unlikely, since a combination of initial sensitivity analyses, modelling guidelines, and basic intuition and reasoning can typically be used to assist in constraining their choice (Beven and Binley, 1992). The parameters and their ranges used in this study (including computational settings) are displayed in Table 1. Information regarding initial, or a priori, parameter distributions is also required, reflecting the modeller's prior knowledge of the parameter space (Beven and Binley, 1992).
In the absence of any information pertaining to a priori parameter probability distributions at Dig Tsho, uniform distributions were used for stochastic sampling of all available parameter spaces.

Model perturbations
In many instances, the specific mechanisms of moraine dam failure are poorly documented, often as a result of a lack of direct observation. Post-event glaciological and geomorphological analysis of a moraine, its parent glacier, and their surroundings may provide clues as to the triggering event that causes moraine failure or a stand-alone GLOF event (e.g. Hubbard et al., 2005;Osti et al., 2011). However, a degree of uncertainty may still surround the establishment of the precise trigger(s) and the nature of dam failure.
Our approach to dam-breach modelling comprises the simulation of three modes of breach initiation, namely (i) a control scenario where lake waters have risen gradually to a point where the overtopping discharge is large enough to trigger sustained down-cutting and breach development; (ii) breach initiation by a series of overtopping waves, such as those resulting from the rapid input of a mass of rock or ice, that traverse the lake and overtop the moraine, thereby initiating its failure; and (iii) instantaneous mechanical failure of the dam face. All three modes are documented triggers for morainedam failure (Vuichard and Zimmerman, 1986;Richardson and Reynolds, 2000;Worni et al., 2012). The trigger mechanism for the failure of the Dig Tsho moraine is generally accepted to be repeated overtopping and down-cutting of the Table 1. Parameter ranges and geometric characteristics of the Dig Tsho moraine dam used for input to the HR BREACH numerical dam-breach model. Input parameter ranges were established from a combination of initial experimentation and parameter sensitivity analysis (E), in situ field observation (F) and data from similar sites stated in the literature (L). All dam geometry data were extracted from the SfM-DTMs of the moraine and floodplain (SfM-DTM), with the exception of downstream Manning's n, which reflects a sedimentological valley-floor characterisation of gravels, cobbles, and rare large boulders (after Chow, 1959 (Chow, 1959) moraine by waves generated from an ice avalanche entering the lake (Vuichard and Zimmerman, 1986). The purpose of including two additional breach initiation scenarios was to explore whether equifinal final breach morphologies would be produced from a variety of breach initiation scenarios. Scenario design is described in the following sections.

Control scenario
A control scenario was formulated in which breach formation was initiated through down-cutting of a predefined spillway. Inclusion of an existing spillway conforms to pre-GLOF observations of the Dig Tsho moraine by Vuichard and Zimmerman (1986), and we note that many extant morainedammed lakes are drained in this manner (e.g. Hambrey et al., 2008). This down-cutting is a result of flow produced by the pressure head associated with our specified initial lake level (which mirrored the reconstructed dam crest elevation of 4356 m a.s.l.). The modelled spillway measured 0.5 m wide and 0.5 deep and extended from the upstream end of the moraine crest to the dam toe. We note that, in the absence of detailed observations of the spillway prior to dam failure (Vuichard and Zimmerman, 1987), these dimensions are hypothetical and do not necessarily replicate the precise spillway conditions prior to dam failure in 1985.
In addition to a control scenario (DT_control), two styles of system perturbation scenario were introduced to the dambreach models to explore and quantify the impact of systemscale perturbations on model output. These perturbations were (i) the introduction of overtopping waves of varying magnitude and (ii) the instantaneous removal of material from the downstream face of the dam immediately prior to breach development

Overtopping waves
The failure of Dig Tsho has been attributed to the overtopping of the terminal moraine by waves produced by an ice avalanche from the receding Langmoche Glacier (Vuichard and Zimmerman, 1987). Presently, most numerical dambreach models, including HR BREACH, are unable to explicitly simulate the dynamic effects of avalanche-lake interactions. This necessitated an inventive yet relatively simple approach to reconstructing overtopping behaviour in HR BREACH. Instead of simulating the passage of a series of gradually attenuating displacement or seiche waves and dynamic interaction with the dam structure, a solution was devised which involved the rapid increase and subsequent decrease of the lake water surface elevation. These artificial water level variations prompted short-lived overtopping of the dam structure at predefined intervals. Temporary increases in lake level were achieved through systematic variation of reservoir inflow to introduce maximum overtopping discharges of 4659, 3171, and 1809 m 3 s −1 , representing initial volumes of the triggering avalanche of 2.0 × 10 5 , 1.5 × 10 5 , and 1.0 × 10 5 m 3 , respectively. This overtopping discharge range reflects the estimated volumetric range provided by Vuichard and Zimmerman (1987). Each initial wave is followed by successive waves of exponentially decaying magnitude. Overtopping wave scenarios were named DT_overtop_min, DT_overtop_mid, and DT_overtop_max, respectively. Richardson and Reynolds (2000) suggested that the initial failure of a moraine dam may be "explosive" in nature. This explosive force is reflected by the mass of individual transported clasts, which can exceed > 100 t (Richardson and Reynolds, 2000). Osti et al. (2011) documented the destabilisation and landsliding of morainic material from the Tam Pokhari moraine dam, Nepal. Seismic activity and excessive rainfall caused oversaturation and subsequent partial failure of a section of the proximal dam structure, which contributed to its eventual failure and the generation of a GLOF.

Instantaneous mass removal
HR BREACH is currently unable to simulate "explosive" or rapid, large-scale rotational mass failures. In order to simulate the instantaneous removal of material from the dam structure, three mass-removal scenarios were developed. HR BREACH requires that the user specify an initial breach spillway in the dam crest and downstream face of the dam. We use these spillway dimensions to represent the mass of morainic material "instantaneously" removed from the crest and distal face of the dam. However, this is far from an ideal numerical realisation of the precise mechanics of catastrophic, near-instantaneous dam failure, and represents a highly simplified approximation.
The default spillway dimensions used in the control and overtopping experiments were 0.5 m wide and 0.5 m deep. Perturbed spillways cross-sectional dimensions were 1, 3, and 5 m 2 , representing total removal volumes of 159, 1435, and 3987 m 3 of material and named DT_instant_1, DT_instant_3, and DT_instant_5, respectively,

Stochastic parameter space sampling
Although a range of stochastic sampling methods are available for use in the GLUE workflow (e.g. Metropolis, 1987;Iman and Helton, 2006), perhaps the most widely used approach is Monte Carlo analysis (e.g. Beven and Binley, 1992;Kuczera and Parent, 1998;Aronica et al., 1998Aronica et al., , 2002Blazkova and Beven, 2004). With the Monte Carlo method, values from individual parameter spaces are sampled in a truly random manner, thereby eliminating any subjectivity which might be introduced at this stage. The method is fully capable of accounting for predefined probability distribution data, making it a simple yet effective tool for rapidly generating the random parameter ensembles required for model input. However, with a poorly defined prior distribution, or a small number of model simulations, point clustering in specific regions of the parameter space may occur. Clustering can be easily avoided by undertaking a straightforward investigation into patterns of histogram convergence of an output variable (or variables), whereby the minimum number of simulations required for the production of an acceptable level of convergence is established. Equally, alternative stochastic sampling methods, including Latin hypercube sampling, may be employed (e.g. Hall et al., 2005;Iman and Helton, 2006). When faced with multi-dimensional problems, the Latin hypercube method can be used to partition the probability distribution of each input parameter, before proceeding to sample from each partition. This method thereby avoids the clustering that is associated with an insufficient number of Monte Carlo samples and ensures that the final parameter ensembles are representative of the full sampling space of each parameter.
In this study, sustained histogram convergence became minimal after the execution of 1000 model runs. This number of simulations was therefore deemed acceptable for stochastic sampling. However, we note that this number of simulations might not be directly applicable to other sites, and we highly recommend that modellers undertake their own pre- liminary sensitivity analyses before establishing a sufficient number of simulations to undertake.

Likelihood measures and functions
Model evaluation is achieved through quantification of how well a parameter ensemble performs at reproducing a series of observable system-state variables, or "likelihood measures". Parameter ensembles that are unable to do so are deemed to be non-behavioural and are assigned a likelihood score of zero. In contrast, ensembles that reproduce these variables within acceptable limits are deemed to be behavioural and accepted for further analysis. It is not uncommon for all ensembles to be rejected (e.g. Parkin et al., 1996;Freer et al., 2002), thereby suggesting that it is the model structure that is incapable of reproducing the observed data, instead of individual parameter combinations (Beven and Binley, 1992). Such a situation may be overcome by widening the limits of acceptability, but at the potential cost of decreasing confidence in any newly accepted ensembles (Beven and Binley, 1992). Behavioural ensembles were assigned positive likelihood values in the range 0-1, where 1 represents an ensemble that is capable of perfectly replicating the observed data. Likelihood functions are specific to each likelihood measure used for model evaluation. Three likelihood measures were used to evaluate model performance: (i) final upstream breach depth (LH1); (ii) the residual sum of squared errors of the final longitudinal elevation profile of the breach (LH2); and (iii) the location of the critical flow constriction (LH3, Fig. 3). These morphological variables are directly quantifiable by comparing final modelled breach geometry with that extracted from the SfM-DTM of the breached moraine dam. Dam breaching is a fully three-dimensional problem that involves progressive backwasting, mass slumping, and downcutting of morainic material by escaping lake waters. The likelihood measures described above are two-dimensional approximations of the breach geometry, and were deemed appropriate descriptors of the breaching process since in combination they describe the breach end states in the vertical and both horizontal dimensions and relate to the modelled and observed lateral, longitudinal, and vertical expansion of the breach.
An observed final upstream breach depth of 16 ± 2 m was used as the first likelihood measure. We note that this value essentially reflects the maximum depth of lake lowering and, by definition, the volume of released lake water. This depth, or lowering, was quantified using an assumed dam freeboard of zero prior to breach initiation, and differs from the previously described maximum breach depth (∼ 60 m) owing to the slope gradient of the breach channel. The error range was designed to account for observed SfM model georegistration errors, which may have resulted in elevation of the lake exit being under-or overestimated. This range also accounts for any post-GLOF lake lowering or seasonal variations in lake level, which remain unquantified. A triangular www.earth-surf-dynam.net/3/171/2015/ Earth Surf. Dynam., 3, 171-199, 2015 likelihood function was used, with an observed breach depth of 16 m and upper and lower limits defined as 14 and 18 m, respectively. The second likelihood measure is a direct comparison between observed (post-GLOF) and modelled elevation profiles of the breach thalweg. This measure represents a distributed method of quantifying the performance of HR BREACH in producing post-GLOF thalweg elevation profiles that replicate that observed in the field. Thalweg elevation data were directly extracted from the SfM-DTM ( Fig. 12a; Westoby et al., 2012). The residual sum of squares (RSS) method was used to quantify the deviation between the observed and modelled data: where Z obs. is the observed elevation (m) of the breach thalweg and Z pred. is the elevation of the modelled thalweg. For the attribution of a scaled likelihood value, an RSS of 0 corresponded to a likelihood score of 1, whilst the lowest RSS for all scenarios (169 221, dimensionless) was used to represent the lower, non-behavioural limit (i.e. a likelihood of zero). A simple linear likelihood function was applied between these upper and lower values. Choice of an observed value for the location of the critical flow constriction was complicated by the asymmetry of the observed breach planform, whereby flow constrictions on either side of the breach are offset by a distance of approximately 40 m. This asymmetry is most likely a function of complex flow hydraulics and patterns of erosion of the moraine during development and expansion of the breach. Specifically, the wide grain-size distribution of the moraine, which comprises material ranging from silts and sands to boulders with intermediate diameters greater than 5 m, combined with its unconsolidated nature causes the breach enlargement process to differ markedly from a systematic, largely uniform style of expansion typically modelled by numerical breach models. Specifically, side-wall detachment and emplacement of large boulders in the breach thalweg may serve to impede or divert breach flow, such that the breach planform adjusts in response to locally altered flow directions and magnitudes. This behaviour may be one explanation for the observed breach asymmetry. Whilst our numerical model accounts for undercutting and mechanical failure of the breach walls, it is unable to resolve flow obstructions caused by individual large clasts.
Use of a trapezoidal likelihood function would have been possible, whereby any modelled values that fell within the observed constriction "offset zone" would be assigned a likelihood of 1. However, such an approach was deemed inappropriate, because it would render a significant proportion of parameter ensembles as absolutely behavioural. Instead, a triangular likelihood function was used. The mid-point of the observed offset was used as the central, observed value, and a range of ±40 m was applied to encompass observed asymmetry in flow constriction location.

Likelihood updating
Where multiple likelihood measures are used, it is necessary to arrive at a final, global likelihood value for each behavioural parameter ensemble. Bayesian updating represents a statistically robust method for combining multiple likelihood values. It is able to account for the influence of prior likelihood values on the generation of updated values as more data become available. In the context of this study, the initial prior likelihood value relates to final breach depth. Likelihood values from the second (breach-elevation profile) and third (location of the critical flow constriction) are subsequently combined with this initial likelihood through implementation of Bayesian updating, which is summarised as (modified from Lamb et al., 1998) where L i |Y 1,2 is the conditional posterior likelihood for a parameter ensemble, i , conditioned on two sets of observations, Y 1 and Y 2 ; L ( i |Y 1 ) is the likelihood of i conditioned on an initial set of observations (Y 1 ); L ( i |Y 2 ) becomes the likelihood conditioned on a second, or additional, set of observations (Y 2 ); and C is a conditional operator, or scaling constant, to ensure the cumulative posterior likelihood equals unity (Lamb et al., 1998;Hunter et al., 2005). The "i" superscript refers to the inverse of the relevant conditional likelihood values.

Cumulative distribution function generation
Final likelihood values associated with each behavioural parameter ensemble reflect the confidence of the modeller in the ability of each ensemble to reproduce an observed data set. Considering the cumulative distribution of these global likelihood values as a probabilistic function facilitates an assessment of the degree of uncertainty associated with the behavioural predictions (Beven and Binley, 1992). These data are referred to as cumulative distribution functions (CDFs). Measure-specific likelihood values for each behavioural ensemble were re-scaled to sum to unity. The final measure can be treated as a surrogate for true probability, but cannot be used for subsequent statistical inference . Weighted and re-scaled ensembles are ranked and plotted as a CDF curve (Fig. 4), from which cumulative prediction limit data can be extracted (Beven and Binley, 1992). The generation of weighted CDFs is unique to the GLUE approach and represents a multivariate, additive method that accounts for ensemble performance.

ISIS 2-D
The shallow-water flow model ISIS 2-D (Halcrow, 2012) was used to simulate GLOF propagation. The model includes alternating direction implicit (ADI) and MacCormack total variation diminishing (TVD) two-dimensional solvers for hydrodynamic simulation. The ADI scheme solves the shallow-water equations (SWEs) over a regular grid of square cells. Water depth is calculated at cell centres, and flow discharges at cell boundaries. The SWEs are solved by subdividing the computation into x and y dimensions. Wa-ter depth (m) and discharge (q x , per unit width) are solved in the x dimension in the first half of the time step. Water depth and discharge (q y ) in the y dimension are solved in the latter half of the time step. Other variables are represented explicitly for each stage. Since water depths are calculated at the cell centre, depths at cell boundaries are interpolated. Different interpolation methods are used, depending on water depth (Liang et al., 2006). Treatment of the friction term is also depth-dependent, such that, below a user-defined threshold, a semi-implicit scheme is used to improve model stability with decreasing flow depth (Halcrow, 2012). The ADI scheme provides accurate solutions for flows where spatial variations are smooth. Sudden changes in water elevation and flow velocity may give rise to numerical oscillations, making it largely unsuitable for the simulation of transcritical flow regimes (Liang et al., 2006). In contrast, the MacCormack-TVD scheme uses predictor and corrector steps to compute depth and discharge for successive time steps. A TVD term, Var(h), is added at the corrector step to suppress numerical oscillations near sharp gradients, making it highly suited for the simulations of rapidly evolving, transcritical, and supercritical flows (Liang et al., 2006), including sudden-onset floods arising from dam breach Liao et al., 2007). Var(h) is calculated as A far smaller time step is required to achieve numerical stability, resulting in extended model run times (Halcrow, 2012).

Hydrodynamic solver comparison
A comparison was carried out to assess which of the twodimensional solvers would be more appropriate for GLOF simulation. The comparison comprised the simulation of a single dam-breach hydrograph across a reconstructed digital terrain model of the Langmoche Khola (Fig. 5). A 4 m 2 grid discretisation was used, in conjunction with a 0.1 and 0.04 s time step for the ADI and TVD schemes, respectively (following guidance in the model documentation). A global Manning's n roughness coefficient of 0.05 was used and reflects a channel bed and margins composed predominantly of pebbles, cobbles, and large boulders, which is characteristic of floodplains in alpine settings, and also reflects our field-based sedimentological observations. The spatial resolution of the valley-floor grid was finer than that used for subsequent probabilistic GLOF mapping, and was intended to represent the finest amount of topographic complexity on the floodplain. An ASTER GDEM-derived digital elevation model was appended to the downstream boundary of the SfM-derived DTM of the Langmoche Khola to allow water to exit the domain of interest unimpeded and eliminate any artificial upstream tailwater effects. Routing of a breach hydrograph of 5 h duration took approximately 2.5 and 0.5 h of simulation time for the TVD and ADI solvers, respectively. The computational burden of the ISIS 2-D TVD solver far exceeds that of the ADI scheme for identical model setups, owing to the finer temporal resolution required by the TVD solver. Depth-based inundation maps for the results of each solver were created in ISIS Mapper ® and exported for display in ArcGIS (Fig. 6). A difference image of inundation was also created. Floodwaters follow the channel thalweg for both solvers during the 1 h. Thereafter, increasing flow stage results in the inundation of a wide reach between 1.1 and 1.7 km. Total inundation of this particular reach is achieved by 01:15 h for the ADI solver, and ∼ 02:00 h for the TVD solver. In the early stages of the GLOF (0-45 min), the travel distance of the flood wave front is slightly greater for the TVD solver. After 1 h, areas of floodplain inundated exclusively by the TVD code are confined to a zone immediately south of the moraine breach. This difference becomes increasingly evident following the onset of floodwater recession (∼ 3 h onwards).
Inundated area is greater for the ADI solver for all time steps (Fig. 6). Maximum difference in inundation extent (122 816 m 2 ) coincides with the upstream flood peak at 2.5 h. Significant dynamic differences also arise between the two solvers. The most prominent of these is the onset of severe oscillatory behaviour from an early stage (∼ 0.5 h) in the ADI data (Fig. 6). These oscillatory behaviours, which are manifested as discrete, arcuate "waves" of excessive flow depth (often > 10 m), are oriented perpendicular to flow, span the entire width of the main channel in places, and occur at regularly spaced intervals. These oscillatory features in fact signify dramatic numerical oscillations in the ADI code and may be the result of numerical "shocks" being aligned with the topographical grid (e.g. Meselhe and Holly, 1997;Venutelli et al., 2002). In contrast, the TVD data reproduce flow channelisation far more clearly, apparently free of any significant oscillatory behaviour.
The results support the use of the two-dimensional TVD solver for GLOF simulation. The lack of any significant numerical instability, otherwise prevalent in the ADI results, is the predominant advantage of the TVD solver. Although processing times are considerably lengthier for the TVD solver, these were deemed acceptable. The solver used in this study simulates only clear-water flows, with no consideration of sediment entrainment, transfer, and depositional dynamics, including any impact on flow rheology. By default, our breach model does not output time-step-specific sediment outflow discharges. It is possible to undertake a crude calculation of time-varying sediment production through the interpolation and differencing of successive breach crosssectional geometries (see Fig. 7 in Westoby et al., 2014a). These data would in theory be suitable for input to multiphase hydrodynamic modelling.  Figure 6. Comparison of GLOF inundation in the upper reaches of the Langmoche Khola using the ISIS 2-D alternating implicit direction (ADI) and total variation diminishing (TVD) solvers at selected time steps. Also shown is a difference image of inundation, where black and grey shading corresponds to areas inundated exclusively by the TVD and ADI solvers, respectively. See Fig. 1a for location. Inset tick marks spaced at 200 m intervals.
In addition, the DTM that was used to represent the floodplain domain immediately downstream of the moraine dam reflects post-GLOF valley-floor topography. As such, derived maps of inundation and flow depth should not be taken as indicative of the passage of the 1985 GLOF.

Results
The following sections present the results of our probabilistic, coupled dam-breach-GLOF simulation experiments. The performance of the dam-breach model at reproducing observed breach geometry is first evaluated (Sect. 6.1.1) before attention turns to issues associated with the extraction of useful, probabilistic breach hydrograph data for use as input to GLOF modelling (Sects. 6.1.2 and 6.1.3). A variety of approaches for translating the impact of dam-breach model parameter uncertainty into probabilistic maps of GLOF inundation and hazard are presented in Sect. 6.2.

Model evaluation
Simulations that were deemed non-behavioural were assigned a likelihood value of 0 and not considered for further analysis ( Table 2) eters except Manning's n (Fig. 7), whereby increasing the roughness coefficient is associated with decreasing peak discharges. However, this inverse relationship is evident only for the first and second likelihood measures. Near-identical relationships were identified for all additional, perturbationfocused scenarios. Behavioural parameter ensembles possessed a low range of peak discharge (Q p ) values relative to the entire range of unconditioned data; maximum non-behavioural Q p values exceeded 18 000 m 3 s −1 for most scenarios, whereas maximum behavioural Q p was found typically to be approximately ∼ 2000 m 3 s −1 (Fig. 10) and is in line with previous empirical estimates for the Dig Tsho failure (Cenderelli and Wohl, 2001), but smaller than other documented Himalayan GLOFs (e.g. Osti and Egashira, 2009). Whilst the number of simulations retained for the overtopping scenarios was broadly similar, the rather more ex-Earth Surf. Dynam., 3, 171-199, 2015 www.earth-surf-dynam.net/3/171/2015/ tended range possessed by the instantaneous mass-failure scenarios reflects an inverse correlation between the initial volume of material removed and number of simulations retained as behavioural. Full hydrographs were obtained for each of the retained simulations. The range of behavioural peak discharges of similar magnitude to previous estimates provided for the Dig Tsho GLOF (Vuichard and Zimmerman, 1987;Cenderelli andWohl, 2001, 2003), and of equivalent or lower magnitude to palaeo-GLOFs reported from other regions (e.g. Clague and Evans, 2000;Huggel et al., 2004;Kershaw et al., 2005;Worni et al., 2012). The maximum and minimum range of behavioural Q p values for all scenarios are strikingly similar (Table 3). These similarities imply that scenarios with considerably different modes of breach initiation are capable of producing a broadly similar range of behavioural peak discharges. This finding suggests that additional factors, including the constraints on breach development imposed by the initial dam geometry and the various parameter sampling ranges, exert an overriding influence on this aspect of breach outflow.
Maximum and minimum behavioural likelihood scores after the data had been conditioned using the RMSE of modelled elevation profile of the breach thalweg varied from 0.97 to 0.014. Within this range, the distribution of likelihood scores was comparable for the control scenario and all overtopping scenarios (0.970, 0.969, 0.969, and 0.970 for DT_control, min, mid, and max overtopping scenarios, respectively). Whilst the maximum likelihood score for DT_instant_1 was almost equally as high (0.821), this value decreases with increasing volume of mass removed (0.819 and 0.744 for DT_instant_3 and DT_instant_5). All of these scenarios possessed minimum likelihood scores that were appreciably lower than the control and overtopping scenarios. Within these ranges, likelihoods were distributed relatively evenly. Accordingly, the instantaneous mass-removal scenarios, particularly DT_instant_3 and DT_instant_5 exhibited the poorest performance against this likelihood function.
Behavioural elevation profiles are displayed in Fig. 8. The well-defined break in slope which is identifiable in the observed data at ∼ 550 m is reproduced by the output from HR BREACH. However, HR BREACH almost consistently underestimates the distance along the breach at which this break occurs. The majority of modelled profiles in the behavioural DT_instant_3 and DT_instant_5 simulations are located ∼ 50-100 m further upstream than the equivalent SfM-derived, observed profile (Fig. 8). It is unclear why this systematic underestimation occurs, although an important consideration is that the observed centreline profile, viewed from above, is not linear. Instead, the breach thalweg meanders along its length, having exploited localised weaknesses in the degrading dam structure and breach-flow dynamics as it developed. Consequently, the observed and modelled profiles are not truly comparable, introducing a fundamental source of error into the resulting likelihood scores (and also accounting for the discrete zone of variance between 300 and 370 m on all plots (Fig. 8).
Modelled planforms are broadly similar in form both between and within each scenario (Fig. 9). Observed and modelled planforms gradually taper towards a flow constriction, beyond which the breach width expands to form a bellshaped exit. Flow-constriction location varies considerably, with a substantial number of parameter ensembles deemed non-behavioural following conditioning using this likelihood measure ( Table 2). The majority of non-behavioural simulations located the flow constriction upstream of the behavioural limit (Fig. 9). However, no discernible relationship between input parameters and flow constriction location was identified in the parameter-specific likelihood data. Only seven DT_control parameter ensembles were retained (0.7 % of the original simulation pool), following conditioning on this likelihood measure, and only two (0.2 %) of the DT_instant_5 simulations remained. Further reductions in the number of retained simulations were imposed for all scenarios (Table 2).

Percentile hydrograph extraction
CDF curves were extracted from behavioural, scenariospecific likelihood data (Fig. 4). Using these data, time-stepspecific percentile discharges were extracted and combined to construct probabilistic breach outflow hydrographs for each scenario (Fig. 10). Similarities between the percentile hydrographs for each scenario are striking, particularly between data conditioned on modelled final upstream breach depth (LH1) and modelled upstream breach depth and breach centreline elevation profile data (LH1 + 2). All 95th percentile hydrographs possess steep rising limbs and comparatively shallower falling limbs, and they generally trace the upper boundary of the behavioural hydrograph envelope. Exceptions to this rule are DT_overtop_min, DT_overtop_max, and DT_instant_1, whose hydrographs dip below this upper bound by a noticeable margin. DT_instant_3 and DT_instant_5 possess shorter-duration hydrographs -a direct consequence of the form of the behavioural envelope from which the CDFs were derived. Shortening of the hydrograph duration results in greater concentration of all percentile hydrographs for these scenarios, with the result being reduced variation in percentile-specific Q p and Q vol (Table 4).
Individual overtopping waves are preserved for each percentile in the relevant scenarios (Fig. 10). Median (50th) percentile hydrographs exhibit slightly more variation, both between scenarios and following conditioning using the second likelihood measure. This conditioning step results in a decrease in median percentile Q p for control and overtopping scenarios, but it has a negligible effect on the mass-removal simulations (Fig. 10). For all scenarios, 5th percentile hydrographs generally trace the lower boundary of the behavioural hydrograph envelope, and appear to be largely unaffected by additional conditioning using the elevation profile of the breach centreline.
The most noticeable impact on percentile hydrograph form is caused by the additional conditioning of the data on the final likelihood measure. Mass-removal scenarios appear to be affected to a far lesser degree than the control and overtopping scenarios. The exception is DT_instant_5, where discharges for 50th and 5th percentile hydrographs are increased in the first ∼ 80 min, after which 95th percentile discharges decrease, resulting in substantial concentration of the percentile hydrographs (Fig. 10). Conditioning on this final likelihood measure causes increased Q p for DT_overtop_min and DT_overtop_mid median percentile hydrographs, whilst the hydrograph for DT_overtop_max is perturbed between 140 and 200 min. Notably, DT_control exhibits increased discharges for the 5th and 50th percentiles, coupled with decreased Q p and steepening of the falling limb of the 95th percentile hydrograph.
Crucially, percentile hydrograph form is dictated by the time-step-specific CDF data. In turn, CDF form is deter-Earth Surf. Dynam., 3, 171-199, 2015 www.earth-surf-dynam.net/3/171/2015/ mined by variations in the likelihood of individual behavioural hydrographs, and the cumulative distribution of their associated discharges (for each time step). The vast number of simulations which, following conditioning on flow constriction location, were subsequently deemed to be non-behavioural (Table 1) serves to alter the form of scenario-specific CDFs. This effect is particularly dramatic for DT_control, where the number of behavioural simulations reduces from 76 to 7 following conditioning on final breach depth, breach centreline elevation profile, and flow constriction location (LH1 + 2 + 3). Similarly, CDF data for DT_instant_5 are derived from just two behavioural hydro-graphs, which explains both the smaller hydrograph envelope and the high concentration of the percentile hydrographs (Fig. 10).

Data clustering
Issues of mass conservation arose with the extraction of behavioural, percentile-derived breach hydrographs. Both 5th and 50th percentile hydrographs for all scenarios consistently under-predicted the volume of lake water released during breach development (∼ 7.3 × 10 6 m 3 ), whereas volumetric integration of 95th percentile hydrograph data produced substantially increased values of Q vol (Table 4). The median   Figure 10. Percentile hydrographs derived from behavioural Dig Tsho simulations, for successive likelihood updating steps ("LH1", "LH1 + 2", "LH1 + 2 + 3"). Small dashes: 5th percentile; long dashes: 50th (median) percentile; solid black: 95th percentile. Behavioural hydrograph envelope is displayed in grey. percentile performs best at replicating total flood volume, under-predicting Q vol by between 0.2 and 0.7 × 10 6 m 3 , 0.4 and 0.7 × 10 6 m 3 , and 0.1 and 1.3 × 10 6 m 3 for LH1, LH1 + 2, and LH1 + 2 + 3, respectively (Table 4). From a hydrodynamic modelling perspective, an additional and equally important observation is the form of the percentile hydrographs, which generally do not mirror the form of any of the behavioural hydrographs used as input. When combined, issues of mass conservation and the unrepresentative form of the percentile breach hydrographs render them largely unsuitable for use as an upstream boundary condition for subsequent hydrodynamic modelling.
In an effort to further refine the behavioural simulations and improve the representativeness of the derived percentile hydrographs, all behavioural data were clustered, regardless of their inclusion of any modelled system perturbations (Fig. 11). Data clustering was undertaken in an effort to characterise "styles" of breaching, such as those characterised by low peak discharge magnitude and lengthy time to peak or high peak discharge and short, sharp rise to peak. Clustering used Q p and T p data as the sole input, and was independent of the various breaching scenarios that were modelled. An automated, subtractive clustering algorithm was applied to the unified Q p and time-to-peak (T p ) data. This method identifies natural cluster centroids in raw point data sets, and quantifies the density of points relative to one another. Each point is assumed to represent a potential cluster centroid, and a measure of likelihood is calculated based on the density of surrounding points. Following initial cluster identification, the density function is revised and subsequent cluster centres identified in the same manner until a sufficient number of natural clusters are deemed to have been obtained (Math-Works, 2012). Point-cluster membership was determined by calculation of the minimum Euclidean distance between each data point and cluster centre. The subtractive method eliminates any subjectivity associated with manual cluster identification.
Clusters are broadly defined by T p range. Cluster 1 contains all simulations with T p of approximately 60-130 min, cluster 2 is defined by T p of ∼ 130-170 min, and cluster 3 possesses T p values in the range ∼ 170-270 min. Ranges of Q p overlap substantially between cluster 1 and cluster 2 (∼ 900-2100 and ∼ 700-1800 m 3 s −1 , respectively), and to a lesser degree between cluster 2 and cluster 3. These clusters may be taken as approximately representing a number of "types" of breach hydrograph, namely (relatively) highmagnitude, short-duration (cluster 1), moderate peak magnitude and mid-range duration (cluster 2) and low-magnitude, extended duration (cluster 3) GLOF hydrographs. Cluster membership is not as clear-cut as might be anticipated (Fig. 11). The clustering results appear to imply that pigeonholing different breaching scenarios by hydrograph type, or style, is virtually impossible. However, the exceptions to this rule are the instantaneous mass-removal scenarios, which almost exclusively produce high-magnitude, short-duration hydrographs. This finding would appear to imply that alternative factors are required to explain the similarity in the range of hydrograph forms that are produced by each scenario.
Percentile hydrographs were also extracted from the clustered data. Deviations between observed and modelled median percentile Q vol data are in the range −41 to −6 %, demonstrating a minimal improvement over unclustered and scenario-specific values of modelled Q vol (Table 4). Clus- Table 4. Total percentile hydrograph outflow volumes for individual scenarios and following clustering.

The optimal hydrograph
Deterministic approaches to flood reconstruction require the identification of the optimal model, and its subsequent use for predictive flood-forecasting. To illustrate the variability between scenario-specific optimal hydrograph routing patterns, the optimal hydrographs for DT_control, DT_overtop_max and DT_instant_5 were used as upstream input for simulation in ISIS 2-D. Maps of inundation extent and flow depth (Fig. 12) reveal prominent inter-scenario differences in the spatial extent of inundation as the respective hydrographs and GLOF floodwaters progress downstream. Variations include the initial downstream transmission of the DT_overtop_max overtopping wave, which triggers rapid inundation of the entire reach. However, initially high-flow stages are not maintained; these only rise once again with increasing breach discharge associated with breach expansion. Use of the DT_instant_5 hydrograph produces spatial and temporal patterns of inundation and wetting front travel time similar to that of DT_overtop_max (Fig. 12).

GLUE-based GLOF reconstruction
Probabilistic maps of inundation extent and flow depth were constructed through the retention and evaluation of scenario-specific and likelihood-weighted breach hydrographs. In the example presented herein, we simulated the propagation of 76 individual moraine-breach hydrographs using the ISIS 2-D TVD solver (with an 8 m topographic grid discretisation and a 0.04 s time step), representing the behavioural DT_control parameter ensembles after conditioning on final breach depth. However, the method is equally applicable to the use of several, hundreds, or thousands of individual simulations. For each time step, per-cell CDF curves of flow depth were assembled, from which percentile flow depths were extracted and plotted (Fig. 13). Given the inherent uncertainty surrounding the precise mode of moraine-dam failure and outflow hydrograph form, these data effectively convey the resulting variability in likelihood-weighted predictions of reconstructed inundation extent, whilst preserving time-stepspecific percentile flow depths. Because of the nature of their construction, these data do not relate to a specific event or hydrograph but instead provide an indication as to the potential uncertainty in GLOF inundation extents and flow depths associated with a range of behavioural breach hydrographs.

Probabilistic hazard mapping
The final output of a GLOF hazard assessment comprises the production of maps of flood hazards, conditioned by one or more directly quantifiable flood-intensity indicators (e.g. Aronica et al., 2012). Whilst inundation depth is arguably the most significant flood-intensity indicator for predicting monetary losses associated with individual flood events (Merz and Thieken, 2004;Vorogushyn et al., 2010Vorogushyn et al., , 2011, its combination with flow velocity is regarded as an improved indicator of hazard to human life (Aronica et al., 2012).
A global hazard index proposed by Aronica et al. (2012) was used to construct maps of GLOF hazard (Fig. 14). Taking probabilistic flow depth and velocity data as input, probabilistic GLOF hazard maps were produced for the DT_control scenario. Four hazard classes are defined and shaded for distinction (H 1 to H 4 , in order of increasing hazard). Hazard distribution is characterised by extensive zones of H 1 for the 5th and 50th percentiles in the first 0.5 h, with higher hazard classes becoming increasingly prevalent in the 95th percentile data (Fig. 14). As breach outflow increases and the maximum inundation extent is approached, the zone immediately adjacent to the channel thalweg is classed as "very high" (H 4 ) hazard, and is associated with flow depths  Figure 13. GLUE-based percentile maps of inundation for DT_control. 5th, 50th, and 95th percentiles maps of inundation represent the water depth that would be exceeded with 95, 50, and 5 % probability, respectively. See Fig. 1a for location. Inset tick marks spaced at 200 m intervals.
GLOF, the purpose of including these data is to highlight the impact of breach model parametric uncertainty on one of the final products in a typical GLOF hazard assessment exercise. The representivity of these data is associated with a number of caveats, namely the routing of the flood across a post-GLOF digital elevation model and the application of water-only hydrodynamic modelling. In reality, the flow may have exhibited transient debris-flow characteristics given the availability of sediment for entrainment and transport from the eroding breach and unconsolidated valley floor (Huggel et al., 2004;Westoby et al., 2014a).  Figure 14. Percentile flood-hazard maps, based on a global hazard index forwarded by Aronica et al. (2012). Such data facilitate a probabilistic evaluation of the evolution of GLOF flood hazard. See Fig. 1a (this paper) for location, and Fig. 1b in Aronica et al. (2012) for description of the hazard index.

Key controls on breach development
We have demonstrated that the propagation, or cascading, of the parametric uncertainty and equifinality through the dam-breach and hydrodynamic modelling components of the GLOF model chain is not only possible but may also be of considerable value to flood-risk practitioners. A key contribution of the research is the demonstration that the predictive limits of numerical models, in this instance applied to the reconstruction of historical moraine breaching, can be quantified through the use of a weighted, probabilistic modelling framework (Fig. 15). Our approach illustrates how para- Figure 15. A unifying, "end-to-end" framework for probabilistic GLOF reconstruction incorporating high-resolution photogrammetry and probabilistic, GLUE-based numerical dam-breach and hydrodynamic modelling approaches. metric uncertainties may be propagated through the GLOF model chain when the output from one model is used as the input to another. The approach can therefore be used to isolate the most sensitive components of a predictive systeme.g. Manning's n -and thus indicate the most important areas for future model development or data collection to support informed model parameterisation and testing.
Our results highlight the primary influence of the material roughness of moraine material indicating HR BREACH parameter ensemble, and therefore breach-hydrograph, per-formance. Specifically, behavioural Manning's n coefficients were found to be in the range 0.020-0.029 m −1/3 s (Fig. 7), representing a significant refinement of the a priori parameter range (0.02-0.05 m −1/3 s). In contrast, no significant refinement of the remaining material characteristic parameters was made following model evaluation. However, posterior Manning's n values appear to be unrealistic; the Dig Tsho moraine is composed predominantly of pebble-, cobble-, and boulder-sized material, which is more likely to be associated with larger roughness-coefficient values (Chow, 1959).
If all additional resistance effects, including form-, spill-, and curvature, are fully accounted for in the model solution, then the resistance term will reduce to account for the particle, or skin, resistance alone. This would account for the retention of such "low" Manning values. However, since 2-D fluid flow is modelled in HR BREACH at comparatively coarse spatial resolutions relative to the scale of the particles that moraine is composed of, it is reasonable to expect that the effective values of n will incorporate the aforementioned macroscopic properties. In such a case we would expect values of n to be larger. For example, in a reconstruction of the Chukhung moraine breach, also in the Khumbu Himal, Westoby et al. (2014b) found behavioural roughness values for overtopping-type failures to be in the range 0.042 to 0.049. Our results suggest that a reduction in n is more likely to reflect parameter compensation or interaction effects elsewhere in the model chain. Specifically, a reduction in n will increase shear stress, which will in turn offset other parameters that serve to reduce bed erodibility. Given the complexity of the model parameter space, it is these effects that account for the low values of n, rather than uniquely the identification of appropriately scaled particle-roughness values.
The observation that broadly similar behavioural peak discharges are associated with different modes of breach initiation is a particularly significant one (Table 3). This finding appears to suggest that additional factors, such as the sampling ranges for the various input parameter ranges and model boundary conditions (dam geometry and lake bathymetry), exert an overriding influence on this breach hydrograph characteristic. However, the significance of breach initiation mode is a primary control over downstream GLOF wetting front travel times and inundation extent. This is illustrated in Fig. 12, and reinforces the need for moraine-dam failure modelling exercises to consider, wherever possible, the complete range of potential breach-initiation scenarios at the model design stage.
The numerical physicality of the dam-breach model used in this research represents a notable improvement over empirical or analytical models. However, many of the geometric and material characteristics of the moraine and lake complex remain highly simplified. This simplification is a necessary compromise related to our use of a dam-breach model, the primary intended application of which is the investigation of (relatively) simple artificial earthen or concrete dam constructions. Alternative dam-breach models, including unstructured mesh-based variants (e.g. Worni et al., 2012), have been demonstrated to perform well at reproducing historical moraine-dam failure dynamics. In combination with stochastic parameter sampling functionality, they would represent a powerful combination that could be easily incorporated into the framework we present herein (Fig. 15).

Wider use of probabilistic GLOF modelling
The model chain presented in this article is parametrically complex, and whilst there remains an imbalance between the parametric degrees of freedom in the model setup and the objective descriptors of post-GLOF geometry used to evaluate model performance, this imbalance is one common to many areas of numerical modelling in the Earth surface sciences. The value of the GLUE method in this respect is its suitability for actively exploring and quantifying the extent of uncertainty and equifinality in dam-breach model output that results from poorly constrained conditioning of the model input parameters. Similarly, this problem is common to virtually all dam-breach modelling efforts, not only moraine dam breaching. This research has provided specific tools, in the form of probability-based maps of GLOF inundation and hazard, to communicate the effects of model uncertainty to potential end users in a manner that is both open and objective.
Whilst the extraction of percentile maps of inundation extent and flow depth is not necessarily an entirely new concept (e.g. Brasington, 2007, 2008;Vorogusyhn et al., 2010Vorogusyhn et al., , 2011, the application to GLOF reconstruction presented here is an original and novel one. This approach represents a significant improvement in the effective communication of the likelihood associated with a range of moraine-dam failure scenarios. Significantly, our probability-based flood hazard maps (Fig. 14) can be visualised in a GIS environment, and provide a clear picture of patterns of flood hazard zonation, at varying levels of confidence and at successive model time steps, that would prove useful to disaster risk managers.
The likelihood of multiple GLOFs occurring from an individual moraine-dammed lake complex is low. In most cases, the breached moraine dam can comfortably accommodate relict lake discharges. Therefore, the identification and use of posterior parameter distribution data for predictive GLOF forecasting is of limited utility if these ranges prove to be highly site-specific. The identification of a suite of universal or region-specific material characteristics and their probabilistic distributions would facilitate their use in predictive GLOF simulation efforts. We believe that the Dig Tsho failure is regionally representative of comparable styles of glacial lake systems, and as such our results are of value for extension of the technique to similar glaciated regions, with the caveat that a first-pass assessment of GLOF hazard should be first undertaken to identify "priority" glacial lakes (Huggel et al., 2004). In identifying the importance of constraining specific parameters for dam-breach model parameterisation (e.g. Manning's n), we have highlighted the importance of identifying model-critical data, which will aid the design of future field campaigns.
Probabilistic approaches have clear advantages over the deterministic approaches traditionally used for GLOF reconstruction such as the use of palaeohydraulic techniques for at-a-point or reach-scale peak discharge estimation (e.g. Cenderelli and Wohl, 2003;Kershaw et al., 2005;Bohorquez and Darby, 2008). Probabilistic methods embrace and attempt to convey the influence of uncertainty and equifinality in model input on subsequent output. It might be argued that their value outweighs the additional processing time required for their implementation, which may involve the execution of hundreds or thousands of individual simulations.
In considering the source of uncertainty in the GLOF modelling process, we have focused on its influence over the moraine dam-failure process. However, numerous additional sources of uncertainty are present at various stages in the workflow, such as the reconstruction of lake-basin bathymetry, and merit further investigation (Westoby et al., 2014a). The logistical impracticalities of identifying and addressing all sources of uncertainty in the GLOF model chain are currently a significant hindrance to applied modelling efforts. However, simple sensitivity analyses remain of value to quantify the impacts of individual sources of uncertainty on numerical dam-breach and hydrodynamic simulation, and might be incorporated straightforwardly into our modelling framework (Fig. 15).

Conclusions
This paper has outlined and presented results from a workflow for cascading uncertainty and equifinality through the glacial lake outburst flood (GLOF) model chain using a combination of advanced, physically based numerical dambreach and hydrodynamic models. Dam material roughness is the dominant influence on outflow hydrograph form. Morphological characteristics of a GLOF breach are appropriate measures for assessing the performance of individual simulations, or parameter ensembles, at reproducing observed breach morphology. Breach morphology is reproducible by parameter ensembles associated with differing breach-initiation scenarios, lending support for the adoption of probabilistic, as opposed to deterministic, methods for dam-breach outburst-flood reconstruction. We also demonstrate an effective approach for cascading dam-breach simulation likelihood data through to the construction of probability-based maps of GLOF inundation extent and flow depth, and the subsequent derivation of event-specific maps of flood hazard.