Kinematic reverse flood routing in natural rivers using stage data

In many developing countries, due to economic constraints, a single station on a river reach is often equipped to record flow variables. On the other hand, hydrographs at the upstream sections may also be needed for especially assessing flooded areas. The upstream flow hydrograph prediction is called the reverse flood routing. There are some reverse flood routing pocedures requiring sophisticated methods together with substantial data requirements. This study proposes a new reverse flood routing procedure, based upon the simple kinematic wave (KW) equation, requiring only easily measurable downstream stage data. The KW equation is first averaged along a channel length at a fixed time, t, assuming that channel width is spatially constant, and then the spatially averaged equation is averaged in time, Δt. The temporally averaged terms are approximated as the arithmetical mean of the corresponding terms evaluated at time t and t + Δt. The Chezy roughness equation is employed for flow velocity, and the upstream flow stage hydrograph is assumed be described by a two parameter gamma distribution (Pearson Type III). The spatially averaged mean flow depth and lateral flow are related to the downstream flow stage. The resulting routing equation is thus obtained as a function of only downstream flow stage, meaning that the method mainly requires measurements of downstream flow stage data besides the mean values of channel length, channel width, roughness coefficient and bed slope. The optimal values of the parameters of reverse flood routing are obtained using the genetic algorithm. The calibration of the model is accomplished by using the measured downstream hydrographs. The validation is performed by comparing the model-generated upstream hydrographs against the measured upstream hydrographs. The proposed model is applied to generate upstream hydrographs at four different river reaches of Tiber River, located in central Italy. The length of river reaches varied from 20 to 65 km. Several upstream hydrographs at different stations on this river are generated using the developed method and compared with the observed hydrographs. The method predicts the time to peak with less than 5% error and peak rates with less than 10% error in the short river reaches of 20 km and 31 km. It also predicts the time to peak and peak rate in other two brances of 45 km and 65 km with less than 15% error. The method satisfactorily generates upstream hydrographs, with an overall mean absolute error (MAE) of 42 m3/s.


Introduction
Flow rate at a river section is required for many purposes including river restoration, navigation, design of flood control structures, reservoir operations, assessment of climate change effects and mitigation of drought impacts. Since flow rate is a product of cross section and velocity, the velocity and topographic cross-section measurements are required. For that purpose, channel section is equipped with hydrometric sensors for flow stage and current meter for flow velocity measurements (Moramarco et al. 2005).
In many cases, especially in developing countries, due to economic constraints, a single station on a river is often equipped to record flow rate. On the other hand, hydrographs at the upstream sections may be needed for especially assessing flooded areas, and evaluating drought and climate change effects (Das 2009). The hydrograph at a gauging station may have quite different characteristics than the one at the upstream section due to number of reasons such as different resistance, lateral inflow, drainage area and channel storage (Doria and Tanda 2012).
The reverse flood routing can be defined as the upstream hydrograph prediction using the channel reach geometric and downstream hydrograph characteristics (Doria and Tanda 2012). There are several reverse flood routing methods such as: (1) the numerical solution of the St.Venant equations (Eli et al. 1974;Szymkiewicz 1996;Dooge and Bruen 2005;Bruen and Dooge 2007;Artichowicz and Szymkiewicz 2009), (2) the employment of the Muskingum method (Das 2009), (3) the Bayesian geostatistical approach (Doria and Tanda 2012) and (4) the basic conservation of mass principle (Zucco et al 2015). When Eli et al. (1974) performed the reverse flood routing at James River in Virginia by employing the St. Venant equations, they observed the numerical instabilities especially at low discharges. Other studies also observed strong oscillations for sharp hydrographs and step waves when performing the reverse flood routing using the St. Venant equations (Szymkiewicz 1996;Dooge and Bruen 2005;Bruen and Dooge 2007). Bruen and Dooge (2007) stated that the parameters related to the discreatization may have influence on the numerical solution causing instabilities. Das (2009) indicated the separate calibration of the Muskingum method before the application for the reverse flood routing. In the case of the reverse flood routing using the Bayesian geostatistical approach, the upstream hydrograph is treated as a random function defined by statistical properties and some sort of degree of continuity and smoothness is imposed on the hydrograph shape (Doria and Tanda 2012).
Application of the reverse flood routing methods of the St. Venant equations and the Muskingum method requires substantial data on river reach geometry and flow variables. In addition, these methods were, so far, applied to regular rivers having subcritical flows (Doria and Tanda 2012). Although the Bayesian geostatistical reverse flood routing method can be applied to more general cases, it requires similar input requirements as the other methods since it routes the upstream hydrograph many times in downstream direction by employing the one-dimensional St. Venant equations to produce the downstream hydrograph as close as possible. The basic continuity equation-based reverse flood routing (Zucco et al 2015) method also requires not only discharge hydrograph but also cross-sectional area measurements.
This study proposes a new method coupling the continuity equation with the kinematic wave approximation of the momentum equation. The method first averages the depthaveraged kinematic wave equation along the channel length and it then performs the temporal averaging of the resulting equation. A two-parameter gamma distribution is assumed to describe the upstream hydrograph. The spatially averaged flow depth and lateral flow are related to the downstream flow stage. The optimal values of parameters are obtained by the genetic algorithm by producing the downstream flow stage hydrograph as close as possible. The flood routing equation developed in this study is just a function of flow stage whose measurement is fairly inexpensive and easy (Moramarco et al 2005;Perumal et al. 2007;Spada et al 2017;). The model requires data of downstream flow hydrograph, mean channel width, mean channel length, mean channel slope and mean roughness for a channel reach. In other words, only downstream stage measurement is needed and other data can be obtained from available satellite observations (Bjerklie et al. 2003(Bjerklie et al. , 2005.

Reverse flood routing model
The reverse flood routing model is based on the continuity equation: and the kinematic wave approximation where A is cross-sectional area, Q is flow rate, q l is unit lateral flow, u is flow velocity, h is flow depth, = C z √ S o , and β = 1.5 if Chezy's equation is used, C z is the Chezy's coefficient, and S o = S f , S o is channel bed slope and S f is friction slope (energy gradient).
Assuming that river cross section is almost rectangular and channel width is constant, equation can be written as Substitution of Eq. (2) into Eq. (3) results in the depthaveraged kinematic wave equation: The kinematic reverse flood routing equation starts with spatially averaging Eq. (4) along a channel reach length (L) at a fixed time (t) as follows: where L u and L d are the upstream and downstream station locations, respectively, i.e., L = L d -L u . Applying the Liebniz integral rule to the first term on the left-hand side of Eq. (5) and recalling that α and β are constant parameters, where h d is flow depth at the downstream section and h u is the flow depth at the upstream station.
Assuming that the channel length (L) does not change spatially at the station locations, i.e., dL hdx is spatially averaged mean flow depth and q l dx is spatially averaged mean unit lateral flow.
The spatially averaged flow depth,h(t) , and the spatially averaged lateral flow, q(t) , are expressed as power functions of flow stage at the downstream station as: where η and δ are coefficients whose optimal values are to be estimated. Similarly, When Eqs. (11)(12)(13)(14) are substituted into Eq. (10), the result is: Equation (7) is now averaged over a time period of Δt, i.e., from t to t + Δt, as follows: This averaging operation would result in the following equation: where the symbol ⟨.⟩ stands for temporally averaged term, i.e., term that is averaged over time Δt. The temporallyaveraged terms in Eq. (9) are approximated by assuming that each term averaged-over time Δt is approximately equal to the arithmetic mean of the corresponding terms at time t and t i m e ( t + Δ t ) , e . g . , ⟨ . This approximation was also employed in storage routing (PULS) and Muskingum flood routing methods (Henderson 1989). Equation (9) can thus be written as follows: Solving Eq. (15) for h d (t + Δt), the result is the following reverse flood routing equation: The two-parameter gamma distribution (Pearson Type III distribution) ) is proposed for describing the upstream flow stage hydrograph: where h b is the flow depth at base flow, h p is the peak depth value of the distribution, t p u is the time to peak when the peak value of the distribution occurs, and γ is the distribution shape parameter, which assumes values in [1-2] ). In Eq. (17), there are three parameters to be estimated, namely h p , t p u and γ. Note that the peak value of the upstream flow depth is the summation of h b and h p , i.e., In a similar fashion, the upstream hydrograph at (t + Δt) can be expressed as: Zucco et al (2015) also employed the same distribution (Pearson Type III) for expressing the upstream flow discharge hydrograph. Equation (16) is the basic reverse flood routing equation, while Eqs. (17-18) are the auxiliary ones. The reverse flood routing equation uses flow stage at the downstream station and mean values of channel length, channel width, channel bed slope and channel roughness coefficient. It has basically five parameters (η, δ, h p , t p u and γ) whose optimal values were estimated by the genetic algorithm. Note that in the basic equation (Eq. 16), = C z √ S o and β = 1.5 and once the flow depth is obtained, then it is easy to generate the upstream and downstream flow disch a rge hy d r o g r a p h s a s Q u

Genetic algorithm (GA)
GA is an evolutionary search algorithm which has three basic units (bit, gene and chromosome) and five operations (initial gene pool generation, chromosome fitness evaluation, chromosome selection, cross-over and mutation). The initial gene pool can be generated using some random functions. The chromosome fitness is evaluated in two steps by first substituting each chromosme into the objective function to find the respective value and then by dividing its value to the total value. The selection procedure is randomly performed by using either the roulette wheel method or by employing the ranking procedure. After the selection stage, the parent chromosome is formed before subjecting them to the crossover operation to generate new off-springs. After the crossover operation, the mutation can also be used to obtain diverse off-springs by simply reversing about 5% of the bits. Details on GA basics and operations can be obtained from Goldberg (1989), Sen (2004), Tayfur (2012Tayfur ( , 2015. For this study, the GA was employed to obtain the optimal values of the reverse flood routing model parameters by minimizing the mean absolute error (MAE) as an objective function: where N is the number of observations, h d model is the modelproduced downstream flow stage, and h d measured is the measured downstream flow stage.

GA model implementation
GA is applied to obtain the optimal values of model parameters by minimizing the mean absolute error (MAE), Eq. (19).
To start the iterations, random values are assigned within the search space for each parameter. The search space for the following parameters is decided based upon the information on the downstream stage hydrograph for each flood event as: where t rd is the time when the stage hydrograph starts rising at the downstream station and t pd is the time when the peak rate occurs at the downstream station. The optimal values of other parameters are searched within [0-1] for δ, η and [1-2] for γ in agreement with .
GA has employed 80 chromosomes in the initial gene pool, 80% cross-over rate and 4% mutation rate and 2000 iterations. The evolver GA solver for Microsoft Excel (Palisade Corporation, 2012) package program, which employs the Recipe Solving Method to minimize the objective function under specified constraints, is employed. It takes a very short CPU time, in the order of seconds, to run the program for thousands of iterations with 80 chromosomes in the gene pool.

Reverse flood routing procedure
The reverse flood routing procedure can be simply summarized as follows: 1. Initial values are randomly assigned for h p , t p u γ, η, δ within the search space for each parameter. 2. Upstream flow stage hydrographs h u (t) and h u (t + Δt) are computed using Eqs. 17 and18. 3. Downstream flow stage hydrograph h d (t + Δt) is computed using Eq. 16. 4. Mean absolute error for model produced downstream flow stage hydrograph and measured downstream flow stage hydrograph is computed using Eq. 19. 5. Current values of parameters ( h p , t p u γ, η, δ) are then updated. 6. Minimization of the error function (Eq. 19) is continued, while trying to reach the optimal values of parameters by performing steps from 2 to 5. 7. Iterations are stopped when the global minimum error is reached and the optimal values of parameters are obtained.
Note that a single iteration in GA involves steps from 2 to 5. At the end of a certain number of iterations (or satisfying predetermined tolerance limit), the downstream flow stage hydrograph is generated as close to the observed downstream flow stage hydrograph as possible. The parameter values that result in the best hydrograph are considered as optimal values which are employed in Eq. (17) to produce the upstream flow stage hydrograph. Once this is achieved, the upstream flow discharge hydrograph is produced by employing Q u (t) = B h u (t) .

Calibration and validation of the model
The calibration of the model is performed by simulating the measured downstream stage hydrographs as close as possible by minimizing the error function (Eq. 19) and consequently finding the optimal values of the model parameters ( h p , t p u γ, η, δ) by the GA. Note that discharge is computed by Q = h u B where h is the flow depth, u is the flow velocity, and B is the channel width. The flow velocity is computed as u = h −1 where = C z √ S o , and β = 1.5, C z is the Chezy's coefficient, and S o is the channel bed slope. Flow depth at the downstream section is computed by Eq. 16. Since, at a river reach, the channel width, the channel slope and the roughness coefficient are known, one can easily compute the flow discharge by Q d (t) = B h d (t) . The observed downstream section flow stage hydrographs are simulated as close as possible by minimizing the error function (Eq. 19) between the observed and model produced stage hydrographs using the steps outlined in the "Reverse flood routing procedure" section above. Therefore, the simulations of downstream stage hydrographs (and consequently the downstream discharge hydrographs), given in the section "Reverse flood routing application" are in fact the simulations as a result of the calibration procedure.
During the calibration procedure, at each iteration, the GA-based model, at the same time, computes the upstream stage hydrograph by Eq. 17. That is, while simulating the observed downstream stage hydrograph as close as possible by minimizing the error function (the calibration procedure), the model at the same time produces the upstream stage hydrograph by Eq. (17), as presented in the "Reverse flood routing procedure" section above. The upstream section flow stage and hydrograph simulations given in the section "Reverse flood routing application" are in fact the simulations as a result of the validation procedure. Note that the model only minimizes the error between observed downstream flow stage and model-produced downstream flow stage hydrographs (Calibration Stage) while generating the upstream stage hydrographs (consequently the upstream discharge hydrographs) (Validation Stage). The comparison between the observed upstream stage hydrographs (and discharge hydrographs) and the model-produced upstream stage hydrographs is in fact the simulations as a result of the model validation procedure.

River reaches
The model was applied to reverse flood routing in four different reaches on Tiber River in central Italy. Figure 1 shows the location of the hydrometric stations while Table 1 summarizes the main characteristics of the selected river reaches. Each gauged section is equipped with a remote ultrasonic water-level gauge, and velocity measurements are carried out by current meter. Several accurate flow measurements were available which allowed the estimation of the

Upstream hydrograph predictions
The measured downstream flow stage hydrograph was generated as closly as possible (i.e., minimizing the global error between the model-produced downstream stage data and measured downstream stage data), while finding the optimal values of the reverse routing model parameters. The figures given next involve simulations of measured upstream and downstream stage and discharge hydrographs. The measured upstream and downstream hydrographs were obtained by using the measured values of flow depth (h), flow velocity (u) and channel width (B) at the gauging stations (i.e., Q = h u B) while the simulated discharge hydrographs were com- and β = 1.5 where, along with computed flow stage, the average values given in Table 1 were used for channel width, slope and roughness. For each river reach, four simulations were performed and related performance measures are given in Tables 2, 3, 4 and 5, while for the sake of brevity only two simulations for each river reach are shown in Figs. 2, 3, 4, 5, 6, 7, 8 and 9.      Table 2 summarizes the performance measures for these plus December 1990 and June 1997 events. For all these events, the downstream stage hydrographs are simulated satisfactorily (Figs. 2 and 3). This is because these are used for the calibration of model parameters. The optimal values of parameters that have resulted in these simulations are summarized in Table 2. The downstream flow discharge hydrographs are also captured satisfactorily. Time to peak, peak rate, rising and recession limbs are all simulated (Figs. 2 and 3). The primary objective of this study is to produce the upstream hydrographs, especially the upstream flow discharge hydrograph. As seen in Figs. 2 and 3, these hydrographs are produced satisfactorily although the flood volumes are underestimated, especially in the case of May 1995 flood event. The rising and recession periods are well captured. The model has produced hydrographs reaching peak values earlier for all the events (Figs. 2 and 3, Table 2). The peak rates are slightly underestimated (Fig. 2, Table 2) or overestimated (Fig. 3, Table 2). Table 2 summarizes the performance measures for these simulations. As it can be seen, the time to reach peak rates is captured within, on average, less than 10.0% error.
Only for the December 1990 event, the time to reach peak rate takes longer with 42% error. Similarly, for this event, MAE = 1.35 m for the upstream flow depth. The average error for the other three events is MAE = 0.38 m. The model  (Table 2).

S. Lucia-P. Nuovo river reach
This reach is 65 km long with S o = 0.0014, B = 45 m, C z = 26 m 0.5 /s and T l = 5.5 h. Four events are simulated and related parameter values and error measures are summarized in Table 3. The model has performed poorly for this reach, may be due to longer wave travel time of 5.5 h and longer reach length of 65 km. As it can be seen in Fig. 4, the peak rate is underestimated at the upstream station, while the peak of June 1997 event (Fig. 5) is overestimated. As summarized in Table 3, there is, on average, 29% error in estimating the upstream peak rate and 15% in estimating the timing of the peak value. The overall, MAE = 41.0 m 3 /s and h u = 0.77 m.

P. Felcino-P. Nuovo river reach
The reach length is 20 km with S o = 0.0012, B = 47 m, C z = 28 m 0.5 /s and T l = 1.5 h. Table 4 summarizes the optimal values of model parameters and error measures for four events simulations. Figures 6 and 7 show simulations of hydrographs of two events. The model has produced the hydrographs satisfactorily, capturing the rising and recession limbs and as well as peak rates and time to reach peak rates although the flood volumes are underestimated, especially in the case of December 1995 flood Only in the case of event April 1993, it has overestimated the upstream peak discharge by 23.7%. On average, for the other three events, it has less than 10% error in predicting the peak rates of upstream hydrograph. It has generally captured the peak rates with no delay or early rise. Overall, it has made less than 5% error when capturing the time to reach the peak values (Figs. 6, 7,

Nuovo-M. Molino river reach
The reach is 31 km long with, on average, bed slope of S o = 0.0009, channel width of B = 50 m, and roughness coefficient of C z = 31 m 0.5 /s, and the wave travel time is 2.5 h. For this reach, four event hydrographs are simulated and the optimal values of model parameters and the related error measures are summarized in Table 5. Figures 8 and 9 show the simulated hydrographs for events that occurred in April 1998 and February 2004. The upstream stage and discharge hydrographs of the April 1998 event are satisfactorily reproduced with slightly early rise and underestimation (Fig. 8, Table 5) of peak rate and flood volume. The hydrographs of February 2004 event are reproduced satisfactorily in terms of not only the rising and recession limbs but also the time to peak and peak rates (Fig. 9, Table 5). The model generally has reached the peak rate early and underestimated the peak rate with low error for all the four events observed in this particular river reach. The model, on average, has 5% error in timing of the peak rate and 10% error in predicting the peak rate at the upstream station. These errors are lower for the downstream station. It has produced, on average, MAE = 54 m 3 /s and MAE = 0.96 m for upstream hydrograph rate and flow depth, respectively.

Discussion of results
The reverse flood routing model has satisfactorily generated the upstream stage and discharge hydrographs. In the P. Felcino-P. Nuovo River reach of 20 km and in the P. Nuovo-M. Molino reach of 31 km, the model has captured the timing of the peak with around 5% error and had, on average, about 10% error in capturing the peak rate. As the reach length becomes longer, the rate of error also increases. Generally, as the reach length becomes longer, the model produces hydrographs tended to reach the peak rate earlier. While the model makes around 10% error in capturing the peak in 20, 31 and 45 km river reaches, this error has increased to 30% in the 65 km river reach. The upstream peak rates of hydrographs observed and simulated in all these river reaches with different lengths ranging from 20 to 65 km are predicted, on average, with less than MAE = 42 m 3 /s. Similarly, the flow depths at the upstream stations of the four gauging sections are also predicted with, on average, MAE = 0.75 m.
The errors are low for predicting the stage at the downstream station since calibration is carried out by minimizing the error between observed stage and model produced stage hydrographs at the downstream station. Accordingly, the model has satisfactory simulated flow hydrographs at the downstream station. It has caused, on average, less than 6% error in capturing the peak and 8% error in capturing the timing of the peak at the downstream stations.
Considering that no information is used for the upstream end, the reverse flow routing model can be considered satisfactory for short (20 km) and as well as long river reach (up to 45 km) with an intermediate drainage area. This is important from a hydrological point of view, considering that mainly the stage observed at the downstream end can suffice for obtaining satisfactory results in terms of reverse flow routing at the upstream end.
The model has performed very poorly when applied to predict the upstream hydrograph of a longer reach of 96 km, which is the reach between Santa Lucia and Monte Molino sections (see Fig. 1).
At river reach of 45 km of S. Lucia and P. Felcino and at 20 km reach of P. Felcino and M. Molino, the model tends to underestimate the flood volumes. It also tends to reach the peak rates earlier. At 65 km river reach of S. Lucia and P. Nuovo, the model's poor performance is observed especially in capturing the peak rates. At 31 km reach of P. Nuovo and M. Molino, the model generally captures and estimates the peak rates with less than 5% error. The model cannot simulate double (or multiple) peak hydrographs. The model cannot also satisfactorily simulate hydrographs having longer equilibrium periods.

Concluding remarks
The reverse flood routing may be necessary for especially estimating flooding areas upstream region of a gauging station. For that reason, many reverse flood routing procedures are developed, requiring sohpisticated methodologies and substantial data requirements. This study developed a simple approach requiring easily measurable stage data.
The reverse flood routing model is based on the continuity equation and kinematic wave appoximation. It spatially averages the depth averaged kinematic wave equation along a channel length (L) at a fixed time, t, and then carries out the temporal averaging of the resulting terms over a time period of t + Δt to obtain the flood routing equation. The model simplifies the natural channel by assuming that the channel length and cross-sectional area are spatially constant and by approximating the temporally averaged terms as the arithmetical mean of the corresponding terms evaluated at time t and t + Δt. The model assumes a two-parameter gama distirubtion to define the upstream stage hydrograph and relates the spatially averaged lateral flow and mean flow depth to the downstream flow stage. The resulting basic routing equation is thus obtained as a function of only the downstream flow stage. Accordingly, the method mainly requires The model is calibrated using the measured dowstream hydrographs while it is validated using the measured upstream hydrographs. The model is applied to generate upstream hydrographs at four different river reaches of different lengths ranging from 20 to 65 km on Tiber River, Italy. The model can generate upstream stage and discharge hydrographs at river sections whose reach length is not longer than 45 km. It can overall capture hydrograph shape, time to peak and peak rate with less than 10% error. As the reach length becomes shorter (such as 20 km), the capturing of time to peak becomes even better. As the reach length becomes longer (such as 65 km), the generated upstream hydrographs tend to rise to peak rate earlier and the model mostly underestimates the peak rate with about 30% error. However, the model tends to underestimate the flood volumes. It cannot generate double peak (or multi-peak) hydrographs. It cannot also capture the prolonged equilibrium sections of the measured hydrographs. An application to a river reach longer than 65 km, such as 96 km between S. Lucia and M. Molino, could not reproduce satisfactorily the upstream hydrographs. Neverthless, requiring only the stage data measurements that are easy and economical to measure has important implications for sites having poor gauging stations where some water resources projects may be developed. The other required information on mean values of channel slope, roughness, width and length can be We can restate that this model can be used for the reverse flood routing for a river reach of up to 45 km. It may also be cautiously used for a river reach up to 65 km, but not longer than 65 km.