Discussion: New and improved four-parameter non-linear Muskingum model

The original paper by Easa (2014) developed a new storage equation for the Muskingum flood routing model to improve the fit to observed outflows. The original research is both useful and interesting. In this discussion, some notes and comments will be presented to extend and uphold the results of the original paper.

In general, the storage of unsteady flow in a natural river depends primarily on the inflow and outflow and on the geometric and hydraulic characteristics of the river and its control features (Akbari and Barati, 2012; Chow, 1959; Viessman and Lewis, 2003). Easa (2014) developed a new four-parameter non-linear Muskingum model by assuming that the flow and storage characteristics at the upstream and downstream cross-sections are constant. If this assumption is not established, the storages at the upstream and downstream sections can be expressed as

where I and Q are simultaneous values of inflow and outflow; S_in and S_out mirror storage at the upstream and downstream sections, respectively; a₁ and n₁ denote the depth-discharge characteristics; b₁ and m₁ denote the mean depth-storage characteristics at the upstream cross-section; and a₂ and n₂, b₂ and m₂ are corresponding values at the downstream cross-section.

By substituting S_in and S_out from the above equation into Equation 16 of the original paper by Easa (2014), the general form of storage–discharge relationship (or storage equation) of the Muskingum model can be expressed as

where β, γ₁=wb₁/(a₁)^m₁/n₁, α₁=m₁/n₁, γ₂=(1−w)b₂/(a₂)^m₂/n₂, and α₂=m₂/n₂ are flood routing parameters; w defines the relative weights given to inflow and outflow for the river reach.

The model has five parameters, γ₁, α₁, γ₂, α₂ and β, which must be calibrated using available flood events. The available storage equations of the Muskingum model including linear and non-linear models along with routing parameters are presented in Table 7. For the linear model, a finite-difference method is available for the numerical solution of the governing equations (i.e. storage equation together with the one-dimensional continuity equation) (Barati, 2013), while for each form of non-linear model, similar to Easa (2014) or Barati (2011), the Euler method can be adopted for the numerical solution. In order to avoid negative discharges in the routing procedure of the non-linear models, the penalty approach which was proposed by Karahan et al. (2013) can be considered.

Interestingly, a better solution can be obtained by using the four-parameter storage equation for example 3 of Easa (2014) as K = 0·4608, w = 0·1665, α = 0·9215 and β = 1·5679, with sum of the square of the deviation (SSQ) equal to 73 379·35.

For a five-parameter model, the flood routing parameters and SSQ values along with the percentage of the improvement over the best results of the four-parameter model for the three case studies of Easa (2014) are presented in Table 8. The results indicated that the five-parameter model yielded better results than the four-parameter one. The significant improvement is achieved for the first and latter case studies, although the improvement for the second case is not considerable.

Different weighted-flow and storage volume relationships in natural rivers may be occurring around the world. By considering the available storage equations of Table 7, it is possible to select the proper storage equation for a given natural river. As discussed by Barati (2014), the proper storage–discharge relationship for each river must be selected by considering the relationship between weighted-flow and storage volume of the flood events. However, it should be stated that the calibration procedure of the routing parameters may be more challenging; also the dependence of the model's results on the calibration data may increase by increasing the number of parameters of the storage equation. Therefore, by considering the degree of non-linearity of the flood event, the engineers must select a storage equation with lower routing parameters, as far as possible, without losing accuracy in the flood routing simulation.

The author thanks the discussers for their interest in the paper. Their extension to include separate flow parameters for the upstream and downstream sections is useful for the cases in which the sections have different characteristics. There are just a few points that are worthy of note.

The discussers' five-parameter model should naturally improve model performance as it has more parameters, compared with the four-parameter model. As the discussers correctly pointed out, a trade-off should be made between the extra parameters and the resulting improvement in model performance. However, a relatively large number of parameters no longer presents an issue given modern computers and superior optimisation software.

The discussers have obtained a slightly improved solution of SSQ = 73 379·35 for example 3 with the following optimal parameters K = 0·0768, w = 0·1665, α = 0·9215 and β = 1·5679. This is indeed a better solution. The author obtained a solution with identical SSQ and similar optimal parameter values, except that the optimal value of K was 0·4609 (now corrected by discussers). It is worth noting that the flood routing optimisation problem is highly non-linear and non-convex. Generally, it is possible to improve the solution by running the model several times, each time starting with the previous ‘optimal' solution. The process is stopped when the improvement in the optimal solution is not large.

Although the discussers' model is theoretically correct, a better model that preserves the physical significance of the original parameters is given by

where K1 and K2 = storage parameters for the upstream and downstream sections, respectively, given by b1 /aα11 and b2 /aα22. For constant characteristics of the upstream and downstream sections, K1 = K2, a1 = a2, and Equation 26 reduces to the four-parameter model of the original paper by Easa (2014).

The discussers presented in Table 7 a summary of the available storage equations of the Muskingum model. There is, however, another storage equation that should be included in this table. It is a three-parameter non-linear Muskingum model with variable exponent parameter (VEP), which takes the form (Easa, 2013)

where β(u_j) = exponent parameter for time interval j and u_j = dimensionless inflow variable (0 to 1) for time interval j = I_j/I_max, where I_j = inflow for time interval j and I_max = maximum inflow during the routing period. The number of exponent parameters depends on the number of specified inflow levels L and equals 2 + L. The exponent parameter was modelled using a step function which is given by (for L = 3 for example)

where A_i = dimensionless inflow interval i (i = 1, 2, 3) and χ_Ai(u_j) = indicator function of A_i, which is defined as

Note that there are three exponent parameters (β₁, β₂ and β₃) that are defined for three dimensionless inflow ranges (0, v₁), (v₁, v₂), and (v₂, 1), where v₁ and v₂ are dimensionless inflow dividing values.

Extending the idea of VEP to the four-parameter model presented in the paper yields the following storage equation

The total number of parameters of this model is 3 + L. For L = 2, for example, the number of parameters is 5. Note that this five-parameter model with VEP is obtained by using two exponent parameters, whereas the discussers' five-parameter model is obtained by using two flow parameters. To allow comparison of the two models, the five-parameter model with VEP was applied to the three examples and the results are shown in Table 9. As noted, the reductions in SSQ, compared with the discussers' model, are 2, 28 and 14% for examples 1, 2, and 3, respectively (for example 1, which already has an almost perfect fit, a continuous function of the exponent parameter was used). In addition, a six-parameter model (L = 3) was evaluated and SSQ was reduced further to 31 and 33%, for examples 2 and 3 respectively. Clearly, the extra parameter has more than doubled SSQ reduction for example 3 (multiple-peak hydrograph), but produced a small improvement for example 2 (non-smooth, single-peak hydrograph). Thus, the number of inflow levels may vary from one situation to another.

Based on these and the discussers' results, it appears that the exponent parameter β has a more significant effect on model performance than the flow parameter α. This was perhaps the main reason most researchers in the past have focused on the non-linear model with β (NL2) than the one with α (NL1), as mentioned in the paper. The encouraging results of the VEP presented in this response have provided the motivation to extend the paper's four-parameter Muskingum model to account for the variability of all four model parameters with flood characteristics. A forthcoming work by the author, which has been submitted for publication, demonstrates that the model would be quite useful for all types of hydrographs.

In conclusion, over the past several decades researchers have attempted to improve the performance of the non-linear Muskingum model by adopting different solution algorithms. This strategy, however, has resulted in a very slight improvement. To promote alternative thinking, the author has recently focused on changing the structure of the model itself. The author was pleased to see the discussers' contribution as it followed the same thinking.