Inverse Flood Routing Using Simplified Flow Equations

The paper considers the problem of inverse flood routing in reservoir operation strategy. The aim of the work is to investigate the possibility of determining the hydrograph at the upstream end based on the hydrograph required at the downstream end using simplified open channel flow models. To accomplish this, the linear kinematic wave equation, the diffusive wave equation and the linear Muskingum equation are considered. To achieve the hydrograph at the upstream end, an inverse solution of the afore mentioned equations with backward integration in the x direction is carried out. The numerical solution of the kinematic wave equation and the Muskingum equation bases on the finite difference scheme. It is shown that both these equations are able to provide satisfying results because of their exceptional properties related to numerical diffusion. In the paper, an alternative approach to solve the inverse routing using the diffusive wave model is also presented. To this end, it is described by a convolution which involves the instantaneous unit hydrograph (IUH) corresponding to the linear diffusive wave equation. Consequently, instead of a solution of partial or ordinary differential equations, the integral equation with Laguerre polynomials, used for the expansion of the upstream hydrograph, is solved. It was shown that the convolution approach is more reliable comparing to the inverse solution of the simplified models in the form of differential equations.


Introduction
To ensure the requested flow conditions in a selected section of the considered river reach below a dam, the appropriate reservoir operating strategy should be applied. Thus, it faces the problem of flood control, typical in water management. This problem is closely related to inverse flood routing. Usually, inverse routing is solved using the optimization technique (Zucco et al. 2015;Azizipour et al. 2021;Badfar et al. 2021) or using the Bayesian geostatistical optimization approach (D'Oria et al. 2014). Sometimes, it is solved via the formulation of an inverse problem for the system of Saint-Venant equations. Usually, this system is integrated in a conventional way, which means that the equations are solved over a channel reach of length L for an increasing time until t = T, i.e., integration is carried out in the following domain: 0 ≤ x ≤ L and 0 ≤ t ≤ T. To ensure the uniqueness of the solution, the appropriate initial and boundary conditions must be imposed. The flood wave imposed as the boundary condition at the upstream end of the river reach allows the resulting hydrographs in the consecutive cross-sections to be calculated until the downstream end.
For the system of Saint-Venant equations, the inverse problem can be formulated as well. The equations can be integrated over the assumed time simulation [0, T], towards the diminishing x only if the flow is subcritical (Szymkiewicz 1993(Szymkiewicz , 1996. Such an approach is applied when the hydrograph at the upstream end ensuring the required flow condition at the downstream end is searched.
To the best knowledge of the authors, the first attempt of reverse flow routing has been undertaken by Eli et al. (1974). The inverse problem for the Saint-Venant equations may be solved either with the method of characteristics or with the finite difference method using the box scheme. The first approach was proposed by Bodley and Wylie (1978), whereas the second was applied by Szymkiewicz (1993Szymkiewicz ( , 1996. However, the application of these approaches to real-life problems is not easy, especially when a considered river has compound cross-sections and the considered simulation time T is large. For the reasons presented above, many authors believe that inverse flood routing can be successfully carried out using the simplified open channel flow equations (see, for instance, Koussis et al. 2012;Shamaa 2019;Badfar et al. 2021). Such an approach seems to be less complex and, consequently, it can be easier for possible implementation in engineering practice. However, it should be borne in mind that simplified models have a limited applicability.
In the paper, the following linear equations are analyzed from the point of view of inverse flood routing: -the kinematic wave model in the form of the linear single transport equation and the linear Muskingum equation related to the kinematic wave equation, -the linear diffusive wave model expressed by a convolution.
The non-linear simplified flood wave equations are beyond the scope of this paper because of the following reasons.
Among possible forms of non-linear kinematic wave equation, only one form satisfies the mass conservation principle (Gąsiorowski 2013) and therefore, only this form could be chosen. Moreover, the non-linear kinematic wave equation seems to be less interesting for inverse flood routing, because while traveling along an open channel axis it becomes steeper and steeper, leading sometimes even to a kinematic shock. Such steep shape of wave is rather not required at the downstream end. As far as the Muskingum equation is considered, many authors tried to improve the linear Muskingum model via the introduction of a non-linear formula relating to storage, inflow and outflow (Gill 1978;Tung 1985;Singh and Scarlatos 1987;among others). Consequently, instead of the linear equation, its non-linear version is obtained. However, Gąsiorowski and Szymkiewicz (2020) showed that usually the non-linear Muskingum equation becomes dimensionally inconsistent. Of course, such a situation is inappropriate and cannot be accepted from the point of view of physics. In spite of this, the non-linear Muskingum equation is also used for inverse flood routing. For instance, the linear and non-linear forms of the Muskingum equation and their application for inverse flood routing were analyzed by Badfar et al. (2021). Koussis et al. (2012) has also presented a discussion on various aspects related to the reverse solution of this equation. However, the results presented by the mentioned authors do not provide a definitive solution to the inverse flood routing problem.
The main objective of the study is to check the possibility of using simplified equations of unsteady flow in open channels to transform flood waves in the reverse direction. The paper shows that the numerical diffusion generated by the method of solution plays a crucial role in ensuring successful inverse flow routing. Moreover, it is shown that an approach basing on a convolution describing the simplified models is a more reliable than direct solution of the differential equations.

The Linear Kinematic Wave Equation
Consider the kinematic wave model consisting of the original differential continuity equation (without the lateral inflow) and the dynamic equation reduced to the Manning formula (Miller 1984;Chow et al. 1988): where t time, x longitudinal coordinate, h water stage, Q flow discharge, A cross-sectional area of flow, B channel width at the water surface, s channel bottom slope, considered as an average for a river reach of length L, R hydraulic radius, n Manning roughness coefficient.
The Eqs. (1) and (2) can be reduced to a single transport equation. Following the wellknown method (Chow et al. 1988), the advection transport equation is obtained. For a wide and rectangular channel with an energy line slope equal to the slope of the channel bottom, it is expressed as follows: with: and where C(Q) kinematic wave celerity, m kinematic wave ratio (= 3/5 when for friction the Manning formula is used), p wetted perimeter.
Equation (3) can be linearized. To this order C(Q) is replaced by the following constant value: where Q average discharge assumed to be constant.
Therefore, the following linear equation is obtained: For the numerical solution of the kinematic wave Eq. (7), the finite difference box scheme is used. It works on the grid point presented in Fig. 1.
All derivatives are approximated at point P, which can move inside the mesh, using the following formulas (Cunge et al. 1980): where X, θ weighting parameters, both ranging from 0 to 1, Δt time step, The approximation of the linear kinematic wave Eq. (7) using the Formulas (8)  In the case of the conventional solution, the nodal values Q j i and Q j i+1 are known either from the initial condition (Q(x,0) = Q in (x) for 0 ≤ x ≤ L) or from the previous time step (Fig. 1), whereas Q j 1 (j = 1,2,3,…, M) are known from the boundary condition imposed at x = 0 (Q(0,t) = Q 0 (t) for 0 ≤ t ≤ T). Consequently, this system of difference equations is particularly easy to solve as it is linear and, additionally, all of its equations can be isolated and solved separately because in each equation only one unknown exists. It is equal to: where (10) (for i = 1, 2, 3, … , N − 1; j = 1, 2, 3, … , M) is the Courant number. However, the order of computations from x = 0 towards x = L should be respected.
Equation (10) can be formally used for inverse flood routing as well. In such a case, the following initial and boundary conditions should be imposed at the boundaries of the solution domain (Fig. 1): -initial conditions: Q(L,t) = Q in (t) for 0 ≤ t ≤ T; -boundary conditions: Q(x,T) = Q T (x) for 0 ≤ x ≤ L.
In this case, the nodal values Q j i+1 and Q j+1 i+1 are known as the required hydrograph at x = L or as a result of the calculations carried out for the previous cross-section. In addition, Q M i is known as the boundary condition given at the time t = T. Therefore, similarly to the direct solution, in Eq. (10) only one unknown Q j i exists, which can be easily isolated: Computation should be carried out for the subsequent nodes along the t axis following from the last time level t = T towards t = 0, i.e., for j = M-1, M-2, M-3,…, 1. Finally, the required hydrograph at the upstream end of the considered channel reach is obtained.
As Eq. (10) approximates the partial differential kinematic wave equation, then numerical errors are introduced into its solution. It is well known that the wave attenuation obtained while solving the kinematic wave equation being of the hyperbolic type, has a numerical root. It results from the artificial diffusion generated by the numerical method applied for the solution. Detailed information on this subject results from an analysis of the modified equation. The modified equation is obtained by performing the inverse operation to the approximation of the solved differential equation. For this purpose, in the equation approximating the kinematic wave equation with a box scheme (Eq. (10)), all the nodal values of the function are replaced with Taylor series expansions. The box scheme modifies the kinematic wave equation as follows: where D n coefficient of numerical diffusion, and E n coefficient of numerical dispersion are given by the following formulas: In Eq. (14) all terms of the Taylor series neglected during the approximation appear on its right side. Note that for θ = 0.50 and X = 0.50, numerical diffusion disappears, whereas if additionally Cr = 1, all terms on the right side disappear so that the exact solution is obtained. It is important that for a smooth function, the sum of all terms on the right side of Eq. (14) is dominated by the first term (Fletcher 1991). Consequently, if the coefficient of numerical diffusion is sufficiently large, other terms are not important, then in the solution an artificial damping of flood wave will be observed. However, if D n = 0, which takes place for θ = 0.50 and X = 0.50, the term representing numerical dispersion becomes dominating. In such a case, spurious oscillations should be expected in the solution. Detailed information on the modified equation analysis proposed by Warming and Hyett (1974) is provided by Fletcher (1991) and Morton and Mayers (2005).
The modified Eq. (14) makes possible comprehensive interpretation of the results provided by the numerical method of solution. An example dealing both direct and inverse solutions is presented below.
Example 1 Numerical tests concerning the kinematic wave equation are performed as follows: -for the given initial and boundary conditions, Eq. (11) is solved by directly providing the hydrograph at the downstream end, -then, the hydrograph computed at the downstream end was treated as the required one, and the hydrograph at the upstream end was determined by inverse flood routing (Eq. (13)).
In this way, it was possible to compare both the imposed and computed hydrographs at the upstream end. Consequently, the effectiveness of the solution method and its accuracy could be estimated.
In the example presented below, the following data was assumed: It is assumed that the hydrograph imposed at the upstream end is given by the following formula: where , β parameter determining the shape of wave,with the following values of parameters: Q in = 5.0 m 3 /s, Q max = 100 m 3 /s, t max = 4.0 h, t 0 = 0. h, and β = 2.0 (curve a in Fig. 2).
The other parameters were assumed as follows: the time step is equal to Δt = 600.0 s, ensuring the Courant number equal to Cr = 0.40. The weighting parameters in Eq. (11) for the direct solution were θ = 0.50 and X = 0.25. The hydrograph computed at the downstream end is presented in Fig. 2 as curve b. As can be seen, its peak, which appeared at the downstream end ca. 12 h later, was damped from 100 m 3 /s to 80 m 3 /s. Then the calculated hydrograph was assumed as expected and inverse flood routing was performed. In this case, the weighting parameters equal to θ = 0.50 and X = 0.50 are assumed. The computation results are also shown in Fig. 2 as curve c.
The hydrograph computed at the upstream end is significantly different than expected. Numerous tests carried out for different values of both parameters θ and X provided similarly unsatisfying results. However, the same set of parameters as used for the direct solution, i.e., θ = 0.50 and X = 0.25, provided an excellent result (curve d), which is consistent with the expected result (curve a). As can be seen in Fig. 2, indeed, both curves are virtually identical.
A very important fact results from numerical experiments. Any set of parameters X and θ ensuring numerical diffusion with the coefficient D n provides an inverse solution completely consistent with the hydrograph given at the upstream end. Therefore, if a hydrograph at the downstream end (x = L) is computed for assumed values of C, Δx and Δt, giving a numerical diffusion coefficient equal to D n , then another set of parameters ensuring the same coefficient of numerical diffusion will give an inverse solution consistent with the hydrograph at the upstream end. In other words, if parameters X and θ give a direct solution, then assuming θ 1 , an exact solution of the inverse problem is obtained, on condition that the parameter X is equal to: where D n is calculated accordingly to Eq. (15). Summarizing, each pair of values of parameters θ 1 and X 1 provide an exact solution of the inverse problem when a constant wave celerity C and constant mesh of nodes are used (Fig. 3). However, the parameters involved in Eq. (18) should give the same value of coefficient D n as that generated in the direct solution.
Returning to the modified Eq. (14), it enables the interpretation of the results of calculations presented in Fig. 2. The attenuation of the hydrograph computed at the downstream end (Fig. 2, curve b) results from numerical diffusion. In this case, for X = 0.25 and θ = 0.50, it is equal to D n = 1050 m 2 /s. If identical weighting parameters are assumed in the inverse problem, the correct solution is provided. On the other hand, if both parameters are assumed to be equal to 0.50, a solution corresponding to a pure translation of the downstream hydrograph (curve c in Fig. 2) is obtained. It is a completely different hydrograph than that expected at the upstream end. Note that according to Eq. (15), for X = 0.50 and θ = 0.50, numerical diffusion is not generated (D n = 0). Consequently, such a set of parameters is outside the range of permissible combinations of the parameter values X 1 and θ 1 leading to the correct value of the diffusion coefficient D n equal to 1050 m 2 /s (see Fig. 3).
The presented inverse solution of the linear kinematic wave equation can provide the results corresponding to the inverse routing of the diffusive wave equation. This idea is similar to that of Cunge (1969)  where D is coefficient of hydraulic diffusivity given as D = Q/(2b⋅s) (Chow et al. 1988), whereas Q is average flow rate considered as constant. Note that the inverse solution of Eq. (19) in a similar way to the inverse solution of the kinematic wave equation is not possible because Eq. (19) is of the parabolic type. However, if the set of the assumed parameters used for numerical solution of the kinematic wave Eq. (7) ensures the numerical diffusion equal to the hydraulic one involved in the linear diffusive wave, i.e. when D n = D, then both solution will coincide. Owing to such an approach inverse solution of the diffusion wave equation is obtained via inverse solution of the kinematic wave equation. However, in this situation, the following question seems justified. Is it possible to use the kinematic and diffusion wave equations interchangeably if they were derived with different assumptions?

The Linear Muskingum Equation
Another possible approach to inverse routing concerns the Muskingum equation. Its original form has been derived from the storage equation obtained by the integration of the continuity equation (Eq. (1)), which next was completed with an additional linear formula relating to storage, inflow and outflow (McCarthy 1938): where Q i inflow at the upstream end of a space interval, Q i+1 outflow at the downstream end of a space interval, K constant parameter expressed in time units, X dimensionless weighting parameter ranging from 0 to 1.
The parameter K is expressed as follows: and it is interpreted as the time needed for the wave to travel the distance Δx between two neighboring cross-sections i and i + 1.
As far as the Muskingum model is considered, the corresponding algebraic equation approximating Eq. (20) can be derived directly from Eq. (10) which approximates the kinematic wave equation. Using the relation (21), the following formula is obtained: where the number of reservoirs is equal to N-1.Usually the Muskingum Eq. (20) is solved using the implicit trapezoidal rule (Cunge 1969). Note that this approach coincides to Eq. (22) for θ = 1/2. However, the Muskingum equation can be integrated using another . It is given as follows: where: Keeping the order of calculation from the end j = 1 towards the end j = N-1, all nodal values of discharge are obtained at the next time level.
The Muskingum Eq. (22) can also be used for inverse routing. In this case, in Eq. (22) the only unknown is the nodal value of discharge Q j i . Its isolation yields: where: The linear Muskingum Eq. (20) is closely related to the kinematic wave Eq. (7). The similarity of both equations becomes obvious when the Muskingum equation, being the ordinary differential equation, is solved using the implicit trapezoidal method, while the kinematic wave equation, being the partial differential equation, is solved using the finite difference box scheme. A comparison of Eqs. (12) and (20) confirms their identity, because taking into account the relation (21), the Courant number (12) can be expressed as: This result is not surprising. It is well known that the Muskingum equation is a semidiscrete form of linear kinematic wave equation (Szymkiewicz 2002). For this reason, it can be expected that the linear kinematic wave Eq. (7) and the linear Muskingum Eq. (20) will give identical results for both conventional and inverse flood routing.
Examples of both direct and inverse solutions concerning the Muskingum equation are presented below in Example 2.
Example 2 At the upstream end, a hydrograph in a more realistic form was assumed (see Fig. 4-curve a). Its shape is similar to hydrographs that occur in practice during the operation of reservoirs. The calculations were performed for N = 30. The other parameters involved in Eq. (22) were assumed as follows: time step Δt = 0.5 h, weighting parameter θ = 0.5, whereas the constant K = 1 h gives the Courant number equal to C r = 0.50. The weighting parameter X for the direct solution was X = 0.25. The hydrograph computed at the downstream end is presented in Fig. 4 as curve b. As can be seen, the imposed wave with steep fronts is remarkably attenuated and smoothed. Next, this hydrograph was imposed at the downstream end and inverse flood routing was performed. The result of computation for the weighting parameter X = 0.35 is presented in Fig. 4 as curve c (dashed line). As can be seen, the obtained hydrograph differs significantly from the one that was given as the inflow hydrograph.
Taking into account the results from the kinematic wave example, in the next calculation, the weighting parameter equal to X = 0.25 is assumed. The computation results are represented by the curve d (dotted line). In this case the hydrograph computed at the upstream end is very similar to the expected one, i.e., imposed as the upstream hydrograph for the direct solution. This means that even for a rather challenging form of hydrograph, the results of the inverse solution can be considered as excellent.
Relating to the results provided by the Muskingum equation, they can be interpreted in the same way as those presented for the kinematic wave equation. Also in this case, a very helpful analysis of the modified equation can be performed. The modified equation corresponding to Eq. (22) is as follows: where the coefficient of numerical diffusion is equal to: As can be seen, in this case, the coefficient of numerical diffusion depends on both weighting parameters θ and X. However, if θ = 1/2 is assumed, as in the standard approach, then instead of Eq. (29), the following equation is obtained: The parameter θ makes possible the introduction of additional numerical diffusion into the solution if necessary. It is also worth adding that the coefficient of numerical diffusion given by Eq. (29) with relation (21) exactly corresponds to Eq. (15) describing the numerical diffusion coefficient in the kinematic wave equation modified by the box scheme.

Inverse Flood Routing Using the Diffusion Wave Equation and Convolution
It is well known that any linear dynamic system, such as an open channel reach, can be described using the convolution (Chow et al. 1988): where t time, p "memory" of the channel reach, Q(t) flow rate at the downstream end, q(t) flow rate at the upstream end, τ time (dummy parameter).
The function h(t) is an impulse response function, which in hydrology is called the Instantaneous Unit Hydrograph (IUH). Typically, for the known IUH and for the hydrograph q(t) given at the upstream end, the hydrograph at the downstream end Q(t) is computed. Of course, both solutions obtained using the linear partial differential equation and the IUH corresponding to this equation are equivalent.
Conversely, if the required hydrograph at the downstream end Q(t) is imposed and the IUH is known, then inverse flood routing will provide the upstream hydrograph q(t). As can be seen, it is another approach than the inverse solution of the partial differential equation. For instance, in the case of the linear diffusive wave (Eq. (19)), the downstream hydrograph can be obtained via the solution of the partial differential equation or using the convolution (31). This is possible because the IUH for an open channel reach of length L and corresponding to the diffusive wave equation, is well known (Hayami 1951): Regarding the Muskingum equation (Eq. (20)), it can be considered as a semi-discrete form of the kinematic wave equation. On the other hand, owing to the numerical diffusion (29) D n = (2 − 1) 4D ⋅ t generated in its numerical solution, it is able to reproduce the solution of the diffusive wave equation. This idea has been used in the Muskingum-Cunge model (Cunge 1969). Consequently, the linear Muskingum model can be expressed also in the form of a convolution which involves the IUH of the diffusive wave equation. Taking into account this fact, Szymkiewicz (2002) proposed an alternative form of the IUH for the Muskingum equation.
If the relations typical for the Muskingum equation are introduced into Eq. (32), it takes the following form: All symbols involved in Eqs. (32) and (33) were used and explained previously. A very important conclusion resulting from the discussion presented above is as follows: inverse flood routing using the diffusive wave equation and the Muskingum equation, can be performed using the convolution. Consequently, another approach than presented in previous sections for solving this very practical problem is possible.
Equation (33) possesses all the required properties formulated with regard to the IUH: -it holds for X ≤ 1/2 only, including the negative values as well; -for X → 1/2, it tends to the Dirac δ − function: -the parameter N can be any positive number, not necessarily an integer (Chow et al. 1988).
The presented IUH is able to reproduce required lag time of upstream and downstream hydrographs. The IUH for selected sets of parameters is displayed in Fig. 5.
An analyses of the function h(t) (Eq. (33)) reveals that it tends asymptotically to zero for t → ∞. Consequently, it has values significantly different from zero in the range 〈0, p since h(t) ≤ ε for t ≥ p (ε is positive, a small number). Thus, in order to compute the value of Q at time t, we need to know the values of the input function q in the preceding interval ⟨t−p, t⟩ . More information and additional details are given by Szymkiewicz (2002). As a matter of fact, the IUH given by Eq. (33) is valid for all previously considered simplified equations, i.e. for the kinematic wave equation, the Muskingum equation and the diffusive wave equation.
The integral of the convolution has the property of symmetry, which allows Eq. (31) to be written in an equivalent form: Conventional flood routing when the functions h(t) and q(t) are given is not a problem. To this end, the integral in Eq. (34) is calculated numerically using the rectangular formula. A more challenging problem occurs when for a known h(t) and for the required downstream hydrograph Q(t), the hydrograph q(t) at the upstream end should be computed. To this end, any decaying function such as q(τ) can be expanded on the basis of polynomials satisfying the condition of orthogonality. In the paper, weighted Laguerre polynomials are used (de Sousa and Matt 2019).
Assume that the unknown hydrograph q(τ) is represented by a series expansion as: Denoting the integral as Equation (36) can be rewritten as: Note that the relationship (38) is a convolution in which Laguerre polynomials appear as integral kernels.
As far as weighted Laguerre polynomials are concerned, they are defined by the following formula (de Sousa and Matt 2019): The first 4 polynomials are given by the following formulas: The weighted Laguerre polynomials are orthogonal functions for t ≥ 0. Usually, a time scale factor should be introduced into the Laguerre polynomials. It has to ensure the appropriate transformation of the time scale of the considered process (de Sousa and Matt 2019). Comparing the convergence of the Laguerre polynomials with the shape of the expanded functions, in this case, one can assume that the time scale factor is equal to 1. Consequently, the polynomials in the form presented above will be used.
Introducing the discrete time, Eq. (38) can be rewritten as: with Eq. (37), in which the integral is expressed using the approximating rectangular formula: On the other hand, the hydrograph required at the downstream end is known and it is denoted as Q(T) . Therefore, the coefficients c m involved in Eq. (44) should be chosen in such a way that the following objective function related to the least squares method will take the least possible value Because E(c 0 ,c 1 ,…,c M ) linearly depends on the coefficients c m they can consequently be determined via the solution of the linear system of algebraic equations resulting from the condition defining the extreme point of the function E(c 0 ,c 1 ,…,c M ): The solution of the system (47), having the dimensions (M + 1) × (M + 1), provides the coefficients c m , which next, according to Eq. (36), are used for the calculation of the hydrograph at the upstream end. An illustration of inverse flood routing using the proposed method is presented in next section.

Inverse Flood Routing Using the Experimental Data
Consider the unsteady flow in a physical model of cascade consisting of two reservoirs installed in laboratory (Geringer 1995). The results recorded during the experiments are shown in Fig. 6. A wave entering the cascade (in Fig. 6 curve "a") after transformation by two reservoirs in series takes the form of downstream hydrograph (in Fig. 6-curve "b"). The flood routing is carried out using the Muskingum equation expressed by the differential equation and by the convolution. Note that the kinematic wave equation derived for the channel cannot be used in the case of flood wave traveling through reservoirs.
Regarding the conventional solution of the Muskingum model in the form of convolution a good fit of the calculated hydrograph (curve c) to the observed hydrograph (curve b) was obtained. Computations were carried out for N = 2, K = 40 s, X = − 0.20 and Δt = 20 s. Memory p of the system constituted by two reservoirs was assumed to be equal to T p = 900 s. The computations performed for different maximal degree of the Laguerre polynomials M showed that this factor seriously influences quality of the provided results. The curve c has (47) The negative value of the weighting parameter X requires some comment. The considered case of flood routing is rather unusual because the flood wave does not travel through the reach of open channel having prismatic cross sections, but through two reservoirs in series. For this reason the lag time of both hydrographs is relatively small while the wave damping is more intense than in open channel. Consequently, to obtain the required consistency of the calculated and observed hydrograph the untypical negative value of the weighting parameter X had to be used. This is because relatively great numerical diffusion ensuring wave attenuation was required. It is worth adding that formally the assumed IUH (Eq. (33)) allows negative values of the parameter X. For the Muskingum-Cunge method this issue was discussed by Szel and Gaspar (2000). One can expect also that a similar problem will occur when the inverse flood routing is performed in a channel having compound cross-sections.
Next, the observed hydrograph (curve b) was assumed to be required at the downstream end. For the value of above listed parameters an inverse flood routing was performed. It appeared that the Muskingum model in the form of convolution with IUH in the form of Eq. (33) provided satisfying results. This fact is confirmed in Fig. 6. The hydrograph computed at the upstream end is displayed in the form of dotted line (curve d). Indeed, curve d in Fig. 6 is very similar to the expected curve a. In this case, the root mean square error equal to RMSE = 0.00130 m 3 /s is slightly greater than one corresponding to conventional solution. It should be added that if instead of the observed hydrograph the calculated one (curve c) is imposed at the downstream end as required one, then perfect agreement with the hydrograph observed at the upstream end is ensured.
As for the Muskingum model in the form of differential equation, it appeared that it was not able to provide the reasonable results. Although the result of the conventional solution obtained for N = 3, K = 40 s, and dt = 20 s (C r = 0.5) was very good, the inverse flood routing was completely unsuccessful. Since in this case the inverse solution of the Muskingum differential equation was not possible, one can suppose that we are dealing with ill posed problem. Insignificantly changes to the input data correspond to drastic changes to the solution. The results of numerous numerical experiments suggest that the Muskingum model in the form of a differential equation becomes unreliable when the flood wave is intensively attenuated, i.e. when numerical diffusion generated in the solution is very high. It should be remembered that the Muskingum model, which is a semi-discrete form of a kinematic wave, can be used for channels with a significant longitudinal slope. This is not the case in the cascade under consideration.
Presented example suggests that the Muskingum model in the form of convolution is more reliable tool than one in the form of differential equation. Moreover, the convolution (34) with the IUH in the form of Eq. (33) can be considered as a general form of all simplified flood routing models. The choice of the model is determined by the assumed parameters involved in Eq. (33). Finally, it seems that rather the convolution approach should be recommended for inverse flood routing using the simplified open channel flow equations.

Conclusions
Simplified flow equations, in the form of the linear kinematic wave equation, the Muskingum equation and the diffusive wave model were analyzed from the point of view of their ability for solving the problem of inverse flood routing. An alternative approach basing on a convolution was also analyzed. It was shown that both types of mathematical models can be successfully used for inverse flood routing on condition that flood wave's attenuation is not very great.
In the first case, the partial differential equation of the kinematic wave equation is integrated backwards. It was shown that the solution accuracy is strongly determined by the numerical diffusion generated by the solution method. It was found that each set of parameters involved in the formula for the coefficient of numerical diffusion ensures perfect solution accuracy on condition that the parameters produce the same numerical diffusion as that generated in the conventional solution. Moreover, it appeared that the Muskingum equation, described in fact by equivalent equation, provides also equivalent results.
In the second approach, the diffusive wave model is described by a convolution. Using an alternative form of the IUH a general model representing all simplified models was obtained. The inverse flood routing using a convolution bases on Laguerre polynomials used for the expansion of the searched upstream hydrograph. As the assumed objective function linearly depends on the coefficients of expansion, they can be determined using the least squares method. Consequently, instead of the numerical solution of differential equations, an integral has to be computed. This quite different approach for inverse flood routing appeared more reliable than an approach basing on the numerical solution of the differential equations even in the case when they are not valid. Due to convolution approach the difficulties of the inverse solution of parabolic equations, which include the diffusive wave equation, can be avoided.
The preliminary results presented in this work seem encouraging. However, in order to better evaluate the practical aspects of the proposed approaches, more hydrological data acquired from field experiments are needed.
Author Contribution All authors contributed to the study's conception and design. Material preparation, data collection and analysis were performed by Romulad Szymkiewicz and Dariusz Gąsiorowski. The manuscript was written by Romulad Szymkiewicz and Dariusz Gąsiorowski. All authors read and approved the final version of manuscript.
Funding The authors declare that no funds, grants, or other support were received during the preparation of this manuscript.

Availability of Data and Material
The authors confirm that some data and codes are available from the corresponding author on request.

Ethics Approval
The authors declare that the manuscript is original and has not been published in any journal.

Consent to Publish
The authors declare their consent to publication of the manuscript in "Water Resources Management" journal.

Competing Interests The authors have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.