# Network flow and flood routing model for water resources optimization

Ayoub Tahiri, Daniel Che, David Ladeveze, Pascale Chiron, Bernard

Archimède

► To cite this version:

Ayoub Tahiri, Daniel Che, David Ladeveze, Pascale Chiron, Bernard Archimède. Network flow and flood routing model for water resources optimization. Scientific Reports, 2022, 12 (1), pp.3937. 10.1038/s41598-022-06075-0 . hal-03956137

HAL Id: hal-03956137

https://hal.science/hal-03956137

Submitted on 11 Apr 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

Check for updates

OPEN

# Network flow and flood routing model for water resources optimization

Ayoub Tahiri<fnref n="1" />,², Daniel Che³, David Ladeveze², Pascale Chiron<fnref n="1" />✉ & Bernard Archimède<fnref n="1" />

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.

One of the most important aspects of minimizing the impacts of flooding is the proper operation of flood control systems itself<fnref n="1" />. The hydraulic system includes rivers and their tributaries, catchments, and natural or artificial hydraulic structures. Operation management of these hydraulic system is very challenging since it involves conflicting objectives and complicated variables. Real-time management consists in maximizing benefits, minimizing costs, satisfying the required flows in the river and storing water in reservoirs, answering to water demands, avoiding floods, and preserving the quality of water. These management requirements cause a need for a river basin optimization model that provides appropriate results.

Modeling choices for monitoring hydraulic system depend on numerous criteria linked to the studied case. Therefore, compromises have to be done among decision scale, precision needs, expected robustness, computing cost and so on. In this paper we focus on the case of hydraulic constraints for large scale and strongly influenced water resource systems. A flood routing approach that can be integrated in the mathematical formulation of an optimization model is proposed. The flood routing model was designed in order to limit the number of iterated simulations and the complexity of the optimization problem to solve.

In this work, a comparison between this model and the Muskingum's model for the Wilson standard test and for a real case study is provided. Some applications of this approach to several French rivers are also presented.

## State of the art

Real-time operation of multireservoir systems involves various operational, hydrologic, and hydraulic considerations². For efficient operation, real-time management model should contain a flow routing procedure to predict the impacts of the observed and/or predicted inflows hydrographs on the downstream parts of the river system³. As flood waves travel from upstream to downstream, they attenuate and get delayed. Flood waves are subject to two principal movements: uniformly progressive flow and reservoir action. A uniformly progressive flow designates a shifting of the wave from upstream to downstream without a change in shape, which would occur only under ideal conditions. Reservoir action designates the modification of a flood wave by reservoir pondage. Flood routing is a technique that determines the flood hydrograph at a section of a river using an upstream hydrograph.

<footnote n="1">¹Université de Toulouse, INP-ENIT, Laboratoire Génie de Production, LGP, Tarbes, France. ²Compagnie d'Aménagement des Coteaux de Gascogne, Tarbes, France. ³Department of Civil Engineering, Ohio University, Ohio, USA. ✉email: pascale.chiron@enit.fr</footnote>
<figure><img src="image_3.png" /><figcaption>Figure 1. Interfacing a simulation model and an optimization algorithm.</figcaption></figure>

The optimization problem for the operation of multireservoir systems under flooding conditions can be stated as follows:

$$ \text{Minimize : } z = f(h, Q) \qquad (1) $$

subject to

$$ \text{Hydraulic constraints : } g(h, Q) = 0 \\ \text{Bounds on discharges : } Q \le Q \le \bar{Q} \qquad (2) $$

$$ \text{Bounds on elevations : } \underline{h} \le h \le \bar{h} $$

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 and depends on the total flood damage. Different optimization methods were proposed in the literature 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⁴ and Yeh⁵ 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⁶.

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⁷. 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⁸. In order to integrate flow routing into the reservoir operation model, many researchers interfaced a simulation model with the optimization algorithm²,⁹⁻¹² 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³ 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¹³. 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¹⁴. The typical structure of evolutionary computation makes them suitable to adapt vaguely defined objective functions and constraints¹⁵. However, optimization-simulation process is time consuming and convergence problems may occur¹⁶. 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¹⁷,¹⁸.

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¹⁹.
$$ \frac{dS}{dt} = I - O \qquad (3) $$

The storage function can take many forms. Storage can be expressed as a function of inflow, outflow or both:

$$ S = f\left(I, \frac{dI}{dt}, \dots, O, \frac{dO}{dt}, \dots\right) \qquad (4) $$

The Muskingum flood routing model is the most used model¹⁹. The absolute storage is a function of both outflow and inflow discharges. The evolution in time of the absolute storage is expressed by:

$$ S_t = K(xI_t + (1-x)O_t) \qquad (5) $$

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:

$$ O_t = C_0 I_{t-TT} + C_1 I_{t-TT-1} + C_2 O_{t-1} \qquad (6) $$

where $C_0$, $C_1$ and $C_2$ are constants that are computed from K and x and $\Delta T$²⁰.

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²¹ formulated a theoretical recursive linear programming model for the operation of flood control, using the Muskingum method for channel routing. Hsu and Wei²² 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.²³ 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²⁴. 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²⁵). The LP is expressed as:

$$ \text{Minimize : } z = \sum_{j=1}^{n} c_j x_j \qquad (7) $$

Subject to:

$$ \sum_{j=1}^{n} a_{ij}x_j \leq b_i \quad \text{for } i = 1, \dots, m \\ x_j \geq 0 \quad \text{for } j = 1, \dots, n \qquad (8) $$

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²⁶, 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³,²⁷. 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:

$$ \text{Minimize : } z = \sum c_{ji}q_{ij} \quad \text{for all arcs} \qquad (9) $$

Subject to:

$$ \sum q_{ij} - \sum \bar{q}_{ji} = 0 \quad \text{for all nodes} \\ \bar{q}_{ij} \leq q_{ij} \leq \bar{q}_{ij} \quad \text{for all arcs} \qquad (10) $$

where $q_{ij}$ is the flow rate in the arc connecting node i to node j; $c_{ij}$ is a cost for $q_{ij}$; $\bar{q}_{ij}$ and $\bar{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²⁸. Braga and Barbosa²⁹ report on inclusion of channel routing into multiple time step optimization using an advanced
<figure><img src="image_4.png" /><figcaption>Figure 2. Residual storage schematic representation.</figcaption></figure>

Network Simplex Solver that can handle non-network side constraints required for inclusion of channel routing. Nonetheless, non-network constraints may disturb the allocation priorities as stated by Ferreira³⁰ and Chou and Wu³¹. To eliminate the aforementioned disadvantages and difficulties of the use of the existing flood routing models in water resource optimization models based on a network approach, the conceptual model developed herein can be proposed. The flood routing model presented here uses the computational properties of the network flow technique, and can be coupled to a network flow structure.

## 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:

$$S'_{t+1} = S'_{t} + I_{t-TT} - O_{t} \quad (11)$$

where $TT$ is the Transit Time. This parameter represents the time from when the inflows $\{I_{t-TT}, I_{t-TT-1}, \dots, I_0\}$ impact the outflow at time t. Hence, inflows $\{I_t, I_{t-1}, \dots, 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}$:

$$O_t = (1 - \alpha)(S'_t + I_{t-TT}) \quad (12)$$

The parameters of the model are the reach's Transit Time ($TT$), the proportionality coefficient ($\alpha \in [0:1]$) and the initial residual storage $S'_0$.

The proportionality coefficient $\alpha$ 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³². For an arc $e_{ij}$, $i$ is the origin node, and $j$ is the end node. Let $\gamma(n)$ and $\gamma^{-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 $\Phi_e$. For each node $n$ of $V$, except for $S$ and $T$, the conservation of flow is satisfied:

$$\forall n \in V \setminus \{S, T\} \sum_{e \in \gamma(n)} \Phi_e = \sum_{e \in \gamma(n)^{-1}} \Phi_e \quad (13)$$

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
<figure><img src="image_5.png" /><figcaption>Figure 3. Network model corresponding to a reach, with $TT = 1$ and $Horizon = 4$.</figcaption></figure>

of the reach at time $t$, respectively. Inflow at instant $t$ is represented by the flow carried by the arc $e_{SU_t}$. Inflow at instant $t$ joins the reservoir at instant $t + TT$ through the arc $e_{U_t R_{t+TT}}$. The reservoir node at instant $t + TT$ preserves a proportional portion of the incoming flows: $S'_{(t+TT)+1} = \alpha(S'_{t+TT} + I_t)$, and releases the left portion: $O_{t+TT} = (1 - \alpha)(S'_{t+TT} + I_t)$.

**Calibration of the model.** For the simulated outflows to be close to the measured downstream flows, the model's parameters of the flood routing model have to be calibrated. The objective function to minimize is stated as follows:

$$ \text{Minimize} \quad \sum_{t=0}^{Horizon-1} (O_t^{\text{Measured}} - O_t^{\text{RSM Model}})^2 \qquad (14) $$

Subject to:

$$ \begin{array}{l} \forall n \in V \setminus \{S, T\} \quad \sum_{e \in \gamma(n)} \Phi_e = \sum_{e \in \gamma(n)^{-1}} \Phi_e \\ S'_{t+1} = \alpha(S'_{t} + I_{t-TT}) \\ \alpha \in [0:1] \\ TT \ge 0 \\ S'_{0} = S'_{Horizon-1} \end{array} \qquad (15) $$

where $O_t^{\text{Measured}}$ and $O_t^{\text{RSM Model}}$ are the measured and computed outflow at time $t$, respectively.

The unknowns of the problem are: the reach's transit time ($TT$), the distribution coefficient ($\alpha$) and the initial residual storage $S'_0$. The final residual storage $S'_{Horizon-1}$ should be equal to the initial one $S'_0$ in order to conserve the volume of the routed hydrograph.

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, \alpha, 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-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.
<table><caption>**Table 1.** Recorded and reconstructed outflow for the Wilson event.</caption><thead><tr><th>T(h)</th><th>Recorded inflow (m³s⁻¹)</th><th>Recorded outflow (m³s⁻¹)</th><th>Simulated outflow (RSM) (m³s⁻¹)</th><th>Simulated outflow (Muskingum) (m³s⁻¹)</th></tr></thead><tbody><tr><td>0</td><td>22</td><td>22</td><td>17.30</td><td>20.95</td></tr><tr><td>6</td><td>23</td><td>21</td><td>18.74</td><td>21.26</td></tr><tr><td>12</td><td>35</td><td>21</td><td>19.75</td><td>22.23</td></tr><tr><td>18</td><td>71</td><td>26</td><td>20.78</td><td>26.86</td></tr><tr><td>24</td><td>103</td><td>34</td><td>24.57</td><td>38.75</td></tr><tr><td>30</td><td>111</td><td>44</td><td>36.33</td><td>54.67</td></tr><tr><td>36</td><td>109</td><td>55</td><td>54.38</td><td>68.60</td></tr><tr><td>42</td><td>100</td><td>66</td><td>71.10</td><td>78.38</td></tr><tr><td>48</td><td>86</td><td>75</td><td>82.79</td><td>83.35</td></tr><tr><td>54</td><td>71</td><td>82</td><td>88.63</td><td>83.59</td></tr><tr><td>60</td><td>59</td><td>85</td><td>88.75</td><td>80.16</td></tr><tr><td>66</td><td>47</td><td>84</td><td>84.35</td><td>74.51</td></tr><tr><td>72</td><td>39</td><td>80</td><td>77.39</td><td>67.45</td></tr><tr><td>78</td><td>32</td><td>73</td><td>68.92</td><td>60.13</td></tr><tr><td>84</td><td>28</td><td>64</td><td>60.29</td><td>53.00</td></tr><tr><td>90</td><td>24</td><td>54</td><td>52.12</td><td>46.61</td></tr><tr><td>96</td><td>22</td><td>44</td><td>44.99</td><td>40.90</td></tr><tr><td>102</td><td>21</td><td>36</td><td>38.84</td><td>36.14</td></tr><tr><td>108</td><td>20</td><td>30</td><td>33.82</td><td>32.31</td></tr><tr><td>114</td><td>19</td><td>25</td><td>29.95</td><td>29.19</td></tr><tr><td>120</td><td>19</td><td>22</td><td>26.97</td><td>26.64</td></tr><tr><td>126</td><td>18</td><td>19</td><td>24.60</td><td>24.67</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Wilson event: recorded inflow and outflow, reconstructed outflow for the RSM and the Muskingum models.</figcaption></figure>

The root-mean-square (RMS) and the error (Error) are computed using Eq. (16):

$$ RMS = \sqrt{\frac{\sum_i (Q_i^{\text{observed}} - Q_i^{\text{model}})^2}{\text{Horizon}}} \\ Error = 100 \times \frac{\sum_i |Q_i^{\text{observed}} - Q_i^{\text{model}}|}{\sum_i Q_i^{\text{observed}}} \quad (16) $$

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₀ and C₁ are negligible compared to C₂ which shows that the outflow at time t + 1 depends
<table><caption>**Table 2.** Calibration results for the Wilson event.</caption><thead><tr><td></td><td>RSM</td><td>Muskingum</td></tr></thead><tbody><tr><td>RMS (m³s⁻¹)</td><td>4.73</td><td>4.48</td></tr><tr><td>Error (%)</td><td>8.37</td><td>9.32</td></tr></tbody></table>

<table><caption>**Table 3.** Parameter of the calibration results for the Wilson event.</caption><thead><tr><td></td><td>TT(h)</td><td>α</td><td>S<sub>0</sub> (m³)</td><td>ΔT (s)</td><td>K</td><td>X</td><td>C<sub>0</sub></td><td>C<sub>1</sub></td><td>C<sub>2</sub></td></tr></thead><tbody><tr><td>RSM</td><td>11</td><td>0.94</td><td>270.13</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Muskingum</td><td>5</td><td></td><td></td><td>10800</td><td>229117.15</td><td>0.023</td><td>12.345E-05</td><td>0.047</td><td>0.952</td></tr></tbody></table>

<figure><img src="image_7.png" /><figcaption>**Figure 5.** Absolute storage vs. outflow (Muskingum model).</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Residual storage vs. outflow (RSM).</figcaption></figure>

primarily on outflow at time *t*. Even if in the RSM model the influence of the output is only depending on the residual storage state and not directly on the output flow as in the Muskingum model, the observed results are similar. Referring to 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
<figure><img src="image_9.png" /><figcaption>**Figure 7.** Synoptic of the Arrats' river (with KP indicating the kilometric point).</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>**Figure 8.** Map of the Arrats' river and the hydrometric stations created with QGIS version 3.16 (www.qgis.org).</figcaption></figure>

highlights the non-linear relationship between the absolute storage and the outflow assumed in the Muskingum model. The plot consists in a loop formed by a forward and a reverse path. The loop reflects the non-symmetrical relationship existing when the reach is storing or emptying. On the other hand, in 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.
<figure><img src="image_11.png" /><figcaption>**Figure 9.** Calibration: Recorded inflow in Bives, recorded outflow in St-Antoine, reconstructed outflow for the RSM and the Muskingum models.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Simulation: Recorded inflow in Bives, recorded outflow in St-Antoine, reconstructed outflow for the RSM and the Muskingum models.</figcaption></figure>

**Calibration.** The hydrograph chosen for the calibration is a 200 hours' event measured in the Bives and St-Antoine hydrometric stations during the period 03/07/2018 – 11/07/2018. The total simulation time was 200 hours starting at $t = 0h$, using a computation interval of one hour. The RSM was calibrated and the parameters $TT, \alpha, \text{and } S_0'$ were found. The calculated parameters were $TT = 11h, \alpha = 0.82, \text{and } S_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 \text{ and } \Delta T = 10580s \approx 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^3s^{-1})$ and $(6.6\%; 0.43m^3s^{-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 incorrectly used. The routing method
<table><caption>**Table 4.** Simulation results for the Arrat's reach.</caption><thead><tr><td></td><td>RSM</td><td>Muskingum</td></tr></thead><tbody><tr><td>RMS (m³s⁻¹)</td><td>0.44</td><td>0.43</td></tr><tr><td>Error (%)</td><td>7.1</td><td>6.6</td></tr></tbody></table>

<table><caption>**Table 5.** Parameter of the simulation results for the Arrat's reach.</caption><thead><tr><td></td><td>TT(h)</td><td>α</td><td>S'₀(m³)</td><td>ΔT (s)</td><td>K</td><td>X</td><td>C₀</td><td>C₁</td><td>C₂</td></tr></thead><tbody><tr><td>RSM</td><td>11</td><td>0.82</td><td>16</td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Muskingum</td><td>5</td><td></td><td></td><td>10580</td><td>58412.7</td><td>0.0717</td><td>0.0185</td><td>0.1593</td><td>0.8222</td></tr></tbody></table>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** Scheme for the optimization strategy integrating hydraulic models.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>**Figure 12.** Hydrographic perimeter of the operational water resource system, created with QGis version 3.16 (www.qgis.org).</figcaption></figure>

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-
<figure><img src="image_15.png" /><figcaption>**Figure 13.** Right tributaries of the Garonne dashboards screenshots.</figcaption></figure>

trated in Tahiri et al.³⁷ 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.

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.
<figure><img src="image_16.png" /><figcaption>**Figure 14.** Left tributaries of the Garonne dashboards screenshots.</figcaption></figure>
<figure><img src="image_17.png" /><figcaption>Figure 15. Adour tributaries dashboards screenshots.</figcaption></figure>

Received: 17 September 2021; Accepted: 21 January 2022

Published online: 10 March 2022

References

1. Wei, C.-C. & Hsu, N.-S. Multireservoir real-time operations for flood control using balanced water level index method. *J. Environ. Manag.* **88**, 1624–1639 (2008).

2. Ahmed, E.-S.M.S. & Mays, L. W. Model for determining real-time optimal dam releases during flooding conditions. *Nat. Hazards* **65**, 1849–1861 (2013).

3. Wurbs, R. A. *Comparative evaluation of generalized river/reservoir system models* (Tech. Rep, Texas Water Resources Institute, 2005).

4. Labadie, J. W. Optimal operation of multireservoir systems: state-of-the-art review. *J. Water Resour. Plan. Manag.* **130**, 93–111 (2004).

5. Yeh, W.W.-G. Reservoir management and operations models: A state-of-the-art review. *Water Resour. Res.* **21**, 1797–1818 (1985).

6. Koch, H. & Grünewald, U. A comparison of modelling systems for the development and revision of water resources management plans. *Water Resour. Manag.* **23**, 1403 (2009).

7. Chow, V. T., Maidment, D. R. & Mays, L. W. *Applied Hydrology* (McGraw-Hill Series in Water Resources and Environmental Engineering (McGraw-Hill, New York, 1988).
8. Barat, R., Akbari, G. H. & Rahimi, S. Flood routing of an unmanaged river basin using muskingum-cunge model; Field application and numerical experiments. *Casp. J. Appl. Sci. Res.* **2**, 2 (2013).

9. Che, D. & Mays, L. W. Development of an optimization/simulation model for real-time flood-control operation of river-reservoirs systems. *Water Resour. Manag.* **29**, 3987–4005 (2015).

10. Yazdi, J. & Neyshabouri, S. S. A simulation-based optimization model for flood management on a watershed scale. Water Resources Management. *Publ. Eur. Water Resour. Assoc. (EWRA)* **26**, 4569–4586. https://doi.org/10.1007/s11269-012-0167-1 (2012).

11. Malekmohammadi, B., Zahraie, B. & Kerachian, R. A real-time operation optimization model for flood management in river-reservoir systems. *Nat. Hazards* **53**, 459–482. https://doi.org/10.1007/s11069-009-9442-8 (2010).

12. Unver, O. I. & Mays, L. W. Model for real-time optimal flood control operation of a reservoir system. *Water Resour. Manag.* **4**, 21–46. https://doi.org/10.1007/BF00429923 (1990).

13. Carson, Y. & Maria, A. Simulation optimization: Methods and applications. In *Winter Simulation Conference Proceedings.*, 118–126 (1997).

14. Ahmad, A., El-Shafie, A., Razali, S. F. M. & Mohamad, Z. S. Reservoir optimization in water resources: A review. *Water Resour. Manag.* **28**, 3391–3405. https://doi.org/10.1007/s11269-014-0700-5 (2014).

15. Rani, D. & Moreira, M. M. Simulation-optimization modeling: A survey and potential application in reservoir systems operation. *Water Resour. Manag.* **24**, 1107–1138 (2010).

16. Dhar, A. & Datta, B. Optimal operation of reservoirs for downstream water quality control using linked simulation optimization. *Hydrol. Process.* **22**, 842–853, https://doi.org/10.1002/hyp.6651 (2008). https://onlinelibrary.wiley.com/doi/pdf/10.1002/hyp.6651.

17. Yin, J. et al. Blending multi-satellite, atmospheric reanalysis and gauge precipitation products to facilitate hydrological modelling. *J. Hydrol.* **593**, 125878. https://doi.org/10.1016/j.jhydrol.2020.125878 (2021).

18. Yin, J. et al. Does the hook structure constrain future flood intensification under anthropogenic climate warming?. *Water Resour. Res.* **57**, e28491. https://doi.org/10.1029/2020WR028491 (2021).

19. Subramanya, K. *Engineering Hydrology* (Tata McGraw-Hill Publishing, 2007).

20. Al-Humoud, J. M. & Esen, I. I. Approximate methods for the estimation of muskingum flood routing parameters. *Water Resour. Manag.* **20**, 979–990 (2006).

21. Windsor, J. S. Optimization model for the operation of flood control systems. *Water Resour. Res.* **9**, 1219–1226. https://doi.org/10.1029/WR009i005p01219 (1973).

22. Hsu, N.-S. & Wei, C.-C. A multipurpose reservoir real-time operation model for flood control during typhoon invasion. *J. Hydrol.* **336**, 282–293. https://doi.org/10.1016/j.jhydrol.2007.01.001 (2007).

23. Kumar, D. N., Baliarsingh, F. & Raju, K. S. Optimal reservoir operation for flood control using folded dynamic programming. *Water Resour. Manag.* **24**, 1045–1064 (2010).

24. Needham, J. T., Watkins, D. W. Jr., Lund, J. R. & Nanda, S. Linear programming for flood control in the Iowa and des Moines rivers. *J. Water Resour. Plan. Manag.* **126**, 118–127 (2000).

25. Ponnambalam, K., Vannelli, A. & Unny, T. An application of karmarkar's interior-point linear programming algorithm for multi-reservoir operations optimization. *Stoch. Hydrol. Hydraul.* **3**, 17–29 (1989).

26. Ahuja, R. K., Magnanti, T. L., Orlin, J. B. & Reddy, M. Applications of network optimization. *Handbooks Oper. Res. Manag. Sci.* **7**, 1–83 (1995).

27. Kuczera, G. Network linear programming codes for water-supply headworks modeling. *J. Water Resour. Plan. Manag.* **119**, 412–417 (1993).

28. Ilich, N. Shortcomings of linear programming in optimizing river basin allocation. *Water Resour. Res.* **44**, 2 (2008).

29. Braga, B. & Barbosa, P. S. Multiobjective real-time reservoir operation with a network flow algorithm. *J. Am. Water Resour. Assoc.* **37**, 837–852 (2001).

30. Ferreira, I. C. L. *Deriving Unit Cost Coefficients for Linear Programming-Driven Priority-Based Simulations* (University of California, 2007).

31. Chou, F.-F. & Wu, C.-W. Determination of cost coefficients of a priority-based water allocation linear programming model-a network flow approach. *Hydrol. Earth Syst. Sci.* **18**, 1857–1872 (2014).

32. Tahiri, A., Ladeveze, D., Chiron, P., Archimede, B. & Lhuissier, L. Reservoir management using a network flow optimization model considering quadratic convex cost functions on arcs. *Water Resour. Manag.* **2**, 1–14 (2018).

33. Wilson, E. *Engineering Hydrology* (Macmillan, 1983).

34. Gill, M. A. Routing of Floods in River Channels. *Hydrology Research* **8**, 163–170, https://doi.org/10.2166/nh.1977.0013 (1977). https://iwaponline.com/hr/article-pdf/8/3/163/9361/163.pdf.

35. Singh, V. P. & McCann, R. C. Some notes on muskingum method of flood routing. *J. Hydrol.* **48**, 343–361 (1980).

36. Yoon, J. & Padmanabhan, G. Parameter estimation of linear and nonlinear muskingum models. *J. Water Resour. Plan. Manag.* **119**, 600–610 (1993).

37. Tahiri, A., Ladeveze, D., Chiron, P. & Archimede, B. Improving the characterization of initial conditions for streamflow prediction using a precipitation reconstruction algorithm. *IFAC-PapersOnLine* **51**, 19–24 (2018).

## Author contributions

A.T., D.C. and D.L. conceived the experiments. A.T. and D.L. developed the algorithms, conducted the experiments and wrote the main manuscript text. The work was proposed from A.T., D.C. and D.L. suggestions and revised by P.C. and B.A. All authors reviewed the manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Correspondence and requests for materials should be addressed to P.C.

Reprints and permissions information is available at www.nature.com/reprints.

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
<figure><img src="image_18.png" /></figure>

**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://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022