The 2015 magnitude 7.8 Gorkha earthquake and its aftershocks weakened mountain slopes in Nepal. Co- and postseismic landsliding and the formation of landslide-dammed lakes along steeply dissected valleys were widespread, among them a landslide that dammed the Kali Gandaki River. Overtopping of the landslide dam resulted in a flash flood downstream, though casualties were prevented because of timely evacuation of low-lying areas. We hindcast the flood using the BREACH physically based dam-break model for upstream hydrograph generation, and compared the resulting maximum flow rate with those resulting from various empirical formulas and a simplified hydrograph based on published observations. Subsequent modeling of downstream flood propagation was compromised by a coarse-resolution digital elevation model with several artifacts. Thus, we used a digital-elevation-model preprocessing technique that combined carving and smoothing to derive topographic data. We then applied the 1-dimensional HEC-RAS model for downstream flood routing, and compared it to the 2-dimensional Delft-FLOW model. Simulations were validated using rectified frames of a video recorded by a resident during the flood in the village of Beni, allowing estimation of maximum flow depth and speed. Results show that hydrological smoothing is necessary when using coarse topographic data (such as SRTM or ASTER), as using raw topography underestimates flow depth and speed and overestimates flood wave arrival lag time. Results also show that the 2-dimensional model produces more accurate results than the 1-dimensional model but the 1-dimensional model generates a more conservative result and can be run in a much shorter time. Therefore, a 2-dimensional model is recommended for hazard assessment and planning, whereas a 1-dimensional model would facilitate real-time warning declaration.
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).
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.
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.Chanson 2004). Integrating Equation 1 over the water column results in ū/umax = 0.86, where ū is the depth-averaged flow speed. Based on this, the particle image velocimetry result corresponded to a depth-averaged flow speed of 7.7 to 8.6 m/s.
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.
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 m3, impounding a lake with a volume of 8 million m3. 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).
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.
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.38S0.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/m1/3 was determined and applied throughout the domain. As a sensitivity analysis, another model run was conducted with roughness n = 0.03 s/m1/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.
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.
The BREACH hydrograph of flow exiting the dam site is shown in Figure 3, with a maximum flow rate of 1696 m3/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 m3/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).
Validation based on observed flow depth and speed at Beni Bridge
At Beni Bridge, the 2-D coarse grid simulations hindcast a maximum flow depth of 4 to 5 m (Figure 5), which is in agreement with the video of the flood discussed above. In this instance, the 1-D simulation is conservative, hindcasting a larger flow depth. All of the simulations using the smoothed/carved topography hindcast a maximum depth-averaged flow speed at Beni Bridge close to the range of 7.7–8.6 m/s measured from the video.
Because the 2-D simulations using raw (unsmoothed) topography showed maximum flow speeds much slower than that observed at Beni Bridge, it is likely that they were losing too much energy because of spuriously large obstacles (pools and riffles). This hypothesis was tested via a sensitivity analysis of the raw topography simulation with an unrealistically small roughness of n = 0.01 s/m1/3. As Figure 5 shows, this also resulted in a flow speed at Beni Bridge smaller than the measured value, lending credence to this hypothesis.
Validation based on flood wave arrival time at Beni Bridge
In the raw topography simulations, the flood wave arrived at Beni Bridge more than one-half hour later than in the smoothed topography simulations. Based on direct observations, this delay also appears to be spurious. Witnesses reported that landslide overtopping began about 5 PM local time (Anonymous 2015b; Collins and Jibson 2015; Marahatta and Parajuli 2015). Because sunset on 24 March 2015 was 7 PM local time, darkness would have prevented local residents from recording videos of the flood arriving at Beni Bridge after 7 PM. However, as Figure 2 shows, videos were recorded even after the bore arrived. The flood wave hindcast by the BREACH hydrograph and routed over raw topography arrived at Beni Bridge over 2 hours after overtopping began, at which time darkness would have prevented locals from recording video. Thus, as with the spuriously slow flow speed discussed above, the simulation using raw topography hindcast a spuriously late flood arrival time. In contrast, in simulations using the Marahatta and Parajuli (2015) hydrograph, both the smoothed and raw topography flood waves arrived at Beni Bridge before sunset, with maximum flow reached 1 hour earlier.
The 1-D simulation hindcast a flood wave arrival one-half hour earlier than the 2-D simulation's hindcast. However, observations showing the absolute time of flood wave arrival were not recorded, so it is not possible to decide which is more accurate based on observations. Nonetheless, the 1-D simulation, as shown above, hindcast a flow depth at Beni Bridge larger than the observed flow depth. Because a deeper flow equates to greater wave front speed, the early arrival of the 1-D front is likely a result of the conservatively large flow depth hindcast by the 1-D model. This holds true for both the BREACH and Marahatta and Parajuli (2015) hydrographs.
Sensitivity to roughness
In order to investigate the effect of Manning's n on the model result, another 2-D simulation was evaluated with n = 0.03 s/m1/3. As expected, this caused an increase in flow speed, a reduction in flow depth, and a shortening of the arrival lag time of the flood wave at Beni Bridge. However, the difference was small, and both depth and speed fell close to the range of observed maximum flow depth and speed at the bridge. Nonetheless, because Jarrett's (1984) formula estimated n = 0.04 s/m1/3, the simulation with n = 0.04 s/m1/3 was used as the base case for the current analysis.
Sensitivity to model grid resolution
The 2-D simulation with fine grid resolution hindcast a deeper, faster flow at Beni Bridge than the base-case 2-D simulation. The arrival time of the flood wave hindcast by the fine-resolution simulation was also earlier than that of the base-case simulation. This is likely because the fine-resolution simulation undergoes less numerical dissipation (Morton and Mayers 1994) than the coarse simulation.
Flood behavior at Maldhunga Bridge
Figure 5 also shows the results for Maldhunga Bridge, further downstream. The simulation using raw topography hindcast flow depth and speed much lower, and a flood wave arrival over 1 hour later, than that hindcast by the simulations based on smoothed topography. The 1-D model hindcast a greater flow depth than the 2-D models but a slower flow speed. Sensitivity to roughness was similar to that at Beni Bridge, whereas the fine-grid simulation hindcast slightly less flow depth and greater flow speed than the coarse-grid model.
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 m3/s; the BREACH result of 1696 m3/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).
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.
This work was supported by a grant from the president of Tohoku University and the director of the International Research Institute of Disaster Science. We thank Phillip Mineart for help with the BREACH model and Fabien Retif for helping us obtain Malpasset dam-break topography data. Topographic data were provided by Alpine Consultancy Pvt Ltd Nepal, and the Department of Survey, Nepal.