A new method for dividing flood period in the variable-parameter Muskingum models

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 subperiods. 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.


INTRODUCTION
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 can't 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 In this paper, a three-parameter nonlinear Muskingum model, whose Muskingum parameters vary during the routing period, is employed as follows: The continuity Equation is written in the following fashion: Considering Equations (1) and (2), the outflow is determined by the following equation: The rate of storage changes with time is calculated by substituting Equation (3) into Equation (2): The storage in the next step is calculated by the following equation: 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: 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: where O p represents the observed peak outflow, and Ô p stands for the calculated peak outflow. T p andT p 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 EQ p , the more accurate the model. Plus, the smaller the value of ET p , 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.

Genetic algorithm
The optimization in GA is performed through crossover and mutation, presented in the following equations: X i t and X iþ1 t form a pair of population members before the crossover, whereas X i tþ1 and X iþ1 tþ1 are a pair after the crossover. c i 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 jv t i j V max , the velocity and position of the particles are updated by Equations (14) and (16) in each iteration (Xu et al. 2019): x ¼ 2 where c 1 and c 2 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 c 1 and c 2 . Rand(·) acts for a random positive number that is calculated by a uniform distribution between 0 and 1. pbest i t is the best position particle 'i' has ever been in. x i t and v i t are the position and velocity of the ith particle in the tth iteration in the search space. After all, gbest i t 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.

RS method
Suppose I i points are to be divided into three categories, the points that are placed in categories 1, 2, and 3 are given the values (k 1 , x 1 , m 1 ), (k 2 , x 2 , m 2 ), and (k 3 , x 3 , m 3 ), respectively. First, the values (k 1 , x 1 , m 1 ), (k 2 , x 2 , m 2 ), and (k 3 , x 3 , m 3 ) are randomly generated in the search space. In the next step, points I i , 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, for each I i , 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 (k 1 , x 1 , m 1 ), (k 2 , x 2 , m 2 ), and (k 3 , x 3 , m 3 ) are performed using the optimization algorithm and assigning I i 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.
where s is the intended variable; m(s) represents the membership degree of the variable s. s 1 is the minimum value, while s 2 and s 3 represent the median and the maximum values, respectively. Lastly, U is the uncertainty value (

RESULTS AND DISCUSSION
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 (k i , x i , m i ) 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 ET p and EQ p represent the time and the quantity of the peak discharge, respectively. The calculated ET p ¼ 0 indicates the fact that the time of the peak discharge has not changed or shifted. In addition, the calculated EQ p ¼ 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.
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 (m) 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 Uncorrected Proof 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 / ¼ 0, the greater the uncertainty values, as it is inferred from width of the membership function, whereas the closer it is to the / ¼ 1, 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  Uncorrected Proof 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 / ¼ 0:1 to the point at which the membership function value is 1. It is determined using Equation (18). The horizontal lines on Figure 3 illustrate / ¼ 0:1. 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 ¼ a o x þ a 1 ) in the scatterplots, the a o and a 1 coefficients are closer to the 1 and 0 with a high determination coefficient value (R 2 ).   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 Uncorrected Proof 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 (k i , x i , m i ) 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 it 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, ET p 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 EQ p 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. ET p 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 (R 2 ) values are high.

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 7. It should be noted that the maximum outflow is 619 in the observed outflow hydrograph, which occurred at nine time points. Indeed, ET p 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 EQ p ¼ 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

CONCLUSIONS
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 constantparameter 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 n k sub-period, for x, it is divided into n x subperiod, and for m, into n m 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.

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