Reza BARATI*, Sajjad RAHIMI, Gholam Hossein AKBARI

1. Faculty of Civil and Environmental Engineering, Tarbiat Modares University, Tehran, Iran

2. Social Development and Health Promotion Research Center, Kermanshah University of Medical Sciences,Kermanshah, Iran

3. Department of Civil Engineering, University of Sistan and Baluchestan, Zahedan, Iran

1 Introduction

The simulation of flood events by flood routing methods commonly follows hydraulic and hydrologic approaches. Although both hydraulic and hydrologic approaches use the principle of conservation of mass, the former approach considers dynamic effects of flow through the momentum equation, whereas the latter approach does not and simply regards the volume of water in a channel reach as a single-valued function of discharge with the storage-continuity equation. In general, hydrologic models such as linear and nonlinear Muskingum models need to determine hydrologic parameters using recorded data in both upstream and downstream sections of rivers and/or by applying robust optimization techniques(Barati 2011a, 2011b). On the other hand, hydraulic models such as dynamic and diffusion wave models require the gathering of a lot of data related to river geometry and morphology and consume a lot of computer resources (Samani and Shamsipour 2004). Based on the available field data and goals of a project, one of these approaches is utilized for the simulation of flooding in rivers and channels.

Among hydraulic approaches, dynamic wave models that include the Saint-Venant equations and that can perform well for most rivers have been widely applied to river flood routing (Zhang and Bao 2012). These equations are nonlinear, and have no general analytical solutions. Therefore, numerical methods such as the method of characteristics, finite element method, finite volume method, and finite difference method have been used to solve the unsteady-flow equations considering all of the terms of the momentum equation: the pressure gradient, inertia, gravity, and flow resistance terms (Chow et al. 1988; Zhang 2005;Moghaddam and Firoozi 2011). On the other hand, some terms of the momentum equation have been omitted to simplify the problem in the kinematic wave model (i.e., without the acceleration and pressure terms), noninertia or diffusion wave model (i.e., without local and convective acceleration terms), and quasi-steady dynamic wave model (i.e., without the local acceleration term)(Tsai 2003; Wang et al. 2003; Moramarco et al. 2008). Other simplified flood routing methods such as the Muskingum-Cunge and variable parameter Muskingum models (Ponce and Lugo 2001; Wang et al. 2006; Perumal and Sahoo 2007; Song et al. 2011)have been developed as semi-distributed models. However, when the magnitudes of different terms of the momentum equation are widely varying, the use of these simplified models may not yield accurate simulations for all natural rivers. For example, excluding the inertia and pressure gradient terms of the momentum equation (i.e., the kinematic waves)may lead to significant errors for a channel with a milder bed slope as well as with a larger value of roughness coefficient (Barati 2010).

The sensitivity analysis plays an important role in developing computer models because it allows environmental professionals, regulatory reviewers, and policy-makers to better understand the results obtained by use of such models and consequently to make more effective decisions. Researchers have considered different aspects of the dynamic wave model for flood routing in rivers and channels (Cunge et al. 1980; Venutelli 2002, 2011; Helmi? 2005;Anderson et al. 2006; Zhang 2005; Kuiry et al. 2010; Akbari et al. 2012; Akbari and Barati 2012). In this study, a weighted four-point implicit finite difference scheme was developed with several subroutines in the environment of the MATLAB software for the simulation of flooding in rivers and channels. Then, the developed model was evaluated by examining the accordance of the simulated results of the flood events with real conditions of several natural rivers. Finally, in order to investigate the effects of errors in parameters of the model structure,parameters of the geometries and bed surface of rivers, parameters of catchment features, and forcing data of floods on the modeling results, sensitivity analyses of the input parameters of the dynamic wave model for flood routing in rivers and channels were performed with consideration of the effects of variations of each parameter on the variations of other parameters. The key point of the present research is that the parameters with high values of the sensitivity index in special situations must be attentively marked during the selection of the input parameters of the model.

2 Model development and verification

2.1 Dynamic wave model

The following assumptions are used in the derivation of the governing equations: (1)the pressure distribution is hydrostatic, (2)the velocity is uniformly distributed over a channel section, (3)the average channel bed slope is small, (4)the flow is homogeneous and incompressible, and (5)there is no lateral flow.

Based on the assumptions, the continuity and momentum equations, the Saint-Venant equations, can be respectively expressed as

where Q is the discharge, A is the cross-sectional area of flow, x is the horizontal coordinate along the channel, t is time, g is the acceleration due to gravity, y is the flow depth, S0is the slope of the bottom of the channel, and Sfis the friction slope.

The equations above have two independent variables, x and t, and two dependent variables, the discharge Q and flow cross-sectional area A. When the Manning formula is used to represent the flow resistance, the friction slope is expressed as

If appropriate initial and boundary conditions are prescribed, the numerical solutions of Eqs. (1)and (2)can be obtained. Implicit finite difference schemes have been proven to be more efficient in the numerical treatment of the one-dimensional unsteady flow in rivers with a free surface than other methods such as the explicit and characteristic methods (Cunge et al.1980; Chow et al. 1988; Chaudhry 1993; Venutelli 2002; Anderson et al. 2006). For example,because of the numerical stability characteristics of the finite difference equations,theoretically, the implicit method does not restrict the size of the time step. Larger values of time steps enable the implicit method to be more computationally efficient than other methods,particularly for long-duration floods.

In the weighted four-point scheme, the time and space derivative and non-derivative terms of the Saint-Venant equations are approximated as follows, respectively:

where D is a generic parameter that represents the two dependent variables, Δt is the time step, Δx is the space step, i is the spatial index, j is the temporal index, and θ is the weighting factor that ranges from 0 to 1.0. When θ = 1.0, a fully implicit scheme is formed;when θ = 0.0, a fully explicit scheme is formed; while the value of θ = 0.5 gives a box scheme(Hassan et al. 2009). This four-point implicit method is unconditionally stable when 0.5≤θ≤1.0 (Akan 2006).

By substituting the finite difference approximations above and the coefficients into the equations of gradually varied unsteady flow, and assigning the initial conditions and two boundary conditions, a set of nonlinear algebraic equations are obtained. These equations can be solved using a functional iteration method such as the Newton-Raphson method.

Initial guesses of the two dependent variables in the implicit scheme are approximated as follows:

For subcritical flow, boundary equations for upstream and downstream boundaries are used, while, for supercritical flow, both boundary equations are for the upstream end. More details on the types of boundary conditions are presented in Vreugdenhil (1994). In this study,for subcritical flow, the two boundary conditions required by the model are the inflow discharge hydrograph at the upstream boundary and the stage-discharge curve at the downstream boundary. It is notable that downstream sections used in the computation are distant from the real downstream section of reach to weaken the effects of the downstream boundary condition on results. The values of the two dependent variables at the beginning of the time step are specified at all the nodes along the channel as initial conditions.

2.2 Model verification

In order to compare the field observations and the results simulated by the developed model, the field data from several natural rivers in the Persian Gulf region were utilized. Over ten flood events with single- and multi-peaked hydrographs were investigated (Barati 2010).For brevity, the results of only two flood events are illustrated in Fig. 1. The results of all of the flood events that were simulated by the developed model were satisfactory in terms of attenuation, lag, and mass conservation (Barati 2010). For example, for the flood event shown in Fig. 1(a), the peak outflow discharge of the routed hydrograph generated by the dynamic wave model and the peak outflow discharge of the observed hydrograph are 1 594 and 1 577 m3/s, respectively. These values occur at 24.75 and 26 h, respectively, in terms of the time to peak. Furthermore, the value of the flood volume estimated by numerical integration of Simpson’s rule shows a slightly greater error than the observed value. In general, the results indicate that the developed model has high accuracy and consistent simulation results with the field data with different sets of input variables from several natural rivers. In other words, the results of the developed model are compatible with natural characteristics addressed in this paper.

Fig. 1 Observed inflow and outflow hydrographs and simulated hydrograph generated by dynamic wave model

3 Performance evaluation criteria (PEC)

For the evaluation of the results of the dynamic wave model, the attenuation of the peak outflow (ε)and the lag of the peak outflow (η)are considered. These criteria, which are dimensionless factors, are presented in Eqs. (8)and (9).

where Qpiand Qpoare the peak discharges for the upstream and downstream hydrographs,respectively; and Tpiand Tpoare the time to peaks related to Qpiand Qpo, respectively.

The attenuation refers to the reduction in the peak and dispersion of the flood hydrograph as it propagates, whereas the lag refers to the deferment in time of the peak discharge at downstream points. The attenuation and lag are the main criteria in flood routing because the peak discharge and the time to peak that consider the attenuation and lag criteria, respectively,define the shape of hydrograph (Perumal and Sahoo 2007; Barati 2011b).

4 Numerical experiments

4.1 Experimental runs

In order to investigate the effects of input parameters of the dynamic wave model on output results, about 800 simulating experiments based on different combinations of the channel characteristics (i.e., the bed slope S0and Manning’s roughness coefficient n), the flood and catchment characteristics (i.e., the time to peak Tpand skewness factor γ), and the model characteristics (i.e., the weighting factor θ, time step Δt, and space step Δx)were performed. The simulated channel has a rectangular cross-section with a bottom width B of 50 m and a length L of 30 km. The base and peak discharges are considered to be 100 and 500 m3/s, respectively. The details of the parameters of the numerical experiments are listed in Table 1. The following synthetic inflow hydrograph of the form of the Pearson Type III distribution was used for the upstream boundary condition:

where Q( t)is the discharge hydrograph at the upstream end of the channel reach, Qbis the base flow, and Qpis the peak flow.

Table 1 Values of channel, model, and flow characteristics in experimental runs

4.2 Model sensitivity analyses

In order to investigate the variation of Manning’s roughness coefficient, bed slope,skewness factor, and time to peak with the attenuation and lag, a total of 360 routing experiments were performed. The dynamic wave model was performed using different combinations of the aforementioned parameters (n, S0, γ, and Tp)that presented in Table 1.

4.2.1 Graphical multi-parametric sensitivity analyses (GMPSA)

Sensitivity analysis consists of investigating how the variation in the output of a model can be qualitatively or quantitatively allocated to different sources of variation, and how the outputs of a given model depend upon the information fed into it (Refsgaard et al. 2007;ASCE 2008). The greater the parameter sensitivity, the greater the effect an error in that parameter will have on the computed results (McCuen 2003; McCuen and Knight 2006;USEPA 2003).

The developed procedure for the sensitivity analysis, that is, the graphical multiparametric sensitivity analysis (GMPSA), is illustrated in Fig. 2. The first step is the selection of parameters to be tested. Fundamentally, computer models may include ill-de fi ned parameters that cannot be measured with a high degree of accuracy in the fi eld or in the laboratory and, therefore, will severely influence the accuracy of any single simulation and increase the dif fi culty in assessing the relative importance of parameters. For this purpose, the skewness factor and time to peak are considered hydrograph shape factors. The effects of the variation of the two parameters on inflow hydrographs are presented in Fig. 3.In general, the variations of the hydrograph shape factors represent the variations of the characteristics of catchments, such as the catchment area, catchment shape, river morphology, lithology, and vegetation. For example, the duration of hydrographs in vast catchments is longer (i.e., higher values of the skewness factor and time to peak)than that in small catchments. On the other hand, the characteristics of flood and rainfall events such as the rainfall intensity and rainfall duration are important for hydrograph shapes. For example, the hydrograph shape in a flash flood with a high intensity is more tapered (i.e., with lower values of the skewness factor and time to peak)than that of rainfall events with a lower intensity. The values of Manning’s roughness coefficient always have some degree of uncertainty. Some important factors used for selection a roughness coefficient are: (1)surface irregularities, (2)variations in channel shape and size,(3)amount of vegetation, (4)obstructions, (5)channel meandering and curvature, (6)change of season, (7)temperature, (8)scour and deposition, and (9)channel alignment. Therefore,selecting a value of Manning’s roughness coefficient for a natural river is not easy (Akan 2006;Kim et al. 2010). The channel slope represents the effects of the gravity. Gravity is the driving force in open channel flow with free surface. Therefore, values of the bed slope are essential in the simulation of flooding in rivers and channels. On the other hand, some parameters such as the river length and peak flow that can essentially be as small or as large as an engineer desires based on the goals of projects are not varied in the numerical experiments.

The sensitivities of simulation results to input parameters need to be evaluated by setting the range of parameters based on the variation of the parameters in the real world (Table 1).After each parameter was varied, the attenuation and lag criteria (PEC)were calculated, and the sensitivity index (SI), that is, the relative change of PEC with the change of each input parameter, was calculated. Because the sensitivity index is a dimensionless factor, it can be used to compare the sensitivities of parameters. It is notable that a negative value of the sensitivity index indicates an inverse relationship between the input and output parameters.These steps must be repeated until PEC and the corresponding sensitivity index for all sets of the input parameters are calculated. Then, the effects of variations of input parameters on the results can be investigated through illustration of the variation of the sensitivity index. Finally,for the comparison of the sensitivities of different input parameters, the mean of the absolute sensitivity index (MASI)for each parameter can be calculated.

Fig. 2 Flowchart of procedure of GMPSA

Fig. 3 Effects of skewness factor and time to peak on inflow hydrographs

4.2.2 Results and discussion of sensitivity analysis

In general, for all the parameters that were changed in the numerical experiments, the values of the sensitivity index in terms of the attenuation criterion are larger than those of the sensitivity index in terms of the lag criterion. In other words, most of the effects of the variations of the parameters are on the attenuation criterion rather than on the lag criterion.The results of the sensitivity indices for Manning’s roughness coefficient and bed slope in terms of the attenuation and lag criteria are presented in Figs. 4 through 7. For Manning’s roughness coefficient, the value of the sensitivity index increases when the skewness factor and time to peak increase in terms of both the attenuation (Fig. 4)and lag (Fig. 5)criteria.However, these values, particularly for larger values of the time to peak, are not significantly sensitive to the skewness factor for the steeper bed slope. Furthermore, when the bed slope increases, the sensitivity index increases until a specific value of bed slope appears, and then it decreases. These variations are particularly sensible in terms of the attenuation criterion. In other words, the effects of an error in the estimation of Manning’s roughness coefficient on the output results are more significant in some cases, such as a bed slope of about 0.000 8 under lower values of the time to peak, and a bed slope of about 0.000 4 under larger values of the time to peak. On the other hand, for the bed slope, like Manning’s roughness coefficient, the absolute value of the sensitivity index increases when the skewness factor and time to peak increase in terms of both the attenuation (Fig. 6)and the lag (Fig. 7)criteria. Moreover, when Manning’s roughness coefficient increases, the absolute value of the sensitivity index decreases for the lower value of time to peak while the variations of the sensitivity index do not show a particular trend for the larger values of time to peak. In conclusion, the sensitivity of the results to the variations of the channel characteristics (i.e., Manning’s roughness coefficient and the bed slope)increases when the flood and catchment characteristics (i.e., the skewness factor and time to peak)increase. The physical concepts of these results are, for example, that the values of the sensitivity index of Manning’s roughness coefficient and the bed slope for vast catchments and/or longer-duration rainfalls are larger than those of small catchments and/or abrupt rainfalls. In other words, errors of the estimated Manning’s roughness coefficient and bed slope parameters have more significant effects on the output results in vast catchments and/or for longer-duration rainfalls. Therefore, a reasonable effort should be made to reduce errors for these situations in the estimation of Manning’s roughness coefficient and bed slope when using the dynamic wave model to simulate a flood event.

Fig. 4 Variation of SI of n in terms of ε

Fig. 5 Variation of SI of n in terms of η

The variation of the sensitivity index has no particular trend for both the skewness factor and the time to peak in terms of the lag criterion. The results of the sensitivity indices for the two parameters in terms of the attenuation criterion are presented in Figs. 8 and 9. For the skewness factor, the absolute value of the sensivity index increases when Manning’s roughness coefficient decreases, and/or when the bed slope and time to peak increase. On the other hand, for the time to peak, the value of sensivity index increases when Manning’s roughness coefficient decreases, and/or when the bed slope and skewness factor increase. In brief, the sensitivity of the results to the variations of the hydrograph shape factors (i.e., the skewness factor and time to peak)increases when Manning’s roughness coefficient decreases,and/or the bed slope increases. It can be concluded that errors of a design hydrograph that can be generated using the synthetic unit hydrograph (SUH)methods, such as Snyder’s method,the Taylor and Schwarz (TS)model, the Soil Conservation Service (SCS)method, and Gray’s method in design projects of the real world (Bhunya et al. 2011), have more significant effects on the output results for channels with a steeper bed slope or with a lower roughness coefficient. Furthermore, the skewness factor and time to peak have a direct relationship with one another in terms of the sensitivity of the results (i.e., the sensitivity of the results to the variations of the skewness factor increases with the time to peak and vice versa).

Fig. 6 Variation of SI of S0 in terms of ε

Fig. 7 Variation of SI of S0 in terms of η

Analyses of the sensitivity of the results to the variations of the channel, flood, and catchment characteristics indicate that the rankings of the parameter importance in terms of MASI are as follows: the skewness factor (396.36%), time to peak (131.68%), Manning’s roughness coefficient (111.82%), and bed slope (96.98%)for the attenuation criterion; and the time to peak (79.35 %), Manning’s roughness coefficient (54.74%), bed slope (40.65%), and skewness factor (17.11%)for the lag criterion. On the other hand, if the attenuation and the lag criteria are simultaneously considered, the rankings are as follows: the skewness factor(206.74%), time to peak (105.52%), Manning’s roughness coefficient (83.28%), and bed slope(68.82%). These results indicate that the effects of the variations of the hydrograph shape factors (i.e., the time to peak and skewness factor)on the results are more significant than the effects of the variations of the parameters of the characteristics of the geometries and bed surface of rivers (i.e., Manning’s roughness coefficient and the bed slope).

Fig. 8 Variation of SI of γ in terms of ε

Fig. 9 Variation of SI of Tp in terms of ε

4.3 Model accuracy analyses

The implicit model can be unconditionally stable, but may not be unconditionally convergent due to the changes of the values of the space and time steps. In order to achieve reasonable accuracy, both the space and time step size should be small. On the other hand, the value of the weighting factor also has an effect on the convergence of the model. In order to investigate the effects of the parameters of the model structure (i.e., space step Δx, time step Δt , and weighting factor θ)on output results, a total of 432 numerical experiments were implemented. The variations of the attenuation criterion for the space and time steps and also the weighting factor are listed in Tables 2 through 4, respectively. Corresponding results for the lag criterion are only described briefly.

The results of experiments for the attenuation and lag criteria show that the variations of the space step only have slight effects on the attenuation criterion, whereas these variations are not significant in terms of the lag criterion. When the space step increases, the attenuation criterion increases for channels with a steeper bed slope while the criterion decreases for channels with a milder bed slope. It is notable that the relationship between the attenuation criterion and the space step is a quadratic curve relationship with a high correlation coefficient in most cases. On the other hand, the results indicate that the variations of the time step have significant effects on both the attenuation and the lag criteria. When the time step increases,the attenuation criterion increases as a linear relationship with high correlation coefficients in all cases, whereas the relationship between the lag criterion and the time step does not have a special trend.

Table 2 Variations of attenuation criterion with space step

Although the variations of the weighting factor are not significant in terms of the lag criterion, these variations have significant effects on the attenuation criterion. The relationship between the attenuation criterion and weighting factor is similar to the relationship between the attenuation criterion and time step (i.e., the criterion increases as a linear relationship with high correlation coefficients in all cases as the weighting factor increases).

Table 3 Variations of attenuation criterion with time step

Table 4 Variations of attenuation criterion with weighting factor

5 Conclusions

In this study, the dynamic wave model based on the implicit finite difference scheme was developed for flood routing in rivers. The field application of the model indicated a good agreement between the simulation results and the observed data. For the evaluation of the results, the attenuation and lag of the peak outflow were used as PEC. Then, sensitivity analyses of the input parameters of the model were performed through numerical experiments.These experiments consist of different combinations of the parameters of the model structure,geometries and bed surface of rivers, and different characteristics of catchments and floods.The variations of the sensitivity of input parameters with the variation of other factors as well as the most influential parameters on the output results were investigated. The most important conclusions can be summarized as follows:

(1)When Manning’s roughness coefficient increased, and/or when the bed slope became milder, the time to peak of the output hydrograph increased (i.e., the lag criterion increased)and the peak discharge decreased (i.e., the attenuation criterion increased).

(2)The attenuation criterion decreased as both the skewness factor and time to peak increased, whereas the lag criterion increased when the skewness factor increased, and/or the time to peak decreased.

(3)When the space step increased, the attenuation criterion increased for channels with steeper bed slopes while the criterion decreased for channels with milder bed slopes, whereas the attenuation criterion increased when the time step and/or weighting factor increased.

(4)Of the hydrograph shape factors (i.e., the skewness factor and time to peak)and the channel characteristics (i.e., Manning’s roughness coefficient and the bed slope), the most influential parameter in regard to the attenuation criterion was the skewness factor, whereas the most influential parameter in regard to the lag criterion was the time to peak.

(5)The characteristics of a design hydrograph had significant effects on the results,particularly for lower values of Manning’s roughness coefficient and/or a steeper bed slope.

(6)The effects of the variation of the channel characteristics on the output results were more significant for larger values of the skewness factor and/or time to peak.

(7)Of the model structure parameters (i.e., the space step, time step, and weighting factor), the most influential parameter was the weighting factor of the dynamic wave model.

Akan, A. O. 2006. Open Channel Hydraulics. Oxford: Elsevier.

Akbari, G. H., and Barati, R. 2012. Comprehensive analysis of flooding in unmanaged catchments.Proceedings of the Institution of Civil Engineers-Water Management, 165(4), 229-238. [doi:10.1680/wama.10.00036]

Akbari, G. H., Nezhad, A. H., and Barati, R. 2012. Developing a model for analysis of uncertainties in prediction of floods. Journal of Advanced Research, 3(1), 73-79. [doi:10.1016/j.jare.2011.04.004]

American Society of Civil Engineers (ASCE). 2008. Comments of the American Society of Civil Engineers on the Proposed Revisions to the Economic and Environmental Principles and Guidelines for Water and Related Land Resources Implementation Studies. Washington, D.C.: ASCE.

Anderson, B. G., Rutherfurd, I. D., and Western, A. W. 2006. An analysis of the influence of riparian vegetation on the propagation of flood waves. Environmental Modelling and Software, 21(9), 1290-1296.[doi:10.1016/j.envsoft.2005.04.027]

Barati, R. 2010. Investigation of Flood Routing Methods in Natural Waterways. M. S. Dissertation. Zahedan:University of Sistan and Baluchestan.

Barati, R. 2011a. Discussion of “Parameter Estimation for Nonlinear Muskingum Model based on Immune Clonal Algorithm”. Journal of Hydrologic Engineering, 16(4), 391-393. [doi:10.1061/(ASCE)HE.1943-5584.0000311]

Barati, R. 2011b. Parameter estimation of nonlinear Muskingum models using Nelder-Mead Simplex algorithm. Journal of Hydrologic Engineering, 16(11), 946-954. [doi:10.1061/(ASCE)HE.1943-5584.0000379]

Bhunya, P. K., Panda, S. N., and Goel, M. K. 2011. Synthetic unit hydrograph methods: A critical review. The Open Hydrology Journal, 5, 1-8.

Chaudhry, M. H. 1993. Open-Channel Flow. Englewood Cliffs: Prentice Hall.

Chow, V. T., Maidment, D. R., and Mays, L. W. 1988. Applied Hydrology. Singapore: McGraw-Hill Science.

Cunge, J. A., Holly, F. M., and Verwey, A. 1980. Practical Aspects of Computational River Hydraulics.London: Pitman Publishing Ltd.

Hassan, A., Norio, T., and Nobuyuki, T. 2009. Distributed water balance with river dynamic-diffusive flow routing model. Journal of Hydrodynamics, 21(4), 564-572. [doi:10.1016/S1001-6058(08)60185-7]

Helmi?, T. 2005. Unsteady 1D flow model of a river with partly vegetated floodplains: Application to the Rhine River. Environmental Modelling and Software, 20(3), 361-375. [doi:10.1016/j.envsoft.2004.02.001]

Kim, J. S., Lee, C. J., Kim, W., and Kim, Y. J. 2010. Roughness coefficient and its uncertainty in gravel-bed river. Water Science and Engineering, 3(2), 217-232. [doi:10.3882/j.issn.1674-2370.2010.02.010]

Kuiry, S. N., Sen, D., and Bates, P. D. 2010. Coupled 1D-Quasi-2D flood inundation model with unstructured grids. Journal of Hydraulic Engineering, 136(8), 493-506. [doi:10.1061/(ASCE)HY.1943-7900.0000211]McCuen, R. H. 2003. Modeling Hydrologic Change: Statistical Methods. Boca Raton: CRC Press.

McCuen, R. H., and Knight, Z. 2006. Fuzzy analysis of slope-area discharge estimates. Journal of Irrigation and Drainage Engineering, 132(1), 64-69. [doi:10.1061/(ASCE)0733-9437(2009)135:2(149)]

Moghaddam, M. A., and Firoozi, B. 2011. Development of dynamic flood wave routing in natural rivers through implicit numerical method. American Journal of Scientific Research, (14), 6-17.

Moramarco, T., Pandolfo, C., and Singh, V. P. 2008. Accuracy of kinematic wave approximation for flood routing, II: Unsteady analysis. Journal of Hydrologic Engineering, 13(11), 1089-1096. [doi:10.1061/(ASCE)1084-0699(2008)13:11(1089)]

Perumal, M., and Sahoo, B. 2007. Volume conservation controversy of the variable parameter Muskingum-Cunge method. Journal of Hydraulic Engineering, 134(4), 475-485. [doi:10.1061/(ASCE)0733-9429(2008)134:4(475)]

Ponce, V. M., and Lugo, A. 2001. Modeling looped ratings in Muskingum-Cunge routing. Journal of Hydrologic Engineering, 6(2), 119-124. [doi:10.1061/(ASCE)1084-0699(2001)6:2(119)]

Refsgaard, J. C., van der Sluijs, J. P., H?jberg, A. L., and Vanrolleghem, P. A. 2007. Uncertainty in the environmental modelling process: A framework and guidance. Environmental Modelling and Software,22(11), 1543-1556. [doi:10.1016/j.envsoft.2007.02.004]

Samani, H. M. V., and Shamsipour, G. A. 2004. Hydrologic flood routing in branched river systems via nonlinear optimization. Journal of Hydraulic Research, 42(1), 55-59. [doi:10.1080/00221686.2004.9641183]

Song, X. M., Kong, F. Z., and Zhu, Z. X. 2011. Application of Muskingum routing method with variable parameters in ungauged basin. Water Science and Engineering, 4(1), 1-12. [doi:10.3882/j.issn.1674-2370.2011.01.001]

Tsai, C. W. 2003. Applicability of kinematic, noninertia, and quasi-steady dynamic wave models to unsteady flow routing. Journal of Hydraulic Engineering, 129(8), 613-627. [doi:10.1061/(ASCE)0733-9429(2003)129:8(613)]

U.S. Environmental Protection Agency (USEPA). 2003. Draft Guidance on the Development, Evaluation, and Applications of Regulatory Environmental Models. Washington, D.C.: The Council for Regulatory Environmental Modeling.

Venutelli, M. 2002. Stability and accuracy of weighted four-point implicit finite difference schemes for open channel flow. Journal of Hydraulic Engineering, 128(3), 281-288. [doi:10.1061/(ASCE)0733-9429(2002)128:3(281)]

Venutelli, M. 2011. Analysis of dynamic wave model for unsteady flow in an open channel. Journal of Hydraulic Engineering, 137(9), 1072-1078. [doi:10.1061/(ASCE)HY.1943-7900.0000405)]

Vreugdenhil, C. B. 1994. Numerical Methods for Shallow-Water Flow. Dordrecht: Kluwer Academic Publishers.

Wang, G. T., Chen, S., Boll, J., and Singh, V. P. 2003. Nonlinear convection-diffusion equation with mixing-cell method for channel flood routing. Journal of Hydrologic Engineering, 8(5), 259-265.[doi:10.1061/(ASCE)1084-0699(2003)8:5(259)]

Wang, G. T., Yao, C., Okoren, C., and Chen, S. 2006. 4-Point FDF of Muskingum method based on the complete St Venant equations. Journal of Hydrology, 324(1-4), 339-349. [doi:10.1016/j.jhydrol.2005.10.010]

Zhang, X. Q., and Bao, W. M. 2012. Modified Saint-Venant equations for flow simulation in tidal rivers.Water Science and Engineering, 5(1), 34-45. [doi:10.3882/j.issn.1674-2370.2012.01.004]

Zhang, Y. 2005. Simulation of open channel network flows using finite element approach. Communications in Nonlinear Science and Numerical Simulation, 10(5), 467-478. [doi:10.1016/j.cnsns.2003.12.006]