Introduction

On 25 April 2015, Nepal was struck by a magnitude 7.8 earthquake that caused over 8000 casualties (GoN National Planning Commission 2015). Many aftershocks followed, including a 7.3 temblor on 12 May. These earthquakes weakened steep hillslopes (Kargel et al 2015), including those of Myagdi district (Anonymous 2015c). Near the village of Baisari (28.40°N, 83.60°E), a damaged slope failed early on 24 May, causing a landslide, which blocked the Kali Gandaki river. Officials quickly ordered evacuations of low-lying areas in the town of Beni and the informal settlement of Maldhunga (Figure 1) downstream of the landslide, preventing casualties when the landslide dam was overtopped and breached later that day (Marahatta and Parajuli 2015).

FIGURE 1 

Study site and stream thalweg profile. (A) Locator map; (B) western/central Nepal; (C) topography of the model domain, in World Geodetic System (WGS) 1984 Universal Transverse Mercator (UTM) Zone 44N horizontal datum with units of meters; (D) stream thalweg profile of the site with raw and smoothed topographies. (Source: C, GoN Department of Survey 1996)

This event is just one example of the frequent flash floods experienced in Nepal. In addition to landslide dam breaching (Weidinger 2006; Awal 2008a, 2008b), flash floods can be caused by intense rainfall events in steep mountainous areas, rockfall induced by glacial melting (Bhandary et al 2012), or glacial lake outburst floods (Bajracharya et al 2007). Rapid prediction of flood arrival time and depth is essential for issuing timely warnings to downstream population centers, yet it is often compromised by lack of the data that would make it possible to run simulations with the increasingly sophisticated hydrodynamic models available today (Schwanghart et al 2016). Using the breach of the Baisari landslide dam as an example, this study investigated options for predicting flood behavior in a data-scarce environment, including the effect of model complexity—1-dimensional (1-D) versus 2-dimensional (2-D)—and digital-elevation-model preprocessing on the results of hydrodynamic simulations.

Field observations by residents

Local residents used smartphone cameras to record videos of the Kali Gandaki flood at Beni Bridge (28.34°N, 83.56°E) at the time of maximum flow depth and speed (Figure 2). Construction drawings position the bridge's low chord (lowest elevation of the girder under the bridge deck) at 6.5 m above the streambed. Based on this, the video indicates a maximum flow depth of 4 to 5 m. Use of image rectification (Bourgault et al 2013) followed by particle image velocimetry to trace the motion of floating debris indicated a maximum flow speed of 9 to 10 m/s at the water surface.

FIGURE 2 

Frames from video of Beni Bridge taken during the peak of the Kali Gandaki flood, with flow speed evidenced by floating debris (indicated by black squares). The first, second, and third frames were recorded at 1-second intervals, but the time interval between the third and fourth frames was 2 seconds. Small circles and crosses are reference points used in image rectification and may be ignored. (Photo source: Anonymous 2015a)

For comparison with depth-averaged river flood models, we converted surface flow speed to depth-averaged flow speed. In typical open channels, the vertical profile of horizontal flow speed can be represented by a power law profile (Chanson 2004) as shown in Equation 1.

Numerical models

To generate a hydrograph for the landslide dam breach itself, the BREACH dam-break model (Fread 1988a) was used. Using this hydrograph, flood routing downstream was evaluated with 2 models: the 1-D unsteady HEC-RAS (USACE 2010; Brunner 2014) and 2-D Delft-FLOW (Deltares 2011).

BREACH dam-break model

BREACH is a physically based mathematical model of dam breach by overtopping or piping (internal erosion). Overtopping was known to be the cause of the flood in the present case (Marahatta and Parajuli 2015; Wang et al 2016). Following the onset of overtopping, the model uses shear-stress-based scour formulas and geotechnical slope-stability formulas based on the material properties of the dam to calculate erosion of the breach into the dam face as well as breach widening by side slope failure.

BREACH is one of several dam-break models used by the US National Weather Service (Fread 1988b), but it is the only one of these that uses a physical basis to calculate the flood hydrograph. The other models are empirical, using observations of past failures of man-made dams to approximate the rates of breach deepening and widening. Few data exist on breaching of landslide dams (O'Connor and Beebee 2009), but the physically based BREACH model has successfully been applied to such cases (Fan et al 2012). Therefore, BREACH was chosen to generate the hydrograph in this study. Other physically based dam-break models include Awal's (2008a, 2008b); Table 1 lists several empirical dam-break models.

TABLE 1 

Empirical dam breach models and their peak flow hindcast.

The Baisari area is remote, so data are sparse. Landslide geometry was estimated based on Collins and Jibson (2015). The landslide was reported to be approximately 200 m long, 50 m wide, and 30 m high, with a total volume of 300,000 m³, impounding a lake with a volume of 8 million m³. Assuming a triangular cross section over a local riverbed slope (horizontal:vertical) of 50:1 (GoN Department of Survey 1996) allowed estimation of the upstream and downstream face slopes of the landslide (Table 2). Together with Collins and Jibson's (2015) description of the landslide, photographs of the landslide (Petley 2015) allowed gross estimation of geotechnical parameters (Table 2). Reservoir surface area as a function of water level was determined using a topographic contour map (GoN Department of Survey 1996).

TABLE 2 

BREACH landslide dam parameters.^a)

Flood-routing models and digital-elevation-model preprocessing

Routing of flash floods is a topic on which research is actively being carried out (Lighthill and Whitham 1955; Mudd 2006; Abderrezzak et al 2009; George 2011). The flood hydrograph shown in Figure 3 was routed downstream by both 1-D and 2-D models, in order to compare the performance of each type of model. Topography used in each model was derived from a contour map of the area (GoN Department of Survey 1996) digitized to 15-m resolution.

FIGURE 3 

Hydrograph of landslide breach with sensitivity analysis to geotechnical parameters.

In addition to flood routing using these raw topography data, other flood-routing models were run using topography data that had been conditioned with hydrological smoothing (Koenker 2005; Schwanghart and Scherler 2014). Unprocessed digital elevation models frequently contain artificial topographic depressions along a river's longitudinal profile; previous approaches have attempted to remove these by either filling or carving (Lindsay 2015). However, both methods tend to generate long zero-gradient sections followed by steep steps. The smoothing approach applied here used curvature-regularized linear interpolation with additional constraints on downslope gradient. This smooths the longitudinal river profile while ensuring that elevations decrease uniformly downstream, thus avoiding both sinks and macroroughness elements absent in reality.

Based on the smoothed profile, we derived a 1-m-deep, 45-m-wide channel carved into the thalweg. Figure 1D shows the effect of hydrological smoothing on the stream thalweg profile. Results were extracted at the locations of Beni Bridge and Maldhunga Bridge (28.25°N, 83.61°E) for comparison between models, with the particle image velocimetry result for maximum flow depth and speed based on the video recorded by a resident.

Two-dimensional Delft-FLOW: A depth-averaged Delft-FLOW model (Deltares 2011) was set up for the domain shown in Figure 1C, using a curvilinear grid with resolution from 20 to 40 m. The model time step was 0.6 second. Horizontal turbulent eddy viscosity was assumed to be a uniform 1 m2/s (Bricker and Nakayama 2007); sensitivity to this parameter was investigated by running the model with an eddy viscosity of 0 m2/s as well, but the result was indistinguishable from the 1 m2/s result. To correctly resolve the flow of the wetting front over dry land, Delft-FLOW's FLOOD advection solver (Stelling and Duinmeijer 2003) was used. The base-case hydrograph shown in Figure 3 was applied to the model as a point discharge into the stream thalweg at the model's upper boundary, with the boundary itself a solid wall. The lower boundary was inconsequential, because it had no impact on the lowest point of interest (Maldhunga Bridge) during the simulation.

Manning's n was specified per the formulation of Jarrett (1984), in which n = 0.38S^0.38H^−0.16, where S is the stream slope and H is the water depth. For an average bed slope of S = 0.007 (from the topographic data) and flow depth H = 5 m (based on trial-and-error model runs), a uniform n = 0.04 s/m^1/3 was determined and applied throughout the domain. As a sensitivity analysis, another model run was conducted with roughness n = 0.03 s/m^1/3. Likewise, to test sensitivity of the model result to grid spacing, a run was conducted on a fine grid with resolution 3 times that of the base-case (coarse) grid, using a time step of 0.2 second. The smaller time step for the fine model run was necessary to preserve model stability because of the Courant number condition (Morton and Mayers 1994).

The Delft-FLOW model was validated for use in flash flood modeling by comparison with the results of Hervouet and Petitjean (1999) for the 1959 Malpasset dam break, as measurements of the wave travel time and flood elevations were reported. Model topography, initial conditions, and roughness were set up as in Hervouet and Petitjean (1999). Uniform rectangular grids at 20- and 10-m resolution as well as coarse (approximately 15-m resolution in the channel) and fine (approximately 7.5-m resolution in the channel) nonuniform stream-following curvilinear grids were evaluated. The model time step was 0.01 second.

Figure 4A shows the flood wave elevation at a transformer that was known to have been knocked out of service by the flood approximately 21 minutes after the dam break. All of the Delft-FLOW configurations tested showed this time to be the onset of a rapid rise in water level; the coarse (20-m mesh) rectangular grid simulation was the most diffusive, whereas the fine curvilinear grid simulation was the least diffusive. As Figure 4B shows, all model configurations reasonably reproduced maximum measured water depths from the physical model described in Hervouet and Petitjean (1999). Because a curvilinear grid allowed higher resolution and less diffusive modeling than a rectangular grid, given the computational resources available, a curvilinear grid was chosen for modeling the Kali Gandaki flood.

FIGURE 4 

Validation of Delft-FLOW 2-D model for the 1959 Malpasset dam break, based on Hervouet and Petitjean (1999). (A) Flood wave arrival time at transformer; (B) maximum flood depths at measurement points in the laboratory model.

One-dimensional HEC-RAS: For comparison between 1-D and 2-D models, the flood was also routed using HEC-RAS (USACE 2010). A 1-D HEC-RAS model of the river reach shown in Figure 1C was generated with both raw and smoothed/carved topography, but only the latter ran in a stable manner. Transects were defined transverse to the stream thalweg with a maximum spacing of 75 m. The upstream boundary condition was the base-case hydrograph shown in Figure 3 at the landslide location, and the downstream boundary condition (far enough downstream to be inconsequential) was normal depth with a slope of 0.007. Because the HEC-RAS model requires water in all transects in order to run correctly, an initial flow rate of 5 m3/s was specified throughout the model domain. HEC-RAS was run in unsteady-flow mode, set to operate in the mixed-flow regime (allowing flow to transition between subcritical and supercritical). The model time step was 10 seconds, and a uniform roughness of n = 0.04 s/m1/3 was applied, as discussed above.

Results

Flood hydrograph

The BREACH hydrograph of flow exiting the dam site is shown in Figure 3, with a maximum flow rate of 1696 m³/s occurring 2.4 hours after the beginning of overtopping. In addition to physically based models, many empirical models have been developed based on historical dam failures; they give maximum flow rate as a function of dam height and impounded volume and are summarized by Brunner (2014) and Khanal et al (2015). As Table 1 shows, these formulas produce widely varying results, with the BREACH result near the lower end. However, the agreement of the BREACH hydrograph with observations of the flood (discussed below) lends it credence.

Dam geometric parameters were specified based on the observations of Collins and Jibson (2015), which estimated a mean grain size of less than 20 cm. Wang et al (2016) estimated a median grain size of 0.02 mm in the upper 40-cm airfall layer of the deposit. Other geotechnical substrate characteristics were not measured. Figure 3 shows an analysis of the hydrograph's sensitivity to variation in soil quality parameters, with the base case as described in Table 2. The peak flow rate varied slightly (about 10%) within a reasonable range of input values, and the shape of the hydrograph remained relatively unchanged. Sensitivity analysis considered only uncertainty in soil quality, not dam geometry or the impounded stage–area relation, as the observations of Collins and Jibson (2015) and map contours were used for the latter.

The hydrographs shown in Figure 3 are at odds with the report by Wang et al (2016) that the time from the onset of overtopping to maximum outflow was only one-half hour, and that the breach was asymmetric because of the cliff face constraining the left-bank erosion (BREACH cannot account for this constraint). Unfortunately, Wang et al (2016) did not provide any quantitative information to use in hydrograph generation. However, Marahatta and Parajuli (2015) estimated that breaching began about 1 hour after the onset of overtopping, that the maximum breach flow rate was 1291 m³/s, and that this flood flow lasted for approximately 30 minutes. A simple hydrograph with a linear rise was added to Figure 3 to represent the breach described by Marahatta and Parajuli (2015), and this hydrograph was routed in 1-D and 2-D over smoothed and raw topography to compare downstream flooding with the BREACH hydrographs (Figure 5).

FIGURE 5 

Time series of water depth and flow speed at Beni Bridge and Maldhunga Bridge. Base case uses smoothed topography, n = 0.04 s/m^1/3, and the BREACH hydrograph.

Discussion

The results described above indicate which models are applicable to the Kali Gandaki landslide flood of 2015, and hint at which models should be applied in similar environments elsewhere. Hindcasting the flood hydrograph at the location of the dam itself showed the most uncertainty. The empirical formulas referenced in Table 1 produced a range of peak flow rates from 87 to 8970 m³/s; the BREACH result of 1696 m³/s was near the lower end of that range. Nonetheless, the BREACH hydrograph, when combined with either 1-D or 2-D routing models, produced a result for flow speed and depth at Beni Bridge in close agreement with observations. If the peak flow had been much larger, as it would be under many of the empirical formulas, the flow at Beni Bridge would also have been deeper, causing overtopping of the bridge. In reality, overtopping did not occur. Likewise, the maximum flow depth hindcast at Maldhunga Bridge in this study was approximately 10 m. Although it is a taller bridge, with a lower chord 16 m above the riverbed, the larger flow rates indicated by the Table 1 formulas could have caused overtopping here as well, but none was reported.

Khanal et al (2015) applied a similar set of empirical dam-break formulas to glacial lake outburst floods in the Himalaya, and concluded that they are likely to overestimate peak flows. Davies et al (2007) modeled a landslide dam-break flood in New Zealand and compared it to field observations of the flood, and noted that “empirical relations produced some wild values” of peak flow rates, whereas the BREACH model slightly underestimated the flow rate. The reason for the BREACH underestimation was that the landslide dam had slopes steeper than 5:1 (horizontal:vertical), which BREACH could not simulate because of limitations in its sediment transport equation. In the case of the Baisari landslide (Table 2), however, slopes based on geometry reported by Collins and Jibson (2015) were less than 5:1 and so did not suffer from the same limitation.

Empirical formulas need to be applied with great care, as they account for only dam height and/or lake volume. Uncertainty of input parameters (eg uncertainty of lake hypsometry and dam characteristics) to a physically based model by Walder and O'Connor (1997), however, also produced distributions of peak discharge that ranged between 3 orders of magnitude (Schwanghart et al 2016). A physically based model like BREACH thus requires well-constrained input data and knowledge about the geometry and geotechnical composition of the dam. When these data are available, physically based models can be more robust when applied within the limits for which they were formulated.

Concerning downstream flood routing, our analysis showed that hydrological smoothing and carving is a necessity in regions where detailed topographic data are unavailable. At Beni Bridge, the flood hindcast with unsmoothed topography arrived later, and its flow speed was slower, than observed in the field (Figure 5). Further downstream at Maldhunga Bridge, the flood based on unsmoothed topography arrived over 1 hour later and was much weaker and shallower than that shown by the other models.

If a model severely underestimating flood hazard were adopted for use in evacuation planning or warning, the result could be catastrophic. For example, in the case of the 2011 Great East Japan Earthquake and Tsunami, an initial evacuation warning was issued based on the first few seconds of recorded seismic data, resulting in an estimated earthquake of magnitude 7.9 (Japan Meteorological Agency 2013). Based on this estimate, predicted tsunami heights at the coast were underestimated, so much of the populace did not feel an urgent need to evacuate. As more seismic data were acquired, the estimated earthquake magnitude was raised to 9.0, allowing accurate hindcasting of tsunami heights, but not in time to issue effective evacuation warnings (Japan Meteorological Agency 2013). Likewise, if a model such as BREACH were used with unsmoothed topography for flood evacuation planning or warning, the underestimation of flood magnitude and overestimation of lag time before flood wave arrival could result in insufficient evacuation and increased casualties (Jonkman et al 2002).

At the opposite extreme, the 1-D model used in this study hindcast a slightly deeper flood and shorter flood wave arrival lag time at both bridges than did the other models (Figure 5). This result is conservative, which is a positive aspect of a warning system if used only once. However, if a conservative model is used in a warning system over a long time period, resulting in multiple false alarms or near misses, the population at risk is likely to stop heeding warnings, which would negate the protective effect of the warning system (Frieser 2004; Barnes et al 2007; Parker et al 2009; Dillon-Merrill et al 2011).

An interesting aspect of this study was the value of crowd-sourced video available on the Internet for assessing flood models. Particularly in countries with few gauging stations, such videos are an invaluable source of information on alpine flash floods, whose extents and dynamics can usually not be derived from satellite-based imagery as is the case for large, lowland rivers. Future work should develop methods that utilize and augment this increasingly available video material for flood analysis. Similar techniques were used for validation of inundation models during the Great East Japan Tsunami (Bricker et al 2015) and Typhoon Haiyan (Roeber and Bricker 2015).

Conclusions

The importance of hydrological smoothing and carving is clearly displayed by the results above, as the simulation using raw topography hindcast lower flow depth and speed, and a longer flood wave arrival lag time, than was observed in the field. In developing countries such as Nepal, where accurate data on river channel topography do not exist, basing a flash flood warning system on raw global topography such as SRTM (NASA 2014) or ASTER (NASA 2012) could result in delayed evacuation and increased casualties.

Using smoothed/carved topography and basing a warning system on either a 1-D or 2-D model would give useful and reliable information to those in the hazard area. In the case studied here, the 1-D model hindcast a more conservative result than the 2-D model. Because 1-D models can be executed much more quickly than 2-D models (on the Intel Core i7 processor used for this study, the 1-D model ran in about 4 minutes, whereas the coarse-grid 2-D model took 4 hours and the fine-grid 2-D model took 24 hours), and because the 1-D model produces a conservative result, the 1-D model is more useful for real-time warning systems. Because of the run time required for a 2-D model, application to a real-time warning system might not be practical, especially in developing countries with insufficient access to powerful computing resources. However, the accuracy of 2-D models makes them useful for predisaster hazard assessment and hazard map generation. Furthermore, they can be used in real-time warning systems if they are executed en masse beforehand to generate a large number of flood scenarios, from which a pattern-recognition system such as a neural network (Campolo et al 1999, 2003) or decision tree can be used to extract the appropriate scenario during an actual flood.

The delicate middle ground between insufficient warning (as results from the model with unsmoothed topography) and excessive warning (as results from the 1-D model) is different for every locale and needs to be investigated before floods occur in order for credible warning systems to be implemented. In the case of the 2015 Kali Gandaki flood, the 2-D model generated a result most in line with field observations. However, this cannot be assumed to hold true universally. For disaster preparedness planning in a given basin, a range of models should be evaluated by hindcasting historical events for which field data or crowd-sourced videos exist.