Flood Hazard Assessment for the Tori Levee Breach of the Indus River Basin, Pakistan

: Levee breaches are some of the most common hazards in the world and cause the loss of lives, livelihoods, and property destruction. During the 2010 ﬂood in Pakistan, the most devastating breach occurred at Tori Levee on the right bank of the Indus River, downstream of the Guddu Barrage, which caused residual ﬂoods in northern Sindh and the adjoining regions of the Balochistan province. In this study, 2D unsteady ﬂow modeling performed for Tori Levee breach computed residual ﬂood inundation by coupling a HEC-RAS (Hydrological Engineering Centre—River Analysis System) 2D hydraulic model with remote sensing and Geographic Information System techniques. The model performance was judged by comparing the observed and simulated water levels (stage) during peak ﬂow at seven different gauging stations located within the Indus River reach and daily ﬂood extents and multi-day composites. The quantitative values for the calibration and validation of the HEC-RAS model showed good performance with a range of difference from 0.13 to − 0.54 m between the simulated and observed water levels (stage), 84% match for the maximum ﬂood inundation area, and 73.2% for the measure of ﬁt. The overall averages of these values for the daily ﬂood comparison were 57.12 and 75%, respectively. Furthermore, the simulated maximum ﬂow passed through the Tori Levee breach, which was found to be 4994.47 cumecs (about 15% of peak ﬂow) with a head water stage of 71.56 m. By using the simulated ﬂows through the Tori Levee breach, the ﬂood risk maps for the 2010 ﬂood identiﬁed hazard zones according to the ﬂood characteristics (depth, velocity, depth times velocity, arrival time, and duration). All the ﬂood risk maps concluded the fact that the active ﬂood plain was uninhabitable under ﬂood conditions.


Introduction
Floods are among the most common natural disasters around the world, claiming a substantial loss of life and property in the vicinity of rivers. Extreme rainfall is a major contributor to flooding along with minor contributions from snow-glacier melt, earthquakes, landslides, glacial lake outburst flooding (GLOF), and dam and levee failure/overtopping [1,2]. Pakistan is one of the most vulnerable countries in the world in terms of exposure to different natural calamities [3]. According to the Long-Term Several researchers [12,[24][25][26][27] have used a 1D-2D coupled approach for modeling levee breach flood inundation. However, this approach is not suitable for large rivers and vast floodplains [28]. The exchange of flow between 1D and 2D elements is also complex and may introduce substantial uncertainties. Hence, the use of 2D models for modeling rivers and floodplains becomes the right choice. 2D models based on the Saint Venant Equations calculate flows in x and y spatial dimensions [29]. Although, the computational burden for 2D models is larger than 1D and 1D-2D models, more accuracy can be achieved. Unlike using cross-section interpolation to represent river topography as in the case of 1D, 2D models provide flexibility to incorporate a physical and topographic representation of terrain to improve accuracy for complex terrains. 2D models with the integration of remote sensing data are frequently used by several researchers [1,[30][31][32] for levee breach and urban flood inundation modeling.
In Pakistan, flood hazard mapping and modeling are mostly performed using 1D models for active flood plain, while the risk of flooding behind levees (inactive flood plain) is often underestimated or ignored. Therefore, for this study, we chose a region which is highly vulnerable to devastating floods and the failure of levees, particularly Tori Levee, historically breached more than 19 times [33]. To the best of our knowledge, not a single published study has been carried out on levee breaching, dam breaching, nor residual flood mapping using 2D unsteady flow modeling in Pakistan. Therefore, this study is an attempt to provide a comprehensive methodology for the research community in Pakistan and globally to analysis flood risk behind levees through the modeling of a devastating flood event which occurred in Pakistan during the 2010 monsoon period. The Tori Levee breach event was modeled using HEC-RAS 2D Model, by comparing observed verses simulated water levels at several gauging stations (within the active floodplain) and daily flood extent of MODIS (MODerate Resolution Imaging Spectroradiometer) with the simulated flood extent for that day (in avulsion area that is inactive floodplain). The discharge passed through the Tori Levee breach was unknown and predicted by using the HEC-RAS model. Subsequently, several flood components including water depth, velocity, product of depth and velocity, flood arrival time and flood duration maps were determined by using a calibrated and validated HEC-RAS 2D model of a Tori Levee breach.

Study Area
The Indus River is the largest river system and a major source of water resources in Pakistan. The study region is mainly influenced by the summer monsoon precipitation pattern with approximately 130 mm mean annual areal precipitation. However, the river reach selected in the current study mostly receives major flows from high-altitude regions of Pakistan due to western distribution precipitation and snow-and glacier-melt. In this study, the Indus River reach, affected by a devastating flood in 2010, was selected with a length of 120 km between two barrages (i.e., the Guddu and Sukkur Barrages). The study area (reach and flood plain area) lies between 28 • 23.08 km (downstream). The channel slope or river gradient (vertical to horizontal ratio) of the selected reach between the Guddu-Sukkur is 1:7000 (V:H) (as per data provided by Sindh Irrigation and Drainage Authority (SIDA)). The study area lies in the flat topography of the lower Indus Basin in the northern part of the Sindh. Historically, the 2010 flood is considered the most devastating disaster with a return period of 20 years as estimated using generalized Pareto (GPA) distribution at site flood frequency analysis performed by [34], resulting in 2000 fatalities and the displacement of more than 20,000,000 people for weeks to several months [9,35]. Sindh and Khyber Pakhtunkhwa (KP) were the most impacted provinces in Pakistan. KP was affected due to flash floods, while embankment breaches were the major cause of flooding in Sindh and Punjab. The most devastating breach of the 2010 flood occurred at Tori Levee, the main embankment on the right bank of Indus downstream of Guddu with the top, bottom elevation and length equal to 72, 69.13 and 2879 m, respectively. The riverbed and flood plain elevation in the vicinity of Tori Levee was approximately was 65 and 68 m, respectively. The Tori Levee on the night of 7th and 8th of August, collapsed when Guddu reached a flood peak of 31,714.8 cumecs [36][37][38]. The Tori Levee breached at three different locations and depths equal to 272, 103 and 629 m with a cumulative bottom width of the breach sections equal to 1004 m, estimated using the method proposed by [14] (Figures S1 and S2).
This breach on the right bank of Indus caused the avulsion of the Indus River and overwhelmed the irrigation and flood protection infrastructure in northern Sindh and adjoining Baluchistan and caused widespread devastation in areas well outside the Kacha area active flood plain [2]. The Tori Levee breach immediately transformed the scale of the disaster. The flow from the breach traveled over 300 km and rejoined the mainstream of the Indus River via Manchar Lake near Shewan ( Figure 1). This residual flood lasted for a long time as the Indus River remained above flood stage at Guddu until 2nd September 2010.

Datasets
The daily discharge, stage data, and the architectural design specification of Guddu and Sukkur Barrages were obtained from the Guddu and Sukkur Barrage authorities, for the flood events of 2010 and 2015. High flood level (HFL) data for the 2010 flood along the right levees were acquired from Sind Irrigation and Drainage Authority (SIDA). The longitudinal geometrical and 3D profiles of all levees were generated using data obtained from Sindh  [14]. In this study, channel bathymetry was also incorporated using cross-section (extracted from ALOS DEM) interpolation and combining the interpolated channel and floodplain DEM [14,16]. The daily flood extent was extracted for the flood period using daily MODIS Terra images at 250 m spatial resolution for the study area. The land use land cover (LULC) classification of Landsat 5 was acquired from the Ministry of Environment

HEC-RAS Model Application
The HEC-RAS 5.0 two-dimensional (2D), a widely used hydrodynamic model developed by the U.S. Army Corps of Engineers, was employed in this study for the hazard mapping of the 2010 flood in Pakistan. The HEC-RAS model simulated the flood using the combination of 2D Saint-Venant and diffusive wave equations [39]. The model can use the 2D Saint-Venant equations for multiple situations such as dynamic flood wave, mixed flow regime, abrupt contraction, wave propagation modeling and super elevation around bends. However, the diffusive wave equation is applicable for floodplain flows with minimal acceleration and turbulence. The HEC-RAS 2D model provides flexibility to calibrate using a manning (n) roughness coefficient based on LULC, the contraction and coefficient of expansion using boundary conditions of hydrographs (stage and flow), normal depth and rating curve [40]. Additionally, in HEC-RAS, the RAS mapper facilitates the creation, and improvement of terrain data (channel bathymetry), creation and definition of user-defined land cover data or importing from web services [41]. Additionally, the RAS mapper can visualize the outputs of HEC-RAS such as the static and dynamic flood hazard maps in form of water depth, velocity, product of depth and velocity, water surface elevation, shear stress, inundation boundary, flood arrival time, duration time and recession time at specific depth levels [40]. The detailed description of HEC-RAS is provided by [40,41]. The detail of the methodology adopted in this study is shown in Figure 2.  Figure S3. Furthermore, the Tori Levee breach was modeled using breach specifications extracted from the previous study [36] (see Figure S2), carried out for the 2010 flood. Overall, the Tori Levee was breached into three sections with different bottom widths (272, 103 and 629 m) with sum of 1004 m (see Figures  S2 and S3), also confirmed by [36]. However, the top and final bottom elevation of the Tori Levee was incorporated as equal to 72 and 69.13 m obtained from SIDA, respectively. The breach section was considered as trapezoidal, and the final specifications/dynamics of the bottom width (b), side slope (S), cross-sectional length (top width of Tori Breach) and height (h) of the breach section were determined to be 1004 m, 80 degrees, 1459 and 2.87 m, respectively. The top width of the Tori Levee breach was measured with the help of the bottom width (i.e., 1004 m) by using the trapezoidal channel design equation, as suggested by [14]. However, the other breach parameters like breach formation time and water elevation triggering the breach (water surface elevation) was set at 168 h (with a oneminute interval) and 71.34 m, respectively, estimated using methodology adopted by [14]. The aforementioned factors and parameters were used for the 2010 flood (calibration flood event). The 2015 flood (validation event) occurred within the river reach without any levee breach, therefore, the validation was carried out by keeping the river reach and structural design factors of the barrages and parameters constant, except Tori Levee breach parameters.
Furthermore, the LULC data were adopted to define Manning's roughness (n) values for the model domain. The built-up clusters from open street map (OSM) and Wikimapia were extracted and used to improve the original LULC, as used by several researchers [42,43]. Furthermore, the calibration of HEC-RAS for the 2010 flood was carried out by the optimization of Manning's roughness (n) values for different LULC classes as per the guidelines provided by the United States Geological Survey (USGS) [44] and some other available sources in the literature [45,46], altering breach dimensions and breach parameters and integrating other flood control structures. The LULC data were spatially bifurcated for the optimization of Manning (n) roughness coefficients in the active flood plain (river channel and floodplain within the levees) and inactive floodplain (avulsion area). The final optimized Manning's (n) roughness values for different classes are shown in Table 1. Finally, after providing complete data inputs, the model was run with computation at one (1) minute interval and 168 h as a breach develop time in levee calculated using a method proposed by [14] as shown in Figure S2, on a core i5 machine with 4 Gb RAM, that successfully computed the simulation in 18 h. The statistical indices were applied for the evaluation of HEC-RAS simulations during calibration and validation in two different ways. Firstly, within levees calibration was performed using an observed high flood level at seven gauging stations and simulated water level data at peak flow. Secondly, in the case of an avulsion area, due to the nonavailability of gauging data, the satellite (UNOSAT (for composite flood extant) and MODIS (daily basis flood propagation)) images were adopted for the calibration. It is important to note that the complete simulation was carried out within the levee water levels and avulsion plain extant in a single HEC-RAS model run. The HEC-RAS model automatically takes the outputs of first step daily flow (cumecs) and stage (m) time series  October) at the Tori Levee during the flood event as an input. The observed and simulated water levels (stage) were evaluated using difference, Nash-Sutcliffe efficiency (NSE), index of agreement (d), relative deviations of Nash-Sutcliffe efficiency (E rel ) as well as relative deviations of index of agreement (d rel ) at seven gauging stations (i.e., Guddu, Begari, Kacha Kharif, Ghora Ghat, Tori, Sukkur Beghari and Sukkur Barrage), as the Equations (1)-(3) are given below: The Nash-Sutcliffe Efficiency (NSE) was computed by using Equation (1), as follows: where Q o = mean of water levels, Q m = simulated water levels Q t o = observed water levels at time t.
The difference was estimated using following established Equation: (2): The index of agreement (d) was computed using Equation (3) given below: where O i = observed water levels, P i = simulated water levels, O = average observed water levels, P = average simulated water levels.
The relative deviations of Nash-Sutcliffe efficiency (E rel ) and relative deviations of index of agreement (d rel ) are as follows: With modifications in the Nash-Sutcliffe efficiency and Index of Agreement, the differences between the observed and simulated water levels can be computed as relative deviations, which can reduce the substantial effect of the absolute differences during peak flows or water levels. Furthermore, the aforementioned efficiency indices are discussed in detail by [16].
Furthermore, the 'percent match' and 'measure of fit' statistics were used for the comparison of daily flood extents simulated by HEC-RAS and extracted from MODIS, and maximum composite flood extent provided by UNOSAT [47]. The best simulation values for the 'percent match' and 'measure of fit' were obtained by tuning the manning roughness coefficient (n) up to the benchmark of the least betterment in statistical descriptors, as final optimized 'n' values for each LULC class are given in Table 1. Finally, after setting breach parameters and channel calibration, the HEC-RAS integrated model was used to compute a stage-flow hydrograph at the Tori Levee during the breach event and for the modeling of flood in the avulsion plain. The modeling of the avulsion plain was comprised of flood inundation/hazard simulations (mapping of flood-depth, -velocity, -depth to velocity product-, -inundation duration for 0.5, 1.0, 1.50 and 2.0 m depth and -arrival up to 96 h or 4 day maps).

Calibration and Validation of HEC-RAS Model
Initially, three different DEMs (ALOS 30 m, ASTER 30 m and SRTM 30 m) were adopted for the river morphology, topography and breach formation dynamics measurements. However, during the analysis, it was observed that the ALOS 30 m DEM performed far better than ASTER and SRTM 30 m DEM. The ASTER and SRTM DEM were unable to show the proper river morphology for the selected river reach. The sensitivity of the model for ALOS, ASTER and SRTM was evaluated for maximum depth (m), velocity (m/s) and flood extent (km 2 ). It was observed that the ALOS DEM 30 m performed fairly well for maximum depth, velocity and flood extent in comparison to ASTER and SRTM (Table 2), as confirmed by [40].
The calibration and validation process of the HEC-RAS model include the water levels (stage) at seven gauging stations within the Indus River levees and flood extent over the avulsion plain. The results for the stage within the Indus River levees at seven different gauging stations produced a difference ranging between 0.13 and −0.54 m during the peak flow discharge (Table 3). Overall, the results showed that the HEC-RAS model performed fairly well to simulate water levels during peak floods in comparison to the observed. The computed statistical descriptors show substantial efficiency for the comparison between observed and HEC-RAS simulated water levels at all gauging stations, including five along the right bank Indus River levees and two on the barrages (Table 3).
However, for the calibration of an avulsion plain, the flood extent simulated by the HEC-RAS model was compared with the maximum flood extant composite provided by UNOSAT and the daily flood extent extracted from MODIS. Overall, a mixed pattern of over-and under-estimations of flood extant was found over the avulsion plain for 2010 flood (Figure 3). The HEC-RAS model performed fairly well for the reproduction of flood extent with 84% (≈12,362 km 2 ) overall matching, 16% (≈2290 km 2 ) over-and under-estimations and 73% measure of fit.   Additionally, the simulated flood propagation and flood extents on a daily basis were compared with cloud free MODIS Terra images (Table S1). The factual values in Table S1 shows that the HEC-RAS model performed fairly well in reproducing daily flood propagation in comparison to the MODIS extent with average, overestimation and underestimation matching values of 75, 37 and 25%, respectively (see Figure S5). Since the model was calibrated for 2010 flood (within the river reach and avulsion plain), therefore, by keeping all modeling parameters constant, the HEC-RAS model was validated for 2015 flood inundation over the avulsion plain, which occurred in the same Indus River reach. The result for the validation with 89% match and 77.22% measure of fit shows HEC-RAS model efficient ( Figure S4). The over-(15%) and under-estimations (11%) were found to be slightly less for the validation in comparison to the calibration of the HEC-RAS model. The over-and under-estimations may be associated with the use of the coarse resolution DEM (i.e., DEM at 30 m resolution), uncertainties involved in modeling and MODIS based flood extent, as confirmed by [30]. Additionally, during the 2010 flood, the local flood control and management efforts, including engineered breaches (intentionally created), can cause of different simulated flood propagation in contrast to actual. The results obtained for the calibration and validation of HEC-RAS are in line with previous studies conducted in different regions of the world [1,17,40].

Tori Levee Breach Outflow Hydrograph
The flow passed through the Tori Levee breach during the flood event was unknown. Therefore, to investigate the flood inundation simulations over the alluvial plain, the Tori Levee breach flow (i.e., stage-flow hydrograph) was obtained from the calibrated HEC-RAS Model (Figure 4). The maximum flow passed through the Tori Levee breach simulated by HEC-RAS and was found equal to 4994.47 cumecs (about 15% of the upstream Guddu peak flow) with a head water (HW) stage of 71.56 m and occurred on 19 August and 10 August, 2010. Subsequently, the flood rejoined the Indus River on the 17 September 2010 after making its way through Manchar Lake (see Figure 1 for locations). The results obtained in this study are in line with a previous study conducted in this region with slight underestimations demonstrated in the current study [33] as reported that the 22% (≈7000 cumecs) Indus River flows which were diverted through Tori Levee at the time of breach. This may be associated with the use of a different simulation method, as [33] used multi-temporal remote sensing images for the assessment in contrast to the current study (HEC-RAS model).

Flood Inundation and Risk Maps
The calibrated HEC-RAS model was further used for flood inundation and risk hazard mapping of the study area by using flows simulated at Tori Levee as shown in Figures 5-9. These maps included maximum flood-depth, -arrival, -duration, -velocity and -depth times velocities. The depth threshold, taken as 0.35 m, is beyond which the water depth can be life-threatening [48]. The floodplain depth (1.0 to 15 m) shown in Figure 5 shows the maximum depth >5 m within the main Indus River reach, Sukkur Barrage and Manchar lake. However, the depth over most of the alluvial areas varied between 0.36 and 2.25 m, with few places in different cities (i.e., Jacobabad, Shuhdadkot, Qambar, Johi and Dadu) ranged between 2.25 and 5.0 m which shows a greater life threat. Additionally, the maximum depth was noticed upstream of the barrages and near Tori Levee ( Figure 5). Furthermore, the flood velocity is another major hazard for human life and property, so extreme velocity maps were simulated using HEC-RAS to highlight the most fatalityprone areas due to high flood velocity. The floodwater velocity of 1.50 m/s can cause instability for an adult [49]. In relation to that, Figure 6 depicts the maximum velocity within the main Indus River course substantially above the threshold, resulting in a hazardous conditions, in comparison to the alluvial and floodplain where the range of velocities (0.01 to 0.80 m/s) was demonstrated by itself. Overall, low velocities were observed in the floodplain which may be associated with a low slope topography and an increase in elevation level on the western flood protection levees, as well as a high level of roads and irrigation canals compared to the natural surface level in the floodplain [33,36]. The slow propagation of the breaching water was dually confirmed by other studies [33] and reported the average velocity of 0.1 to 0.3 m/s over the avulsion plain down-valley, some propagation slowdown due to the high elevation of roads, irrigation canals and topographic gradient on the west side of the floodplain, as confirmed by [20,33]. However, in the southern part, the flood propagation was driven by the topographic gradient without the contribution of major breach flows from the Indus River. Conversely, the flood velocity within the main course of the Indus River channel travelled at three times this rate [33]. Additionally, the combined effect of the depth and velocity of water is a more realistic measure for the flood hazard analysis. The individual 1.5 m/s of velocity and 0.35 m of depth are life threatening, so the product of both can provide a combined threshold value (0.52 m 2 /s) for the quantification of human hazard [48]. Figure 7 depicts that most of the floodplain was at a <0.5 m 2 /s velocity to depth product, which is close to the threshold value (i.e., 0.52 m 2 /s). However, several spots within the floodplain were also found to be extremely dangerous to life in the flood situation with values of 0.51-0.60 m 2 /s, which make the active floodplain uninhabitable, particularly close to Jacobabad and Jhatpat city. Furthermore, it was observed that the main route of the flood in the floodplain starts from the Tori Levee breach and includes major populated cities (i.e., Jacobabad, Jhatpat, Gari Khairo, Shahdadkot, Qambar and Johri) where between 0.31 and 0.50 m 2 /s values resulted in greater life threat (Figure 7).  Overall, extensive damage over the avulsion plain was not only associated with the flood wave impacting inadequate and weak levee structures but also due to the lack of enough storage capacities of two major reservoirs (Mangla and Tarbela dams), barrages located upstream of Tori Levee, and flood diversion channels along the Indus River even at the locations where breaches occurred in the past. According to previous studies, both reservoirs and all barrages located upstream of Tori Levee were at peak storage capacity during the 2010 flood event [9,33], which created substantial pressure on flood protection bunds/levees constructed along the Indus River. Furthermore, reportedly, the Tori Levee breach event leads to a devastating flood towards the northern avulsion plain and particularly the 2.7 km breach at the location of Tori Levee. The Tori Levee was already historically in poor condition due to continuous failures and had lost 1.7 m in design height prior to the 2010 flood event due to soil erosion and a lack of maintenance. Although this extreme event damaged the region, it was not a rare event occurring along this river. Several major floods have occurred due to a breach of different levee structures and spilling of water, aforementioned facts confirmed by [33]. Overall, this study can conclude that the Indus River and its tributaries are highly vulnerable to precipitation extremes and structurally inadequate flood protection levees along the reach. Therefore, both non-structural flood management strategies such as flood forecasting and warning systems [50][51][52][53][54] and structural flood management strategies such as additional storage reservoirs [55][56][57], flood diversion structures, and the widening of the river channel at the most critical locations are needed [58].This study is an attempt to help the future planning for improving flood management by identifying potential flood diversion channels location and management options and by considering the flood extent and inundation maps along the Indus River.

Conclusions
This study carried out the modeling of the devastating 2010 flood caused by the Tori Levee breach in Pakistan by using a well-known hydrodynamic model (i.e., HEC-RAS 2D). In this study, the HEC-RAS model was calibrated and validated for the peak water levels at different gauging stations and flood extent for the years 2010 and 2015. Additionally, the flood inundation risk analysis using the HEC-RAS model provided maps for the spatial variations in flood depth, velocity of water, depth times velocity product, flood arrival, and flood retention duration over the Northern Sindh Province of Pakistan (study area).
This study showed that the HEC-RAS 2D model is capable of modeling the water levels (stage) at seven different gauging stations during the peak flood with a sum of square errors value of 0.90. The HEC-RAS 2D model performed well in calibrating and validating the 2010 and 2015 flood extents with overall 84 and 89% matching, respectively, between the extents simulated by HEC-RAS and the UNOSAT and MODIS images. Additionally, the calibrated and validated HEC-RAS model simulated how the flows (4994.47 cumecs) passed through the breached Tori Levee, resulting in the devastating flood. The inundation map of flood depth over the study area showed maximum water depths over several spots in floodplain. However, the velocity of flooded water was maximum within the main course of the Indus River. The flood poses very little treat to life in floodplain with the values below the life-threatening threshold except within the vicinity of the breach. However, the product of the velocity and depth map showed that the active floodplain is uninhabitable due to life threat with several spots far beyond the threshold (0.52 m 2 /s) values over the floodplain of the study area. Furthermore, the flood duration map at different depth levels showed the depth of flood water retained at particular spots for particular times. The results showed that the 0.35 m depth of water was retained for approximately 70 days at most of the floodplain after the initiation of the flood event. Several spots were under the 2.0 m depth of water retention over floodplain.
Overall, the study showed that apart from the Ghauspur Town and nearby settlements which were instantly flooded after the breach, the majority of the areas were completely flooded after the four days of Tori Levee breach, resulting in human life loss and economic loss. These losses could be avoidable by providing residual flood risk maps for the avulsion area in advance, resulting in quick evacuation from the active floodplain to safe havens located in nearby regions. The results and methodology adopted in this study can be used for potential hazarding mapping and management for the regions frequently facing floods. Additionally, since, the Tori Levee is a hotspot for breaching and was already breached several times before the 2010 flood, therefore, the results of this study can be useful for the local government authorities, federal flood commission (FFC) and National Disaster Management (NDMA) to evacuate and shift the population under threat during expected flooding conditions in the future.

Informed Consent Statement:
It is stated that the informed consent was obtained from all subjects involved in the study. Data Availability Statement: Data sharing not applicable.