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.

  • 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.

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.

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.

Most river basin models that incorporate channel routing usually use the Muskingum channel routing procedure due to the ease of its implementation, since it involves fitting of the three routing coefficients:
formula
(1)
where the current channel outflow Ot is related to the current inflow into the channel It and the previous time step inflows and outflows It−1 and Ot−1. The principal issue with the Muskingum approach is that the routing coefficients Ci 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:
formula
(2)
By subtracting both sides of the above equation with Ot−1, multiplying by t/(Ot-Ot−1) and by letting ΔS/(Ot-Ot−1)=TS, the above equation becomes:
formula
(3)
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 Ot from the reach at the end of the current time step t, implying that Ts corresponds to the time that would be required to ‘empty’ the current storage in the river reach assuming constant outflow Ot. Expression (3) can also be converted to the following explicit form:
formula
(4)

The above routing model is non-linear, so it typically requires a few numerical iterations where Ot is recalculated and the travel Ts is updated using the interpolated from the table of Ts 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 Ci 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.

The effects of routing are particularly important for larger systems with multiple reservoirs. Fixed rate storage releases are transformed before they reach the end of each river reach, as shown in Figure 1 where three downstream river reaches (A, B, and C) are assumed, and they may also be altered by the effects of additional tributaries, diversions or other changes due to river regulation. The ultimate goal is to find reservoir releases that will not exceed the full bank capacity at the critically affected downstream areas (highlighted in orange in Figure 1), by taking into account the hydrological transformation of flow through the river reaches and downstream influences that may affect the flows. All these effects should be properly included in the modeling process.
Figure 1

Routing of fixed storage release through downstream reaches.

Figure 1

Routing of fixed storage release through downstream reaches.

Close modal

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.

As originally proposed by Windsor (1973) and later confirmed by Needham 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).
Figure 2

Multiple time step solution network for three consecutive time steps.

Figure 2

Multiple time step solution network for three consecutive time steps.

Close modal
Typically, flood condition is indicated by flows that exceed the full bank capacity. Flood damage function can then be created for several benchmark flood levels above the full bank capacity, representing increasing cost of flooding as a function of the magnitude of flood. This would internally cause the model to split a river reach into parallel channels to enable piece-wise linearization, with each having a unique pricing vector applicable per unit of flow, as originally proposed by Windsor (1973). If the optimization problem is posed as maximization of benefits, the pricing vector Pi 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:
formula
(5)
Subject to:
formula
(6)
formula
(7)

Equation (6) denotes the mass balance at network node n, with Ai 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 Ui,t 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 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 Pi is easy with early LP-based applications (especially those based on NFA solvers), since the only requirement is to select the Pi 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 Pi 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:

  • (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 Pi are provided in the ‘case study’ section. At this point, it should be noted that Pi 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.

The issue of automated search for the best values of the pricing vector Pi 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.

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.

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.

Option (b) with the hydrological channel routing is particularly useful for modeling floods, where the parameter n represents the length of the forecast time horizon given as the number of simulated time steps that are solved simultaneously. For example, STO-0 implies a solution that assumed only one time step forecast, without any knowledge of inflows and outflows beyond the length of the time step. This is equivalent to simulation used by most NFA-based models such as the MODSIM or REALM. The time step length can be any multiple of days (or a fraction of days for individual events). The model can utilize daily, weekly, monthly or custom time step lengths, where the custom lengths are defined by the user. For multi-year simulations, users should define custom time step lengths such that they sum up to 365 days for regular years and to 366 days for leap years. Any simulation with time steps that are shorter than a day can also be specified by using the custom time step length and specifying the length of time step as a fraction of a day. For example, the time step length of 0.25 would correspond to a 6-h time step since it is a quarter of 1 day. No matter what the time step length, the user can always select to run the model using the STO-n option. Figure 3 demonstrates the use of the STO-2 option.
Figure 3

Example of STO-2 solution mode.

Figure 3

Example of STO-2 solution mode.

Close modal

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.

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.

Input data are based on daily natural flow series that were recalculated using the historic records of storage levels and outflows, as well as the recorded flows at Durgapur Barrage. The model was run by using a daily time step with 2 days of runoff forecast. Although two daily time steps were solved simultaneously, only the solution for the first day is taken as final before proceeding to the next 2-day time horizon. This ‘rolling’ solution allows evaluation of all 11 severe floods that have happened historically from 1961 to 2017. The operating policy is therefore modified during floods by the results of optimization that are based on the runoff foresight of 2 days. This makes it possible for the model to decide how to arrange storage releases and balancing of all five reservoirs such that the negative impacts of floods are minimized, while simultaneously taking into account channel routing and local runoff downstream of Durgapur Barrage. Figure 4 shows the system modeling schematic.
Figure 4

WEB.BM modelling schematic of the Damodar River Basin.

Figure 4

WEB.BM modelling schematic of the Damodar River Basin.

Close modal

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.

The downstream channel at Durgapur Barrage is subdivided into four flow zones: below 3,000 (no cost); between 3,000 and 3,200 (moderate cost); between 3,200 and 3,500 m3/s (high cost); and above 3,500 m3/s (excessive cost). Historically, the full bank capacity used to be around 3,700 m3/s in the past has been reduced to 2,850 m3/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 (Qfb) to 3,000 m3/s as shown in Figure 5. It should be noted that the daily maximum flows have exceeded 5,000 m3/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 m3/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).
Figure 5

Flood damage cost functions at Durgapur and Jamalpur.

Figure 5

Flood damage cost functions at Durgapur and Jamalpur.

Close modal
Figure 6

Simulated channel flows at Durgapur Barrage and below Jamalpur.

Figure 6

Simulated channel flows at Durgapur Barrage and below Jamalpur.

Close modal
The responses of the reservoirs have resulted in the flows at Durgapur and Jamalpur that are shown in Figure 7 below. Note that the optimization solver adjusted reservoir releases in each day such that the sum of the local runoff downstream of both reservoirs with reservoir outflows did not exceed 3,000 m3/s at Durgapur, and its routed component at Jamalpur did not exceed 2,850 m3/s.
Figure 7

Simulated flows at Durgapur and Jamalpur for 2009 flood.

Figure 7

Simulated flows at Durgapur and Jamalpur for 2009 flood.

Close modal
One can also observe the effects of routing, where the red line represents the downstream flows at Jamalpur, which are associated with reduced peaks and time lag compared to the flows at Durgapur (blue line). The area underneath each line represents the total flow volume, and these volumes should be the same for inflow and outflow hydrograph over each hydrological event. Nine out of eleven severe floods that have happened in the 1961–2017 period were handled in a manner similar to the solution for the 2009 flood. Only two of them (2006 and 2017) caused overbank spills within the model, which were still kept within 3,500 m3/s at Durgapur and within 3,200 m3/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.
Figure 8

Routing effects on the Damodar River above the Panchet Dam.

Figure 8

Routing effects on the Damodar River above the Panchet Dam.

Close modal
The flood in 2006 was selected to demonstrate the effect of the length of the available time horizon, which were tested for 1-, 2-, 3- and 4-day lead times by running a daily simulation for the 2006/2007 year and by comparing the solutions in terms of the downstream flows at Durgapur and Jamalpur along with the reservoir levels. Earlier start of pre-flood drawdowns can be observed for solutions with longer time leads, which is to be expected, resulting in reduction of the peak flow at the downstream control channels, as depicted in Figures 9, 10, 11 and 12.
Figure 9

Flows at Durgapur for various forecast (lead time) lengths.

Figure 9

Flows at Durgapur for various forecast (lead time) lengths.

Close modal
Figure 10

Flows at Jamalpur for various forecast (lead time) lengths.

Figure 10

Flows at Jamalpur for various forecast (lead time) lengths.

Close modal
Figure 11

Maithon reservoir levels for various forecast (lead time) lengths.

Figure 11

Maithon reservoir levels for various forecast (lead time) lengths.

Close modal
Figure 12

Panchet reservoir levels for various forecast (lead time) lengths.

Figure 12

Panchet reservoir levels for various forecast (lead time) lengths.

Close modal

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.

Table 1

Key reservoir elevations (m)

TilaiyaKonarTenughatMaithonPanchet
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 
TilaiyaKonarTenughatMaithonPanchet
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 

Maximum flow through the turbines is usually governed by the amount of available storage, as shown in Figure 13. Once the designed power is achieved (typically at the FSL or in its vicinity), any increase in storage levels that may happen during floods will not generate more power. Rather, it will require reduction of the maximum flow since the installed power capacity should not be exceeded, while power is a product of net head and flow, so the higher net head (reservoir elevation) should result in reduced flow.
Figure 13

Typical relationship between the maximum turbine flow and reservoir water level.

Figure 13

Typical relationship between the maximum turbine flow and reservoir water level.

Close modal

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 m3/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 m3/s); 2000 (3,263.5 m3/s); 2006 (3,494.5 m3/s) and 2017 (3,888.7 m3/s). The objective function has a progressive penalty of 9,000 for each unit of flow between 3,000 and 3,200 m3/s, 25,000 for each unit of flow between 3,200 and 3,500 m3/s, and 150,000 for each unit of flow above 3,500 m3/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.

Table 2

Summary of the tested scenarios

ScenarioFlood damage ($)Irrigation deficits (%)Mean annual generated energy from hydropower (GWh)
2-day lead time 89.3 * 106 288.03 
4-day lead time 297.54 
6-day lead time 304.28 
10-day lead time 316.81 
≥20-day lead time 329.67 
ScenarioFlood damage ($)Irrigation deficits (%)Mean annual generated energy from hydropower (GWh)
2-day lead time 89.3 * 106 288.03 
4-day lead time 297.54 
6-day lead time 304.28 
10-day lead time 316.81 
≥20-day lead time 329.67 

Runoff forecasts of 10 or 20 days ahead do not have a practical value, but they are shown here for the purpose of demonstrating the effects of optimization. Having a perfect forecast of incoming high flows for 2 or 3 weeks in advance would initiate the reservoir drawdown earlier using lower outflow rates that can fit the turbine capacity, thus minimizing spills bypassing turbines both during the pre-flood drawdown period as well as during the flood routing period. It has always been known that perfect forecast was the key to deriving perfect solution, but the development of various scenarios which are only differentiated by the length of the assumed perfect forecast can help establish the ultimate length of the forecast that provides the best possible solution, as shown in Figure 14 below, where the maximum generated energy from Panchet and Maithon hydro power plants is expressed as the length of the forecast horizon (the hydro power plant at Tilaiya reservoir was ignored due to its capacity of only 2 MW).
Figure 14

Mean annual generated energy as a function of the length of the forecast horizon.

Figure 14

Mean annual generated energy as a function of the length of the forecast horizon.

Close modal
Historically, reservoirs have been operated based on the rule curves that were developed in the 1970s using spreadsheet calculation and simulation models. Figure 15 shows compares the historic levels of Panchet reservoir which has the largest hydro power plant with installed capacity of 80 MW with the simulated storage levels for the scenario with 4 days forecast horizon. It is evident that the rule curves used historically often resulted in lowering the storage levels much more than necessary.
Figure 15

Comparison of historical and simulated Panchet reservoir levels.

Figure 15

Comparison of historical and simulated Panchet reservoir levels.

Close modal

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 m3. 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.

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.

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.

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 cannot be made publicly available; readers should contact the corresponding author for details.

The authors declare there is no conflict.

Andreu
J.
,
Capilla
J.
&
Sanchís
E.
1996
Aquatool, a generalized decision support system for water-resources planning and operational management
.
Journal of Hydrology
177
(
3
),
269
291
.
doi:10.1016/0022-1694(95)02963-X
.
CALVIN project
2023
Available from: https://calvin.ucdavis.edu/calvin-project-overview (accessed 2 October 2023)
.
Cunge
J. A.
1969
On the subject of a flood propagation method (Muskingum method)
.
Journal of Hydraulic Research
7
(
2
),
205
230
.
DHI 2023 Mike Hydro Basin User Guide
, https://manuals.mikepoweredbydhi.help/2021/Water_Resources/MIKEHydro_Basin_UserGuide.pdf
last (accessed November 2023).
Dobson
B.
,
Wagener
T.
&
Pianosi
F.
2019
An argument-driven classification and comparison of reservoir operation optimization methods
.
Advances in Water Resources
128
,
74
86
.
Elsevier Ltd. doi:10.1016/j.advwatres.2019.04.012
.
EWater Academy
2022
EWater Model. Available from: https://ewater.org.au/ewater-academy/ (accessed October 2022)
.
Giuliani, M., Castelletti, A., Pianosi, F., Mason, E. & Reed, P.
2016
Curses, tradeoffs, and scalable management: Advancing evolutionary multiobjective direct policy search to improve water reservoir operations
.
Journal of Water Resources Planning and Management
142 (2), 04015050
.
Ilich
N.
2023
Use of mathematical optimization to construct optimal reservoir operating rules—A case study of the Barna reservoir in Narmada Basin, India
.
Natural Resources Conservation and Research
Ilich
N.
&
Manglik
N. K.
2022
The battle of river basin models: The Narmada River Basin challenge
.
Hydrological Sciences Journal
. .
Israel
N.
&
Lund
J.
1999
Priority preserving unit penalties in network flow modelling
.
Journal of Water Resources Planning and Management
125
(
4
),
205
214
.
Labadie
J. W.
2004
Optimal operation of multireservoir systems: A state of-the-art review
.
Journal of Water Resources Planning and Management
130
(
2
),
93
111
.
Labadie
J. W.
,
Bode
D. A.
&
Pineda
A. M.
1986
Network model for decision-support in municipal raw water supply
.
Water Resources Bulletin
22
(
6
),
927
940
.
Needham
J. T.
,
Watkins Jr
D. W.
&
Lund
J.
2000
Linear programming for flood control in the Iowa and Des Moines rivers
.
Journal of Water Resources Planning and Management
126
(
3
),
118
127
.
Perumal
M.
&
Sahoo
B.
2008
Volume conservation controversy of the variable parameter Muskingum–Cunge method
.
Journal of Hydraulic Engineering
134
(
4
),
475
485
.
Randall
D.
,
Cleland
L.
,
Kuehne
C. S.
,
Link
G. W.
&
Sheer
D. P.
1997
Water supply planning simulation model using mixed-integer linear programming ‘engine’
.
Journal of Water Resources Planning and Management
123
(
2
),
116
124
.
Rani
D.
&
Moreira
M.
2010
Simulation-optimization modelling: a survey and potential application in reservoir systems operation
.
Water Resources Management
24
,
1107
1138
.
doi:10.1007/s11269-009-9488-0.
RTC-Tools
2022
DeltaRes. Available from: https://oss.deltares.nl/web/rtc-tools/features (accessed October 2022)
.
Tahiri
A.
,
Che
D.
,
Ladeveze
D.
,
Chiron
P.
&
Archimède
B.
2022
Network flow and flood routing model for water resources optimization
.
Scientific Reports
12
,
3937
.
https://doi.org/10.1038/s41598-022-06075-0
.
US Army Corps of Engineers
2021
Program Description and User Manual for SSARR – Streamflow Synthesis and Reservoir Regulation
.
US Corps of Engineers
1998
HEC-5 Simulation of Flood Control and Conservation Systems. User's Manual, Version 8.0
.
Victoria State Government
2022
Resources Allocation Model (REALM)
. .
Williams
J.
1969
Flood routing with variable travel time or variable storage coefficients
.
Transactions of the American Society of Agricultural Engineers
12
(
1
),
100
103
.
doi:10.13031/2013.38772
.
Windsor
J. S.
1973
Optimization model for operation of flood control systems
.
Water Resources Research
9
(
5
),
1219
1226
.
Yates
D.
,
Seiber
J.
,
Purkey
D.
&
Huber-Lee
A.
2005
WEAP21 – a demand-, priority-, and preference-driven water planning model: Part 1
.
Water International
30
(
4
),
501
512
.
doi:10.1080/02508060508691894
.
Yeh
W. W.-G.
1985
Reservoir management and operations models: A state-of-the art review
.
Water Resources Research
21
(
12
),
1797
1818
.
Zagona
E. A.
,
Fulp
T. J.
,
Shane
R.
,
Magee
T.
&
Goranflo
H. M.
2001
Riverware: A generalized tool for complex reservoir system modeling
.
The Journal of the American Water Resources Association
37
(
4
),
913
929
.
This is an Open Access article distributed under the terms of the Creative Commons Attribution Licence (CC BY 4.0), which permits copying, adaptation and redistribution, provided the original work is properly cited (http://creativecommons.org/licenses/by/4.0/).