The Modified One-Dimensional Hydrodynamic Model Based on the Extended Chezy Formula

Abstract

Although steady uniform friction formulas have been introduced to the framework of a one-dimensional (1D) hydrodynamic model for centuries, the error of friction calculation inevitably undermines the performance of flood routing. Based on successful results of unsteady channel friction research studies, a newly proposed unsteady friction model is introduced to establish a modified 1D hydrodynamic model (namely, the modified SVN model). With the help of a carefully designed parameter calibration method, the performance of the modified SVN model was compared with that of the original SVN model in a simulation test for a hydraulic experiment. This study’s results revealed that compared with the original SVN model, the modified SVN model could achieve a better simulation in both the flow depth and the sectional averaged velocity simulations. Furthermore, it could reduce the peak value error and the time-at-peak error as well, indicating that the use of an unsteady friction model can efficiently improve the performance of a 1D hydrodynamic model.

1. Introduction

Based on Saint-Venant equations, the one-dimensional (1D) hydrodynamic model is a fundamental tool for the simulation of 1D flow movement in the case of river or open channel flows [1,2]. It has been widely used in the fields of hydraulic simulation like flood forecasting and waterway engineering [3,4]. The partial differential equations for rivers and channels can help us better understand the unsteady flow behavior and prevent the potential damage of flood events. However, the performance of the 1D hydrodynamic model critically depends on the frictional calculation accuracy. In the case of unsteady flow, flow unsteadiness and flow non-uniformity make the frictional resistance become more complex. As a result, friction factors (e.g., the Chezy coefficient or the Manning roughness factor) become very hard to identify and the error of friction calculation inevitably undermines the performance of flood routing [5].

Two types of approaches have been developed to address this problem. First, many advanced optimal algorithms have been applied to calibrate the parameter of friction term. For example, the influence coefficient algorithm was introduced by Becker and Yeh for the parameter identification problem in unsteady open channel flow [6]. The genetic algorithm (GA) was introduced by Giustolisi to determine the Chezy coefficient in corrugated channels [7]. A numerical method based on the adjoint equation algorithm was proposed by Ding and Wang to identify Manning’s n in river networks [8]. Further information about identifying friction factors can be found in Liu [9], Nocedal and Wright [10], Ding et al. [11], Zou et al. [12], etc. Despite the fact that parameter calibration algorithms take advantage of optimal control theories, it seems that a calibrating friction factor as a constant is not enough to compensate all the discrepancies between model calculations and observations [13]. Second, many real-time parameter updating methods have been used to dynamically update the friction factor in the engineering field. For example, Hsu et al. proposed an iteration Manning roughness factor updating technique for real-time flood forecasting [14]. Bao and Zhao applied the Kalman filter to deal with the roughness coefficient updating problem [15]. In addition, Liu et al. presented a multi-segmental roughness coefficient updating method to consider both spatial distribution and a time-varying process of roughness coefficient in the channel from Cuntan to Wanxian in Three Gorges of the Yangtze River [16]. These dynamic routing models with real-time roughness updating methods give good results for hydraulic simulations. However, because the updated friction factors of the present time step are usually quite different from those of the previous time step, the results deteriorate markedly when the updated parameter is used for hydraulic forecasting [17].

Neither carefully designed parameter calibration algorithms nor dynamically updating methods can significantly improve the accuracy of friction calculation, indicating that the classic friction models like the Chezy formula and the Manning formula only work for steady uniform flows but not for unsteady flows. If one neglects the impact of flow unsteadiness on frictional resistance, the friction coefficients become grey-box parameters [18]. Rouse [19], Yen [20], and Ding et al. [11] argued that in the case of unsteady flow, friction parameters such as the Chezy coefficient and the Manning roughness factor actually include the components of surface friction resistance, form resistance, wave resistance, and resistance due to flow unsteadiness. As a result, these coefficients are very unstable and unpredictable in an unsteady flow regime [21]. Many studies in the literature are devoted to developing new models for calculating 1D unsteady channel friction, like Tu and Graf [21], Rowiński et al. [22], Mrokowska et al. [23], and Bao et al. [24]. The proposed unsteady friction models contain the terms of time and space derivatives of hydraulic parameters, which better demonstrate the dependence of unsteady friction on flow unsteadiness and flow non-uniformity. Nevertheless, on one hand, most studies considered the problem of unsteady friction as the inverse problem of flood routing. Hence, the unsteady friction formulas proposed in them are derived from the momentum equation of Saint-Venant equations. These equations cannot be used to calculate the friction term of the 1D hydrodynamic model. On the other hand, although some of the proposed formulas can be used to compute the friction slope in the 1D hydrodynamic model, these formulas are more complex than the Chezy and Manning formulas. For example, Bao et al. proposed a 1D unsteady friction formula based on dimensional analysis instead of the inverse analysis of momentum equation [24]. However, in their proposal, more empirical parameters were introduced to describe the dependence of unsteady friction on spatio-temporal derivatives of hydraulic parameters. It is known that the determination of empirical parameters is very difficult due to the high nonlinearity of the 1D hydrodynamic model. Hence, the estimation of friction parameters would be more difficult and time-consuming when the proposed unsteady friction model is used as a substitute for the Chezy or Manning formula. It is evident that the advantages of using these unsteady friction models in the 1D hydrodynamic model have not yet been discussed.

To uncover the potential advantages of unsteady friction formulas, a modified 1D hydrodynamic model is proposed in this paper, referred to as the modified SVN model, based on the unsteady friction formula presented by Bao et al. [24]. Since the steady uniform component of the unsteady friction model is the same as the Chezy formula, a classic 1D hydrodynamic model with the friction calculator of the Chezy formula is used for comparison, referred to as the original SVN model. The performances of the modified and the original SVN models are compared using a dataset of hydraulic experiments, where a genetic algorithm-based parameter calibration method is used to overcome the difficulties of empirical parameter estimation.

2. Flume Experiments and Data

The hydraulic experiments were conducted in a 20.0 m long, 0.3 m wide, and 0.4 m high glass-made flume with a recirculating water system in the hydraulic laboratory of Hydraulic and Hydro-Power Engineering, Hohai University. The experiment setup is shown in Figure 1. The water recirculation system consists of the glass-made flume, a recycle pond, a large water tank, and two pumps. The pump between the flume and the large water tank was used to pump water into the flume. In order to relieve the vibration to the flume, a flexible tube was used for connecting the water pump and the flume. A small water tank was placed ahead of the flume to straighten the incoming flow and eliminate large eddies, which made the flow state across the experimental range more stable and without significant ripples. The used water passed through the tail gate and then entered the cycle pond. It entered the large tank through another pump for reutilization. In addition, gravel was paved evenly on the bottom to simulate a rough riverbed [25]. The flume bed remained flat throughout the experiments. Three observation cross sections were placed along the flume at an interval of 5 m, referred to as cross sections #1, #2, and #3. The incoming flow was controlled using a hydraulic valve and its discharge was measured using a pipe ultrasonic flow meter with a precision of 0.5% [26].

In order to obtain the cross-sectional averaged velocity, a NORTEK Acoustic Doppler Velocimeter (ADV) (Nortek, Vangkroken 2, 1351 Rud, Norway) was used to measure instantaneously the vertical distributions of point velocities at each cross section [13]. The sampling rate and number of the ADV were set to 50 Hz and 5000, respectively, according to Niño et al. [27]. The relationships of the cross-sectional averaged velocity versus a fixed-point velocity were carefully calibrated at cross sections #1, #2, and #3, referred to as V−V* relationships. They are expressed using Equations (1), (2), and (3), respectively:

$$V=0.8373\times V^{*}+0.1211$$

(1)

$$V=0.7274\times V^{*}+0.2185$$

(2)

$$V=0.7405\times V^{*}+0.2832$$

(3)

where V* is the averaged velocity over a time which is detected by the ADV at a fixed point and V is the cross-sectional averaged velocity. Fitting curves are displayed in Figure 2, where the coefficients of determination are 0.9990, 0.9695, and 0.9937, respectively.

In the unsteady flow experiments, the discharge of incoming flow which was controlled by the hydraulic valve gradually rose from the base value, reached its peak value, and then fell back to the original one, simulating flood events of single-peak hydrographs. In all, 10 different man-made flood events were conducted. Each flood event was performed within 3–4 min. The discharge of incoming flow in the No. 1 experiments (1-1, 1-2, 1-3, 1-4, 1-5) was in the range of 60–110 (m³/h), and the discharge in the No. 2 experiments (2-1, 2-2, 2-3, 2-4, 2-5) was in the range of 60–100 (m³/h). During the unsteady flow experiments, both the cross-sectional averaged velocity V and the flow depth D were recorded every 5 s and the observation precision was evaluated using a posteriori error inspection method. According to the principle of conservation of mass, variables at each cross section should be compatible with following water balance equation:

$${\displaystyle \sum _{t=0}^T{Q}_{in}\left(t\right)\cdot \Delta t}=B\cdot {\displaystyle \sum _{t=0}^{T}D\left(t\right)\cdot V\left(t\right)\cdot \Delta t}+B\cdot L\cdot [\overline{D}\left(T\right)-\overline{D}\left(0\right)]$$

(4)

where the left side equation is the accumulative discharge of inflow (i.e., ∑Q_in), B is the width of the flume, L is the length of the part ahead of the cross section, D(t) is the flow depth series, V(t) is the average velocity series, $\overline{D}$(T) is the average stage of the part ahead of the cross section at the end of experiment, and $\overline{D}$(0) is average stage of the part ahead of the cross section at the beginning of the experiment. According to Equation (4), the ∑Q_in can be calculated based on the time series of D(t) and V(t), referred to as calculated accumulative discharge. Moreover, ∑Q_in was measured using the pipe ultrasonic flow meter as well. Therefore, the relative error of ∑Q_in in each experimental flood event can be obtained, as shown in

Table 1. Relative error between the observed and calculated accumulative discharges at each cross section in man-made flood events.

No.	#1	#2	#3	No.	#1	#2	#3
1-1	0.1%	1.4%	1.4%	2-1	0.4%	−0.5%	1.9%
1-2	1.2%	0.5%	1.2%	2-2	0.3%	1.6%	0.1%
1-3	−2.1%	−3.1%	−1.0%	2-3	0.0%	0.8%	−1.0%
1-4	−2.2%	−1.2%	−1.1%	2-4	−0.3%	−0.4%	1.0%
1-5	−1.3%	−0.8%	−2.7%	2-5	−0.3%	0.6%	1.3%

.

Table 1 shows that the relative error between the calculated accumulative discharge and the observed one is less than 4%, mostly controlled within 2% (86.7%). This illustrates that the observations of D(t) and V(t) at each cross section are reasonable. With the dataset of hydraulic experiments, the characteristic of unsteady flow is simply illustrated. This study’s experiments showed that the wave crests passed through cross sections #1, #2, and #3 in sequence. As an example, the flow depth hydrographs of the No. 1-1 flood event are shown in Figure 3.

In addition, the relationships between the cross-sectional averaged velocity and the flow depth displayed an anticlockwise loop-shaped curve in the experiments (see Figure 4).

3. Methods

3.1. Extended Chezy Formula

In engineering practice, empirical steady uniform friction formulas are widely used to estimate the unsteady friction slope of open channels and rivers [28]. For example, the Chezy formula can be written as:

$$S_c=\frac{V^2}{C^{2}R}$$

(5)

where S_c is the friction slope computed with the Chezy formula, V is the cross-sectional averaged velocity, R is the hydraulic radius, and C is the Chezy coefficient. The relationships among Chezy’s C, Manning’s n, and the Fanning friction factor f can be written as follows:

$$C=\frac{1}{n}{R}^{1/6}=\sqrt{\frac{2g}{f}}$$

(6)

where n is the Manning roughness coefficient, f is the Fanning friction factor, and g is the gravitational acceleration.

However, the frictional resistance of unsteady flow is more complex than that of steady uniform flow. A symbolic function which was deduced by Rouse [19] reported that the Fanning friction factor f is associated with several dimensionless hydraulic parameters, as follows:

$$f=F\left(\mathrm{Re},\mathrm{K},\mathsf{\eta },\mathrm{N},\mathrm{Fr},\mathrm{U}\right)$$

(7)

where Re is the Reynolds number, K is the relative roughness, η is a dimensionless parameter that describes the effect of the cross-sectional geometric shape, N is the dimensionless flow non-uniformity factor, Fr is Froude number, and U is the dimensionless flow unsteadiness factor.

Equation (7) is a general equation for both steady uniform and unsteady frictional resistances. In Equation (7), some dimensionless parameters such as Re, Fr, K, and η are also the factors which affect the frictional resistance of steady uniform flow. On the other hand, the other parameters like U and N which represent the impacts of flow unsteadiness and flow non-uniformity are the additional factors that emerge in the case of unsteady flow. Therefore, they can be reorganized as the following expression:

$$f=F_1\left(\mathrm{Re},\mathrm{Fr},\mathrm{K},\mathsf{\eta }\right)\times F_2\left(\mathrm{U},\mathrm{N}\right)$$

(8)

where f is the Fanning friction factor, F₁ is a function of the factors which affect both steady uniform and unsteady frictions, and F₂ contains the factors which only affect unsteady friction, which depends on time and space derivatives of hydraulic parameters, such as ∂D/∂x, ∂D/∂t, ∂V/∂x, and ∂V/∂t. Based on dimensional analysis, a modified friction model for unsteady flow has been proposed by Bao et al. [24]. The modified model states the following:

$$S_{\mathrm{m}}=c_0\frac{V^2}{gR}+c_1\frac{V^2}{gR}\frac{\partial D}{\partial x}+c_2\frac{V}{gR}\frac{\partial D}{\partial t}+c_3\frac{1}{g}\frac{\partial V}{\partial t}$$

(9)

where D is the flow depth, V is the cross-sectional averaged velocity, and S_m is the friction slope calculated by the modified friction model, in which the steady uniform component of the modified model is the same as the Chezy formula. Moreover, it considers the effect of flow non-uniformity and flow unsteadiness by additional components related to ∂D/∂x, ∂D/∂t, and ∂V/∂t. Thus, the modified friction model can be regarded as an extension of the Chezy formula, namely, the extended Chezy formula. Since it is not deduced from the momentum equation of the Saint-Venant equations, the extended Chezy formula can be recognized as a substitute for the Chezy formula or Manning formula to modify the friction calculation of a 1D hydrodynamic model.

3.2. Modified SVN Model

Saint-Venant equations as one-dimensional shallow water equations are derived from the conservation of mass and the dynamic relation for momentum conservation principles, which can be applied to unsteady flow appropriately [29]. Its expressions can be written as follows:

$$D\frac{\partial V}{\partial x}+V\frac{\partial D}{\partial x}+\frac{\partial D}{\partial t}=0$$

(10)

$$\frac{\partial D}{\partial x}-S_0+\frac{1}{g}\frac{\partial V}{\partial t}+\frac{V}{g}\frac{\partial V}{\partial x}+S_f=0$$

(11)

where D is the flow depth, V is the cross-sectional averaged velocity, S₀ is the bed slope, S_f is the friction slope, and g is the gravitational acceleration. In this paper, the extended Chezy formula is used to compute the friction term in the momentum equation of the Saint-Venant equations. That is, the S_f is computed using Equation (9). Therefore, a new version of a 1D hydrodynamic model is proposed, which is referred to as the modified SVN model.

The numerical solution of the modified SVN model can be obtained at given initial and boundary conditions. In this paper, the Preissmann implicit difference method was used for discretization, in which the arithmetic average method and the weighted average method are applied in spatial partial derivative approximation and temporal partial derivative approximation, respectively [30,31,32]:

$$f\left(x,t\right)=\frac{\theta }{2}\left(\Delta f_{j+1}+\Delta f_j\right)+\frac{1}{2}\left(f_{j+1}^n+f_j^n\right)$$

(12)

$$\frac{\partial f}{\partial x}=\theta \frac{\Delta f_{j+1}-\Delta f_j}{\Delta x}+\frac{f_{j+1}^n-f_j^n}{\Delta x}$$

(13)

$$\frac{\partial f}{\partial t}=\frac{\Delta f_{j+1}+\Delta f_j}{2\Delta t}$$

(14)

where f represents a certain variable such as flow depth and cross-sectional averaged velocity. θ is the difference coefficient. Its range is 0.5 ≤ θ ≤ 1 [33]. The superscripts represent time t. n and n + 1 are the current moment and the next moment, respectively. The subscripts represent location x. j and j + 1 are the two adjoining cross sections, in which jth cross section locates upstream of the j + 1th cross section. $\Delta f_j=f_j^{n+1}-f_j^n$ and $\Delta f_{j+1}=f_{j+1}^{n+1}-f_{j+1}^n$.

Friction is of great importance in hydrodynamic modelling. Since it is one of the major sources of uncertainties in flood routing, every finding on friction formula modification may help modelers to cope with this uncertainty. Hence, the extended Chezy formula is employed to compute the friction slope in Equation (16), establishing a modified SVN model. In the process of discretization, four hydraulic variables of the two adjacent cross sections are organized into two equations [14,34]:

$$p{a}_{1j}\times \Delta V_j+p{b}_{1j}\times \Delta D_j+p{k}_{1j}\times \Delta V_{j+1}+p{d}_{1j}\times \Delta D_{j+1}=p{e}_{1j}$$

(15)

$$p{a}_{2j}\times \Delta V_j+p{b}_{2j}\times \Delta D_j+p{k}_{2j}\times \Delta V_{j+1}+p{d}_{2j}\times \Delta D_{j+1}=p{e}_{2j}$$

(16)

where Equations (15) and (16) are the two discretized equations derived from the continuity equation and the momentum equation between the jth and the j + 1th cross sections, respectively. Subscript j and j + 1 represent the index of sections. pa_1j, pb_1j, pk_1j, pd_1j, pe_1j, pa_2j, pb_2j, pk_2j, pd_2j, and pe_2j are coefficients. They are determined by the current state of the hydraulic parameters and the four parameters involved in the extended Chezy formula, (pa_1j… pe_1j, pa_2j… pe_2j) = f (q, V_j, D_j, V_j+1, D_j+1, c₀, c₁, c₂, c₃). In addition, ∆V_j, ∆D_j, ∆V_j+1, and ∆D_j+1 are the increment of the next time relative to the current moment.

3.3. Original SVN Model

It is possible to substitute the Chezy formula into the Saint-Venant equations, set 1/C² = c₀, and then obtain the classic 1D hydrodynamic model. In this paper, it is referred to as the original SVN model. Similarly, the numerical solution of the original SVN model can be derived using the Preissmann differencing method to scatter its governing equations. Obviously, the discretized equation of the continuity equation is the same as Equation (15). On the other hand, the discretized equation of the momentum equation is different from Equation (16) due to the different friction formula, which can be expressed as follows:

$$p{a}_{3j}\times \Delta V_j+p{b}_{3j}\times \Delta D_j+p{k}_{3j}\times \Delta V_{j+1}+p{d}_{3j}\times \Delta D_{j+1}=p{e}_{3j}$$

(17)

where, pa_3j, pb_3j, pk_3j, pd_3j, and pe_3j are coefficients. Their calculation formulas are the same as pa_2j, pb_2j, pk_2j, pd_2j, and pe_2j in the case of c₁ = c₂ = c₃ = 0.

For a river or open channel which has N cross sections, both the original and modified SVN models are (2N − 2) equations with 2N unknown variables. Therefore, at the given initial conditions, two boundary conditions are needed for solving equations. In this study, the initial conditions were obtained by the linear interpolation method using the observations of the upstream and downstream cross sections (i.e., cross sections #1 and #3). In addition, the cross-sectional averaged velocity of cross section #1 and the flow depth of cross section #3 served as boundary conditions.

3.4. A Genetic Algorithm (GA)-Based Parameter Calibration Method

It is known that the determination of empirical parameters is very difficult due to the high nonlinearity of 1D hydrodynamic models. Therefore, the estimation of friction parameters would be more difficult and time-consuming since the extended Chezy friction model in the modified SVN model has more empirical parameters that must be calibrated compared with the Chezy formula. In order to overcome the difficulty of parameter estimation, a GA-based parameter calibration method is proposed in this paper. GA is an adaptive optimized algorithm of probability search which simulates heredity and evolution of biology in the wild. According to the natural evolution rules of “Survival of the fittest” and genetic choice theory, the algorithm begins from a random population and evolves to an increasingly good approximation. GA-based methods have been widely used in the parameter calibration of complex and nonlinear problems [35]. For example, Lei et al. [36] and Tang et al. [37] successfully calibrated numerical models using a GA algorithm in complex river networks.

For both the original and modified SVN models, when they are used to simulate man-made flood events, there must be an optimal group of parameters which make the model calculation best fit the experimental observation. Hence, the calibration problem is recognized as an optimization problem. In this study, based on the given initial and boundary conditions, the cross-sectional averaged velocity and the flow depth of cross section #2 were calculated. The objective function can be written as follows:

$$\min E\left(c_0,\dots ,c_k\right)={\displaystyle \sum _{j=1}^m\left[{\alpha }_j\times {\displaystyle \sum _{i=1}^n{\left({D^{\prime }}_{2,i}-D_{2,i}\right)}^2+\left(1-{\alpha }_j\right)}\times {\displaystyle \sum _{i=1}^n{\left({V^{\prime }}_{2,i}-V_{2,i}\right)}^2}\right]}$$

(18)

$$c_k^a\le c_k\le c_k^b\left(k=0,1,\dots ,p-1\right)$$

(19)

where m is the number of experimental flood events (for this study, m = 10 and all 10 flood events were used for parameter calibrations); n is the sequence length of the jth experimental flood event; V_2,i and D_2,i are the observed cross-sectional averaged velocity and flow depth at cross section #2 during the jth experiment, whereas $V_{2,i}^{\prime }$ and $D_{2,i}^{\prime }$ are the computed ones; and α_j is a weighing coefficient between the cross-sectional averaged velocity and flow depth, which is determined by the relative magnitude of the cross-sectional averaged velocity and flow depth [37]. It is given by Equation (20). c_k is the kth parameter to be calibrated and $c_k^a$ and $c_k^b$ are the lower bound and the upper bound, respectively. The frictional slope of the modified SVN model is computed using the extended Chezy formula and four parameters must be calibrated (p = 4), whereas the friction slope of the original SVN model is computed using the Chezy formula and only one parameter must be calibrated (p = 1):

$${\alpha }_j={\displaystyle \sum _{i=1}^n{V}_{2,i}^2}/\left({\displaystyle \sum _{i=1}^n{D}_{2,i}^2}+{\displaystyle \sum _{i=1}^n{V}_{2,i}^2}\right)$$

(20)

where V_2,i and D_2,i are the observed cross-sectional averaged velocity and flow depth at cross section #2.

In this paper, the GA-based parameter calibration method was designed based on the University of Sheffield’s Genetic Algorithm Toolbox [38]. The binary system coding mode was adopted to encode parameters during the parameter calibration. The length of coding that corresponded to each of the parameters was 20. As there were four parameters in the modified SVN model, the length of the individual was 80 (L = 80), whereas in the original SVN model, there was only one parameter and the length of the individual was 20 (L = 20). The process of applying a genetic algorithm for parameter calibration is summarized as follows, with the parameters of GA presented in

Table 2. Parameters involved in the genetic algorithm (GA)-based parameter calibration method.

Parameters of GA	The Original SVN Model	The Modified SVN Model
L	20	80
N	40	40
Px	0.7	0.7
Pm	0.01	0.01
GGAP	0.95	0.95
M	20	80

.

• Generating the initial community: at the beginning of the algorithm, a population of size N is randomly generated and enters the loop.
• Calculating the fitness of parent: after decoding the genetic code of each individual, the parameters are obtained. Then, flood events are simulated using the modified or original SVN model with the parameters to get the value of the objective function. Finally, the mapped method is adopted to calculate the fitness according to the objective function result.
• Generating offspring by crossing and mutation operation: genes are selected by using the roulette method in which the selection probability is proportional to the fitness of parent individuals. In addition, a certain size of offspring population is generated by crossing and mutation, where the simplex crossover operator is used to perform the crossing operation, with the crossover probability being P_x, and the mutation operator is used based on the discrete mutation operator, with the mutation probability being P_m. The ratio of the offspring population to parent population is GGAP, where GGAP = 0.95. Hence, the number of the offspring population is GGAP × N.
• Competition restructuring operation between parent and offspring: according to the principle of survival of the fittest, GGAP × N low fitness parent individuals will be eliminated by offspring ones, that is, 5% excellent individuals of parent and all of offspring will be recombined to form a new population.
• Judging whether the termination conditions are satisfied: if the hereditary algebra is greater than the maximum number of generations M, choose the best individual decoding output in the current population as the optimal approximation solution of the original problem; otherwise, substitute the current population as a parent into the next round of cycle operation, and return to step (2).

The technical flowchart is shown in Figure 5.

3.5. Efficiency Criteria

The Nash–Sutcliffe model efficiency coefficient (NSE), the relative peak value error (RPE) and the relative time-at-peak error (RTE) are used to assess the model accuracy [34,39]. The formulas are as follows:

$$\mathrm{NSE}=1-{\displaystyle \sum _{i=1}^n{\left[X_{\mathrm{c}}\left(i\right)-X_{\mathrm{o}}\left(i\right)\right]}^2}/{\displaystyle \sum _{i=1}^n{\left[X_{\mathrm{c}}\left(i\right)-\overline{X_{\mathrm{o}}}\right]}^2}$$

(21)

$$\mathrm{RPE}=\left(P{X}_{\mathrm{c}}-P{X}_{\mathrm{o}}\right)/P{X}_{\mathrm{o}}\times 100\%$$

(22)

$$\mathrm{RTE}=\left(P{T}_{\mathrm{c}}-P{T}_{\mathrm{o}}\right)/T\times 100\%$$

(23)

where X_c(i) and X_o(i) represent the ith calculated value and the observed value, respectively, $\overline{X_{\mathrm{o}}}$ is the average observed value, and n is the length of X series. NSE indicates how well the observed values match the simulated ones and can range from −∞ to 1; larger NSE values indicate a better model performance. PX_c is the calculated peak value, PX_o is the observed peak value, and RPE > 0 indicates that the calculated peak value is greater than the observed one, and vice versa. PT_c is the time of calculated peak occurrence, PT_o is the time of observed peak occurrence, RTE is the relative time-at-peak error, and RTE > 0 indicates that the calculated peak lags behind the observed peak, and vice versa.

4. Results and Discussion

4.1. Model Calibration

The parameters of both the original and modified SVN models were carefully identified using the proposed GA-based parameter calibration method. Figure 6 illustrates the relationship between the best objective function value and the generation during the parameter calibration for both the original and modified SVN models. It can be seen that the best objective function value decreases with the increasing level of generation. Moreover, this gradual decline trend of the best objective function value ends with a stable state. Therefore, it can be concluded that the proposed parameter calibration method is convergent.

The optimal values of the empirical parameters are listed in

Table 3. Parameter calibration details for the modified and original models.

Empirical Parameter	The Original SVN Model			The Modified SVN Model
	Value	Lower Bound	Upper Bound	Value	Lower Bound	Upper Bound
c0	0.00130	0	1	0.00977	0	1
c1	-	-	-	3.42	3	4
c2	-	-	-	1.12	1	2
c3	-	-	-	−0.724	−1	0
E(c0,…,ck)	0.0196	-	-	0.0042	-	-

, where the lower and upper bounds of c₀ are given after using the Chezy formula. However, the lower and upper bounds of other empirical parameters in the extended Chezy formula could not be estimated due to a lack of similar applications. In this paper, they were determined using the trial and error method.

Moreover, because this paper compared the performance of the original and modified SVN models according to the numerical experiment, the accuracy of the solution was crucial. Mass balance is a good measure for the quality of the computation. In this paper, the relative error (REM) was computed to evaluate the conservation of mass for the measuring segment from cross section #1 to cross section #3. REM is defined as the ratio of the non-conservative discharge and inflow discharge, whereas the concept non-conservative discharge consists of the net outflow and the rate of change in storage [40]. This study’s results showed that REM ranged from 0.3‰ to 0.8‰ in the original SVN model and from 0.4‰ to 0.8‰ in the modified SVN model. This indicated that the non-conservative mass error introduced from the numerical methods was negligible.

4.2. Performance of the Original SVN Model

The optimal parameter values were then inserted into the numerical model (i.e., the original and modified SVN models), and the man-made flood events in the experiment were simulated. The model performance was evaluated by computing the NSE, RPE and RTE between the model calculations and observations at cross section #2.

Table 4. The Nash–Sutcliffe model efficiency coefficients (NSEs), relative peak value errors (RPEs), and relative time-at-peak errors (RTEs) of the original SVN model.

No.	Flow Depth			Sectional Averaged Velocity
	NSE	RPE (%)	RTE (%)	NSE	RPE (%)	RTE (%)
1-1	0.86	−2.2	4.7	0.62	−3.6	−9.3
1-2	0.88	−3.8	2.3	0.44	−7.2	0.0
1-3	0.88	−3.5	8.3	0.49	−5.2	−16.7
1-4	0.89	−3.4	4.3	0.47	−4.1	−8.5
1-5	0.84	−4.5	4.7	0.63	−3.8	−9.3
2-1	0.80	−4.2	11.1	0.62	−2.6	−2.8
2-2	0.63	−5.7	2.2	0.46	−3.6	4.4
2-3	0.65	−4.0	7.5	0.42	−3.8	5.0
2-4	0.77	−5.1	4.5	0.68	−1.6	0.0
2-5	0.86	−4.0	6.4	0.40	−4.2	8.5
AAV	0.81	4.0	5.6	0.52	4.0	6.5

Note: AAV represents averaged absolute value.

shows the NSEs, RPEs, and RTEs of the original SVN model for each flood event. According to Table 4, the averaged absolute value (AAV) of the NSEs for flow depth simulating was 0.81. In addition, the AAVs of the RPEs and RTEs were 4.0% and 5.6%, respectively. Generally, although a good agreement was obtained between observed and simulated flow depths in the original SVN model, there was some pronounced error in the simulations for the peak of flow depth. Specifically, the RPEs between the calculated and observed flow depth were negative at each flood event, whereas the RTEs were positive. This indicated that the peak of the simulated flow depth in the original SVN model lagged behind the observed one and that the simulated maximum flow depth was smaller than that of the observed one.

The AVV of the NSEs for the sectional averaged velocity simulation was 0.52, much lower than that for the flow depth simulation. This indicated that the original SVN model could not simulate the sectional averaged velocity accurately. Since the AAVs of the RPEs and RTEs were 4.0% and 6.5%, respectively, there was some pronounced error in the peak simulation as well. Specifically, the RPEs for the sectional averaged velocity were negative at each flood event, indicating that the simulated maximum sectional averaged velocity was smaller than the observed one. Furthermore, randomly distributed RTEs showed that different flood events had different time-at-peak errors in the sectional averaged velocity simulation.

In order to investigate the mechanism behind the poor performance of the original SVN model, the authors re-ran the original SVN model with three different parameters: c₀ = 0.0013 (the objective function value = 0.0196), c₀ = 0.0018 (the objective function value = 0.0441), and c₀ = 0.0008 (the objective function value = 0.0448). Then, the performances of the original SVN model with different parameters were compared (see Figure 7). According to Figure 7, it can be concluded that although c₀ = 0.0013 could achieve the best overall model performance according to the objective function value, it was not the best parameter for each situation. For example, c₀ = 0.0013 was good for flow depth simulating, whereas c₀ = 0.0018 was better for the sectional averaged velocity simulation. Moreover, c₀ = 0.0013 was a better choice for the Nos. 2-1 and 2-2 flood events, whereas for the Nos. 1-1 and 1-2 flood events, c₀ = 0.0018 was the preferred parameter because it could improve the performance of the sectional averaged velocity simulation without hurting the agreement between observed and simulated flow depths. Hence, not only the flow depth simulation and the sectional averaged velocity simulation had preferred different parameters, but also different flood events had preferred different parameters. In such a case, it is very difficult to obtain an ideal parameter which can support good performance for both the flow depth simulation and the sectional averaged velocity simulation and obtain a sound agreement for every flood event.

4.3. Performance of the Modified SVN Model

The performance of the modified SVN model was evaluated using NSE, RPE, and RTE as well, as is shown in

Table 5. The NSEs, RPEs, and RTEs of the modified SVN model.

No.	Flow Depth			Sectional Averaged Velocity
	NSE	RPE (%)	RTE (%)	NSE	RPE (%)	RTE (%)
1-1	0.97	0.6	0.0	0.89	0.78	4.7
1-2	0.99	−1.7	0.0	0.84	−1.59	4.5
1-3	0.99	−0.1	6.3	0.89	−0.53	2.1
1-4	0.99	−0.5	0.0	0.86	−0.73	−6.4
1-5	0.98	−0.6	4.7	0.85	0.40	0.0
2-1	0.99	0.0	2.8	0.73	0.54	8.3
2-2	0.91	−4.0	0.0	0.93	−0.01	8.9
2-3	0.93	−1.7	0.0	0.86	−0.13	7.5
2-4	0.88	−4.3	−2.3	0.72	2.69	4.5
2-5	0.95	−2.8	0.0	0.88	−0.31	−2.1
AAV	0.96	1.65	1.60	0.84	0.77	4.91

Note: AAV represents averaged absolute value.

. By comparing the results of the modified SVN model in Table 5 with those of the original SVN model in Table 4, one notes that the AAV of NSEs for the flow depth simulation was 0.81 in the original SVN model results and that it was increased to 0.96 in the modified SVN model. The maximum increase in NSE for flow depth simulation was from 0.63 to 0.91 in the No. 2-2 flood event. Moreover, the AAV of NSEs in the sectional averaged velocity simulation has improved from 0.52 to 0.84. The maximum increase was from 0.4 to 0.88 in the No. 2-5 flood event. Therefore, the NSE values were increased in both the flow depth simulation and the sectional averaged velocity simulation. As a result, the agreement between observations and simulations in the modified SVN model is better than those in the original SVN model.

Furthermore, compared with the original SVN model, the modified SVN model reduced the AAV of RPEs from 4.0% to 1.65% for the flow depth simulation and from 4.0% to 0.77% for the sectional averaged velocity simulation. In addition, the AAV of RTEs in the original SVN model was 5.6% for the flow depth simulation and 6.5% for the sectional average velocity simulation, but they were reduced to 1.60% and 4.91%, respectively, in the results of the modified SVN model. The modified SVN model reduced the peak value error and the time-at-peak error in both the flow depth simulation and the sectional averaged velocity simulation.

From the above results, the modified SVN model had higher NSEs and lower RPEs and RTEs than the original SVN model in both simulations of the flow depth and the sectional averaged velocity. It can be concluded that the modified SVN model performs better than the original SVN model. It also indicates that the extended Chezy formula is reasonable.

Figure 8 shows the results of the modified SVN model compared with the original model and the observations. On one hand, Figure 8 confirms that by employing the extended Chezy formula, the modified SVN model can provide a better performance. On the other hand, Figure 8 reveals that some errors remain in the simulation result. For example, for both the modified and original SVN models, there is a decrease in the calculated velocity in the first time steps. It is known that in the numerical experiment, except for the friction calculation, errors will be introduced by the initial and boundary conditions as well. In this paper, this situation was analyzed as follows.

In the numerical experiment, the initial conditions (V and D and time t = 0) were obtained by the linear interpolation method using data from the upstream and downstream cross sections (i.e., cross sections #1 and #3). However, according to Figure 8, the initial V of cross section #2 which was obtained by the linear interpolation method was greater than the observed one, whereas the initial D of cross section #2 which was obtained by the linear interpolation method was slightly less than the observed one (see Figure 8(c1,c2), for example). To counteract the effects of initial condition errors, the calculated velocity decreased and the calculated water depth increased in the first time steps. However, the initial condition error did not result in the failure of the numerical experiment. As the calculation continued, the effect of the initial condition error became less and less important. According to Figure 8, the impact of the initial condition error lasted about four time steps. By comparing the residuals of D and V in the first four time steps and that of the remaining time steps, the additional residual which was introduced by the initial condition error was obtained. Generally, it had an impact on the objective function in the range of 3%–5%. Following the same logic, except for the initial condition error, the result of numerical experiments was also affected by the boundary condition error. It is known that measuring cross-sectional averaged velocity is more difficult than observing water depth. In the numerical experiment, some fluctuations in the calculated averaged velocity hydrographs occurred due to the error in the averaged velocity measurement of cross section #3 (downstream boundary condition). In addition, because the friction slope computed by the extended Chezy formula was based not only on V and D but also on ∂D/∂x, ∂D/∂t and ∂V/∂t, the modified SVN model which employed the extended Chezy formula to compute the friction slope was more sensitive to the error in the averaged velocity measurement of cross section #3 (see Figure 8(a2), for example). In a further study, more data assimilation techniques should be introduced to eliminate the influence of the observation error.

5. Conclusions

In this paper, a newly proposed 1D unsteady friction formula was introduced to the framework of a 1D hydrodynamic model to establish the modified SVN model. With the help of a GA-based parameter calibration algorithm, the performance of the modified SVN model was compared with that of the original SVN model in a simulation test for a hydraulic experiment. The results revealed the following:

The GA-based parameter calibration method was convergent and could identify the empirical parameters efficiently. However, in the original SVN model, not only the flow depth simulation and the sectional averaged velocity simulation had preferred different parameters, but also different flood events had preferred different parameters. In such a case, it is very difficult to obtain an ideal parameter which can support good performance for both the flow depth simulation and the sectional averaged velocity simulation and obtain a sound agreement for every flood event. As a result, although the optimal parameter c₀ = 0.0013 can achieve the best overall model performance, it is not the best parameter for each situation.

Compared with the original SVN model, the modified SVN model can achieve a better simulation in both the flow depth and the sectional averaged velocity simulations. Furthermore, it can reduce the peak value error and the time-at-peak error as well. The modified SVN model is better than the original SVN model. It also indicates that the extended Chezy formula can be used as a substitute for steady uniform friction formulas in the framework of a 1D hydrodynamic model.

This paper provides a simple test verifying the advantages of using the extended Chezy formula in a 1D hydrodynamic model. Despite the fact that similar flume experiments have been widely used to investigate unsteady channel flow behavior, a laboratorial channel is relatively short. For this reason, the flow unsteadiness and flow non-uniformity are limited. Further investigation should concentrate more on flood simulations for field measurements.