Network flow and flood routing model for water resources optimization

Real-time management of hydraulic systems composed of multi-reservoir involves conflicting objectives. Its representation requires complex variables to consider all the systems dynamics. Interfacing simulation model with optimization algorithm permits to integrate flow routing into reservoir operation decisions and consists in solving separately hydraulic and operational constraints, but it requires that the water resource management model is based on an evolutionary algorithm. Considering channel routing in optimization algorithm can be done using conceptual models such as the Muskingum model. However, the structure of algorithms based on a network flow approach, inhibits the integration of the Muskingum model in the approach formulation. In this work, a flood routing model, corresponding to a singular form of the Muskingum model, constructed as a network flow is proposed and integrated into the water management optimization. A genetic algorithm is involved for the calibration of the model. The proposed flood routing model was applied on the standard Wilson test and on a 40 km reach of the Arrats river (southwest of France). The results were compared with the results of the Muskingum model. Finally, operational results for a water resource management system including this model are illustrated on a rainfall event.

where Q is the flow rate, and h is the water surface elevation. Bars above and below a variable denote the upper and lower bounds for that variable, respectively. g is a function that describes the flow in the different components of the hydraulic system. The objective function z is to be minimized ans depends on the total flood damage. Different optimization methods were proposed in the litterature for the optimization of complex water resource systems. An optimization method suited to the case to be addressed, i. e. to solve Eq. (1), depends on the nature of the objective function and the constraints. Labadie 4 and Yeh 5 reviewed the state of the art in optimization of multireservoir systems management and operations .
The constraints of the model (Eq. 2) can be divided into two groups: the hydraulic constraints and the operational constraints. The consideration of hydraulic constraints is crucial for the quality of the operations, since release decisions are made upstream and the target areas are usually downstream. Hence, every water resource optimization model should take into consideration the attenuation of the wave during the transfer 6 .
Many models exist to represent the flood routing in a reach. Basic routing approaches may be classified into two main families: Hydraulic and conceptual approaches.
The hydraulic approach applies the governing Barre de Saint Venant equations represented by the continuity equation and the momentum equation, respectively 7 . The integration of the Barre de Saint Venant equations in the formulation of water optimization algorithms presents several difficulties since these partial differential equations are nonlinear, and their numerical resolution requires a large amount of calculation 8 . In order to integrate flow routing into the reservoir operation model, many researchers interfaced a simulation model with the optimization algorithm 2,9-12 discussed the combination of simulation and optimization for real time flood operation for reservoir system and suggested the coupling of nonlinear optimization and simulation to close the gap between theory and practice. Wurbs 3 discussed simulation, optimization and combined simulation-optimization modeling approaches and presented the strengths and weakness of the reviewed models.
The methodology of interfacing a simulation model to an optimization algorithm consists in solving separately the hydraulic and operational constraints: the optimization algorithm generates the optimal reservoir operating decisions, and the simulation model appropriately simulates the propagation of the flow for a given flood hydrograph and a set of operating decisions. Figure 1 presents the process of interfacing a simulation model and an optimization algorithm 13 . Simulation models produce outputs that are used by the optimization strategy to find an optimal solution. Interfacing a simulation model with an optimization algorithm is only possible when the algorithm is based on an evolutionary computation approach 14 . The typical structure of evolutionary computation makes them suitable to adapt vaguely defined objective functions and constraints 15 . However, optimization-simulation process is time consuming and convergence problems may occur 16 . It is also to be noted that simulation models require: geometric data, initial condition, boundary condition, and hydraulic parameters which are not always available and bathymetric survey campaigns are very expensive. In addition, modeling biases can be observed that question the parameters of the model, and more so in the context of climate change 17,18 . In the literature, the consideration of channel routing in optimization problem for the operation of multireservoir systems was also performed using conceptual routing procedures. Conceptual models are characterized by the fact that one does not seek to understand in detail the physical phenomena that occurs within the flow, but consider the network in its entirety; in other words, as a simple input-output transformer. The calibration of the model using input and output values allows fixing the parameters of the model. These models reflect only the consequences of the phenomena occurring in the system and therefore get over the difficulties of the hydraulic complexity. Most conceptual models are reservoir models; that is, the functioning of each reach is assimilated to the operation of one or more reservoirs in series or in parallel. Conceptual models are based on the continuity equation (the variation of the reach storage (S) corresponds to the difference between the inflow (I) and the outflow (O), see Eq. 3); and a second empirical relation (storage function, see Eq. 4) that connects the reach storage and the outflow rate 19 .  www.nature.com/scientificreports/ The storage function can take many forms. Storage can be expressed as a function of inflow, outflow or both: The Muskingum flood routing model is the most used model 19 . The absolute storage is a function of both outflow and inflow discharges. The evolution in time of the absolute storage is expressed by: where S t is the absolute reach storage at time t; I t and O t are the rates of inflow and outflow at time t, respectively; K is a coefficient with time dimension that represents the storage time of the reach or the travel time of flood waves through the channel reach; and x is a weighting dimensionless coefficient (between 0 and 0.5) that modulates the influence of the inflow and the outflow. Expressing Eqs. (3) and (5) in finite difference form for a time interval, while considering a transit time (TT) as additional pure delay, leads to: where C 0 , C 1 and C 2 are constants that are computed from K and x and T 20 . Muskingum flood routing model is widely used in optimization problem for the operation of multireservoir systems to represent channel routing, mainly because its formulation is linear and does not increase the complexity of the problem. Muskingum model has always been the first choice to model channel routing when the optimization problem is in a linear form. Windsor 21 formulated a theoretical recursive linear programming model for the operation of flood control, using the Muskingum method for channel routing. Hsu and Wei 22 developed a reservoir real-time operation model for determining the optimal real-time release during a typhoon. The formulation of the problem is based on linear programming and the streamflow routing along a reach is modeled by using the Muskingum method of linear channel routing. Kumar et al. 23 adopted Folded Dynamic Programming for developing optimal reservoir operation policies for flood control with channel routing based on Muskingum model imbedded within the algorithm.
Linear programming (LP) objective functions and constraints are restricted to summations of linear terms. It is the optimization technique that is most often applied in modeling reservoir/river systems as well as flood management 24 . The main advantages of linear programming are its ability to optimize large problems, its convergence towards the global optimal and the availability of efficient software packages under free license 25 ). The LP is expressed as: where z is the objective function, x j are the decision variables, c j , a ij , and b i are constants, n is the number of decision variables, and m is the number of constraints.
As a particular form of linear programming, network flow programming 26 , because of its intuitive formulation and short resolution time, is often used in water management applications, and is suitable for solving large-scale allocation problems of multi-reservoirs and multi-periods 3,27 . In order to applied graph theory algorithms, the hydraulic system is modeled as a directed graph where convergence and diversion points, demand locations and water sources are represented by nodes, and reservoir releases, channel flows, carryover storage and withdrawals are represented by arcs. Network flow optimization problem is expressed as: where q ij is the flow rate in the arc connecting node i to node j; c ij is a cost for q ij ; q ij and q ij are lower and upper bound on q ij , respectively. The only constraints allowed are the ones in the form of a "mass balance" equation.
The special structure of network flow programming inhibits the integration of the Muskingum flood routing model in the network flow problem formulation, since Eq. (6) violates the form of Eq. (10) for nodes. Considerable efforts are made to include proper modeling of hydrologic channel routing into network flow models 28 . Braga and Barbosa 29 report on inclusion of channel routing into multiple time step optimization using an advanced

Flood routing model
Mathematical formulation of the flood routing model. In the following, we will refer to the flood routing model developed herein as the Residual Storage Model (RSM). Although inflows during immediate moments have a marginal impact on the outflow, all reservoir routing models consider that the outflow at time t is a function of the absolute reach storage (total inflow minus total outflow) at the same time. Herein, the portion of the absolute reach storage that impacts the outflow significantly is defined as the residual storage of a reach. Let S ′ t denotes the residual storage at time t. Figure 2 illustrates the upstream and downstream of a reach, and the residual storage.
The conservation equation can be written as: where TT is the Transit Time. This parameter represents the time from when the inflows {I t−TT , I t−TT−1 , . . . , I 0 } impact the outflow at time t. Hence, inflows {I t , I t−1 , . . . , I t−TT+1 } , are not used in the computation of O t and are represented by the white area in Fig. 2. The RSM considers that outflow at time t is proportional to S ′ t and I t−TT : The parameters of the model are the reach's Transit Time (TT), the proportionality coefficient ( α ∈ [0 : 1] ) and the initial residual storage S ′ 0 . The proportionality coefficient α physically represents the proportion of the residual storage that stays in the reach. It is between 0 and 1. The initial residual storage S ′ 0 , represents the reach storage during a steady flow. The RSM form is a singular form of Muskingum model (see Eq. (5) with x = 0 ) and is suited for long reach where the downstream flow has a limited influence on the upstream flow.
Network flow model. In a water management problem based on a network representation, the RSM can be easily integrated if it is described as a network. Let G = (V ; E) be a directed single source network, with node set V and arc set E. Let S and T be the source and the sink fictive nodes of the network, respectively. The source node supplies the upstream of the system and the sink node collects the downstream flows. Convergence or diversion points and reservoirs are represented by nodes and water transfer by arcs. In fact, flows do not cross the network instantly, thus, in order to account for the dynamics of the flows, the nodes are duplicated at each time step over the duration of the simulation. In the arcs connecting those copies, the transit times and flows are implicit 32 . For an arc e ij , i is the origin node, and j is the end node. Let γ (n) and γ −1 (n) respectively denote the sets of the outgoing and incoming arcs of a node n. Each arc e is associated with a positive flow e . For each node n of V, except for S and T, the conservation of flow is satisfied: Considering the reach's transfer time TT, the outgoing flow from an upstream node at time t is connected to the reservoir node at time t + TT . In order to model a reach's flood routing, an intermediate node denoted as reservoir node is introduced between every upstream and downstream node. The reservoir node separates the incoming flow into 2 flows: a flow corresponding the volume released from the reservoir at time t and a flow corresponding to the residual storage remaining in the reservoir. Figure 3 presents an example of a network model corresponding to channel routing over one reach. To simplify the example, we considered that the transit time is equal to one-time step, and we only provided the corresponding network for a horizon of 4 times step. Nodes U t and D t represent the upstream and the downstream In order to account for the complex objective functions involved in the calibration step, a genetic algorithm is used instead of traditional optimization algorithms. The Genetic Algorithm (GA) solver in MS Excel is coupled to a network model. Unlike classical optimization search methods, such as the simplex method and gradient-based methods, the genetic algorithm does not necessarily require well-defined functions or derivatives of functions. A genetic algorithm is a metaheuristic based on three bio-inspired operators: selection, crossover, and mutation. In the optimization model, the population in GA is the vector ( TT, α, S ′ 0 ), and the constraints are defined through the network model. Validation of the model. In order to evaluate the validity of the proposed model, it was tested on a hydrograph proposed by Wilson 33 , which is a standard test event that has been extensively studied by other researchers 20,[34][35][36] . The inflows and outflows of the Wilson event, the model routed, and Muskingum outflows are given in Table 1.
The results of the RSM and Muskingum model are illustrated in Fig. 4.  For the two reconstructions of the outflow hydrograph, the root-mean-square (RMS), the error, and the model's parameters are listed in Tables 2 and 3. The determination of K and x is performed by employing a least-square technique. Figure 4 shows a good agreement between the RSM reconstructed outflow and the observed ones. For the Muskingum model, C 0 and C 1 are negligible compared to C 2 which shows that the outflow at time t + 1 depends   Table 2, the errors for the residual storage model and the Muskingum model is 4.73 and 4.48 respectively; and the RMS is 8.37% and 9.32% respectively. With the recorded outflow peak time of t = 60hr , the RSM model was able to yield the same peak time with minor error. Whereas the Muskingum incorrectly predicted the peak time. The RSM also yields a very small RMS from practicing hydrologist point of view. Figures 5 and 6 represent the absolute storage of the reach plotted against the outflow for the Muskingum method, and the residual storage of the reach plotted against the outflow for the RSM, respectively. Figure 5 Table 3. Parameter of the calibration results for the Wilson event.   Fig. 6, the loop is nearly closed and a highly linear relationship between the outflow and the residual storage can be observed. Figures 5 and 6 emphasize that the outflow is more linearly correlated to the residual storage than the absolute storage of a reach.

Application
Application on the Arrats river. The RSM developed herein has been applied on a reach of the Arrats river. The Arrats river is an affluent of the Garonne river in the south west of France. The river is equipped with five hydrometric stations: Astarac (S1), Isle Arne (S2), Mauvezin (S3), Bives (S4), and St-Antoine (S5). Figure 7 presents a synoptic of the river and the distances between the hydrometric stations, and Fig. 8 presents a map of the river and the hydrometric stations. In this case of study, we will focus on the last reach between Bives and St-Antoine, with a length of 40km. The aim of this study is to compare the outflow routed by the RSM and the measured outflow.  hours starting at t = 0h , using a computation interval of one hour. The RSM was calibrated and the parameters TT, α, andS ′ 0 were found. The calculated parameters were TT = 11h, α = 0.82, andS ′ 0 = 6.23m 3 . In order to compare the RSM and the Muskingum model, the latter was also calibrated, and the calculated parameters were K = 58412.7, x = 0.07 and T = 10580s ≈ 3h. Figure 9 shows that the calibration of the RSM and Muskingum model indicates good performance. In fact, the errors for RSM and Muskingum model are only of: 5.07% and 4.48% , respectively.
Simulation. Once the RSM and Muskingum model were calibrated on the last reach of the Arras river, they were tested on other hydrological events in order to examine the constancy of the methodology and the parameters found in the calibration stage. The hydrograph chosen for the simulation was a 320 hours' event measured in the Bives and St-Antoine hydrometric stations during the period 06/05/2018 − 20/05/2018 . The parameters used to simulate the outflow at St-Antoine, were the ones found in the calibration stage. Recorded inflow, recorded outflow and simulated outflows with RSM and Muskingum model are illustrated in Fig. 10. Figure 10 shows that the outflow simulated with the RSM is close to the measured one. It also shows that RSM results are similar to Muskingum's results. The errors and root-mean-square for RSM and Muskingum model are: (7.1%; 0.44m 3 s −1 ) and (6.6%; 0.43m 3 s −1 ) , respectively, see Table 4 and Table 5. The simulation results analysis shows that the proposed routing model can be readily calibrated and provides pertinent results. The Muskingum routing model is convenient in deriving the outflow with given inflow hydrograph, however the result can be greatly affected by inappropriate K and X. There is no general rule of thumb of how the Muskingum K and X should be determined, thus uncertainties remain if these parameters are incorrected used. The routing method  www.nature.com/scientificreports/ incorporated with network flow presented in this paper is to minimize deviations between the observed and derived hydrograph before applying for outflow hydrograph derivation.
Operative results. The RSM model has been integrated into the operational management tool Rio used by the -CACG following the functional principle described in the Fig. 11. Feedback on measurement are operated directly by the optimization model. The feedback on progress loop previously presented on the Fig. 1, is reduced to non-hydraulic models such as streamflow or withdrawal models. This principle is more detailed and illus-  Table 5. Parameter of the simulation results for the Arrat's reach.    37 where it was implemented in order to improve initial conditions for streamflow prediction using rainfall reconstruction algorithms. RSM model is actually in operation on about thirty managed rivers in the south of France presented in Fig. 12. The Neste system is composed by the rivers in the yellow zone which form a singular system alimented by a single canal. All these rivers are left tributaries of the Garonne river except the most westerly river that provides another river named Adour (tributaries represented in orange in Fig. 12). The green and purple rivers in Fig. 12 constitute part of the right tributaries of the Garonne river.

RSM
Rio integrates models to managed dams and canal intakes: • 194 watersheds are modeled by rainfall-runoff model based on 3 interconnected reservoirs.
• 84 RSM models are used to represent the flow dynamics.
• Withdrawals are partially measured and forecasts are statistically interpolated.
Mid-term weather forecasts are obtained from GFS model (10 days). For the current day French model AROME rainfall forecasts are used. For both, 8 representative weather stations are consulted.
The rainfall event presented herein (Figs. 13, 14, 15) for 20 of the managed rivers occurred between the 22/10/2019 and the 25/10/2019. During this period, weather forecast was particularly uncertain due to the stormy conditions. On the screenshot herein, discharge measurements are colored, full and bold lines. The model outputs are represented by a surface of the same color as that used for the measurement station. Vertical dotted green line corresponds to the time of the screen capture. In order to be accurate, the model outputs have to match the measured flows before this line and the forecast flows after. For each river, three screenshots from 22/10/2019 at 6PM, 3/10/2019 at 6PM and 25/10/2019 at 10AM are displayed.

Conclusion
Real-time operation of multireservoir systems requires flood routing in river reaches. This issue has usually been solved by coupling a simulation model to the optimization model when it is based on an evolutionary approach. Conceptual models for channel routing like the Muskingum model can also be imbedded within the mathematical formulation of the optimization algorithm when it is linear. Nonetheless, even if network flow programming is a special form of linear programming, its special structure inhibits the integration of the Muskingum flood routing model in the network flow problem formulation. A new flood routing model destined primarily to be coupled with a network flow structure is developed. The routing model is considered as a network flow and the parameters are calibrated using a genetic algorithm. The model was approved on the Wilson standard test, studied by other researchers and known to present a nonlinear relationship between weighted discharge and storage. The methodology was tested and provided satisfying results on reaches of many rivers in the south west of France. This flood routing model was integrated to various operational systems since two years on several rivers and more recently experimented on the whole system for which calibration and weather forecast accuracy improvements are needed and are actually under studies. Thus, sensibility analysis hasn't been studied in this paper, because of its importance, it is considered as the main perspective of this work.   www.nature.com/scientificreports/ 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/.