## Abstract

This paper presents the benefits of using mathematical optimization for reservoir operation based on the assumed availability of short-term runoff forecasts. The novelty is the inclusion of the SSARR hydrological routing as optimization constraints in multiple time step optimization, where the routing coefficients are adjusted dynamically as functions of the channel flows. The paper shows significant reduction of downstream peak flows in flood-prone areas even with a forecast horizon of only 2 days, and it also includes the results of testing the effects of different lengths of forecasting horizons on model results. The case study is conducted on the Damodar River Basin in the Indian State of West Bengal, where basin development started in the 1950s, with flood protection of the downstream river valley as the highest management priority, in addition to water supply and hydro power. The solution methodology and the model results presented in this paper pave the way for eventual introduction to automated management of reservoir outflows that could shape water resources industry in much the same way that auto-pilot and driverless cars are revolutionizing the transportation industry, assuming that runoff forecasting capabilities continue to improve.

## HIGHLIGHTS

Use of linear programming guarantees finding global optimums for large multi-reservoir systems.

Hydrological channel routing embedded in the optimization solver adjusts the routing parameters based on optimized channel flows.

The paper demonstrates the capabilities of the WEB. BM model which can be accessed and tested by any interested practitioner free of charge.

## INTRODUCTION

The use of computer models is a common component in the process of creating river basin management plans. Early review papers of the available models and their capabilities have been provided by Yeh (1985) and Labadie (2004), while more recent reviews are available from Rani & Moreira (2010) and Dobson *et al.* (2019), which focus on the advantages of optimization models over those based on simulation. In spite of the long history of diverse optimization model development, their actual use in real-time reservoir operation is almost non-existent. Some of the principal reasons for this are (a) the lack of reliable runoff forecast and (b) the inability of most optimization models to combine reservoir release decisions with their routed effects at downstream control points. The effects of channel routing should be combined with the effects of additional tributaries and downstream hydraulic structures. Also, the propagation of reservoir releases to the downstream control points may take longer than the short length of the calculation time steps that are required for adequate modeling of floods, thus requiring simultaneous solutions over multiple calculation time steps – a feature that is not available to simulation models that have traditionally been used to model flood events. This paper deals with the issue (b), i.e., the focus of the paper is not the runoff forecasting capabilities, but rather the required technical capabilities of river basin management models, where reservoir operation takes the center stage as the principal mechanism for taming the otherwise often hostile natural flow regime. Since reservoir operation is an essential part of river basin management, the terms ‘reservoir operation models’ and ‘river basin management models’ in this paper are assumed to refer to the same class of mathematical tools.

The first reservoir operation models such as the HEC-5 model (US Corps of Engineers 1998) appeared in the literature more than 50 years ago, and most of them were initially based on the concept of using simulation tools that allowed for inclusion of ‘what-if’ rules. Although some simulation models such as the Mike Hydro Basin (DHI 2023) and RIBASIM (DeltaRes 2023) are still used, their use has become increasingly difficult with the growth of the complexity of river basin management issues. These models often fail to simultaneously apply the zoning management concept on multiple reservoirs and multiple water users, so as to include not only the storage management policies, but also water rationing policies and deficit sharing in dry years as a function of the available storage and hydrologic conditions in the basin. Efforts to offset these difficulties led to the development of models based on the Network Flow Algorithms (NFAs), which are also based on the theory of linear programming (LP). The first NFA-based model was launched by the Texas Water Development Board (SIMYIELD) in the 1970s, followed by a number of models that have become popular over time, such as the MODSIM (Labadie *et al.* 1986), Aquatool (Andreu *et al.* 1996), E-Water (EWater Academy 2022), or the REALM model (Victoria State Government 2022). These models all have the capability to apply the zoning concept with reservoir rule curves in combination with the zoning concept on demands, such that demand management is also coupled with the management of water supply. Typically, the purpose of a reservoir operation model is to allocate water based on the allocation policies that are provided by the model users. While this approach may be helpful for managing storage during droughts, its usefulness is questionable during floods, where several upstream reservoirs can be managed concurrently, while inflow hydrographs into each reservoir always have unique temporal and spatial distributions. An additional limitation of NFAs is that most of them do not allow losses or gains of flow along a flow link in the network, which significantly complicates the inclusion of hydrological channel routing in the solution process.

Some river basin models, such as the WEAP model (Yates *et al.* 2005), also include an internal rainfall-runoff module that can generate runoff estimates. However, when it comes to river basin management models, this feature is more of an exception rather than the rule. Most river basin management models do not generate runoff estimates based on rainfall data, but rather rely on a previously prepared runoff time series for a period of multiple years or for individual hydrologic events. This makes the evaluation of the reservoir operation models easier, since the model performance can be separated from the ability to accurately predict runoff forecasts. The rainfall-runoff models should be calibrated and validated independently, and they are not discussed further in this paper. The principal question that this paper tries to address is the ability to improve real-time flood management assuming that an accurate short-term (3 days) runoff forecast is available, and to demonstrate the fact that an optimization-based computer model can be a powerful tool for real-time decision-making if coupled with a reliable runoff forecasting tool.

## PREVIOUS DEVELOPMENTS AND THE IMPORTANCE OF HYDROLOGICAL ROUTING

The use of optimization algorithms for reservoir operation during floods was first postulated as a conceptual LP model by Windsor (1973), and later refined by Needham *et al.* (2000). These applications were based on the development of flood damage functions, which were piece-wise linearized and provided to the LP optimization routines that are part of the model. Both papers detail the ability to derive solutions over multiple time steps by an LP solver, which assumes the known inflow hydrographs. However, the issue of channel routing has been limited to the use of classical Muskingum approach, where the routing coefficients are constant and calibrated on the basis of a particular historic runoff event, without discussing the effect of using these constant coefficients on alternate regulated flows that can be significantly different from the hydrologic events used during calibration of routing coefficients. Larger basins with multiple reservoirs offer more possibilities for finding better solutions using joint reservoir operation, but they require inclusion of hydrologic channel routing as constraints to optimization. Some recent efforts to address this have been reported by Tahiri *et al.* (2022), although the routing constraints were associated with iterations imposed within an NFA-based model, which is not an ideal tool to represent the dynamics of channel flow routing since flows along flow links must remain the same, requiring modeling of hydrological transformation of flood ways as a sequence of additions and subtractions to otherwise constant flows in flow links. Other recent attempts to provide comprehensive modeling solution platforms are designed by DetaRes (RTC-Tools 2022) available as a Python package, although the particular details on the implementation of the hydrologic routing schemes are dealt with in a cursory manner requiring more information for independent evaluation.

*O*is related to the current inflow into the channel

_{t}*I*and the previous time step inflows and outflows

_{t}*I*and

_{t−1}*O*. The principal issue with the Muskingum approach is that the routing coefficients

_{t−1}*C*are calibrated for a specific hydrological event, while for longer term modeling that involves transition from medium flows to high flows these coefficients change, since they are a function of the travel time over the length of the river reach, which is dependent on the magnitude of flow. The Muskingum-Cunge (Cunge 1969) method introduces channel geometry (cross-sections), Manning's n and the channel gradient as additional information that is required to enable variation of the routing coefficients as a function of the travel time along the reach. However, most of the required data are usually not available or subject to rough estimates. There are also reports in the literature that claim that the Muskingum-Cunge method is unable to preserve volume conservation in a satisfactory manner (Perumal & Sahoo 2008). One of the most elegant variants of the Muskingum routing method has been implemented in the Stream Synthesis And Reservoir Regulation (SSARR) model (US Army Corps of Engineers 2021), based on the routing approach originally proposed by Williams (1969). As with the other river routing methods, the governing equation is related to channel storage change over a time step, which is a function of average inflow and outflow:

_{i}*O*, multiplying by

_{t−1}*t/(O*and by letting

_{t}-O_{t−1})*ΔS/(O*

_{t}-O_{t−1})*=*

*TS*, the above equation becomes:where the term

*TS*represents the average travel time along a river reach for given flow conditions, evaluated either by reading from a previously generated table of

*TS*vs. flow values for a given river reach. In his original publication of this method, Williams (1969) defined the travel time along a reach on the basis of the updated outflow

*O*from the reach at the end of the current time step

_{t}*t*, implying that

*T*corresponds to the time that would be required to ‘empty’ the current storage in the river reach assuming constant outflow

_{s}*O*. Expression (3) can also be converted to the following explicit form:

_{t}The above routing model is non-linear, so it typically requires a few numerical iterations where *O _{t}* is recalculated and the travel

*T*is updated using the interpolated from the table of

_{s}*T*vs. flows, thus allowing the update of routing coefficients in Equation (4). The process can be repeated a few times until convergence based on a desired tolerance criterion is achieved. The structural form of Equation (4) is identical in form to the well-known Muskingum linear routing form shown in Equation (1), while the coefficients are calculated using a different approach. It should be noted that the SSARR routing coefficients

_{s}*C*listed in Equation (4) sum up to 1.0, which is also the condition for the Muskingum routing coefficients, and which confirms the conservation of volume. The advantage of the SSARR routing method is that the routing coefficients are a function of flow. As such, they are internally updated as the model simulation moves from low flow season to high flow season. Also, the required travel time vs. flow information for a river reach is much easier to obtain, since it requires the average cross-sectional velocity (which is determined at each flow measurement) and the length of a river reach. For very high flows observed during flood peaks, the hourly timing of the peak is observed for the same hydrograph as it propagates downstream between the two observational points, which provides estimates of travel times during high flows when flow measurements are not conducted due to the associated safety risks. Travel times vs. flow can also be obtained from HEC-RAS modeling by assuming various flow rates, or by obtaining the mean cross-sectional velocities by using the Manning's equation are applying them to the length of the entire river reach.

_{i}Note that a typical distance the outflows from the reservoir will travel within 1 day is around 85 km (assuming close to 1 m/s of mean cross-sectional stream velocity), which limits the distance where the storage releases will have effects on downstream flows within a day. Most river basins are longer than 85 km, which implies that storage releases made on day 1 will be felt at downstream locations in subsequent days, depending on the distance. These considerations regarding the travel times are valid both for steady-state model runs as well as for modeling that includes hydrological channel routing. In this context, deriving model releases for individual time steps may or may not meet the downstream operational goals in subsequent time steps, since each time step is an independent model solution. The concept of deriving simultaneous solutions to multiple time steps was introduced to overcome this limitation.

## MULTIPLE TIME STEP OPTIMIZATION METHOD

*et al.*(2000), proper solution of flood routing optimization is based on the assumed knowledge of inflow hydrographs over a specified time horizon. There are no perfect forecasts in reality, although the forecasting tools have been improving recently. The principal approach in testing the reservoir optimization models is to use historical or synthetically generated hydrographs which are assumed to be known. This removes the issues of uncertainties of inflows, thus allowing for an easier comparison of various optimization techniques. Some models such as the RiverWare (Zagona

*et al.*2001) and OASIS (Randall

*et al.*1997) are capable of solving optimization for multiple time steps simultaneously, using the dynamic approach shown in Figure 2, but none of these models are able to include hydrological channel routing as a constraint in the multiple time step solution process the way it was demonstrated in the case of the WEB.BM model (Ilich 2021). The dashed lines in Figure 2 represent the carry over storage between consecutive time steps, the triangle is a symbolic representation of a reservoir while the green square is a symbolic representation of consumptive water use (e.g. an irrigated area).

*P*associated with flooding would have a negative sign for any flow that exceeds the full bank capacity, indicating that benefits would be maximized if the reservoirs could be operated such that the overbank flow is minimized, or completely avoided if possible. The objective function would be applicable for all time steps that are solved simultaneously, and for all river reaches that may be associated with the flood damage, which explains the double summation over both time and space:

_{i}Equation (6) denotes the mass balance at network node *n*, with *A _{i}* representing the incidence matrix coefficient that is typically equal to either −1, 0 or 1, with −1 and 1 denoting the incoming and outgoing flows for network arcs that are associated with a given node

*n*. In multiple time step optimization, the mass balance Equation (6) is applied to each node

*n*in the network and to each time step

*t*over the solution horizon of

*m*time steps. Expression (7) represents the upper bound on flows, which can be either set to constants

*U*representing for example the maximum storage or canal capacity that should not be exceeded, or it can be defined as a function of flow in some other network component

_{i,t}*k*, such as for example the return flow from irrigation being a linear function of consumptive use or the maximum outflow through an outlet structure such as an orifice or hydro power plant being a piece-wise linearized function of the available storage for a simulated time step.

A complete list of components that can be included in the WEB.BM model is available in previous publications (Ilich 2021). It should be noted that both reservoir and hydrological channel routing are included directly in the optimization process, with only the channel routing constraints requiring iterative calls of the optimization solver, as explained in the ‘Implementation Procedure’ section.

Due to significant reliance on the use of reservoir rule curves, many practitioners consider them to be a useful mechanism for reducing flood damage. However, rule curves were primarily introduced to reduce the chance of premature depletion of storage during dry years. They typically have a fixed shape for each year, and they are unable to offer necessary flexibility for managing floods in large multi-reservoir systems, where each incoming flood has a different spatial distribution and different timing of the peaks, thus requiring a unique response from the system that cannot be obtained from a fixed rule curve shape. The purpose of using optimization is to find the best response within a limited time horizon that was 2 days long in this study, although the time horizons of 1, 3, and 4 days were also tested to evaluate their impact on the quality of model solutions. Hence, the ultimate goal of this study was not to improve the reservoir rule curves, but rather rely on the model to develop short-term operational strategies based on the established allocation priorities and assumed availability of a short-term runoff forecast. The reasons for using LP within the WEB.BM model are listed below:

LP is the only optimization solver that guarantees finding the global optimum, due to the convexity of it search space that does not allow the existence of local optima;

There is no need to discretize decision variables into a finite number of possible states as is the case with DP applications, since LP works with decision variables that are double precision real numbers; and

There is no curse of dimensionality associated with LP. The WEB.BM can solve heavily constrained problems with over 500,000 variables in less than 11 min, while other researchers reported solving problems with close to 5 million variables within hours, as was achieved in a study aimed to set up water management plan for the entire State of California (CALVIN 2023).

It is common wisdom among Operations Research (OR) practitioners that if a problem can be successfully linearized, LP is the best choice for an optimization solver. It is for the above reasons that LP solutions are used as benchmarks to many new experimental solvers in the OR industry. This is unfortunately not the case among water resources practitioners, where respectable benchmark problems and their solutions have yet to be established.

Setting up the pricing vector *P _{i}* is easy with early LP-based applications (especially those based on NFA solvers), since the only requirement is to select the

*P*values with a relative difference between water use allocated to various stakeholders. For example, giving stakeholders A, B, and C the priority factors of 1, 2, and 3 would yield the same water allocation as would the priority factors 7, 11, and 32, since in both instances the case of maximization (or minimization of deficits) would allocate water first to stakeholder C, followed by stakeholder B and then stakeholder A. Complications started to emerge with the introduction of constraints, as originally discussed by Israel & Lund (1999), and they were further escalated by introduction of full LP solvers and various piece-wise constraints, along with multiple time step solution mode. In spite of this, selection of Pricing vector

_{i}*P*that will deliver a desirable solution that follows a prescribed operating policy is not a major problem, and it is still considered as a relatively straight forward task for most systems where operational priorities are clearly defined. Although most dams are built for multi-purpose water use, their original design usually had the prime and secondary purpose, i.e., either irrigation was the prime purpose, with hydro power being generated as a byproduct of irrigation releases, or the other way around (considering the year-round operational priorities). Flood events typically last up to one week, and they are usually associated with operation of hydro power plants at full capacity and absence of any shortage of irrigation supply flows, as is the case in this study, where three operational goals have been included in the model in the order in which they are listed:

_{i}- (1)
Flood control;

- (2)
Water supply for irrigation and industrial demands; and

- (3)
Hydropower generation.

With major irrigation canals are located downstream of all reservoirs, it is obvious that the only mechanisms that should cause lowering the storage level below the designed full supply level (FSL) are irrigation demands in the low flow season when storage inflows are lower than the downstream demands and pre-flood drawdown during the monsoon season that would minimize the flows bypassing turbines while simultaneously assisting with keeping downstream flows within the full bank capacity. The selected values of pricing vector *P _{i}* are provided in the ‘case study’ section. At this point, it should be noted that

*P*values have been determined arbitrarily by the user using several trial model runs to ensure that the model operates as expected. In addition to minimizing flood damage, the model was set up to meet the existing water supply demands and maximize generated power, implying a multi-reservoir multi-purpose optimization.

_{i}The issue of automated search for the best values of the pricing vector *P _{i}* has been addressed by many researchers in recent decades. This was precipitated by new emerging water demands that were not originally formulated when the reservoirs were initially commissioned, such as for example the need to maintain the environmental flow targets, which emerged as an issue in the past few decades. To address this, researchers introduced multi-objective optimization, delivering results in the so-called Pareto-optimal front with possibly hundreds of long-term model solutions that are considered ‘equally optimal’ with respect to two or more management goals which some researchers eventually incorporated into a more general solution procedure such as for example the Evolutionary Multi-Objective Direct Policy Search (EMODPS) proposed by Giuliani

*et al.*(2016). In contrast, reservoir operators seek simple guidelines to determine releases from storage in the short-term, and their allocation priorities in most instances are usually well defined, especially in North America, where water allocation priorities are based on the age of respective water licences – a system that is not likely to change easily due to legal constraints. Quinn

*et al.*(2019) listed several reasons for the existing gap between academia and practitioners, one of which is the lack of ability on the part of reservoir operators to fully understand the selected mathematical utility functions and the black-box tools, due to the lack of understanding the processes involved in their creation. Solution procedures such as EMODPS attempt to guide storage releases based on revising the operating rules. The work presented here is different in the sense that it removes the need to use the ‘upper rule curve’ for any storage reservoirs, but rather dynamically creates operating rules based on the results of optimization for each forecast horizon based on the starting storage levels, runoff forecasts and current level of water demands. The lower rule is not an issue here due to sufficient amount of water to meet the current irrigation demands using the Standard Operating Policy (although this situation may change in future as water demands continue to grow). Lower rule curves are pertinent for managing droughts, and system optimization over the entire operation season is a valuable tool for constructing both the operating zones that correspond to various levels of water rationing policies (Ilich 2023), but this is not applicable to the case study presented here. Hence, the proposed modeling framework mimics the operators' decision that takes into account the existing priorities, starting storage, runoff forecasts and demands, and moves from one day to another by sequentially updating the decision made for a single day based on the available runoff forecast horizon. The WEB.BM model could replicate the same decision-making process for time steps that are shorter than a day, going down to 1 h, but the selection of shorter time steps has a bearing on the length of the downstream river reaches to enable meaningful channel routing calculation. For example, if 1 h is selected as a desired calculation time step, the length of each river reach should be defined such that the average channel flow requires 2 h to traverse the entire river reach, which would typically limit the length of river reach no more than 8 km, thus resulting in creation of numerous river reaches for systems that are more than 400 km in length.

## IMPLEMENTATION PROCEDURE

The process of implementing hydrological channel routing into an LP optimization framework requires repetitive calls of the optimization solver until a desired level of convergence is achieved, or until the previously set maximum number of iterations has been reached. Each call of the optimization solver requires that a set of routing coefficients be precalculated based on the most recent channel flows. The routing coefficients are updated be recalculating travel times using the channel flows obtained from the most recent call of the solver, before initiating the next iterative call of the solver. Typically, the maximum number of 5–6 iterations is sufficient to achieve the accuracy of the assumed and recalculated travel times of less than 1%. User have the option to set the tolerance limit (TL) on accuracy and the maximum number of iterations. If the desired accuracy is not achieved within the maximum number of iterations, the model proceeds with the solution from the last iteration and prints out a warning message. This process is presented by using a pseudo code presented below, where an iteration refers to a single call of the optimization solver routine.

For each time horizon solution {

set the tolerance limit TL to a desired threshold; set the maximum number of iterations Imax; set the initial travel times Ts using initial solution flows for a given time of the year;

for all iterations i < Imax {call the optimization solver to get optimal channel flows and reservoir levels;

calculate travel times Ts’ for all channels based on the optimal channel flows;

if[abs(Ts-Ts’) < TL)] // check the compliance with the tolerance limit TL

return to main program with the final solution for a time horizon;

else

set Ts = Ts’ // update the assumed travel times along each channel;

Ci = F(Ts) // recalculate coefficients Ci based on updated the travel times Ts;

i = i + 1; // update the iterations counter;

}

adopt the best optimization solution achieved after the Imax iterations;

}

It should be noted that optimization solver provides flows through the hydro power plants, reservoir spills and diversion canals for the entire system. In essence, these are all reservoir releases (except for the large diversions at Durgapur Barrage), but they are made through different conduits and for different purposes, and as such they are modeled as independent variables.

## WEB.BM MODEL

The WEB.BM model is a web-based application that uses COIN-OR library of public domain solvers Ilich (2021). The model has three solution modes:

- (a)
Solutions for a single time step forecast (STO-0);

- (b)
Simultaneous solutions for Multiple Time Step Solutions for n time step forecast (STO-n); and

- (c)
Simultaneous multiple time step solutions for the entire year (MTO year by year) or for all years.

While options (a) and (b) can be executed with or without the use of hydrological channel routing, option (c) is usually conducted using weekly or longer time steps on a steady-state basis, without the use of hydrological channel routing.

While the solution involves three simultaneous time steps and assumed foreknowledge of flows over a time horizon of 3 consecutive days (Figure 2), the final adopted flows for the first solution are only applicable to day 1, shaded in black (Figure 3). The ending storage for each reservoir on day 1 is used as the starting storage for day 2, and the solution proceeds over the next 3 days. Once this solution is finalized, reservoir operation for day 2 is adopted as final, its storage at the end of day 2 is assumed to be the starting storage for day 3, and the solution proceeds of the next three-day horizon. In this way, the model mimics the decisions of an operator who would take into account a short-term runoff forecast and operational objectives at the start of each day. In addition to the STO-n solution model, for time steps that are multiples of 1 day, the model can be run using two more modes: MTO for year individual year (i.e., all time steps within one year are solved simultaneously) and MTO for all years, which assumes foreknowledge of runoff conditions through the entire period for which input data are available. This can be useful when studying carry over storage requirements from year to year. These two options are typically used for planning studies based on steady-state solutions with weekly or 10-daily time steps. The default model output is shown in the units of flow for channels and net evaporation on reservoirs, elevation (m) for reservoirs, while hydro power output provides both flow through the hydro power plant and the mean power in MW achieved over a simulated time step. The output can optionally be generated in the units of volume, which converts the storage, net evaporation and all channel flows to volumes, thus making it easier to conduct mass balance checks for reservoir nodes. The Google Maps layer shown in Figure 12 is created by the model user with the drawing tools that are internally available in the WEB.BM, which also provides charting tools for both time series and exceedance plots, as well as the statistical analyses of deficits related to the targets for both in-stream and consumptive use targets.

## DAMODAR RIVER BASIN: CASE STUDY

Located in in West Bengal, the Damodar River Basin is one of the oldest basins in India in terms of its history of early water resources developments that date back to the 1950s. There are currently five multi-purpose reservoirs in the Damodar River Basin, built to provide water supply for industry and irrigation, with three installed hydro power plants, but one of main functions assigned to the reservoir operators is protection. Panchet and Maithon are the largest dams in Damodar basin, with total respective installed hydropower plant capacities of 80 and 60 MW. Four out of five reservoirs are currently operated by the Damodar Valley Corporation, and one by the State of Jharkhand. In this paper, the assumption was that all five reservoirs are operated as a system so as to best manage the historic flood events based on the 2-day forecast horizon and the original reservoir design specifications, although some of the reservoirs are currently operated below the originally intended range due to unresolved resettlement and land acquisition issues.

Given a short-term runoff forecast of only 2 days, the optimization program can quickly search through all the available options that also include pre-flood drawdowns. Effective use of the flood storage zones combined with the effects of channel routing and local runoff downstream of Durgapur Barrage can be significant. The principal advantage of LP is that it guarantees finding the global optimum, which is not the case with the heuristic search engines or non-LP solvers.

^{3}/s (high cost); and above 3,500 m

^{3}/s (excessive cost). Historically, the full bank capacity used to be around 3,700 m

^{3}/s in the past has been reduced to 2,850 m

^{3}/s as a result of silt deposits in the river channel below Durgapur Barrage that have built up over the years, revising the full bank capacity flows (

*Q*

_{fb}) to 3,000 m

^{3}/s as shown in Figure 5. It should be noted that the daily maximum flows have exceeded 5,000 m

^{3}/s ten times in the 1961–2017 period. Some of those violations are shown in Figure 6. Modelling results show that if the runoff forecast was available for a 2-day time horizon, daily storage releases derived by the model would have kept all flows downstream of Durgapur Barrage below 3,000 m

^{3}/s for all but two years out of the total of 57 simulated years. Figure 6 shows the results for most of the historic years that have been affected by floods with a comparison of historic flows (dark dashed line) at Durgapur Barrage, and simulated flows at Durgapur Barrage (blue line) and Jamalpur (red line).

^{3}/s at Durgapur, and its routed component at Jamalpur did not exceed 2,850 m

^{3}/s.

^{3}/s at Durgapur and within 3,200 m

^{3}/s at Jamalpur. The effects of including the hydrologic channel routing into the optimization process is also demonstrated for two sequential Damodar River reaches above the Panchet Dam and shown in Figure 8.

## MODELLING OF HYDRO POWER, WATER SUPPLY AND DISCUSSIONS OF OPTIMALITY CRITERIA

Based on the established priorities in the Damodar Valley Corporation which is in charge of managing the river basin, hydro power is to be generated as a byproduct of irrigation releases. Had there been no irrigation and no floods, the obvious optimal solution would be to maintain the FSL at all times, in an effort to route all storage inflows through the turbines while maintaining the design operating head. Only one of the following conditions should cause deviations from the FSL:

- (a)
Model derived drawdown to manage the incoming flows that exceed the hydro power plant capacities and route them all through the turbines by minimizing spills bypassing turbines;

- (b)
A drawdown caused by irrigation demands being higher than storage inflows; and

- (c)
A pre-flood drawdown followed by the use of flood storage zone in an effort to maintain the flows in designated downstream river reaches below the full bank capacity.

An obvious optimal model solution is the one that removes flood damage from all historical floods and meets all irrigation water requirements, while maintaining the storage levels as high as possible while minimizing unnecessary spills so as to maximize generated hydro power subject to the above operational constraints. Optimality of model solutions is often hard to establish if arbitrarily selected user-defined utility functions are optimized. When comparing outputs from different models, the relevant optimality criteria should include the number and magnitude of violations of the stated operating goals. Table 1 shows the key reservoir levels used in this study.

. | Tilaiya . | Konar . | Tenughat . | Maithon . | Panchet . |
---|---|---|---|---|---|

Max WL | 374.50 | 428.00 | 267.00 | 152.00 | 135.00 |

FSL | 371.00 | 425.81 | 263.66 | 149.00 | 130.00 |

Min WL | 363.32 | 410.57 | 249.00 | 132.59 | 119.48 |

. | Tilaiya . | Konar . | Tenughat . | Maithon . | Panchet . |
---|---|---|---|---|---|

Max WL | 374.50 | 428.00 | 267.00 | 152.00 | 135.00 |

FSL | 371.00 | 425.81 | 263.66 | 149.00 | 130.00 |

Min WL | 363.32 | 410.57 | 249.00 | 132.59 | 119.48 |

The diagram in Figure 13 should be related to the net head which is defined as the difference between the head and tail water elevations, rather than the reservoir water level, but in many instances where tail water elevation is modeled using a constant average for median flows (as is the case in this study), the net head elevation translates to the storage levels of the upstream reservoir. It is obvious that low storage levels cannot produce maximum power, with the typical reduction being in the order of 20–30% of the installed power capacity. The WEB.BM can model this constraint easily using piece-wise linearization of the maximum flow function which can be derived from the net head-flow-efficiency diagrams that are provided by turbine manufacturers. It should be noted that numerous recent publications that detail the use of various heuristic solvers for modeling hydro power typically ignore this constraint.

Flood damage occurs when the flow exceeds 3,000 m^{3}/s at Durgapur Barrage. It is encountered only in the optimization run with a 2-day lead time, and it happens in the following years 1978 (3,092.8 m^{3}/s); 2000 (3,263.5 m^{3}/s); 2006 (3,494.5 m^{3}/s) and 2017 (3,888.7 m^{3}/s). The objective function has a progressive penalty of 9,000 for each unit of flow between 3,000 and 3,200 m^{3}/s, 25,000 for each unit of flow between 3,200 and 3,500 m^{3}/s, and 150,000 for each unit of flow above 3,500 m^{3}/s. Table 2 summarizes the results of four scenarios that differ only by the length of the forecast horizon which was set to 2, 4, 6, and 10 days.

Scenario . | Flood damage ($) . | Irrigation deficits (%) . | Mean annual generated energy from hydropower (GWh) . |
---|---|---|---|

2-day lead time | 89.3 * 10^{6} | 0 | 288.03 |

4-day lead time | 0 | 0 | 297.54 |

6-day lead time | 0 | 0 | 304.28 |

10-day lead time | 0 | 0 | 316.81 |

≥20-day lead time | 0 | 0 | 329.67 |

Scenario . | Flood damage ($) . | Irrigation deficits (%) . | Mean annual generated energy from hydropower (GWh) . |
---|---|---|---|

2-day lead time | 89.3 * 10^{6} | 0 | 288.03 |

4-day lead time | 0 | 0 | 297.54 |

6-day lead time | 0 | 0 | 304.28 |

10-day lead time | 0 | 0 | 316.81 |

≥20-day lead time | 0 | 0 | 329.67 |

It is also evident from Table 2 that scenarios with 4, 6, and 10 days forecast horizon provide solutions with no flood damage and no deficits to irrigation, in spite of the total annual requirement of 3.2 billion m^{3}. The amount of generated power increases with the increase in the forecasting horizon, until the maximum is found with the required forecasting horizon of 20 days. Increasing the length of the forecasting horizon beyond 20 days does not increase the generated power any more, and therefore represents the best possible solution, as it confirms optimality of all scenarios subject to the stated constraints and the length of the forecasting horizon. These results demonstrate that the combined use of optimization and reliable short-term forecasts can effectively replace the current practice of using the reservoir rule curves during the high flow seasons.

## CONCLUSIONS AND RECOMMENDATIONS

This paper demonstrates that managing reservoirs using optimization can deliver outstanding results provided that accurate runoff forecasts are available. These results indicate that future efforts should focus on improving the runoff forecasting models that should rely on better weather forecasts, remote sensing and possibly the use of artificial intelligence for improved transformation of rainfall to runoff. The synergy of improved forecasts and mathematical optimization for reservoir operation holds out a promise of delivering significant improvements that would save lives and reduce property damage resulting from severe floods. It should also be noted that the WEB.BM can be used both as an operational or a planning model, with a range of time step lengths, and a combination of different multiple objectives that could simultaneously involve maximization of hydro power or minimization of water supply deficits.

Multi-purpose multi-reservoir optimization presented in this paper is easy to understand and verify. It would be useful to compare other similar applications in river basins where all relevant input data could be posted in public repositories, thus enabling cross comparison of results that were obtained with different tools while the input data, constraints and stated objectives remain the same. Obtaining equally good or better results using the EMODPS or other similar modeling approaches compared to multi-purpose multi-reservoir optimization models can help gain traction and interest among practitioners for alternative technologies. Since most modeling approaches claim to provide optimal guidance to reservoir operators for a given set of state variables and short-term runoff forecasts, it would make sense to compare their results. The need for establishing rigorous benchmarks in the water resources modelling sector has been identified in the past (Ilich & Manglik 2022) as a useful guideline to water resources practitioners, especially for water resources systems with moderate complexity.

## ACKNOWLEDGEMENTS

This project was part of a consultancy agreement funded jointly by the Government of India, Ministry of Jal Shakti and the World Bank under the National Hydrology Project office in New Delhi. The WEB.BM was created by Optimal Solutions Ltd from Calgary, Canada.

## DATA AND MODEL AVAILABILITY

The WEB.BM model is accessible free of charge at www.riverbasinmanagement.com requiring only an email address and a password for registration and access, resulting in the creation of a personal workspace within the application for each registered user. Signing on to the platform provides instant access to instructional videos and the User's Manual. Data used in this study belong to the Ganga River Basin, and as such they are restricted from public sharing. However, the WEB.BM User's Manual contains sample data from other basins that can be used to test the capabilities presented in this paper.

## DATA AVAILABILITY STATEMENT

Data cannot be made publicly available; readers should contact the corresponding author for details.

## CONFLICT OF INTEREST

The authors declare there is no conflict.

## REFERENCES

**142**(2), 04015050

**6**(2). EnPress Publisher, https://systems.enpress-publisher.com/index.php/NRCR/article/view/2256.