The Muskingum routing model is favored by water engineers owing to its simplicity and accuracy. A large amount of research is done to improve the accuracy of the model. One way to do so is to consider variable hydrological parameters during the flood routing period. In this study, the random selection (RS) method was proposed to divide the flood period of the nonlinear Muskingum model into three sub-periods. The proposed method was based on RS of members in each sub-region. It was applied to rout three flood hydrographs, and the objective function was the sum of squared errors. Comparing the results from the three variable-parameter nonlinear Muskingum model with those from the variable-parameter nonlinear Muskingum models in previous studies, the proposed model optimized the objective function in these hydrographs up to 61%. The uncertainty analysis of Muskingum parameters for Wilson's hydrograph was performed by the fuzzy alpha cut method, and it was found that the uncertainty of the parameter x is greater than the uncertainty of the parameters k and m.

  • A new random selection method was used to increase the accuracy of the variable-parameter Muskingum model.

  • The PSO–GA algorithm was used to estimate the Muskingum parameters.

  • The proposed model was applied to three examples, and the results showed that this model has estimated the outflows more accurately in all three examples than the previous methods.

  • The fuzzy alpha cut method was used to quantify the uncertainty of Muskingum parameters.

Flood routing is a procedure to determine the time and magnitude of a flow. The flow hydrograph at a point on a watercourse is determined using known or assumed hydrographs at one or more points upstream. Flood routing is an effective step toward flood prediction and control and can reduce its harmful consequences. Flood routing is classically done in two main ways: hydraulic flood routing and hydrologic flood routing (Barati 2013). Among hydrologic routing methods, the Muskingum method is very popular with researchers and water engineers on account of its ease of use and simplicity. In addition, the Muskingum model makes accurate predictions with fewer data types than common hydraulic methods (Barati 2013). This method can be used to determine flood hydrographs in places with similar morphological aspects (Bazargan & Norouzi 2018), whereas it cannot be employed in cases where significant backwater effects are observed (Chow et al. 1988). Some researchers have considered constant parameters for the Muskingum model, while some have used the variable ones. Clearly, the Muskingum model is called ‘constant-parameter’ as long as the parameters are considered constant during the routing period, while, if these parameters are considered variable during the routing period, the Muskingum model is called a ‘variable-parameter Muskingum model’.

The determination of the parameters of the Muskingum model has been a topic of research for many years. Gill (1978) implemented the segmented curve method for nonlinear flood routing, in which the least squared error was taken advantage of so as to estimate the parameters. Among optimization methods, the evolutionary algorithms have been used to ascertain the parameters of the Muskingum model owing to their high accuracy, as well as their capability of solving multi-variable problems. In this regard, Mohan (1997) employed a genetic algorithm (GA), and Kim et al. (2001) used a harmony search (HS) method to estimate the parameters of the Muskingum model. Chu & Chang (2009) determined the parameters of the model using the particle swarm optimization (PSO) method, whereas Luo & Xie (2010) took advantage of the immune clonal selection algorithm to determine the parameters. The hybrid model, combining HS and PSF methods (parameter-setting-free HS), was introduced by Geem (2011) to deal with complex systems. It is noteworthy to emphasize that, this way, 5,000 functions are optimized in only 1 s. Karahan et al. (2012) proposed a hybrid model of HS and Broyden–Fletcher–Goldfarb–Shanno algorithm methods and compared the results from the proposed model with those from eleven other methods. Hamedi et al. (2016) developed a nonlinear Muskingum model by introducing a parametric initial storage condition and applied the weed optimization algorithm to route the flood. Niazkar & Afzali (2017b) proposed a hybrid model of the modified honey bee mating and generalized reduced gradient (GRG) algorithms to develop a nine constant-parameter Muskingum model. Qiang et al. (2020) and Yuan et al. (2021) utilized the whale optimization algorithm with elite opposition-based learning and Polak–Ribière–Polyak methods, respectively, to route the flood in the Haihe River. Norouzi & Bazargan (2020, 2021) applied the PSO method to optimize the Muskingum parameters of the model. Okkan & Kirdemir (2020) routed four flood hydrographs using PSO to optimize the Levenberg–Marquardt model. Khalifeh et al. (2020b) managed to route the flood in the Kardeh River by the grasshopper optimization algorithm and compared the results from the model with the ones from GA and HS algorithms. Not only have the researchers employed a variety of algorithms, but they have also suggested different structures for the Muskingum model in order to improve the precision of the model in estimating the Muskingum parameters. Khalifeh et al. (2020a) implemented a Muskingum model with seven variables for flood routing in Karun River, in which the parameters were determined by the symbiotic organisms search algorithm. Bozorg-Haddad et al. (2020) utilized a 15-parameter Muskingum equation and employed the Excel Solver in order to optimize the outflow coefficient. They assumed the coefficient as a factor of the inflow rate. The results showed that although this adds to the complexity of the model, it ends in much more precise results.

The variable-parameter Muskingum model has been introduced as an approach to improve the results. The nonlinear variable-parameter Muskingum model was first used in the case of three flood hydrographs by Easa (2013), in which the exponential parameter was set to be variable. Niazkar & Afzali (2017a) and Afzali (2016) split the flood period into two and three sub-regions, where a unique parameter was studied for each. Zhang et al. (2017) used an improved real-coded adaptive GA (RAGA) with an elite strategy for precise parameter estimation of the nonlinear Muskingum model. They divided the flood period into five sub-regions. Kang & Zhou (2018) and Akbari et al. (2020) specified all Muskingum parameters, associated with the flood period, for three sub-regions and showed that the sum squared deviation (SSQ) significantly decreased using the introduced model. Bazargan & Norouzi (2018) determined k, x, and t parameters for the first hydrograph and then calculated the outflows for the second hydrograph based on the parameters calculated for the first one. They divided the flood period into three sub-regions and considered Muskingum parameters variable.

The present study aims to increase the accuracy of the Muskingum model. In this regard, a new method for dividing the flood period into a number of sub-regions in the variable-parameter Muskingum model is proposed. Different methods have so far been used to divide the flood hydrograph period into sub-regions based on the time, or the inflow hydrograph, and the peak flows. In this research, the random selection (RS) method has been utilized to divide the flood period. The efficiency of the RS method has been evaluated by routing three flood hydrographs, including Wilson's hydrograph, O'Donnell hydrograph, and a flood hydrograph in a river with the same morphological aspect as the Karun River, located in Iran. In order to assess the predictive ability of the present model, evaluative measures were used for these hydrographs. The performance of the present model was compared with that of the models from the literature. The values of objective function have been decreased by 5% in the first case study, compared to the nonlinear Muskingum model with eight constant parameters. The values decreased by 61% in the second case study, compared to the nonlinear Muskingum model with four-variable parameters. In the third case study, the objective was optimized by 96%, compared to the Muskingum model with three constant parameters. The Muskingum equations related to all models have been presented in Tables 3, 5 and 7. This does indicate an improvement in the accuracy and efficiency of the proposed model and proves that the RS method can be an attractive alternative. The uncertainty of the Muskingum parameters obtained using the RS method was examined in Wilson's hydrograph by the fuzzy alpha cut (FAC) method. It was found that the uncertainty of parameter x is greater than that of parameters k and m.

The Muskingum model is formed on account of the relationship between the storage, and weighted inflow and outflow. Based on the nonlinear relationship between these two factors, different forms of the nonlinear Muskingum model have been created, which are presented in Niazkar & Afzali (2016).

In this paper, a three-parameter nonlinear Muskingum model, whose Muskingum parameters vary during the routing period, is employed as follows:
(1)
The continuity Equation is written in the following fashion:
(2)
Considering Equations (1) and (2), the outflow is determined by the following equation:
(3)
The rate of storage changes with time is calculated by substituting Equation (3) into Equation (2):
(4)
The storage in the next step is calculated by the following equation:
(5)

The outflow can be determined using Equation (3) in the next step (Kang & Zhou 2018).

Moreover, in this research, the objective is set to be minimizing the sum of the squared errors:
(6)
In addition, the other performance indicators of SAD, VarexQ, EQp, ETp, and MARE were taken advantage of as follows in order to compare the results from the proposed method to the ones from the methods in other works of research:
(7)
(8)
(9)
(10)
(11)
where Op represents the observed peak outflow, and Ôp stands for the calculated peak outflow. Tp and constitute the time of the observed peak outflow and the time of calculated peak outflow, respectively. Moreover, Ōt represents the mean of the observed outflow. The lower the value of EQp, the more accurate the model. Plus, the smaller the value of ETp, the more accurate the prediction of peak discharge. The proximity of the shape and the size of the hydrograph are measured using Equation (11) (Moghaddam et al. 2016).

The choice of the objective function has a direct effect on flood routing results. Many researchers have chosen SSQ as the objective function (Akbari et al. 2020), whereas some others have gone for the SAD (Norouzi & Bazargan 2021). However, both have been chosen as the objective function in some cases (e.g., Orouji et al. 2013), making it a two-objective optimization problem.

Optimization algorithm

Genetic algorithm

The optimization in GA is performed through crossover and mutation, presented in the following equations:
(12)
(13)
and form a pair of population members before the crossover, whereas and are a pair after the crossover. ci is a random number between 0 and 1 (Wu et al. 2015).

The steps of exchanging genetic information are repeated until the algorithm terminates.

PSO algorithm

The PSO is a swarm intelligence algorithm, which is based on the social behavior of birds and fish to find food. First, particles with different velocities and positions are randomly generated in the search space. Each particle, then, updates its position and velocity based on the particle's own experience and the best experience of the population. On condition that , the velocity and position of the particles are updated by Equations (14) and (16) in each iteration (Xu et al. 2019):
(14)
(15)
(16)
where c1 and c2 represent the social parameters 1 and 2, and are both considered 2.05 in this research. In addition, x is the constriction coefficient, which is a function of c1 and c2. Rand(·) acts for a random positive number that is calculated by a uniform distribution between 0 and 1. is the best position particle ‘i’ has ever been in. and are the position and velocity of the ith particle in the tth iteration in the search space. After all, is the best position of the whole population.

PSO–GA algorithm

Gaining the advantages of both algorithms, while avoiding the disadvantages of each, the theory of a new hybrid algorithm was developed. The PSO and GA algorithms have their own advantages and disadvantages, which are used in combination with each other, can lead to increasing the advantages and lowering the disadvantages. In the previous section, the steps of PSO and GA algorithms were separately described. The PSO–GA algorithm optimizes Muskingum parameters with a high convergence rate and a significant accuracy. In this method, Muskingum parameters or so-called ‘the agents’ are randomly generated in the search space at first, and then the steps of the PSO algorithm are executed, through which the velocity and position of the particles are updated. In the following step, the next generation is created by comparing the population's best experience with the particle's best experience. Now, the crossover and mutation are applied to the new population and the members of the new population are sorted according to the objective function. The succeeding generation is produced by comparing the population and particle's best experience. Provided that the termination condition is satisfied, the results are reported as outcomes of the algorithm. Otherwise, the steps of the algorithm continue until the termination condition of the algorithm is reached (Garg 2016). The steps of the hybrid PSO–GA algorithm are illustrated in Figure 1.

Figure 1

Hybrid PSO–GA algorithm (Akbari et al. 2020).

Figure 1

Hybrid PSO–GA algorithm (Akbari et al. 2020).

Close modal
Figure 2

Division of the flood period into three sub-regions and the specification of hydrological parameters for each category.

Figure 2

Division of the flood period into three sub-regions and the specification of hydrological parameters for each category.

Close modal

RS method

Suppose Ii points are to be divided into three categories, the points that are placed in categories 1, 2, and 3 are given the values , , and , respectively. First, the values , , and are randomly generated in the search space. In the next step, points Ii, i= 1, …, J are randomly assigned to each of the categories 1, 2, and 3. At this stage, the number of members in each category is not a fixed number, but the sum of the members of the three categories needs to equal J. Afterwards, the flow routing operation is performed. Now as shown in Figure 2, for each Ii, there is a specific value (k, x, m) depending on which of the categories 1, 2, and 3 is located in, and then the outflows are calculated using it, and the SSQ is obtained. Generating values of , , and are performed using the optimization algorithm and assigning Ii points to categories 1, 2, and 3. This way, the SSQ is calculated. Likewise, the whole process continues until the optimization algorithm termination condition is satisfied.

Uncertainty

Uncertainty analysis depends on many elements, such as knowing the origins of uncertainty and the complexity of the model. In order to determine the uncertainty of the models, the possibility theory is used in the case of exact and quantitative measurements, while the probability theory is taken advantage of for qualitative and rough measurements. Due to the simplification of equations and/or the use of random variables, the quantitative estimation of parameters is always accompanied by uncertainty. Quantifying the uncertainty is one of the ways to determine the uncertainty. It uses the triangular fuzzy membership numbers. In this method, the membership value of the maximum and minimum is assigned 0, while 1 is assigned to the median so as to lower the impact of the scattered data.
(17)
(18)
where s is the intended variable; represents the membership degree of the variable s. s1 is the minimum value, while s2 and s3 represent the median and the maximum values, respectively. Lastly, U is the uncertainty value (Raina & Thomas 2012; Najafi & Hessami 2017).

As the basic concepts of fuzzy logic are so simple, yet very vast, one can return to Varón-Gaviria et al. (2017) for the purpose of further studying.

The routing process was performed on three hydrographs, which were chosen to assess the performance of the RS method.

Case Study 1

Wilson's hydrograph is chosen in this example. Wilson's is a single-peak hydrograph, which has been routed by many researchers. This hydrograph consists of 22 inflow points. In this research, for all models, the flood period is divided into three sub-regions. Each of the inflows is randomly selected and placed in a separate category. However, this RS must be such that the SSQ is minimized. Therefore, the SSQ is determined in each step and is then compared with the SSQ calculated in the previous step. The lower SSQ is reported as the outcome of the model. Table 1 provides the outflows obtained from the method, the values of k, x, and m optimized by the PSO–GA algorithm, as well as the category number, to which each point belongs. The first, second, and third categories include 6, 12, and 4 members, respectively. The values of the objective function and those of the performance indicators are presented in Table 2, in which the results from 10 different methods are compared with each other. Different methods are sorted in lines 1–10 of this table, respecting the SSQ values, i.e., the Easa method has the highest SSQ value, while the lowest value of the SSQ is related to the RS method. The values of are (0.7443, 0.6377, 0.9352), (0.1767, 0.2913, 0.5000), and (1.8232, 1.8254, 1.7805). In addition, the value of SSQ obtained from the proposed RS method is 0.6148. The lowest SSQ value ever obtained is 0.65 using Excel Solver by Bozorg-Haddad et al. (2020). As it turns out, the RS method has reduced the SSQ value by 5%. It is noteworthy to mention that the structure of the Muskingum model by Bozorg-Haddad et al. (2020) consists of eight Muskingum parameters, while in the proposed RS model, the model is made up of three Muskingum parameters. In fact, in order to make a more accurate comparison, it is necessary to compare the results of the proposed model with those from the model by Kang & Zhou (2018), which is similarly formed on three Muskingum parameters. In this model, the SSQ is 2.272, which demonstrates that the RS method has improved the value of the objective function by 73%. This implies the fact that, in both comparisons, the SSQ obtained through the RS method is the lower one. The values of five performance criteria of SAD, EQp, ETp, MARE, and VarexQ were calculated for the RS method and are also presented and compared with the performance indices of nine other methods in Table 2. As it can be inferred from the comparison, the RS method has the highest VarexQ. This index is linked with the closeness of shape and size of hydrograph. The ETp and EQp represent the time and the quantity of the peak discharge, respectively. The calculated ETp = 0 indicates the fact that the time of the peak discharge has not changed or shifted. In addition, the calculated EQp = 0.00085 shows that the quantity of the peak discharge has been forecasted with very little difference. The magnitude of the peak discharge equals 84.5 cms in this model, although it is 85 cms on the observational hydrograph. The MARE value of this model is 0.0034, which is lower than that of the five models listed in Table 2.

Table 1

Comparison of routed outflows obtained from various methods for application of case study 1

Observed data (cms)
Computed outflow (cms)Category number
jTime (h)IjOjRS
22 22 22.0 
23 21 20.8 
12 35 21 21.0 
18 71 26 26.0 
24 103 34 33.8 
30 111 44 44.0 
36 109 55 55.0 
42 100 66 65.7 
48 86 75 75.4 
10 54 71 82 82.0 
11 60 59 85 84.9 
12 66 47 84 83.7 
13 72 39 80 79.9 
14 78 32 73 73.0 
15 84 28 64 64.2 
16 90 24 54 54.1 
17 96 22 44 43.9 
18 102 21 36 36.2 
19 108 20 30 29.7 
20 114 19 25 25.0 
21 120 19 22 21.8 
22 126 18 19 19.2 
k 0.7443 0.6377 0.9352   
x 0.1767 0.2913 0.5000   
m 1.8232 1.8254 1.7805   
Observed data (cms)
Computed outflow (cms)Category number
jTime (h)IjOjRS
22 22 22.0 
23 21 20.8 
12 35 21 21.0 
18 71 26 26.0 
24 103 34 33.8 
30 111 44 44.0 
36 109 55 55.0 
42 100 66 65.7 
48 86 75 75.4 
10 54 71 82 82.0 
11 60 59 85 84.9 
12 66 47 84 83.7 
13 72 39 80 79.9 
14 78 32 73 73.0 
15 84 28 64 64.2 
16 90 24 54 54.1 
17 96 22 44 43.9 
18 102 21 36 36.2 
19 108 20 30 29.7 
20 114 19 25 25.0 
21 120 19 22 21.8 
22 126 18 19 19.2 
k 0.7443 0.6377 0.9352   
x 0.1767 0.2913 0.5000   
m 1.8232 1.8254 1.7805   
Table 2

Performance criteria for the application of case study 1

ModelAlgorithmSSQSADEQpETpMAREVarexQ
Easa (2013) a GA + GRG 24.881 20.71 0.0067 0.0248 99.80 
Moghaddam et al. (2016) b PSO 8.820 9.77 0.0001 0.015 99.93 
Zhang et al. (2017) c RAGA 5.73 8.58 0.0012 0.012 99.95 
Farahani et al. (2018) d Kidney 5.65 7.140 0.0004 0.0241 – 
Farahani et al. (2019) e Shark 5.124 8.1124 0.0002 0.0251 99.96 
Niazkar & Afzali (2017a) f MHBMO 4.043 5.967 0.0004 0.01 99.97 
Kang & Zhou (2018) g Excel Solver 2.272 5.179 NR NR NR – 
Akbari et al. (2020) h PSO–GA 1.0921 3.669 0.0003 0.0018 99.99 
Bozorg-Haddad et al. (2020) k Excel Solver 0.65 2.66 0.0015 NR NR 
RS (This study) PSO–GA 0.6148 2.8205 0.00085 0.0034 99.995 
ModelAlgorithmSSQSADEQpETpMAREVarexQ
Easa (2013) a GA + GRG 24.881 20.71 0.0067 0.0248 99.80 
Moghaddam et al. (2016) b PSO 8.820 9.77 0.0001 0.015 99.93 
Zhang et al. (2017) c RAGA 5.73 8.58 0.0012 0.012 99.95 
Farahani et al. (2018) d Kidney 5.65 7.140 0.0004 0.0241 – 
Farahani et al. (2019) e Shark 5.124 8.1124 0.0002 0.0251 99.96 
Niazkar & Afzali (2017a) f MHBMO 4.043 5.967 0.0004 0.01 99.97 
Kang & Zhou (2018) g Excel Solver 2.272 5.179 NR NR NR – 
Akbari et al. (2020) h PSO–GA 1.0921 3.669 0.0003 0.0018 99.99 
Bozorg-Haddad et al. (2020) k Excel Solver 0.65 2.66 0.0015 NR NR 
RS (This study) PSO–GA 0.6148 2.8205 0.00085 0.0034 99.995 

a.

b,d,e.

c.

g,f.

h.

k.

The uncertainty of hydrologic parameters in this case study is determined by the FAC method. Among past studies, which have been presented in Table 2, Easa (2013), Niazkar & Afzali (2017a), Kang & Zhou (2018), and the present study have employed the three-parameter nonlinear Muskingum model. The membership values were calculated based on Equation (17) for m, x, and k parameters, individually and presented in Table 3. The maximum and minimum values of each parameter were assigned 0, while the value 1 was assigned to the median. Triangular membership function graphs have been drawn for each parameter in Figure 3. The highest membership value is 1 and marks the lowest uncertainty. However, the lowest membership value, which is 0, is related to the highest uncertainty. The membership degree of the values between the minimum and the median, and the median and the maximum changes linearly from 0 to 1. This proves that the uncertainty of these values ranges from 0 to 100%. In the triangular membership function graph, the closer the supporting surface to , the greater the uncertainty values, as it is inferred from width of the membership function, whereas the closer it is to the , the lower the uncertainty of the data. The maximum values of parameters k, x, and m are 2.021, 1.8254, and 1.906, respectively, while the minimum values of parameters k, x, and m are 0.1767, 0.266, and 0.5. In addition, 1.1402, 0.301, and 1.677 are the values of the median of k, x, and m, respectively. The uncertainty value in the FAC method equals the ratio of the width of the supporting surface at section to the point at which the membership function value is 1. It is determined using Equation (18). The horizontal lines on Figure 3 illustrate . The uncertainty of the parameters x, k, and m have been determined to be 4.66, 1.46, and 0.75, respectively. Therefore, it is concluded that the uncertainty of the parameter x in flood routing is greater than that of the parameters m and k among all methods listed in Table 3. The scatter plots of observed and routed outflows are given in Figure 4. As can be seen from the fit line equations (assume that the equation is y = aox + a1) in the scatterplots, the ao and a1 coefficients are closer to the 1 and 0 with a high determination coefficient value (R2).

Figure 3

Triangle membership functions for hydrologic variables.

Figure 3

Triangle membership functions for hydrologic variables.

Close modal
Figure 4

Scatter plot of observed and routed outflows in three case studies.

Figure 4

Scatter plot of observed and routed outflows in three case studies.

Close modal

Case study 2

The second selected hydrograph in this paper is the one named O'Donnell (1985), which is a non-smooth hydrograph. This hydrograph belongs to a river in England. The time step in this hydrograph is 6 h and the flood period is divided into three sub-regions, and the classification criterion is the minimum value of the objective function. Table 4 shows the routed flows and the values of Muskingum parameters obtained from the RS method. In this table, the category of the inflow is given in the sixth column. In addition, each inflow can be assigned to one of three categories. The first category has six members, and each of the second and third categories has 14 members. The values of are (1.4589, 0.4568, 1.4620), (1.068, 0.5350, 1.4694), and (1.5888, −0.9978, 1.2681), while the value of SSQ is determined to be 3733.5. In order to make a better comparison, the SSQ obtained from the other eight methods and corresponding performance indices are presented in Table 5. As can be seen, the obtained SSQ in this paper is the lowest among all these methods. It is also inferred that the RS method has reduced SSQ, whose lowest value is 9654.5, by 61.3%. Moreover, employing the RS method ends in minimum SSQ, ETp and VarexQ values among all the models in Table 5, proving the proposed superiority of the model to them in that respect. The value of the peak discharge has been determined to be 964 cms in this method. However, the observational value of the peak discharge is 969 cms. As it is shown in Table 5. the EQp is 0.0048, and it is the lowest among all methods mentioned in this table. As a result, it is clear that the RS method predicts the peak discharge very accurately. In the present study, MARE, which equals 0.0391, is the least of all of the methods listed in Table 5. ETp was determined to be 0, which proves that the model predicted the time of the peak discharge flawlessly. In Figure 4, the routed outflows are compared with the observed outflows. The coefficient of determination (R2) values are high.

Table 3

Membership number for non-deterministic variables

Methodkxm
Easa (2013)  0.483 0.32 0.266 1.906 
0.483 0.32 0.266 1.883 0.10 
0.483 0.32 0.266 1.884 0.1 
0.483 0.32 0.266 1.889 0.07 
0.483 0.32 0.266 1.880 0.11 
Niazkar & Afzali (2017a)  1.0804 0.25 0.503 0.87 1.661 0.99 
1.686 0.38 0.382 0.95 1.627 0.96 
1.536 0.55 0.300 0.97 1.676 
Kang & Zhou (2018)  2.021 0.495 0.87 1.678 
1.907 0.13 0.396 0.94 1.632 0.96 
1.743 0.32 0.302 1.646 0.97 
This study 0.7443 0.59 0.6377 0.78 0.9352 0.37 
0.1767 0.2913 0.72 0.5 
1.8232 0.22 1.8254 1.7805 0.55 
Minimum 0.1767 – 0.2660 – 0.5 – 
Maximum 2.0210 – 1.8254 – 1.9060 – 
Median 1.1402 – 0.3010 – 1.677 – 
Uncertainty 1.4558 – 4.6627 – 0.7546 – 
Methodkxm
Easa (2013)  0.483 0.32 0.266 1.906 
0.483 0.32 0.266 1.883 0.10 
0.483 0.32 0.266 1.884 0.1 
0.483 0.32 0.266 1.889 0.07 
0.483 0.32 0.266 1.880 0.11 
Niazkar & Afzali (2017a)  1.0804 0.25 0.503 0.87 1.661 0.99 
1.686 0.38 0.382 0.95 1.627 0.96 
1.536 0.55 0.300 0.97 1.676 
Kang & Zhou (2018)  2.021 0.495 0.87 1.678 
1.907 0.13 0.396 0.94 1.632 0.96 
1.743 0.32 0.302 1.646 0.97 
This study 0.7443 0.59 0.6377 0.78 0.9352 0.37 
0.1767 0.2913 0.72 0.5 
1.8232 0.22 1.8254 1.7805 0.55 
Minimum 0.1767 – 0.2660 – 0.5 – 
Maximum 2.0210 – 1.8254 – 1.9060 – 
Median 1.1402 – 0.3010 – 1.677 – 
Uncertainty 1.4558 – 4.6627 – 0.7546 – 
Table 4

Comparison of routed outflows obtained from various methods for the application of case study 2

Observed data (cms)
Computed outflow (cms)Category number
jTime (h)IjOjRS
154 102 102 
150 140 143 
12 219 169 165 
18 182 190 185 
24 182 209 189 
30 192 218 229 
36 165 210 209 
42 150 194 184 
48 128 172 162 
10 54 168 149 138 
11 60 260 136 130 
12 66 471 228 230 
13 72 717 303 286 
14 78 1,092 366 345 
15 84 1,145 456 478 
16 90 600 615 615 
17 96 365 830 827 
18 102 277 969 964 
19 108 277 665 657 
20 114 187 519 539 
21 120 161 444 453 
22 126 143 321 295 
23 132 126 208 207 
24 138 115 176 170 
25 144 102 148 142 
26 150 93 125 117 
27 156 88 114 113 
28 162 82 106 100 
29 168 76 97 89 
30 174 73 89 81 
31 180 70 81 76 
32 186 67 76 72 
33 192 63 71 68 
34 198 59 66 66 
k 1.4589 1.068 1.5888   
x 0.4568 0.5350 −0.9978   
m 1.4620 1.4694 1.2681   
Observed data (cms)
Computed outflow (cms)Category number
jTime (h)IjOjRS
154 102 102 
150 140 143 
12 219 169 165 
18 182 190 185 
24 182 209 189 
30 192 218 229 
36 165 210 209 
42 150 194 184 
48 128 172 162 
10 54 168 149 138 
11 60 260 136 130 
12 66 471 228 230 
13 72 717 303 286 
14 78 1,092 366 345 
15 84 1,145 456 478 
16 90 600 615 615 
17 96 365 830 827 
18 102 277 969 964 
19 108 277 665 657 
20 114 187 519 539 
21 120 161 444 453 
22 126 143 321 295 
23 132 126 208 207 
24 138 115 176 170 
25 144 102 148 142 
26 150 93 125 117 
27 156 88 114 113 
28 162 82 106 100 
29 168 76 97 89 
30 174 73 89 81 
31 180 70 81 76 
32 186 67 76 72 
33 192 63 71 68 
34 198 59 66 66 
k 1.4589 1.068 1.5888   
x 0.4568 0.5350 −0.9978   
m 1.4620 1.4694 1.2681   

Case study 3

The inflow and outflow hydrographs, which are selected in this example, are related to the flood of the period between February 26, 2012 and March 1, 2012. The selected watershed is very similar to the Karun River watershed in terms of the morphology of the two rivers and was taken advantage of by Bazargan & Norouzi (2018) to perform the flood routing. The recorded data are extracted from the article by Bazargan & Norouzi (2018) and are given in Table 6. This hydrograph is also a single-peak hydrograph. The flood period was divided into three sub-regions, where 121 members were assigned to the categories so as to minimize the value of SSQ. Table 6 illustrates how the dataset is distributed throughout categories. The values of are (6.5029, 6.9192, 2.7281), (0.5015, 0.4797, 0.5480), and (1.0966, 1.0896, 1.2373), and the SSQ corresponding to these Muskingum parameters is equal to 2192.8. The value of SSQ from the model by Bazargan & Norouzi (2018) is 51,301, where the objective function was set to be SAD. The values of performance indicators calculated using the RS method are SAD = 420.6, EQp = 0.0081, ETp = 4, MARE = 0.0079, and VarexQ = 99.81, which are given in Table 7. It should be noted that the maximum outflow is 619 in the observed outflow hydrograph, which occurred at nine time points. Indeed, ETp being equal to 4 and the dislocation of the maximum by four units will not be seen as a weakness of the model owing to the fact that the flood has a maximum value at nine time points and, above all, the peak value has not been shifted. The calculated EQp = 0.0081 shows that the quantity of the peak discharge is close to the observed peak discharge. The scatter plots of observed outflows and routed outflows are shown in Figure 4. In order to study more about the number of sub-regions, the flood period is also divided into two and four sub-regions. For L = 2 and L = 4, the values of the objective function are equal to 3565.8 and 2012.2, respectively. The results showed that increasing the number of sub-regions from 2 to 4 decreases the value of the objective function.

Table 5

Performance criteria for the application of case study 2

ModelAlgorithmSSQSADEQpETpMAREVarexQ
Easa (2013) a GA + GRG 35064 NR NR NR NR NR 
Easa (2014) b GA + GRG 32299.2 743.32 0.078 0.10 98.05 
Moghaddam et al. (2016) c PSO 31099.5 695.77 0.090 0.09 98.12 
Farahani et al. (2018) d Kidney 16121.2 120 0.002 0.01 NR 
Farahani et al. (2019) e Shark 17121.2 123 0.002 0.01 99.11 
Kang & Zhou (2018) f Excel Solver 10368 NR NR NR NR NR 
Akbari et al. (2020) g PSO–GA 9654.5 392.47 0.0183 0.06 99.4 
Bozorg-Haddad et al. (2020) h Excel Solver 19953 621 0.071 NR NR 
RS (This study) PSO–GA 3733.5 269.87 0.0048 0.0391 99.77 
ModelAlgorithmSSQSADEQpETpMAREVarexQ
Easa (2013) a GA + GRG 35064 NR NR NR NR NR 
Easa (2014) b GA + GRG 32299.2 743.32 0.078 0.10 98.05 
Moghaddam et al. (2016) c PSO 31099.5 695.77 0.090 0.09 98.12 
Farahani et al. (2018) d Kidney 16121.2 120 0.002 0.01 NR 
Farahani et al. (2019) e Shark 17121.2 123 0.002 0.01 99.11 
Kang & Zhou (2018) f Excel Solver 10368 NR NR NR NR NR 
Akbari et al. (2020) g PSO–GA 9654.5 392.47 0.0183 0.06 99.4 
Bozorg-Haddad et al. (2020) h Excel Solver 19953 621 0.071 NR NR 
RS (This study) PSO–GA 3733.5 269.87 0.0048 0.0391 99.77 

a.

b,c,d,e.

f.

g.

h.

Table 6

Comparison of routed outflows obtained from various methods for the application of case study 3

Observed data (cms)
Computed outflow (cms)Category number
jTime (h)IjOjRS
376 328 328 
381 329 333.7 
386 329 336.7 
391 329 323.0 
396 333 330.1 
401 336 336.7 
406 340 343.1 
411 340 339.7 
416 348 346.9 
10 429 352 351.1 
11 10 443 356 362.7 
12 11 456 363 364.9 
13 12 469 367 369.3 
14 13 482 371 374.9 
15 14 495 379 381.1 
16 15 508 387 389.4 
17 16 521 396 397.9 
18 17 529 404 408.7 
19 18 537 412 411.6 
20 19 545 421 422.8 
21 20 554 429 433.1 
22 21 562 443 443.1 
23 22 570 451 453.1 
24 23 578 465 462.8 
25 24 586 474 472.2 
26 25 594 488 493.0 
27 26 602 497 501.9 
28 27 610 512 510.7 
29 28 619 521 519.0 
30 29 627 536 539.7 
31 30 635 546 546.7 
32 31 643 556 553.9 
33 32 651 566 561.2 
34 33 649 577 572.1 
35 34 647 582 585.3 
36 35 644 592 596.5 
37 36 642 598 605.7 
38 37 636 603 599.0 
39 38 629 608 609.0 
40 39 623 614 616.5 
41 40 616 614 622.1 
42 41 610 614 609.2 
43 42 604 619 615.1 
44 43 598 619 619.4 
45 44 592 619 622.1 
46 45 586 619 623.6 
47 46 580 619 624.0 
48 47 574 619 623.5 
49 48 568 619 622.2 
50 49 562 619 620.2 
51 50 556 619 617.6 
52 51 550 614 614.4 
53 52 544 614 611.3 
54 53 537 614 607.4 
55 54 531 614 618.8 
56 55 525 608 610.1 
57 56 519 608 601.8 
58 57 510 608 594.8 
59 58 501 592 588.8. 
60 59 491 577 576.5 
61 60 482 561 559.5 
62 61 475 556 554.0 
63 62 469 546 547.2 
64 63 462 536 540.4 
65 64 455 526 533.9 
66 65 452 526 526.1 
67 66 450 521 517.1 
68 67 447 517 525.3 
69 68 444 512 515.5 
70 69 441 507 506.7 
71 70 439 502 505.5 
72 71 436 497 497.7 
73 72 433 493 491.1 
74 73 430 487 485.1 
75 74 427 482 479.5 
76 75 424 476 474.2 
77 76 422 470 469.0 
78 77 419 464 464.0 
79 78 416 459 459.9 
80 79 413 453 455.9 
81 80 410 447 452.1 
82 81 407 443 448.4 
83 82 403 438 445.2 
84 83 400 434 442.0 
85 84 396 429 438.7 
86 85 392 424 419.5 
87 86 388 419 418.2 
88 87 384 413 416.5 
89 88 380 408 414.5 
90 89 379 406 411.1 
91 90 378 403 406.8 
92 91 377 401 403.0 
93 92 376 398 399.6 
94 93 374 396 397.0 
95 94 373 393 394.6 
96 95 372 391 391.9 
97 96 371 388 389.5 
98 97 370 385 387.4 
99 98 369 383 385.4 
100 99 368 380 383.5 
101 100 367 378 381.8 
102 101 365 375 380.6 
103 102 364 372 379.3 
104 103 363 370 377.6 
105 104 362 367 364.3 
106 105 361 366 364.6 
107 106 359 365 365.1 
108 107 358 364 365.3 
109 108 356 363 365.2 
110 109 354 362 365.2 
111 110 353 361 364.7 
112 111 351 360 363.9 
113 112 349 359 363.4 
114 113 349 357 362.0 
115 114 349 355 360.0 
116 115 349 353 358.4 
117 116 350 352 356.7 
118 117 350 350 355.3 
119 118 350 348 354.5 
120 119 350 346 353.8 
121 120 350 344 353.2 
k 6.5029 0.5015 1.0966   
x 6.9192 0.4797 1.0896   
m 2.7281 0.5480 1.2373   
Observed data (cms)
Computed outflow (cms)Category number
jTime (h)IjOjRS
376 328 328 
381 329 333.7 
386 329 336.7 
391 329 323.0 
396 333 330.1 
401 336 336.7 
406 340 343.1 
411 340 339.7 
416 348 346.9 
10 429 352 351.1 
11 10 443 356 362.7 
12 11 456 363 364.9 
13 12 469 367 369.3 
14 13 482 371 374.9 
15 14 495 379 381.1 
16 15 508 387 389.4 
17 16 521 396 397.9 
18 17 529 404 408.7 
19 18 537 412 411.6 
20 19 545 421 422.8 
21 20 554 429 433.1 
22 21 562 443 443.1 
23 22 570 451 453.1 
24 23 578 465 462.8 
25 24 586 474 472.2 
26 25 594 488 493.0 
27 26 602 497 501.9 
28 27 610 512 510.7 
29 28 619 521 519.0 
30 29 627 536 539.7 
31 30 635 546 546.7 
32 31 643 556 553.9 
33 32 651 566 561.2 
34 33 649 577 572.1 
35 34 647 582 585.3 
36 35 644 592 596.5 
37 36 642 598 605.7 
38 37 636 603 599.0 
39 38 629 608 609.0 
40 39 623 614 616.5 
41 40 616 614 622.1 
42 41 610 614 609.2 
43 42 604 619 615.1 
44 43 598 619 619.4 
45 44 592 619 622.1 
46 45 586 619 623.6 
47 46 580 619 624.0 
48 47 574 619 623.5 
49 48 568 619 622.2 
50 49 562 619 620.2 
51 50 556 619 617.6 
52 51 550 614 614.4 
53 52 544 614 611.3 
54 53 537 614 607.4 
55 54 531 614 618.8 
56 55 525 608 610.1 
57 56 519 608 601.8 
58 57 510 608 594.8 
59 58 501 592 588.8. 
60 59 491 577 576.5 
61 60 482 561 559.5 
62 61 475 556 554.0 
63 62 469 546 547.2 
64 63 462 536 540.4 
65 64 455 526 533.9 
66 65 452 526 526.1 
67 66 450 521 517.1 
68 67 447 517 525.3 
69 68 444 512 515.5 
70 69 441 507 506.7 
71 70 439 502 505.5 
72 71 436 497 497.7 
73 72 433 493 491.1 
74 73 430 487 485.1 
75 74 427 482 479.5 
76 75 424 476 474.2 
77 76 422 470 469.0 
78 77 419 464 464.0 
79 78 416 459 459.9 
80 79 413 453 455.9 
81 80 410 447 452.1 
82 81 407 443 448.4 
83 82 403 438 445.2 
84 83 400 434 442.0 
85 84 396 429 438.7 
86 85 392 424 419.5 
87 86 388 419 418.2 
88 87 384 413 416.5 
89 88 380 408 414.5 
90 89 379 406 411.1 
91 90 378 403 406.8 
92 91 377 401 403.0 
93 92 376 398 399.6 
94 93 374 396 397.0 
95 94 373 393 394.6 
96 95 372 391 391.9 
97 96 371 388 389.5 
98 97 370 385 387.4 
99 98 369 383 385.4 
100 99 368 380 383.5 
101 100 367 378 381.8 
102 101 365 375 380.6 
103 102 364 372 379.3 
104 103 363 370 377.6 
105 104 362 367 364.3 
106 105 361 366 364.6 
107 106 359 365 365.1 
108 107 358 364 365.3 
109 108 356 363 365.2 
110 109 354 362 365.2 
111 110 353 361 364.7 
112 111 351 360 363.9 
113 112 349 359 363.4 
114 113 349 357 362.0 
115 114 349 355 360.0 
116 115 349 353 358.4 
117 116 350 352 356.7 
118 117 350 350 355.3 
119 118 350 348 354.5 
120 119 350 346 353.8 
121 120 350 344 353.2 
k 6.5029 0.5015 1.0966   
x 6.9192 0.4797 1.0896   
m 2.7281 0.5480 1.2373   
Table 7

Performance criteria for the application of case study 3

ModelAlgorithmSSQSADEQpETpMAREVarexQ
Bazargan & Norouzi (2018) a PSO 51301 2069.4 0.0141 0.04 95.7 
RS (This study)b PSO–GA 2192.8 420.6 0.0081 0.0079 99.81 
ModelAlgorithmSSQSADEQpETpMAREVarexQ
Bazargan & Norouzi (2018) a PSO 51301 2069.4 0.0141 0.04 95.7 
RS (This study)b PSO–GA 2192.8 420.6 0.0081 0.0079 99.81 

a.

b.

Muskingum parameters can be considered constant or variable during the flood routing process. Although the complexity of the model increases in variable-parameter flood routing, the accuracy of flood prediction increases compared to the constant-parameter routing methods. In this study, three variable parameters were considered for flood routing in the Muskingum method. The PSO–GA optimization algorithm was employed to estimate the parameters of the model. A new technique has been proposed for the classification of the sub-regions in order to increase the accuracy and improve the results of the variable-parameter Muskingum model. This technique is formed on the RS of Muskingum parameters. This way, the flood period is divided into L= 3 sub-regions. The sub-regions are divided in such a way that the value of the objective function is minimized. As a result, parameter selection is made randomly, and each of the inflows is placed in these categories so that the objective function is minimized. It is noteworthy that the division intervals of k, x, and m should be chosen differently to improve the results. In other words, dividing method used for sub-periods and also the number of sub-periods may differ for Muskingum parameters. For example, for k, the flood period is divided into nk sub-period, for x, it is divided into nx sub-period, and for m, into nm sub-period. It is a better way to increase the accuracy of the method. The RS method was successfully tested on three models, in which the results of the RS flood routing method are compared with those of previous studies and the performance indicators are calculated. In addition, the objective function is the SSQ in this paper. The results show that the value of the objective function in the first, second, and the third case studies is improved by 5, 61, and 96% compared to the eight-parameter Muskingum model, four-variable-parameter Muskingum model, and three-parameter Muskingum model, respectively.

The number of previous studies indicates that increasing the accuracy of the Muskingum model has been the concern of numerous researchers. Since this method is easy to apply and also increases the accuracy of variable-parameter Muskingum models, this method can be used in all variable-parameter Muskingum models to increase their accuracy. In conclusion, the RS method is recommended in order to divide the flood period. The FAC method was used in order to investigate the uncertainty of the Muskingum parameters. This method has been used to quantify the uncertainty of the Muskingum parameters of Wilson's hydrograph. In this method, the Muskingum parameters of the three-parameter Muskingum model have been taken advantage of. The results show that the uncertainty of the parameter x is greater than that of parameters k and m.

All relevant data are included in the paper or its Supplementary Information.

Afzali
S. H.
2016
Variable-parameter Muskingum model
.
Iranian Journal of Science and Technology – Transactions of Civil Engineering
40
(
1
),
59
68
.
Bozorg-Haddad
O.
,
Mohammad-Azari
S.
,
Hamedi
F.
,
Pazoki
M.
&
Loáiciga
H. A.
2020
Application of a new hybrid non-linear Muskingum model to flood routing
. In
Proceedings of the Institution of Civil Engineers – Water Management
.
Thomas Telford Ltd
. Vol.
173
, No.
3
, pp.
109
120
.
Chow
V. T.
,
Maidment
D. R.
&
Mays
L. W.
1988
Applied Hydrology
.
McGraw-Hill
,
New York
, p.
572
.
Chu
H. J.
&
Chang
L. C.
2009
Applying particle swarm optimization to parameter estimation of the nonlinear Muskingum model
.
Journal of Hydrologic Engineering
14
(
9
),
1024
1027
.
Easa
S. M.
2013
Improved nonlinear Muskingum model with variable exponent parameter
.
Journal of Hydrologic Engineering – ASCE
18
(
22
),
1790
1794
.
Easa
S. M.
2014
New and improved four-parameter non-linear Muskingum model
.
Proceedings of the Institution of Civil Engineers
167
(
5
),
288
298
.
Farahani
N. N.
,
Farzin
S.
&
Karami
H.
2018
Flood routing by Kidney algorithm and Muskingum model
.
Natural Hazards
1
19
.
doi:10.1007/s11069-018-3482-x
.
Farahani
N.
,
Karami
H.
,
Farzin
S.
,
Ehteram
M.
,
Kisi
O.
&
El Shafie
A.
2019
A new method for flood routing utilizing four-parameter nonlinear Muskingum and shark algorithm
.
Water Resources Management
33
(
14
),
4879
4893
.
Garg
H.
2016
A hybrid PSO-GA algorithm for constrained optimization problems
.
Applied Mathematics and Computation
274
,
292
305
.
Gill
M. A.
1978
Flood routing by the Muskingum method
.
Journal of Hydrology
36
(
3–4
),
353
363
.
Hamedi
F.
,
Bozorg-Haddad
O.
,
Pazoki
M.
,
Asgari
H. R.
,
Parsa
M.
&
Loaiciga
H. A.
2016
Parameter estimation of extended nonlinear Muskingum models with the weed optimization algorithm
.
Journal of Irrigation and Drainage Engineering
142
(
12
),
1
11
.
Kang
L.
&
Zhou
L.
2018
Parameter estimation of variable-parameter nonlinear Muskingum model using excel solver
.
IOP Conference Series: Earth and Environmental Science
121
(
5
),
052047
.
Karahan
H.
,
Gurarslan
G.
&
Geem
Z. W.
2012
Parameter estimation of the nonlinear Muskingum flood-routing model using a hybrid harmony search algorithm
.
Journal of Hydrologic Engineering
18
(
3
),
352
360
.
Khalifeh
S.
,
Esmaili
K.
,
Khodashenas
S. R.
&
Khalifeh
V.
2020a
Estimation of nonlinear parameters of the type 5 Muskingum model using SOS algorithm
.
MethodsX
7
,
101040
.
Khalifeh
S.
,
Esmaili
K.
,
Khodashenas
S.
&
Akbarifard
S.
2020b
Data on optimization of the non-linear Muskingum flood routing in Kardeh River using Goa algorithm
.
Data in Brief
30
,
105398
.
Kim
J. H.
,
Geem
Z. W.
&
Kim
E. S.
2001
Parameter estimation of the nonlinear Muskingum model using harmony search
.
JAWRA – Journal of the American Water Resources Association
37
(
5
),
1131
1138
.
Moghaddam
A.
,
Behmanesh
J.
&
Farsijani
A.
2016
Parameters estimation for the new four-parameter nonlinear Muskingum model using the particle swarm optimization
.
Water Resources Management
30
(
7
),
2143
2160
.
Mohan
S.
1997
Parameter estimation of nonlinear Muskingum models using genetic algorithm
.
Journal of Hydraulic Engineering
123
(
2
),
137
142
.
Niazkar
M.
&
Afzali
S. H.
2017a
New nonlinear variable-parameter Muskingum models
.
KSCE Journal of Civil Engineering
21
(
7
),
2958
2967
.
Norouzi
H.
&
Bazargan
J.
2020
Investigation of effect of optimal time interval on the linear Muskingum method using particle swarm optimization algorithm
.
Journal of Applied Research in Water and Wastewater
7
(
2
),
152
156
.
Norouzi
H.
&
Bazargan
J.
2021
Effects of uncertainty in determining the parameters of the linear Muskingum method using the particle swarm optimization (PSO) algorithm
.
Journal of Water and Climate Change
12
(
5
),
2055
2067
.
O'Donnell
T.
1985
A direct three-parameter Muskingum procedure incorporating lateral inflow
.
Hydrological Sciences Journal
30
(
4
),
479
496
.
Orouji
H.
,
Bozorg-Haddad
O.
,
Fallah-Mehdipour
E.
&
Marino
M. A.
2013
Estimation of Muskingum parameter by meta-heuristic algorithms
.
Proceedings of the Institution of Civil Engineers
166
(
6
),
315
324
.
Qiang
Z.
,
Qiaoping
F.
,
Xingjun
H.
&
Jun
L.
2020
Parameter estimation of Muskingum model based on whale optimization algorithm with elite opposition-based learning
. In
IOP Conference Series: Materials Science and Engineering
.
IOP Publishing
. Vol.
780
, No.
2
, p.
022013
.
Varón-Gaviria
C. A.
,
Barbosa-Fontecha
J. L.
&
Figueroa-García
J. C.
2017
Fuzzy uncertainty in random variable generation: an α-cut approach
. In
International Conference on Intelligent Computing
.
Springer
,
Cham
, pp.
264
273
.
Wilson
E. M.
1974
Engineering Hydrology
.
Macmillan Education LTD
,
Hampshire
,
UK
.
Xu
G.
,
Cui
Q.
,
Shi
X.
,
Ge
H.
,
Zhan
Z. H.
,
Lee
H. P.
&
Wu
C.
2019
Particle swarm optimization based on dimensional learning strategy
.
Swarm and Evolutionary Computation
45
,
33
51
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).