# Two-Dimensional Hydrodynamic Modelling of Koshi River and Prediction of Inundation Parameters

^{*}

**Corresponding Author:**Mukesh Raj Kafle, PhD Scholar, Department of Civil Engineering, Institute of Engineering, Pulchowk Campus, Nepal, Tel: 9779851041278, Email: [email protected]

*
Received Date: May 08, 2018 /
Accepted Date: May 28, 2018 /
Published Date: Jun 05, 2018 *

### Abstract

This paper presents the results of modelling study of Koshi River. The modelling approach is based on twodimensional hydrodynamic model. The simulation is carried out with model software Nays 2DH. The study analyses the inundation parameters, hazard assessment criteria, flood inundation extent delineation and identification of hazardous areas in different discharge scenarios of 25, 50 and 100 years return periods flow. Based on goodnessof- fit tests and fitting parameters, generalized extreme event (GEV) distribution method is adopted for flood frequency analysis. The model is calibrated with measured water surface elevations and simulation results. The root means square errors (RMSE) and correlation R2 between measured values and simulation results are 0.95 m and 0.98 respectively. The study results conclude that within the stretch of around 50 km from Chatara to Koshi barrage the flood will not overtop embankments and left overbank is under low danger zone. However, islands within the embankments namely Shukrabare, Rajabas, Khairatol, Shivchowk and Galphadiya are vulnerable to inundation. The modelling approach proposed in this study is an attractive option for modelling exceptional flood events when limited data and resources are available.

**Keywords:**
Hydrodynamic model; Return periods; Inundation; Flood prediction; Simulation; Two-dimensional

#### Introduction

Floodplain inundation is a major environmental hazard in both the developed and developing countries [1]. Yet no consensus exists concerning the level of model and data complexity required to achieve a useful prediction and inundation extent and a number of techniques present themselves for the prediction of inundation extent resulting from fluvial flood event [1]. Till the late seventies, the methods used for solving the Saint-Venant equations appeared to be satisfactory with mathematical models found to be adequate for large number of applications [2]. Bates et al. mentioned the further development of two-dimensional finite element models of river flood flow. They applied the two-dimensional finite element model to the Missouri river, Nebraskan with integration of hydraulic modelling and remote sensing [3]. Han et al. and Chang et al. have also reported onedimensional (1-D), two-dimensional (2-D) coupled modelling of river flood plain flow. The models have used a full dynamic equation for the channel flow and for the two-dimensional flood plain flow; a diffusion wave approximation is utilized [4,5]. Anderson, Robayo et al. and Knebl et al. discussed that flood inundation modelling involves hydrologic modelling to estimate peak flows from storm events, hydraulic modelling to estimate water surface elevations, and terrain analysis to estimate the inundation area [6-8]. Wright et al. presented a methodology for using remotely sensed data to both generate and evaluate a hydraulic model of floodplain inundation for a rural case study in the United Kingdom-Upton-upon-Severn [9]. Zheng et al. developed a distributed model for simulating flood inundation integrating with rainfall-runoff processes using SRTM-DEM data and some remote sensing data sets in the environment of GIS for Maruyama River basin, Japan [10]. The increasing use of two dimensional models over one dimensional models during the past decade has been partly driven by developments in digital elevation modelling (DEM’s), especially from airborne LiDAR data [11]. Thus, as the capability has developed it has been necessary to better understand the effects of moving to two dimension's different applications. Comparisons between one dimensional, two dimensional and coupled one-two-dimensional river modelling approaches (Horritt and Bates) have highlighted conceptual problems with the one-dimensional approach applied to overbank flows when compared to the sometimescomplex flow pathways simulated by two dimensional models [12-15]. The past few decades have witnessed exceptional progress in the development of computer models to simulate the propagation of floods, much of which has been driven by the need for tools to evaluate flood risk. Most of these methods are based on the shallow water equations, a nonlinear hyperbolic system of equations describing the conservation of mass and momentum of water in two horizontal dimensions [16]. Numerical flood models are generally used to assess the consequences of a potential failure and produce flood extent maps. They also provide information on flow depth, velocity and derived parameters such as flow force and intensity that can be used to assess potential damage to structures [17]. Bermudez et al. discussed modelling approaches in practical applications [18]. Unsteady openchannel flow equations govern the basis of using either onedimensional (1-D) or two-dimensional (2-D) hydrodynamic models. The physical characteristics of the river channel and flood plain are governing factors in the choice of model approach. Around river confluences and extensive flat flood plains the flow is strongly twodimensional. In such scenario, the choice of one-dimensional approach is not suitable. Having incremental computational power and capability of solving the shallow water equations, two-dimensional (2- D) models have become a good choice in recent years [17].

The modelling reach of the Koshi River can be described as a braided channel containing many islands both large and small. Man made embankments and structures along both left and right banks of the portion have altered the natural state of river significantly. The resulting non-uniform geometry of the river creates strong crosscurrents in some areas. Due to the complex geometry, large flood plain and flow conditions present, a two-dimensional hydrodynamic model is required. A two-dimensional hydrodynamic model would be able to capture the effects of the crosscurrents and the complex configuration of the river. There have been several software developed including TELEMAC 2D [19], RMA2 [20], River 2D [21], Mike 21 [22], TUFLOW [23], JFLOW [24], HYDRO AS-2D [25], InfoWorks 2D [26] and Nays2DH [27] to simulate two-dimensional flow. Among available software, selection of Nays 2DH is based on its numerous worldwide applications and particular applications in fluvial channel. It has capability to incorporate effects of complex boundaries and river bed shapes, confluences of main channel and tributaries, bottom friction evaluation using Manning’s roughness parameter, vegetation effects using drag force and flow field around obstacle.

This study aims to suggest an attractive option of modeling approach for modelling exceptional flood events in large river basin and flood plain when limited data and resources are available.

The rest of this paper is structured in the following manner. A description of the study area and data is presented in Section 2. In Section 3, a brief review of methodology with flow chart is provided. Section 4 summaries the computational results. Analysis of computed results and discussion are dealt in section 5. Finally, the conclusion is reported in Section 6.

#### Study Area and Data

The study area (**Figures 1** and **2**) is the Koshi River reach from Chatara (Hydrological station) to 50 km downstream (Koshi Barrage). The Koshi river (also called Saptakoshi) known from its seven tributaries is a trans-boundary river flowing through Tibet (China), Nepal and India. It is one of the largest tributaries of the Ganges River. The entire Koshi river basin has a catchment area of 69,300 km^{2} up to its confluence with Ganges in India, out of which 29,400 km^{2} lies in China, 30,700 km^{2} in Nepal and 9,200 km^{2} in India. The Koshi basin occupies eastern part of Nepal.

Koshi River in Nepal has seven major tributaries: Sunkoshi, Tamakoshi, Dudhkoshi, Indrawati, Likhu, Arun and Tamor. At Barahkshetra in Nepal it emerges from mountains and become the Koshi River. After flowing another 58 km it crosses into Bihar, India near Bhimnagar and after another 260 km joins the Ganges near Kursela. The river has a total length of 729 km. It is a snow fed river which originates from the Himalayan range, crosses Mahabharat range and debouches into plains at about near Chatara (approximate 50 km upstream of Koshi Barrage) and enters into Indian territory.

The data requirements for proposed model are (a) topographic data of the channel and floodplain to act as model bathymetry, (b) time series of bulk flow rates (c) roughness coefficients for channel and floodplain, which may be spatially distributed, and (d) data for model calibration. The acquired data and sources of data are presented in **Table 1**.

**Meteorology and hydrology data**

Koshi river system receives rainfall from the south-west monsoon, which generally sets in the first or second week of June and withdraws in the first week of October. The spatial variation of annual rainfall is also very high in the basin ranging from 809.2 mm in Nepalthok (Station No 1115) to 4377.8 mm in Nam (Station No: 1301). The hydrological data from Chatara gauging site, located at 26°52ˈ00ˈˈ N and 87°09ˈ30ˈˈE and elevation of 140 m, is taken for the study and analysis of peak flood. The drainage area of the Koshi basin at Chatra hydrological gauging station is 54100 km^{2}. The average monthly variation of flow at this station shows that the peak discharge occurs in August. The mean annual flow of Chatara Station is 1620 m^{3}s^{-1} (**Figures 3** and **4**). The extreme flow in the Koshi River at Chatara station was 25879 m^{3}s^{-1} in 1968 followed by 24240 m^{3}s^{-1} in 1954 and 24000 m^{3}s^{-1} in 1980. The flood flow on 18 August 2008 at Chatara was only about 4250 m^{3}s^{-1} when the flood disaster in Koshi River was initiated (Source: Department of Hydrology and Meteorology, Government of Nepal).

**Geographic data**

The river topographical survey is essential for generating digital elevation model (DEM) required for model simulation. Ground survey and bathymetry survey were carried out by total station and ecosounder. The survey provides information on rivers and structures profile, including cross-sections of the riverbed, the riverbanks and any structures in the river channel. The survey also provides the details of ground elevation close to the riverbank, which is vulnerable to inundation. River channel cross-section and longitudinal survey at 42 locations are carried out between Chatra and Koshi Barrage. Topographical survey data are processed for use in two-dimensional hydrodynamic model. Cross sections were spaced at varying distances throughout the study reach approximately at 1 km interval. Each cross section was given a station number in river km that corresponded to its location on the river (**Figure 5**).

**Bed material (Grain size)**

The riverbed of the Koshi River consists of granular, sandy material. A relative density and a characteristic grain size characterize the sediment. The sediment density is unknown and therefore set at a standard value of 2650 kg/m^{3}. A constant value for the entire model is selected. The graphical representation of characteristics of grain size distribution based on the results of the sieve analysis is presented in **Figure 6**.

**Riverbed roughness**

The friction of riverbed is set using Manning’s roughness n. The choice of methods of estimation of hydraulic roughness depends on the availability of field data (**Figures 7-9**). When field data are not available, the traditional approach is to use handbook methods or analytical methods to predict the hydraulic roughness value. Previous works developed equations based on field data to estimate hydraulic roughness [28-32]. Among of them Strickler type equation estimates Manning’s roughness considering median grain size diameter (D_{50}). Moreover, this equation is applicable to wide-shallow channels having width-depth ratio greater than 10 where the hydraulic radius can be replaced by the mean depth. So, in this study, estimation of Manning’s roughness by Manning-Strickler equation is appropriate. The equation is

− *Eq*(1)

Where n=Manning’s roughness, D_{50} is the median size diameter of bed particle. The summary of results shows that average value of Manning’s ‘n’ for upper reach resulted 0.05 and for lower reach 0. 04 (**Table 2**).

River reach from Chatara to d/s (Km) | D50 trendline (mm) | Manning’s | Average ‘n’ |
---|---|---|---|

4 | 2.2 | 0.054 | 0.05 |

5.8 | 1.6 | 0.51 | |

7 | 1.4 | 0.05 | |

8.4 | 1.2 | 0.049 | |

10 | 1 | 0.047 | |

11.6 | 0.9 | 0.046 | |

14.6 | 0.7 | 0.045 | |

17 | 0.6 | 0.044 | 0. 04 |

18 | 0.6 | 0.044 | |

19 | 0.6 | 0.043 | |

20.8 | 0.5 | 0.043 | |

22.4 | 0.5 | 0.042 | |

23 | 0.5 | 0.042 | |

24 | 0.5 | 0.042 | |

28 | 0.4 | 0.041 | |

30 | 0.4 | 0.041 | |

32 | 0.4 | 0.04 | |

34 | 0.4 | 0.04 | |

35 | 0.3 | 0.04 | |

39.4 | 0.3 | 0.039 |

**Table 2:** Estimation of Manning’s n by Strickler’s method.

For the estimation of Manning’s roughness, in total 209 soil samples are collected and grain size distribution analysis is performed to understand the nature of riverbed. Samples of river bed materials ensuring all representative areas under water, adjacent dry bed, flood plain and embankment material in various points in the study reach are collected at every 2 km from Chatara to barrage and at every 5 km along the banks. Samples from the armour layer (upper 15 km of river reach) that extends from the surface down to a depth of 2 times more than the existing large particles are also collected. Sand beds are sampled by obtaining bulk samples. The results show that most of the samples from upper 15 km stretch consists of gravel with sand whereas the samples downstream of 15 km are composed of poorly graded sand. The difficulties on getting more representative value for the whole cross-section the sieve curves are merged to obtain a combined and more representative sample and D_{50} and D_{90} of this combined sample are determined.

**Floodplain vegetation roughness**

In the model, resistance exerted by vegetation is set with the drag coefficient of vegetation C_{D}, the area of interception by vegetation per unit volume A_{S} and the height of vegetation. The area of interception by vegetation per unit volume as can be set at each computational cell. The area of interception by vegetation per unit volume A_{S} is calculated by the equation proposed by Shimizu et al. as follow:

− *Eq*(2)

Where N_{S} is the number of vegetation, D_{S} is the average diameter of trunk and S_{S} is the sampling grid width.

The land use pattern shows that the study reach is covered with grasses, islands, bushes and forests. So, to estimate drag coefficient (C_{D}), submerged and non-submerged conditions are taken into account. Petryk and Bosmajian approach verified by Baptist is used for the roughness estimation of non-submerged vegetation [33,34]. Whereas, an approach developed by Klopstra et al. and verified and simplified by Baptist is applied to estimate submerged vegetation roughness [35].

The information related to land use pattern including vegetation type, height of vegetation, vegetation density, diameter of stem is collected from satellite images and field survey.

#### Methodology

**Numerical simulation**

This study is structured on numerical simulation. The modelling approach is based on two-dimensional hydrodynamic model. Nays 2DH model software developed by i-RIC [27] is used for simulation. Nays 2DH is a computational model for simulating horizontal twodimensional (2D) flow, sediment transport, morphological changes of bed and banks in rivers. In this model, a general curvilinear coordinate system is adopted allowing direct consideration of complex boundaries and riverbed shapes. Either upwind difference scheme (first order) or Cubical Interpolation Pseudo-Particle (CIP) method are available options for the finite differencing to the advection terms in momentum equations. Numerical computation method is high-order Godunov scheme referred to as the CIP method. This method splits the integration of the momentum equations of flow into a non-advection and pure advection phase. The solution of non-advection phase is cubically interpolated and then advected to the solution of grid points. The CIP method has been shown to solve the problem of boundedness while introducing little numerical diffusion thus enabling highprecision local interpolation and algorithm implementation is more straightforward than other high-order upwind schemes [27]. In order to achieve numerical stability, a minimal amount of turbulence is required. In horizontal two -dimensional flow calculations, since threedimensional flow structures are completely neglected, unrealistic horizontal large vortices tend to be generated. In this model, adjustment of user-defined parameters for controlling such vortices in fluvial channel is difficult. Eddy viscosity coefficient is the apparent kinetic viscosity coefficient of a flow in turbulent state. It specifies the amount of turbulent diffusion that occurs at the intersection of elements with differing velocities. Based on the river and floodplain characteristics constant eddy viscosity is used for turbulent field calculation. Constant eddy viscosity coefficient (v_{t}) in Equations (6) and (7) being taken as 10^{-6} (m^{2}/s) [27].

The basic equations in Orthogonal Coordinate System (x,y) are:

Continuity equation

− *Eq*(3)

Momentum equations in x- and y- directions

− *Eq*(4)

In which

− *Eq*(5)

− *Eq*(6)

− *Eq*(7)

− *Eq*(8)

Where, h=water depth, t=time, u, v=depth averaged velocities in x and y directions, g =gravitational accelerations, H=water depth, τ_{x} τ_{y}=components of shear stress in river bed in x- and y- directions, F_{x}, F_{y}=components of drag force by vegetation in x- and y- directions, C_{f} =drag coefficient of the bed shear stress,ν_{t}=eddy viscosity coefficient, C_{D}=drag coefficient of vegetation,α_{s}=area of interception by vegetation per unit volume, h_{ν}=minimum value of water depth and height of vegetation.

Here, bed shear stress τ_{x} and τ_{y} are expressed by using coefficient of bed shear force C_{f}. the coefficient of bed shear force C_{f} is estimated by Manning’s roughness parameter n as follow:

− *Eq*(9)

Where, n is the Manning’s roughness parameter estimated by Eq (1), g is the acceleration due to gravity.

**Model setup**

The model domain extends up to 50 km upstream from the Koshi barrage to Chatara. Model grids are created with river survey data and discretised into a structured computational mesh with a variable size (from 5 m to 20 m) to ensure correct definition of river channel and flood plain. An inflow condition is defined at upstream boundary. Koshi barrage is considered as the downstream boundary. Riverbed roughness is selected based on particle size distribution analysis and Manning-strickler equations (Eq. 1). Size of bed material is set as median size diameter (D_{50}) based on results of particle size distribution analysis. Manning’s ‘n’ for upper 15 km reach is set 0.05 and for lower reaches 0.04. For flood plain vegetation roughness, drag coefficient is selected based on submerged and non-submerged conditions. Simulations are run for 3 hours with output results of every 60 seconds.

The methodology of the proposed model is presented by the flowchart below:

**Results**

**Flood prediction**

Three basic approaches namely hydro-meteorological approach, frequency analysis and empirical formulae are in practice to the estimation of design flood. Application of hydro-meteorological approach is ruled out in the present case due to the non-availability of the required data. Many empirical formulae have been devised for the purpose of estimating peak flows. These formulae can be safely applied to the areas for which they have been specifically developed. No particular formula will give precise results for all the sites. Moreover, these formulae must be used with great prudence, and must never be used unless their origin has been investigated. Use of empirical formulae for estimation of design flood is, therefore, not recommended.

The frequency analysis of hydrological extremes such as floods and droughts, aims at estimating the design values corresponding to a given probability of occurrence [36]. Flood frequency analysis (FFA) is a statistical technique, which fits a probability distribution to historical stream flow data series observed at a given location within a catchment. The commonly used probability distributions for fitting the flood series are Gamma, Gamma (3P), generalized extreme value (GEV), Gumbel’s extreme value (EV-I), log-Pearson type 3, Pearson type 3 or Weibull distribution. Based on the guidelines, some countries generally adopt a particular distribution to model floods, viz., in the USA the log Pearson 3 (LPE3) is used for extreme floods; generalized logistic (GLO) and Pearson type III (PE3) are normally recommended in the UK and China, respectively; while in Australia, the log-Pearson 3, generalized extreme value, and generalized Pareto distributions are recommended probability distributions [37]. However, selection of a single probability distribution for all sites in a country is not driven by specific theoretical considerations, although different goodness of fit tests can be adopted for selection of most suitable distribution.

For the selection of a particular distribution representing the flood series three different most widely used goodness of fit tests, viz., Kolmogorov-Smirnov (K-S), Anderson-Darling (A-D) and Chi-square (χ^{2}) tests, either based on probability density functions f(x) or cumulative distribution functions f(x), carried out to identify the best fit model. The said tests calculate test- statistics, which are used to ascertain the suitability of given distribution to fit the data [38].

The goodness-of-fit tests are carried out using EasyFit software. The results of Anderson-Darling, Kolmogorov-Smirnov and Chi-square test for Chatara gauging station are summarized in **Tables 3** and **4**. The test statistics are computed, and the probability distributions are ranked based on the lowest values of the test statistics. The summarized parameters and test statistics clearly indicate that Generalized Extreme Value (GEV) distribution best fits the peak discharge series at Chatara from the ranks of goodness of fit tests. Log Pearson 3 is ranked second and Gamma (3P) is ranked third. Flood quantiles for different return periods for generalized extreme value and log Pearson 3 distribution are calculated using cumfreq software and results from best fit distribution from goodness of fit tests are used for further analysis (**Table 5**).

S.N | Distrubution | Kolmogorov Smirnov | Anderson Darling | Chi-Square | Average Rank | |||
---|---|---|---|---|---|---|---|---|

Statistic | Rank | Statistic | Rank | Statistic | Rank | |||

1 | Gamma | 0.22354 | 6 | 3.4179 | 5 | 8.2681 | 6 | 6 |

2 | Gamma(3P) | 0.11616 | 3 | 2.5907 | 3 | 5.0973 | 3 | 3 |

3 | Gen.Extreme Value | 0.07348 | 1 | 0.2846 | 1 | 1.6625 | 1 | 1 |

4 | Log pearson 3 | 0.08253 | 2 | 0.3556 | 2 | 2.8424 | 2 | 2 |

5 | Weibull | 0.19827 | 5 | 4.3573 | 6 | 7.9007 | 5 | 5 |

6 | Weibull(3P) | 0.12641 | 4 | 2.9648 | 4 | 7.5101 | 4 | 4 |

**Table 3: **Goodness of fit-summary.

S.N | Distrubution | Parameters |
---|---|---|

1 | Gamma | α=3.3832, β=2584.6 |

2 | Gamma (3P) | α=0.8198, β=4405.8, γ=4970 |

3 | Gen.Extreme Value | K=0.43778, s=1650.9, µ=6546.8 |

4 | Log Pearson 3 | α=1.7589, β=2997.9, γ=8.4551 |

5 | Weibull | α=2.7458, β=9470.5 |

6 | Weibull (3P) | α=0.87421, β=4094.8, γ=4970 |

**Table 4: **Fitting results.

Return periods (years) | Log Pearson 3 | Gen. extreme value |
---|---|---|

5 | 11416 | 13703 |

10 | 14670 | 17019 |

25 | 20143 | 21208 |

50 | 25456 | 24316 |

100 | 32083 | 27401 |

**Table 5: **Comparison of Discharges (cumecs).

**Performance evaluation of model**

The performance of model has been evaluated by model calibration. Mainly two methods namely traditional methods and remote sensing techniques are in practice for hydrodynamic model calibration. Previous studies discussed traditional methods of model calibration that rely on field measurements (water depth, velocity and flooded area). Such measures represent only discrete information of the flow conditions for selected sections. Therefore, an accurate calibration requires as many observations as possible. Remote sensing techniques tackles the problem-using aerial and satellite imagery. Previous studies discussed hydrodynamic models calibration with topographic information obtained from airborne laser altimetry, Horritt et al. from satellite synthetic aperture radar (SAR) sensors and Dung et al. from inundation maps. In this case, because of unavailability of required imagery, maps, sensors remote sensing techniques are avoided and traditional methods of model calibration relying on field measurements is adopted.

Several statistical indices root mean square error (RMSE), normalized root mean square error (NRMSE), Nash-Sutcliffe coefficient (E) and Pearson correlation coefficient (R) have been used to evaluate the predictions obtained by numerical models. Normalized root mean square error (NRMSE), non-dimensional forms of the RMSE are used to compare RMSE with different units. The Nash- Sutcliffe model efficiency coefficient (E) is commonly used to assess the predictive power of hydrological discharge models. The RMSE error is interpreted as a deviation of the simulated results from the measurements and used to measure of the difference between values predicted by a model and the values actually observed from the environment that is being modeled. As the root mean square error (RMSE) is a dimensional variable it is convenient to seal it to a reference value, which, in most cases is the mean value of measurements. The square of the Pearson correlation coefficient (R^{2}) known as coefficient of determination describes how much of the variance between the two variables is described by the linear fit. Since the available observed data is water surface elevations statistical indices root mean square error (RMSE) and correlation (R^{2}) are used for calibration. The calibration is performed through the comparison between measured and predicted water surface elevations for 74 points spatially distributed throughout the main channel of the reach at a known discharge of 398 m^{3}s^{-1}. The observed water surface elevations are compared with computed results adjusting different relaxation coefficients 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8 and 0.9. The comparison shows that relaxation coefficient of 0.4 yields the better results. However, none of the simulations are good match at all locations. The Root Mean Square Error (RMSE) and Correlation (R^{2}) are 0.95 m and 0.98 respectively (**Figures 10a** and **10b**) that is considered acceptable for further processing of model (**Table 6**).

Profile no. | Simulated energy gradient elevation (m) | Embankment Elevation (m) |
||
---|---|---|---|---|

25 years return period | 50 years return period | 100 years return period | ||

38 | 102.77 | 103.05 | 103.31 | 109.277 |

36 | 101.4 | 101.71 | 102.03 | 105.468 |

35 | 100.97 | 101.25 | 101.54 | 104.179 |

33 | 99.15 | 99.39 | 99.62 | 100.657 |

32 | 98.42 | 98.66 | 98.88 | 99.206 |

29 | 95.86 | 96.1 | 96.32 | 98.109 |

28 | 95.19 | 95.43 | 95.66 | 97.303 |

27 | 94.53 | 94.76 | 94.98 | 96.695 |

26 | 93.71 | 93.93 | 94.13 | 95.634 |

25 | 92.24 | 92.46 | 92.68 | 94.931 |

24 | 91.26 | 92.46 | 91.7 | 93.869 |

23 | 90.54 | 90.76 | 90.97 | 92.984 |

22 | 89.84 | 90.06 | 90.28 | 92.254 |

21 | 89.31 | 89.53 | 89.75 | 91.586 |

20 | 88.69 | 88.92 | 89.15 | 90.259 |

18 | 87.49 | 87.74 | 87.97 | 89.025 |

16 | 86.24 | 86.48 | 86.69 | 87.991 |

14 | 85.18 | 85.41 | 85.63 | 86.666 |

13 | 84.67 | 84.91 | 85.11 | 87.247 |

7 | 80.82 | 81.1 | 81.35 | 83.136 |

6 | 80.55 | 80.82 | 81.22 | 82.836 |

4 | 79.73 | 79.98 | 80.22 | 82.245 |

3 | 79.1 | 79.34 | 79.57 | 80.19 |

**Table 6: **Comparison of water surface profiles.

**Simulation results**

**Water surface elevation profiles:** For inundation extent delineation, a comparison of water surface elevation with the elevation of embankment in different scenarios is must. The simulated results are compared with river survey data for three discharge scenarios of 25, 50 and 100 years return period flow. Comparison of computed results with river survey whose embankment elevation are available show that the flood will not overtop the embankment for all three-discharge scenarios. The numbering of profiles is based on cross section survey lines starting from 42 at upstream end (Chatara) to 1 at downstream end (Koshi barrage).

**Hydraulic parameters and contour maps: **The two-dimensional (2- D) flow simulation is carried out for different discharge scenarios i.e., annual average flow and 5, 10, 25, 50 and 100 years return period floods. The scalar values for maximum depth of water, shear stress, velocity of flow and Froude number for all discharge conditions are computed (**Table 7**). Contour maps of water depth and flow velocity are developed for three discharge scenarios 25, 50 and 100 years return periods (**Figure 11**). The computed values and generated contour maps are compared, discussed and analysed.

**Figure 11:** Simulated contour map (a) water depth for 25 years return period flow (b) velocity for 25 years return period flow (c) water depth for 50 years return period flow (d) velocity for 50 years return period flow (e) water depth for 100 years return period flow (f) velocity for 100 years return period flow.

Discharge (m^{3}/s) |
Max.Depth (m) | Max.Velocity (m/s) | Max.Shear Stress (N/m ^{2}) |
Max.Froude Number | Remarks |
---|---|---|---|---|---|

1620 | 2.48 | 1.75 | 42.1 | 1.23 | Average annual discharge |

13703 | 5.74 | 3.21 | 85.7 | 1.01 | 5 years return period |

17019 | 6.28 | 3.38 | 92.6 | 0.955 | 10 years return period |

21208 | 7.15 | 3.42 | 106 | 0.812 | 25 years return period |

24316 | 7.45 | 3.59 | 127 | 1.78 | 50 years return period |

27401 | 7.76 | 3.86 | 124 | 0.952 | 100 years return period |

**Table 7: **Simulated results.

#### Discussion

**Inundation hazard assessment criteria**

Inundation extent is not normally a direct prediction of flood routing models; but is more often estimated indirectly through predicted water surface elevations, depth of water, velocity of flow etc. Inundation hazard risk may be assessed with the known values of water depth and flow velocity. Based on water depth-flow velocity relationship, Richardson et al. classify the flood risk zones in three categories namely low-danger zone, Judgment zone and high-danger zone. If depth of flow (y) ≤ 0.3 m and flow velocity (v) ≤ 5 ms^{-1} then the flood risk classification is categorized as low-danger zone. For 0.3 m ≤ y ≤ 1 m and flow velocity (v) <2.5 ms^{-1} the risk classification is also under low danger zone. If depth of flow (y) <1 m and flow velocity (v) >2.5 ms^{-1}, then it falls under the category of Judgment zone. Similarly, for 1m ≤ y ≤ 2 m, the risk classification is categorized as judgement zone for all magnitudes of flow-velocity. If the depth of flow (y) >2 m then flood risk is classified as high-danger zone irrespective with magnitude of flow velocity [39].

**Inundation extent delineation**

The flood depth map is a widely distributed instrument. The values of water level (depth) can be derived from flow models (2-D and 1-D) for river flooding, from statistical analyses or from observations. There is a wide range of applications of such maps. The flood depth map provides information about the water depth in a particular location for a given recurrence interval (or probability) of flood. This map is used to serve as a basic product to establish danger and flood damage maps, city and village planning and risk management (evacuation).

Flood velocity distribution maps represent the flow velocity. Flow velocity information is much more difficult to get than water depth information. The flood hazard in a particular location is represented by the velocity of the flowing water (or sediment in case of debris flow) or by the velocity of the flood propagation. Flow velocity map is used for planning of flood defense measures or any structure within the flood area. The best tool for technicians. It can also be used as planning tool for emergency response, evacuation schemes, and implementation of temporal flood protection measures.

Inundation extent for three discharge scenarios of 25, 50 and 100 years return periods are computed and compared with depth of water and velocity of flow distribution maps (**Figure 11**) and values (**Table 7**). In the case of 25 years return period flow, the maximum depth of water 7. 15 m and flow velocity 3.42 ms^{-1} are predicted. For 50 years return period flow, the maximum predicted depth of water and flow velocity are 7.45 m and 3.59 ms^{-1} respectively. Similarly, in the scenario of 100 years return period flow, the maximum predicted water depth and flow velocity reaches to 7.76 m and 3.86 ms^{-1} respectively. The generated contour maps classify the depth of water and velocity of flow distribution with different color bands. In depth contour map, dark blue ranges from 0 to 1.13 m, light blue from 1.13 m to 2.24 m, green from 2.24 m to 5.55 m, yellow from 6.55 m to 6.65 m and red from 6.65 m to 7.76 m (**Figures 11a, 11c** and **11e**). Besides, in velocity distribution map, dark green ranges from 0 to 0.55 ms^{-1}, light green from 0.55 ms^{-1} to 1.10 ms^{-1}, pink from 1.10 ms-1 to 1.65 ms^{-1} and red from 1.65 ms^{-1} to 3.86 ms^{-1} (**Figures 11b, 11d** and **11f**). The generated flood-inundated maps portray that islands within the embankments are vulnerable and the left over bank is under low-danger zone. The higher water depth cells (marked light blue, green, yellow, red) and velocity of flow grids (marked light green, pink, red) close to the settlement specify high-danger zones. This is the indication of vulnerability to the settlements (**Figure 11**).

**Identification of hazardous areas**

The study reach is embanked, and flow is confined between the embankments. So, it is very difficult to produce hazard maps behind the embankment with associated probabilities and to predict embankment breach location owing to lack of sufficient geometric data behind the embankment. The model results and comparison of simulated water surface elevation for different scenarios with the elevation of embankment show that the flood will not overtop the embankments for all three discharge scenarios of 25, 50 and 100 years return period flow. The simulated water surface elevations are validated with historical records of flood events. The reach from Chatara to Koshi barrage is embanked since 1953 after the construction of koshi barrage. Since 1953, flood events of 25879 m^{3}/s, 24240 m^{3}/s and 24000 m^{3}/s are recorded at October 5, 1968, August 24, 1954 and June 25, 1980 respectively. No overtopping of embankments is reported on those flood events. However, breach of left embankment at 10 km upstream of Koshi barrage triggered even at low discharge scenario of 4729 m^{3}/s on August 18, 2008 resulting massive loss of lives and property both in Nepal and India. That means inundation to settlements due to overtopping of embankment will not occur with the embankment in place. Based on generated inundation contour maps, the islands within the embankments namely Shukrabare, Rajabas, Khairatol and Shivchowk are identified as inundation hazardous areas for 25 and 50 years return period flood. Moreover, for 100 years return period flood, the island Galphadiya is added to the aforementioned identified inundation hazardous areas (**Figure 11**).

**Sensitivity analysis**

Sensitivity of model results are investigated with change of input parameters Manning’s n, relaxation coefficient and drag coefficient.

The simulated maximum depth of flow decreases by approximately 32% on the decrease value of Manning’s n by 50%. In contrast, maximum depth increases by around 9% on the increment of Manning’s n by 50%. On the other hand, simulated velocity magnitude increases by 41% on the decrease of Manning’s n by 50%. In reverse, the increment of Manning’s n by 50% results decrease in maximum velocity magnitude by 26%. The effects of changes values of input parameters relaxation coefficient and drag coefficients are negligible in simulated depth of flow. However, velocity of flow is affected noticeably with change of these parameters. Decrease in relaxation coefficient by 50% results increase of velocity magnitude by 3%. Approximately 5% decrease in velocity magnitude is observed in the increment of relaxation coefficient by 50%. Moreover, decrease of drag coefficient by 50% results increment of flow velocity magnitude by 10%. Disparately, increment of drag coefficient by 50% results decrease of velocity magnitude by 9%. Form the results of sensitivity analysis Manning’s roughness coefficient glances the most sensitive input parameter of the model (**Figure 12**).

The uncertainty that correlates with all data and input parameters determines the accuracy of model. The uncertainties of topographic data are depending on the accuracy of the data acquired during survey. With increasing accuracy and resolution of topographic data, the uncertainties get lower. Hydrological data is used for boundary and initial conditions. They give major information about the physics. Data for these values are obtained, in this case, from hydrometric measurements of the catchment area. Gauges installed at the upstream end of the study area record flow. Roughness parameterization is a very important parameter of hydraulic models. They represent land use types, which, themselves indicate conditions for flow velocity and runoff behavior. The roughness coefficient should also be distinguished between the riverbed and floodplains. And of course, scenarios like dam breaches are also possibilities that contribute to the long list of uncertainties. It is unlikely to eliminate all uncertainties. But by performing calibration, the influence of uncertainties is minimized.

**Model applicability**

The approaches discussed in this paper establish the best practice case study on prediction of inundation parameters, inundation extent delineation and consequent risk identification. Sophisticated twodimensional hydrodynamic models support the methodology and analysis. The extent of data availability, physical characteristics of the river systems and previous academic practices have governed the choice of models. The model software Nays 2DH have better computational efficiency and are freely available. Facilities of settingup, running and output visualization are incorporated as an interface in the model. Computational time of 3 h is applied. Similar results can be expected with the application of same methodology in other twodimensional hydrodynamic models based on depth averaged shallow water equations instead of Nays 2DH.

#### Conclusions

This study analysed the inundation parameters and consequences of flood events of Koshi River. The study was based on two-dimensional hydrodynamic modelling approach. The modelling results showed that inundation to settlements owing to overtopping of embankments would not occur with the embankment in place. The islands within the embankments were vulnerable and the left overbank was under low hazard. The vulnerable islands identified from generated inundation maps were Shukrabare, Rajabas, Khairatol, Shivchowk and Galphadiya. The hazardous areas should be viewed with special attention and given utmost importance. The velocity distribution maps showed that the higher velocities occurred near Rajabas area where the flow was concentrated. So, this area would be susceptible to bank erosion and thus the more vulnerable.

The results have uncertainties associated with the accuracy of DEM, hydrological, sediment, vegetation data and the modelling approach. Nevertheless, they can be utilized as base line information for the policy making to combat against induced disaster, planning flood preparedness and show the possibilities for extending worldwide to other rivers having similar physical characteristics and flood plain. The modelling approach proposed in this study is an attractive option for modelling exceptional flood events when limited data and resources are available.

#### Acknowledgments

This work is self-funded by the authors. It is a part of PhD program of correspondent author.

#### Author Contributions

Mukesh Raj Kafle was responsible for this current research article in the framework of his PhD program and initially wrote the manuscript. Narendra Man Shakya directed the study by helping to interpret the results and improving the quality of manuscript.

#### References

- Horritt MS, Bates PD (2001) Predicting floodplain inundation: Raster-based modelling versus the finite-element approach. Hydrol Process 15: 825-42.
- Priessmann A (1976) Use of mathematical models, Unsteady flow in open channels, BHRA, pp: 23-28.
- Bates PD, Horrit MS, Smith CN, Dason D (1997) Integrating remote sensing observations of flood hydrology and hydraulic modeling. Hydrological Process 11: 1777-1795
- Chang TJ, Hsu MK, Huang CJ (2000) A GIS distributed watershed model for simulating flooding and inundation. Journal of the American Water Resources Association 36: 975-988.
- Han KY, Lee JT, Park JK (1998) Flood inundation analysis resulting from levee break. Journal of Hydraulic Research 36: 747-759.
- Anderson DJ (2000) GIS-based hydrologic and hydraulic modeling for flood delineation at highway river crossings. Master’s Thesis, Univ. of Texas at Austin, Austin, Texas.
- Robayo O, Whiteaker T, Maidment D (2004) Converting a NEXRAD map to a floodplain map, Proc., AWRA GIS and Water Resources III Conf., Nashville, Tenn.
- Knebl MR, Yang ZL, Hutchison K, Maidment DR (2005) Regional scale flood modeling using NEXRAD rainfall, GIS, and HEC-HMS/RAS: A case study for the San Antonio River Basin summer 2002 storm event. J Environ Manage 75: 325-336.
- Wright NG, Villanueva I, Bates PD, Mason DC, Wilson MD, et al. (2008) Case study of the use of remotely sensed data for modelling flood inundation on the river Severn, U. K. Journal of Hydraulic Engineering ASCE, pp: 533-540
- Zheng N, Tachikawa Y, Takara K (2008) A distributed flood inundation model integrating with rainfall-runoff processes using GIS and Remote Sensing data. The International Archives of the Photogrammetry, Remote Sensing and Spatial Information Sciences. Beijing 37: 1513-1518.
- Baldassarre GD (2012) Floods in changing climate:Inundation Modelling. International Hydrology Series. Cambridge University Press. 1st edn November 30, 2012.
- Neal J, Villanueva I, Wright N, Willis T, Fewtrell T, et al. (2012) How much physical complexity is needed to model inundation? Hydrol Process 26: 2264-2282.
- Horritt MS, Bates PD (2002) Evaluation of 1D and 2D numerical models for predicting river flood inundation. Journal of Hydrology 268: 87-99.
- Werner MGF (2004) A comparison of flood extent modelling approaches through constraining uncertainties on gauge data. Hydrology and Earth System Sciences 8: 1141 1152.
- Tayefi V, Lane SN, Hardy RJ, Yu D (2007) A comparison of one and two dimensional 794 approaches to modelling flood inundation over complex upland floodplains. Hydrological Processes 21: 3190-3202.
- De Almeida GAM, Bates P (2013) Applicability of the local inertial approximation of the shallow water equations to flood modeling. Water Resour Res 49: 4833-4844.
- Álvarez M, Puertas J, Peña E, Bermúdez M (2017) Two-dimensional dam-break flood analysis in data-scarce regions: The case study of Chipembe dam, Mozambique Water, Switzerland.
- Bermudez M, Neal JC, Bates P, Coxon G, Freer J, et al. (2017) Quantifying local rainfall dynamics and uncertain boundary conditions into a nested regional-local flood modeling system. Water Resources Research 53: 2770-2785.
- Lang P, Desombre J, Ata R, Goeury C, Hervouet JM (2014) Telemac-2D Software. User Manual: EDF-R&D.
- Barbara PD, Joseph LV, McAnally WH (2011) Users Guide for RMA2 Version 4. 5.
- Steffler P, Blackburn J (2002) Introduction to Depth Averaged Modeling and User's Manual. University of Alberta.
- MIKE Powered by DHI 2015 MIKE 21 and MIKE 3 Flow Model FM, Hydrodynamic Module, Short Description.
- BMT WBM (2016) TUFLOW User Manual. BMT WBM.
- JBA (2014) JFlow Award-winning 2D hyraulic model.
- Nujić M (2014) HYDRO-AS-2D. 2D-Flow model for water management practice. Aachen: Hydrotec engineering company for water and environment mbH.
- Innovyze (2011) InfoWorks 2D, Applying the 2D module to Collection Systems. Technical Review.
- Shimizu Y, Takebayashi H (2014) Nays2DH Solver Manual. iRIC Project.
- Strickler A (1923) Contributions to the question of the velocity formula and the roughness numbers for currents, channels and closed. Bern, Switzerland, Communications of the Federal Office for Water Management, no. 16.
- Henderson PM (1966) Open channel flow. New York, Macmillan, p: 522.
- Limerinos JT (1970) Determination of the Manning coefficient from measured bed roughness in natural channels. U. S Geological Survey Water-Supply Paper 1898-B, p: 47.
- Bray DI (1979) Estimating average velocity in gravel-bed rivers. American Society of Civil Engineers, Journal of the Hydraulics Division 105: 1103-1122.
- Jarrett RD (1984) Hydraulics of high-gradient streams. American Society of Civil Engineers. Journal of hydraulic Engineering 110: 1519-1539.
- Petryk S, Bosmajian G (1975) Analysis of flow through vegetation. Journal of the Hydraulics Division 101: 871-884.
- Baptist MJ (2005) Modelling floodplain biogeomorphology. PhD Thesis, Delft University of Technology. Delft, The Netherlands.
- Klopstra D, Barneveld HJ, Van JNM, Velzen VEH (1997) Analytical model for hydraulic roughness of submerged vegetation. Proceedings of the 27th IAHR Congress theme A Managing Water: Coping with Scarcity and Abundance, San Francisco, pp: 775-780.
- Stedinger JR, Vogel RM, Foufula-Georgiou E (1992) Frequency analysis of extreme events, Handbook of Hydrology.
- Rahman AS, Rahman A, Zaman MA, Haddad K, Ahsan A, et al. (2013) A study on selection of probability distributions for at-site flood frequency analysis in Australia. Natural Hazards 69: 1803-1813.
- Solaiman TA (2011) Uncertainty estimation of extreme precipitations under climate change: A non-parametric approach. Dissertation, Department of Civil and Environmental Engineering, The University of Western Ontario, Canada.
- Woo W, Julien HS, Richardson EV (1988) Hazard classification guidelines, USBR. ACER Technical Memorandum No 11.

Citation: Kafle MR, Shakya NM (2018) Two-Dimensional Hydrodynamic Modelling of Koshi River and Prediction of Inundation Parameters. Hydrol Current Res 9:298. DOI: 10.4172/2157-7587.1000298

Copyright: © 2018 Kafle MR, et al. This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.