1 Introduction

If a channel reach is not influenced by existing hydraulic structures, and its bottom slope is relatively steep, then simplified models can be used for flood routing (Ponce et al. 1978; Singh 1996). Among the known existing simplified models, the linear Muskingum model (McCarthy 1938) still attracts the attention of hydrologists. In this approach, a channel reach of length L is divided into N space intervals of length:

$$ \Delta x=L/N. $$
(1)

For each interval, the following relation has been proposed:

$$ \frac{d}{dt}\left[K\left(X\cdotp {Q}_{j-1}+\left(1-X\right){Q}_j\right)\right]={Q}_{j-1}-{Q}_j $$
(2)

where:

t:

time

Qj–1(t):

inflow at the upstream end of a space interval j,

Qj(t):

outflow at the downstream end of a space interval j,

X:

constant parameter,

K = Δx/C:

constant parameter considered as the time needed to cover the space interval of length Δx by a flood wave which propagates with celerity C.

Equation (2), applied for all space intervals – reservoirs (j = 1,2, ..., N), constitutes the system of ordinary differential equations describing the dynamics of a cascade of reservoirs. This cascade is represented by a set of nodes located along the x axis coinciding with the channel axis (Fig. 1).

Fig. 1
figure 1

Discretization of a channel reach for the Muskingum equation

Full size image

Comparing Fig. 1 with Eq. (2), it can be noticed that the expression in brackets at the left-hand side represents the flow rate at point P, calculated using the linear interpolation between the nodes j - 1 and j:

$$ {Q}_P=X\cdotp {Q}_{j-1}+\left(1-X\right){Q}_j. $$
(3)

The position of point P is determined by the weighting parameter X, which ranges from 0 to 1.

For a long time, the real nature of the Muskingum model remained unrecognized. The kinematic character of the Muskingum equation was deduced for the first time by Cunge (1969) via a comparison of the approximated pure advection equation with the Muskingum equation, solved numerically using the implicit trapezoidal rule. Nowadays, this fact can be confirmed formally in different ways. For instance, the Muskingum equation can be derived similarly to the kinematic wave model, i.e. directly from the differential continuity equation and the simplified dynamic equation in the form of the Manning or Chézy formula, which means that both models have the same origin.

The solution of the Muskingum Eq. (2) can be carried out using the following two-level difference scheme (Szymkiewicz 2010):

$$ {y}_j={y}_{j-1}+h\cdotp \left(\left(1-\theta \right)y{\hbox{'}}_{j-1}+\theta \cdotp y{\hbox{'}}_j\right) $$
(4)

where:

yj:

unknown value of function y at node j,

y'j:

value of derivative of function y at node j,

h:

step of integration,

θ:

weighting parameter ranging from 0 to 1.

The application of formula (4) for the solution of Eq. (2) yields the following linear algebraic equation:

$$ {\displaystyle \begin{array}{l}K\left(X\cdotp {Q}_{j-1}^{n+1}+\left(1-X\right){Q}_j^{n+1}\right)=K\left(X\cdotp {Q}_{j-1}^n+\left(1-X\right){Q}_j^n\right)+\\ {}\kern2em +\Delta t\left[\left(1-\theta \right)\left({Q}_{j-1}^n-{Q}_j^n\right)+\theta \left({Q}_{j-1}^{n+1}-{Q}_j^{n+1}\right)\right]\ \end{array}}\kern0.5em \left(j=1,2,\dots, N\right) $$
(5)

where:

n:

index of the time level,

Δt:

time step.

Note that the parameter θ determines the method used for time integration. For θ = 1/2, the implicit trapezoidal rule is obtained, commonly used for the solution of the linear Muskingum Eq. (2). Using the initial condition \( {Q}_j^{n=0}={Q}_{in} \) for j = 1, 2, …, N (where Qin is the given initial flow rate) and the hydrograph Q0(t) imposed at the upstream end, the outflows \( {Q}_j^{n+1} \) (j = 1, 2, …, N) are calculated directly from Eq. (5). Note that the calculation is very simple because in Eq. (5) \( {Q}_j^{n+1} \) is the only one unknown. This simplicity of the solution resulted in the wide use of Eq. (5) for hydrological forecasting. However, the linear Muskingum equation very often provides results representing poor accuracy. For this reason, hydrologists have tried to improve it. Consequently, nonlinear forms of the Muskingum equation were developed (Gill 1978; Tung 1985; Singh and Scarlatos 1987; Hirpurkar and Ghare 2015; Kang et al. 2017; Bozorg-Haddad et al. 2019). Instead of the linear formula relating to storage, inflow and outflow proposed by McCarthy (1938), the storage equation has been completed by a priori assuming the nonlinear relation leading to the following equation:

$$ \frac{d}{dt}\left[k{\left(X\cdotp {Q}_{j-1}+\left(1-X\right){Q}_j\right)}^r\right]={Q}_{j-1}-{Q}_j\kern2em \left(j=1,2,\dots, N\right) $$
(6)

where k, X and r are constant parameters. For r = 1.0, Eq. (6) becomes the linear Muskingum Eq. (2). In such a case, the parameter k is expressed in time units and it coincides with the parameter K involved in Eq. (2).

It is interesting that the introduced nonlinear relation, i.e. the expression in brackets at the left side of Eq. (6), has no physical interpretations. However, this causes that Eq. (6) is written in a divergent form. Consequently, the numerical solution of Eq. (6) does not generate a mass balance error (Gąsiorowski and Szymkiewicz 2007; Gąsiorowski 2013).

The parameters k, X and r occurring in Eq. (6) are considered as free ones and their values are estimated via a comparison of the numerical solution with the field data. To solve this problem, optimization methods should be used. The optimal values of these parameters should minimize the assumed objective function, which determines the discrepancy between the results obtained from the model and those provided by the measurements. In the context of the estimated parameters for the nonlinear Muskingum equation, many optimization methods have been proposed (the Hook-Jeeves algorithm (Tung 1985), the Nelder–Mead simplex method (Barati 2011), the Broyden-Fletcher-Goldfarb-Shanno algorithm (Geem 2006) among others).

The aforementioned methods are effective and lead to a unique optimal solution when the assumed objective function possesses a unimodal property; in other words, when the objective function has one extreme point in a space of parameters. However, in the case of the nonlinear Muskingum Eq. (6), the objective function possesses more extreme points. In such cases, the mentioned methods usually converge towards a local extreme point. In order to avoid this problem, meta-heuristic optimization algorithms have been developed (Glover and Kochenberger 2003). For parameter identification in the nonlinear Muskingum equation, such methods were used, for instance, by Mohan (1997), Chu and Chang (2009), Luo and Xie (2010), Xu et al. (2012), Niazkar and Afzali (2015), Moghaddam et al. (2016), Yuan et al. (2016), Hamedi et al. (2016) and Farzin et al. (2018).

A review of the presented papers shows that their authors mainly focused attention on techniques of parameter identification. Physical backgrounds of the a priori assumed formula, relating to storage, inflow and outflow and the formal interpretation of the governing Eq. (6) have not been discussed. However, comprehensive recognition of the nonlinear Muskingum equation properties seems to be very important. This knowledge allows the removal of a serious formal disadvantage related to the dimensional inconsistency of Eq. (6) and the improvement of its effectiveness and flexibility, which are especially important during model calibration.

Taking into account the above issues, this paper presents some important aspects with regard to the calibration process when the values of parameters involved in the nonlinear Muskingum equation are identified. First of all, it is shown that the best value of the exponent parameter r is closely related to the slope of the channel bed, which determines the character of the open channel flow. When the bed slope increases, the open channel flow tends to a kinematic one. In such a case, the value of r tends to its limit value 3/5 when the Manning equation is used or to 2/3 when the Chézy equation is applied. It is worth adding that these values of the parameter r resulting from the kinematic wave theory ensure the dimensional consistency of the nonlinear Muskingum Eq. (6) (Gąsiorowski and Szymkiewicz 2018). Moreover, it is shown that the shape of the computed flood wave is determined not only by the parameters involved in the Muskingum Eq. (6), as it is commonly assumed. The accuracy analysis carried out using the modified equation approach showed that the solution accuracy depends also on the space interval Δx (or on the number of reservoirs N due to Eq. (1)). Therefore, this must also be considered as a free parameter. The possible impact of the time step Δt on the numerical solution is also considered in the paper. If, for time integration, the method generating additional numerical diffusion is used, i.e. when θ > 1/2 is assumed, then the flood wave attenuation will be increased. In such a case, the time step Δt and the weighting parameter θ should also be considered as free parameters. Consequently, the assumed objective function should depend on all the parameters of the mentioned model and numerical parameters.

All the presented hypotheses are confirmed by the results of numerical tests.

2 Numerical Origin of Flood Wave Attenuation by the Muskingum Equation

When numerical techniques are applied, then the solution accuracy should be comprehensively examined. In the case of hyperbolic problems, “the modified equation approach” (Warming and Hyett 1974; Fletcher 1991) is very useful, which provides full information on the solution accuracy. This technique, applicable for linear equations, delivers very useful information also for the nonlinear Muskingum equation. Some aspects of solution accuracy for its nonlinear form have been invoked by Vatankhah (2010), Das (2010), Karahan (2014) and Easa (2015). However, the discussions undertaken in all the mentioned publications do not contain any detailed accuracy analysis of the applied method nor the resulting consequences for the parameter identification process.

The modified equation approach performed for Eq. (5) leads to the following equation, valid in each node (Gąsiorowski and Szymkiewicz 2018):

$$ \frac{\partial Q}{\partial t}+C\frac{\partial Q}{\partial x}=\left({V}_{n1}+{V}_{n2}\right)\frac{\partial^2Q}{{\partial x}^2}+{\varepsilon}_n\frac{\partial^3Q}{{\partial x}^3}{+}_{\dots } $$
(7)

with

$$ {V}_{n1}=\frac{{\Delta x}^2}{2K}\left(1-2X\right), $$
(8)
$$ {V}_{n2}=\frac{{\Delta x}^2\cdotp \Delta t}{2{K}^2}\left(2\theta -1\right) $$
(9)
$$ {\varepsilon}_n=\frac{\varDelta {x}^3}{6K}\left(\left(2-3\theta \right){\left(\frac{\varDelta t}{K}\right)}^2+3\left(X+\theta -1\right)\frac{\varDelta t}{K}+\left(1-3X\right)\right) $$
(10)

where νn1 and νn2 are the coefficients of numerical diffusion, whereas εn is the coefficient of numerical dispersion. Equation (7) represents the governing equation modified by the applied numerical solution method. The left side of this equation corresponds to the kinematic wave equation (Gąsiorowski and Szymkiewicz 2018), whereas the right side contains all terms of the Taylor series expansion, neglected during approximation. This fact confirms that, indeed, the Muskingum equation is a particular form of kinematic wave equation. Note that for X = 1/2, θ = 1/2 and Δt = K all terms at the right-hand side of Eq. (7) are cancelled. Consequently, the difference Eq. (5) ensures a pure translation of the flood wave, being also the exact solution of the linear kinematic wave equation.

The terms of numerical diffusion occurring in Eq. (7) are responsible for unphysical wave amplitude attenuation. It is worth adding that the coefficient νn1 (Eq. (8)) coincides with the formula derived by Cunge (1969). As Cunge (1969) applied the implicit trapezoidal rule (θ = 1/2) for time integration, in his approach the coefficient νn2 did not occur.

As far as numerical dispersion is concerned, it changes the wave celerity. Consequently, unphysical oscillations occur in the numerical solution. Such an effect is frequently provided by the Muskingum equation. Spurious oscillations are insignificant in the case when numerical diffusion is strong enough.

The abovementioned brief analysis should draw attention to the well-known fact that the attenuation of the flood wave observed in the numerical solution provided by the Muskingum equation has no physical background, but conversely, has a purely numerical origin. Eqs. (7, 8, 9) show that the intensity of flood wave attenuation is determined by the following parameters: X, Δx (or N), Δt and θ.

Since the Muskingum equation is a particular form of kinematic wave (Szymkiewicz 2002), then its exact solution should ensure wave translation without attenuation. Indeed, as it was previously presented for X = 1/2, Δt = K and θ = 1/2, the numerical solution of the linear Muskingum Eq. (5) provides an exact solution of the kinematic wave, because all terms at the right-hand side of Eq. (7) are cancelled. However, hydrologists require significant damping of the flood wave, as observed in rivers. Summarizing, a satisfying agreement between the numerical solution of the Muskingum equation and the measurements can be obtained owing to the numerical error generated in the solution, i.e. when the numerical solution represents poor accuracy. Conversely, a solution of high accuracy usually disagrees with the observations. Thus, to achieve the required wave damping, identification must also be carried out with regard to the numerical parameters involved in determining the solution accuracy.

3 Dimensionally Consistent Nonlinear Muskingum Equation

As far as the parameters X and r involved in Eq. (6) are concerned, they are dimensionless, whereas the units of parameter k are a priori unknown because they are related to the value of parameter r according to the following rule:

$$ {\mathrm{L}}^{3\left(1-r\right)}\cdotp {\mathrm{T}}^r $$

where L denotes units of length while T denotes units of time. It means that for each set of data used for identification, the value of the obtained parameter r will determine different units of the parameter k. Finally, it can be stated that the nonlinear Muskingum Eq. (6) with the parameter r considered as free is generally dimensionally inconsistent (Gąsiorowski and Szymkiewicz 2018).

The problem of the dimensionally inappropriate form of the nonlinear Muskingum equation has not been discussed with regard to the identification process. As noted by Easa (2015), the application of different optimization techniques causes only slight progress in obtaining better results during parameter identification. To eliminate this disadvantage of Eq. (6), Gąsiorowski and Szymkiewicz (2018) proposed its alternative form with the parameter r having a fixed and constant value, which ensures its dimensional consistency. To derive such a version of the nonlinear Muskingum equation, the differential continuity equation and the simplified dynamic equation, expressed by the Manning or Chézy formula, were used as the basic relations. After transformation with no additional assumptions, the final formula:

$$ \frac{d}{dt}\left[k{\left(X\cdotp {Q}_{j-1}+\left(1-X\right){Q}_j\right)}^m\right]={Q}_{j-1}-{Q}_j $$
(11)

valid for all space intervals – reservoirs, i.e. for j = 1, 2, …, N, was derived (Gąsiorowski and Szymkiewicz 2018). A comparison of this equation with Eq. (6) shows that the only difference is with regard to the exponent r. In Eq. (6) it is considered as a free parameter, whereas in Eq. (11) it takes the constant value r = m. The nonlinear Eq. (11) is written in a dimensionally consistent form, in which the exponent is equal to the constant value m = 3/5 for the Manning formula or equal to m = 2/3 for the Chézy formula.

The parameter k involved in Eq. (11) can be related to the parameter K occurring in the linear Muskingum Eq. (2). Since K represents the time needed to cover a channel reach of length Δx by the propagating flood wave, then the related formula is as follows:

$$ K=k\cdotp m{\left({Q}_P\right)}^{m-1} $$
(12)

As can be noticed, for the linear Muskingum equation in which m = 1, the parameter k coincides with the parameter K.

4 Solution of the Nonlinear Muskingum Equation and Parameter Identification

As far as the numerical solution of the nonlinear Muskingum equation is concerned, to this order, the previously applied two-level method (Eq. (4)) was used for both of its versions. The application of this method for the solution of Eq. (6) yields the following nonlinear algebraic equation for each reservoir:

$$ {\displaystyle \begin{array}{l}{\left(X\cdotp {Q}_{j-1}^{n+1}+\left(1-X\right){Q}_j^{n+1}\right)}^r={\left(X\cdotp {Q}_{j-1}^n+\left(1-X\right){Q}_j^n\right)}^r+\\ {}\kern2em +\frac{\Delta t}{k}\left[\left(1-\theta \right)\left({Q}_{j-1}^n-{Q}_j^n\right)+\theta \left({Q}_{j-1}^{n+1}-{Q}_j^{n+1}\right)\right].\left(j=1,2,\dots, N\right)\end{array}} $$
(13)

Since at each node j = 1, 2, …, N Eq. (13) is the nonlinear algebraic equation with only one unknown \( {Q}_j^{n+1} \), it should be solved for the consecutive cross-sections using an iterative method, such as the Picard method or the Newton method.

The best values of the parameters involved in the Muskingum equation must minimize the following assumed objective function:

$$ F\left({P}_i\right)=\underset{t=0}{\overset{T}{\int }}{\left({Q}_d(t)-{Q}_N(t)\right)}^2 dt\to \mathit{\min} $$
(14)

where:

F(pi):

objective function with respect to the parameters pi (p1 = X, p2 = k, p3 = r, p4 = N),

Qd(t):

hydrograph observed at the downstream end of a channel reach,

QN(t):

hydrograph computed at the downstream end of a channel reach,

T:

total simulation time of flow.

In order to find the best values of parameters, the Powell’s algorithm is applied (Powell 1964). Its description is also given, for instance, by Press et al. (1992). Iterations are continued until the following convergence condition is satisfied:

$$ \frac{2\cdotp \left|{F}^{\left(k+1\right)}-{F}^{(k)}\right|}{\left|{F}^{\left(k+1\right)}\right|+\left|{F}^{(k)}\right|}\le \varepsilon $$
(15)

where:

F:

value of the objective function given by Eq. (14),

k:

iteration index,

ε:

specified tolerance.

In all tests presented below, a rectangular channel of the length L = 100 km and of the width B = 50 m is considered. The channel reach is divided into N subintervals (reservoirs). At t = 0 the flow rate along the channel is constant and equal to a given value equal to q0. At the upstream end, the following flood wave enters the channel reach:

$$ {Q}_0(t)={q}_0+\left({q}_{max}-{q}_0\right){\left(\frac{t}{t_{max}}\right)}^{\beta}\mathit{\exp}\left(1-{\left(\frac{t}{t_{max}}\right)}^{\beta}\right) $$
(16)

where:

q0:

base flow discharge of the inflow,

qmax:

peak discharge of the inflow,

tmax:

time of the peak flow,

β:

parameter determining the shape of the wave.

The values of parameters involved in Eq. (16) have been arbitrarily assumed as follows: q0 = 20 m3/s, qmax = 100 m3/s, tmax = 30 h, β = 2.0. The numerical solution of the Muskingum equation is compared with the solution of the system of Saint Venant equations. Owing to this assumption, a discussion dealing with the properties of the Muskingum equation seems to be easier than in the case when field measurements are used. This is because very often the field data do not satisfy the mass balance.

In the dynamic equation of the Saint Venant model, the friction term is expressed using the Manning formula with the roughness coefficient n equal to 0.03 m-1/3s and the bottom slope of a channel s equal to 0.00008. The Saint Venant equations are solved using the implicit Preissmann scheme (Cunge et al. 1980) with the following mesh dimensions: space interval Δx = 500 m, time step Δt = 500 s.

The numerical tests showed that when the parameter r is considered as free, the final effect of the identification is very sensitive on the assumed initial values of the parameters. In such a case, the objective function possesses many local extreme points. When for a fixed value X the variability with regard to r and k is considered, then, indeed, the objective function F(r,k) has an unusual shape, as presented in Fig. 2.

Fig. 2
figure 2

An example contour map of the objective function F(r,k) (Eq. (14)) obtained using the equation in a dimensionally inconsistent form (Eq. (6)) for the tolerance ε = 10−2 (X = 0.255, Δt = 6h, θ = 0.5)

Full size image

Typically, a single global optimum and numerous local extreme points can be distinguished. Moreover, the values of F(r,k) at these points differ insignificantly from each other, as well as from the global minimum. To overcome these computational difficulties, the accuracy of the optimization algorithm should be significantly increased. To this regard, a very small value of tolerance (equal to ε = 10−5) was assumed. However, it should be noted that the higher accuracy of the identification process implies that the average number of iterations significantly increases. For instance, the number of iterations equal to 240 (corresponding to the original tolerance ε = 10−2) was increased to about 2500 for ε = 10−5. An increase in the number of iterations is also associated with an extension of the computational time. In this case, the CPU time was increased from about 0.3 s to over 4 s, respectively.

The untypical form of the contour map of the objective function F(r,k) is the cause of another problem. Since very similar minimal values of the objective function can be obtained for various combinations of the values X, k and r, then in consequence, very similar numerical solutions in terms of the downstream hydrograph can be expected. This situation is obviously unfavorable, because, in fact, Eq. (6) can deliver a non-unique solution when the identification process is carried out with low accuracy.

To avoid the disadvantages listed above, the nonlinear Muskingum equation with a constant value of exponent r is considered. Therefore, in Eq. (13) r = m = 3/5 is assumed, whereas others parameters, i.e. X and k, are considered as free ones. Similarly, the number of reservoirs N as well as the weighting parameter θ and the time step Δt, both involved in the time integration procedure, are considered also as free parameters.

5 Influence of the Number of Reservoirs N on the Objective Function F(pi)

In this test, the parameters X, k and N are identified, whereas other parameters are assumed to be constant: r = m = 3/5, θ = 1/2 (the implicit trapezoidal rule is used), the time step is equal to Δt = 6 h, the channel bottom slope is equal to s = 0.00008, the tolerance for the convergence condition in the optimization algorithm (Eq. (15)) is equal to ε = 10−2. Since the number of reservoirs N is an integer, then the identification is carried out with regard to the two variables X and k only, but for consecutively assumed different values of N. The results of computations are summarized in Table 1.

Table 1 The objective function F computed using Eq. (11) for a different number of reservoirs N with r = m = 3/5 and θ = 1/2
Full size table

Comparing the values of the objective function obtained for a different number of reservoirs N, it can be noticed that, indeed, it significantly influences the results of parameter identification. The best result is obtained for N = 3, for which the extreme value of the objective function is equal to F = 31.66. It is important that the same results are obtained independently from the assumed initial values of the searched parameters X0 and k0. The explanation of this property results from the contour map of the objective function F(X, k) (Eq. (14)) displayed in Fig. 3. It is apparent that in this case, the objective function possesses only one extreme point. This means that it can be reached using even a very simple method of optimization. Depending on the starting set of parameters, the total number of required iterations for ε = 10−2 ranges from 183 to 271. An analysis of the computed best values of the objective function, presented in Table 1, shows that the frequently-used one reservoir only (N = 1) is not acceptable because this makes it impossible to satisfy the matching of the computational results and the observations. In this case, the best value of the objective function is equal to F = 1423.30, while for N = 3, it is equal to F = 31.66 (Table 1).

Fig. 3
figure 3

Contour map of the objective function F(X,k) computed using Eq. (11) for N = 3

Full size image

The identified values of both the parameters X = 0.260 and k = 97.03, together with the assumed data characterizing the considered channel and flood wave, were used for the solution of the nonlinear Muskingum Eq. (11). In Fig. 4, the obtained results are compared with the corresponding solution of the Saint Venant equations. As can be seen, both solutions agree nearly perfectly.

Fig. 4
figure 4

Solution of the Saint Venant equations compared with the solution of the nonlinear Muskingum equation (Eq. (11)) obtained for N = 3, X = 0.260 and k = 97.03

Full size image

6 Influence of the Weighting Parameter θ and the Time Step Δt

In this section, the impact of the accuracy of the time integration method used for the solution of the nonlinear Muskingum Eq. (11) with r = 3/5 is analyzed. As can be seen from Eqs. (7, 8, 9) for θ = 1/2, the value of the time step Δt does not influence the numerical solution. This is because in such a case, any additional numerical diffusion is not produced (νn2 = 0). However, if θ > 1/2 is assumed, then additional numerical diffusion is generated. Consequently, the wave attenuation is increased. Its intensity depends also on the value of the time step Δt (see Eq. (9)). The impact of numerical diffusion is illustrated by the computations carried out for the same data as used previously, but for different values of the weighting parameter θ and for various values of the time step Δt. The bottom slope s = 0.00008 is assumed, whereas the number of reservoirs is equal to N = 3. The obtained results are displayed in Table 2.

Table 2 The objective function F(X, k) computed using Eq. (11) for θ = 1/2 and for the different values of the time step Δt with N = 3
Full size table

It is apparent that when the implicit trapezoidal method is used, then the best values of parameters X and k obtained for various values of the time step Δt are practically the same (X = 0.260, k = 97.03÷98.20). Consequently, the values of the computed objective function F(X, k) are also very similar. Therefore, it can be expected that the computed hydrographs at the downstream end should also be very similar to one another. Indeed, as can be seen in Fig. 5, the curves QN(t) computed for the time steps listed in Table 2 with the corresponding values of parameters X, k and N differ insignificantly. Moreover, regardless of the assumed values of the time step, they are in perfect agreement with the solution of the Saint Venant equations, considered as the benchmark.

Fig. 5
figure 5

Solution of the nonlinear Muskingum equation (Eq. (11) with m = 3/5) obtained using the implicit trapezoidal method (θ = 1/2) for various time steps Δt with X = 0.260, k = 97.96, N = 3

Full size image

The presented results seem to be obvious. The required wave damping is ensured by the numerical diffusion caused by the approximation of the 1st order space derivative only. The value of the coefficient of numerical diffusion, roughly estimated for X = 0.260 and k = 96.03 (Table 2) using Eq. (8), is equal to ca νn1 = 7000 m2/s.

The opposite is true for θ > 1/2. In this case, the identified values of parameters are changed depending on the assumed value of Δt (Table 3). The wave attenuation caused by the approximation of the time derivative takes the greatest value for θ = 1, i.e. when the implicit Euler scheme is used. Moreover, attenuation increases with an increase in the time step. This experiment confirms the results of the accuracy analysis previously presented. According to Eq. (9), for θ = 1.0 the coefficient of numerical diffusion νn2 takes the largest possible value, which, additionally, depends also on the time step Δt.

Table 3 The objective function F(X, k) computed for N = 3 with θ > 0.5 and for various values of the time step Δt
Full size table

Summarizing, when the method of the 1st order of accuracy (i.e. with θ > 1/2) is used for the solution, then the optimal values of parameters k, X and N will also be influenced by the assumed value of the time step Δt. Consequently, when the identification process is carried out for some assumed value of the time step, then the same time step has to be used during the validation process. Conversely, if the calculations are carried out with a greater or smaller time step than that used during identification, then one can expect to obtain a solution affected by variable numerical diffusion. For illustration, consider the following example. The optimal values of parameters obtained during identification with the time step equal to Δt = 3 h are subsequently used for the solution of the nonlinear Muskingum equation with other values of the time step Δt. The results of calculations provided by the implicit Euler method (θ = 1), and for X = 0.381 and k = 93.86, previously obtained for Δt = 3 h (Table 3), are presented in Fig. 6.

Fig. 6
figure 6

Solution of the Muskingum equation (Eq. (11) with N = 3 obtained using θ = 1 with X = 0.381 and k = 98.36 for various time steps Δt

Full size image

As can be seen, in this case, the solution of the nonlinear Muskingum equation significantly deviates from the corresponding solution obtained by means of the system of Saint Venant equations. For Δt = 6 h, an excessive attenuation and smoothing of the flood wave is observed, while for Δt = 1.5 h, an overestimation of the peak discharge is generated. Only in the case of the time step equal to 3 h was the correct solution obtained because for such a value the identification of parameters had been carried out.

7 Influence of the Channel Bed Slope s on the Exponent r

It is well known that the longitudinal slope of a channel determines the character of the occurring flow (Ponce et al. 1978). When the channel bottom slope s increases from small to middle values, the character of the propagating wave varies systematically from a diffusive one towards a kinematic one. Below, the results provided by numerical experiments illustrate the influence of the channel bed slope s on the identification process.

The nonlinear Muskingum Eq. (6) is solved by the implicit trapezoidal rule (θ = 1/2) and the value of the time step is equal to Δt = 6 h. To ensure the variable character of the propagating wave from a diffusive one with a high attenuation of amplitude to a kinematic one without attenuation, it is assumed that computations using the Saint Venant equations with the Manning formula used in the friction term are carried out for different values of the channel bed slope s ranging from 0.00003 to 0.0008. The parameters X, k, r and the number of reservoirs N are considered as free ones and they are then subjected to identification. The results of identification obtained for the assumed tolerance equal to ε = 10−5 are displayed in Table 4.

Table 4 Optimal values of the parameters in Eq. (6) and minimal values of the objective function F computed for various bottom slopes s and for the Manning formula used in the Saint Venant equations
Full size table

As can be seen, an increase in the value of the channel bed slope s causes the corresponding values of the parameter r to also increase. This tendency is more clearly seen in Fig. 7 showing the impact of the channel bed slope s on the value of parameter r. For increasing channel bed slope s the value of parameter r approaches a constant limit value equal to 0.60 (=3/5). In other words, when the character of the propagating wave becomes more kinematic, then the parameter r tends to the value which corresponds to the derived kinematic wave exponent m (Eq. (11)).

Fig. 7
figure 7

The parameter r vs. the channel bed slope s in the case when the Manning formula is used in the Saint Venant equations

Full size image

The impact of the variable flow character on the flood routing process is illustrated in Fig. 8. For the selected values of the channel bed slope s and for the corresponding optimal values of the identified parameters taken from Table 4, the hydrographs at the downstream end are computed and compared with the respective solutions of the Saint Venant equations.

Fig. 8
figure 8

Comparison of the solution of the Saint Venant equations and the nonlinear Muskingum Eq. (6) obtained for selected values of the channel bottom slope: s = 0.4·10−4 (a); s = 1·10−4 (b); s = 8·10−4 (c)

Full size image

As can be seen in Fig. 8, a satisfactory agreement is obtained for all considered cases representing either the diffusive (a) or kinematic (c) character of the open channel flow. Note that both curves corresponding to the case (c) represent nearly a pure translation of the considered flood wave, which coincides with the exact solution of the kinematic wave equation. Taking into account the presented results, it can be concluded that even the calibration process applied for the identification of the nonlinear Muskingum equation parameters confirms its kinematic origin.

The same property of the nonlinear Muskingum equation is observed when, instead of the Manning formula (r = 3/5), the Chézy formula (r = 2/3) is used. Assuming, in the Chézy formula, a roughness coefficient equal to c = 30 m0.5/s, the obtained results of calculations are presented in Fig. 9.

Fig. 9
figure 9

The parameter r vs. the channel bed slope s in the case when the Chézy formula is used in the Saint Venant equations

Full size image

As can be seen, the relationship between the values of parameter r and the channel bed slope s is the same as that previously presented in Fig. 7. However, in this case, the values of the identified parameter r tend to the value 0.666 (=2/3), corresponding to the Chézy formula applied in the Saint Venant equations. This property confirms that the alternative dimensionally consistent nonlinear Muskingum Eq. (11), in which the parameter m takes a fixed value depending on the assumed form of the steady uniform flow equation, is correctly derived.

8 Summary and Conclusions

In the paper, two forms of the nonlinear Muskingum equation have been discussed. Both versions differ in their interpretation of the exponent parameter. While the first version, derived via a priori, assuming an additional formula relating to inflow, outflow and storage has no physical interpretation, the second one, derived directly from the kinematic wave model, is physically correct. The first version, with the exponent parameter considered as free, is dimensionally inconsistent, whereas the second one is dimensionally consistent. The nonlinear Muskingum equation with a constant value of the exponent parameter ensures that the assumed objective function has one extreme point only. For this reason, the dimensionally consistent form always ensures optimal values of parameters regardless of their initial values. In addition, such a solution can be obtained even using a direct search optimization algorithm such as the Powell’s algorithm with low accuracy of the identification process. Conversely, the nonlinear Muskingum equation with a variable exponent parameter provides the objective function with numerous very similar local extreme points.

Taking into account the kinematic nature of the Muskingum equation, as well as the numerical origin of wave attenuation, it was shown that apart from the parameters usually identified in the nonlinear Muskingum model, additional parameters such as the number of reservoirs, the weighting parameter and the time step should also be considered as free parameters. The numerical experiments have shown that the frequently-used one reservoir only seems to be unacceptable because this can make impossible the correct adjustment of the numerical solution and the observations. The best compliance was achieved for the number of reservoirs in the range of 2 to 5.

It can be stated that the most important conclusion is that in the case when the open channel flow has a strong kinematic character, then for an increase in the channel bed slope, the value of the exponent parameter, considered as free, tends to a constant value equal to 3/5 for the Manning formula or 2/3 for the Chézy formula.