Representation and improved parameterization of reservoir operation in hydrological and land-surface models

Reservoirs significantly affect flow regimes in watershed systems by changing the magnitude and timing of streamflows. Failure to represent these effects limits the performance of hydrological and land-surface models (H-LSMs) in the many highly regulated basins across the globe and limits the applicability of such models to investigate the futures of watershed systems through scenario analysis (e.g., scenarios of climate, land use, or reservoir regulation changes). An adequate representation of reservoirs and their operation in an H-LSM is therefore essential for a realistic representation of the downstream flow regime. In this paper, we present a general parametric reservoir operation model based on piecewise-linear relationships between reservoir storage, inflow, and release to approximate actual reservoir operations. For the identification of the model parameters, we propose two strategies: (a) a “generalized” parameterization that requires a relatively limited amount of data and (b) direct calibration via multi-objective optimization when more data on historical storage and release are available. We use data from 37 reservoir case studies located in several regions across the globe for developing and testing the model. We further build this reservoir operation model into the MESH (Modélisation Environmentale-Surface et Hydrologie) modeling system, which is a large-scale H-LSM. Our results across the case studies show that the proposed reservoir model with both parameter-identification strategies leads to improved simulation accuracy compared with the other widely used approaches for reservoir operation simulation. We further show the significance of enabling MESH with this reservoir model and discuss the interdependent effects of the simulation accuracy of natural processes and that of reservoir operations on the overall model performance. The reservoir operation model is generic and can be integrated into any H-LSM.

1 Introduction

1.1 Background and motivation

Human interventions in natural hydrologic systems, through damming and storing water, diversion, surface and groundwater abstraction, irrigation, and land use change, have significantly altered the natural river flow regimes and the terrestrial water cycle of many river basins (Vörösmarty et al., 1997, 2003; Oki and Kanae, 2006; Wisser et al., 2010; Haddeland et al., 2014; Biemans et al., 2011). These interventions are to fulfill different types of demands such as domestic, industrial, irrigation, and hydropower demands and to meet other needs, such as flood control and conservation of aquatic habitats. With a total storage volume of more than 8000 km³ (ICOLD, 2003; Vörösmarty et al., 2003; Hanasaki et al., 2006), more than 50 000 dams have been constructed globally to regulate more than half of the world's large river systems (Nilsson et al., 2005). The aggregate storage volume of these dams is greater than 20 % of the global mean annual runoff (Vörösmarty et al., 1997) and is 3 times the annual average water storage in world's river channels (Hanasaki et al., 2006).

Despite the benefits in terms of enhancing water availability in support of food security, power supply, etc., dams result in several negative environmental and social consequences. Adverse environmental effects include changes in natural river dynamics in terms of water temperature, sediment and nutrient transport, etc., and the fragmentation and loss of biodiversity (Vörösmarty et al., 2010). Reservoirs can also intensify evaporation by increasing the surface area of water exposed to direct sunlight and air and through water supply for irrigation (de Rosnay, 2003; Pokhrel et al., 2012). Other environmental impacts of dams include the alteration of landscape due to dam construction and changes to land–atmosphere interaction that can have a profound impact on local and regional climate (Hossain et al., 2012; Degu et al., 2011). Adverse social effects include the displacement of people living near the dam site, changes to fishing patterns, and downstream erosion (Strobl and Strobl, 2011, p. 449). There are research gaps remaining in evaluating both positive and negative social impacts of dams (Kirchherr et al., 2016). Such gaps have been the subject of many studies in both academia and industry for years and, recently, have led to the formalization of the study area of “socio-hydrology” (Sivapalan et al., 2012; Sivakumar, 2012).

Dams and reservoirs change the natural flow regimes in rivers both in terms of magnitude and timing of flows. As a result, for rivers that contain large or small dams and reservoirs, flow regimes are a combination of natural and managed flows. Various modeling communities manage this mix of natural and managed flows differently. Archfield et al. (2015) compare three families of models that can be used at continental scales: catchment models (CMs), global water security models (GWSMs), and land-surface models (LSMs). CMs generally ignore water management and focus on unmanaged headwater catchments. GWSMs have been utilized in global-scale streamflow simulations and generally focus on large-scale water management issues, which are hindered by a lack of data on large-scale water management and operational decisions. LSMs have traditionally focused on providing lower boundary conditions for atmospheric models but are increasingly being used for hydrological applications in which they are referred to as hydrologic land-surface models (H-LSMs). LSMs generally ignore water management (Clark et al., 2015; Davison et al., 2016), with a few exceptions (e.g., Voisin et al., 2013a, b). A fourth family of water models, which is relevant to the work presented here, are water management models (WMMs; Labadie, 1995; Yates et al., 2005). Water modelers who know how the water is managed within their basins of interest generally use WMMs (Lund and Guzman, 1999; Labadie, 2004; Kasprzyk et al., 2013). These models contain very detailed representations of water management decisions but often consider natural flow processes in a much more rudimentary fashion than CMs.

Modeling the many managed basins around the world using the current generation of CMs or LSMs can result in models with limited fidelity, questions the credibility of their predictions of future water resources in basins with dams and reservoirs. Therefore, there is a pressing need for better characterization and integration of the operation of dams and reservoirs into hydrological modeling frameworks using CMs and LSMs (Nazemi and Wheater, 2015a, b; Pokhrel et al., 2016; Wada et al., 2017). This need motivated the objectives of this study, described in Sect. 1.2, and some previous research, outlined in Sect. 2. The integration of reservoir regulation into hydrological modeling frameworks will improve our ability to simulate highly regulated basins around the globe, leading to better understanding of historical conditions of water resource systems and improved assessment and prediction of their future vulnerability to climate and environmental change.

1.2 Objectives

Building upon previous research, this study aims to do the following:

• Develop and test an improved reservoir operation model that can be integrated into any CM and LSM at any scale but in particular at large scales. Of interest is a simple but effective parameterization that can be adjusted to varying levels of data availability.

• Integrate the developed reservoir operation model into an LSM and evaluate its performance when working in combination of other processes in the model. Also of interest is assessing the potential conceptual and technical issues in this integration.

Another potentially very fruitful but largely unexplored approach would be to couple CMs and LSMs with WMMs, but that approach is not examined here due to the fact that WMMs generally require extensive information on how water is managed within a basin, whereas we are particularly interested in the more generic case when this information is likely to be limited or unavailable.

The organization of the remainder of the paper is as follows. Section 2 reviews different existing approaches in the literature for the representation of reservoir operations in hydrologic models. Section 3 presents the proposed reservoir operation model and the metrics used to evaluate it in comparison with other existing models. Section 4 provides a description of the reservoir dataset used for the developments and testing. Section 5 presents the assessment results and comparisons. Section 6 ends the paper with a summary of the main findings and conclusions.

2 Existing reservoir models in catchment models and land-surface models

An adequate representation of human interventions in Earth system models is a major challenge. Systematic approaches towards full integration are needed, as outlined in the recent studies of Nazemi and Wheater (2015a, b), Wada et al. (2017), and Pokhrel et al. (2016). In this work, our focus is on the representation of dam and reservoir operations in CMs and LSMs, particularly when used at large scales. While there has been tremendous progress in the last decades in modeling the operation and management of reservoir systems at local to regional scales (e.g., Castelletti et al., 2010; Chang et al., 2010; Fraternali et al., 2012; Razavi et al., 2012; Asadzadeh et al., 2014; Guo et al., 2013), a gap still exists between the methodologies applied for local- and regional-scale reservoir operations and management and the representation of reservoir operation in Earth system models, particularly in LSMs. This gap is due to a twofold challenge. First, the upscaling of methodologies used at smaller scales to larger scales is non-trivial; second, the availability of data on reservoir operation and water use is often limited in many parts of the world. For example, the reservoir purpose and operational details are not always known, and large reservoirs typically serve several purposes (Wisser et al., 2010). As a result, most current hydrological modeling activities with CMs and LSMs, if not all, offer only a limited capability in simulating reservoir operations, whereas reservoir operation in practice involves a complex set of human-driven processes and decisions.

The existing reservoir operation methods in hydrologic models can be categorized roughly into three groups based on their level of complexity in representing flow regulation: (1) natural lake methods, (2) inflow- and demand-based methods, (3) artificial neural-network (NN) techniques, and (4) target storage-and-release-based methods.

2.1 Natural lake methods

The most primitive methods use formulations developed for the simulation of natural lakes or uncontrolled reservoirs. In these methods, the downstream release is calculated as a function of reservoir storage characterized by some empirical parameters (Meigh et al., 1999; Döll et al., 2003; Pietroniro et al., 2007; Rost et al., 2008). For instance, Meigh et al. (1999) calculate the release by $Q_t=S_t^{\mathrm{1.5}}$, where Q_t and S_t are release and reservoir storage, respectively. Their method was later modified by Döll et al. (2003) such that $Q_t=b_{\mathrm{1}}(S_t-S_{\min }){\left(\frac{S_t-S_{\min }}{S_{\max }-S_{\min }}\right)}^{b_{\mathrm{2}}}$, where b₁ and b₂ are release coefficients, and S_min and S_max are minimum and maximum allowable reservoir storages. The advantage of this method, as shown in Döll et al. (2003), is its minimal data requirement, which supports its global applicability to model lakes, reservoirs, and wetlands. However, it has limited functionality in adequately representing managed reservoirs due to not accounting for reservoir operation policies to constrain or increase releases at different phases of reservoir storage dynamics. Such simplistic methods ignore the fact that the operation of a reservoir depends on the reservoir purpose and the seasonal pattern of the mismatch between the demands it supports and the inflow it receives.

2.2 Inflow-and-demand-based methods

The inflow- and demand-based methods include reservoir water balance models that determine reservoir release using a function that accounts for inflow or a combination of inflow and demands. The simplified method in this group is the method used in Wisser et al. (2010); it estimates the release as a function of mean annual inflow and a set of empirical parameters that can be calibrated in the absence of information on the actual operation of a reservoir.

Hanasaki et al. (2006) pioneered the development of inflow and demand reservoir models and laid the foundation for many subsequent developments. The method of Hanasaki et al. (2006) simulates reservoir release at a monthly time step within a global routing model and accounts for water withdrawals for reservoirs categorized as irrigation reservoirs. They grouped reservoirs serving all other purposes as non-irrigation reservoirs. This approach first estimates a provisional total annual release at the beginning of the water year based on the long-term mean annual inflow adjusted by an annual release coefficient. Then, a monthly provisional release is estimated based on the purpose of the reservoir (irrigation or non-irrigation). Downstream demands are accounted for in irrigation reservoirs only. The provisional monthly release for large reservoirs is then modified by the annual release coefficient to calculate the actual monthly release, and the provisional monthly release for small reservoirs is additionally adjusted based on the monthly inflow to calculate the actual monthly release. The release coefficient is estimated as a function of the reservoir storage at the beginning of the operational year and the reservoir capacity (the formulation of Hanasaki et al., 2006, is briefly explained in Sect. 3.4). The release coefficient reduces the current year release if the storage at the beginning is low and vice versa. Thus, the release coefficient accounts for inter-annual variability and facilitates the representation of strategies to overcome reservoir depletion in dry years and flood overtopping in wet years.

The implementation of the release coefficient is one of the limitations of Hanasaki et al. (2006) because it depends only on the year's initial storage and does not account for the actual inflow of the current operational year; i.e. it does not use foresight. The initial storage reflects the recent past of operation of the reservoir, while the actual inflow could be considerably different than the long-term mean annual inflow. For instance a sequence of low-flow years would result in a low initial storage, while the current year inflow (which is not known yet) could be high, and vice versa. Additionally the simplification of complex reservoir operation in Hanasaki et al. (2006) by using the mean annual inflow and a release-constraining coefficient produces errors. However, the method is generic and has low data requirements, which are advantageous. The results showed that the reservoir algorithm improved monthly discharge simulation compared to the natural lake method (Hanasaki et al., 2006). The approach is effective and has found wide applicability in several global hydrological and land-surface models.

The original Hanasaki et al. (2006) reservoir model has been modified in subsequent studies to address some of its limitations. For example, it has been modified for water extraction and other reservoir functions, such as fulfilling environmental flows (Hanasaki et al., 2008a, b; Pokhrel et al., 2012), and has been adjusted to address direct precipitation over, and evaporation from, the reservoir (Döll et al., 2009).

Biemans et al. (2011) added new functionalities to the Hanasaki et al. (2006) reservoir model related to irrigation water demand and supply distribution and ran it at a daily time step. Their contributions include (1) modifying irrigation withdrawals to account for conveyance losses and irrigation efficiency, (2) adjusting the minimum release to 10 % of the mean monthly inflow, (3) prioritizing irrigation over flood control, (4) using regulated flow instead of natural flow to estimate mean annual inflow, and (5) storing the “flow to be released” for 5 d in the reservoir – to mimic the storage within the conveyance system – before it is released to the river. Voisin et al. (2013a) further modified the reservoir model of Hanasaki et al. (2006) to include multipurpose functionalities (irrigation and flood control) by changing the operation to release more before the onset of snowmelt-flood season so that there will be enough room to store flood waters from snowmelt in the reservoir. The modification requires the specification of a flood control period. Voisin et al. (2013a) have also evaluated the uncertainty of reservoir simulation by comparing withdrawal vs. consumptive demands and natural vs. regulated flow for configuring operating rules. The results of Voisin et al. (2013a) demonstrated that adding flood control in reservoir operation, along with a parameterization using mean annual natural inflow and mean monthly withdrawals, improves the reservoir storage and flow simulation.

Haddeland et al. (2006) developed another pioneering generic reservoir model that has been implemented in a routing model at a daily time step to study the impact of reservoir and irrigation water withdrawals on continental surface water fluxes. The model is retrospective; i.e., it assumes full knowledge of the upcoming operation-year reservoir inflow. The reservoir operation is conducted using an optimization scheme to determine the optimal release to satisfy different sectoral demands and targets that are defined in the form of objective functions. In the case of a multipurpose reservoir, the model gives priority to irrigation demand, followed by flood control and hydropower production. Minimum flow is estimated using natural flow based on 7 d consecutive low flows with a 10-year recurrence period. The flood protection objective function is minimizing reservoir release above the bankfull discharge, which is estimated using the long-term mean of annual maximum discharge. Irrigation is optimized to satisfy downstream irrigation demand, while hydropower is optimized to increase power production. Predicting inflows for the current operational year, if possible, would allow the method to optimize the release while accounting for the whole operational year; otherwise optimizing day-to-day release without accounting for the remaining operational year would require several constraints. The maximum daily release is set based on the reservoir water balance that sets the storage at the end of the operational year to vary from 60 % to 80 % of the maximum capacity.

Similar to Hanasaki et al. (2006), the model of Haddeland et al. (2006) is favorable due to its generic formulation and capability to operate multipurpose reservoirs and to extract water for irrigation from the reservoir. These make the model applicable for large-scale hydrologic models, when data on operational policies are limited (Adam et al., 2007; van Beek et al., 2011). One limitation of Haddeland et al. (2006) is that it requires knowledge of the future inflow for each reservoir so that the optimization can be conducted to determine the optimal release. Another limitation is that the release can deviate from the actual value because of simplifications of the objective function and errors from irrigation demand calculation. The algorithm does not represent reservoirs with multi-year operational policies (Adam et al., 2007) and also requires running the model many times to optimize the reservoir release.

Adam et al. (2007) modified the Haddeland et al. (2006) reservoir model parameterization to include (1) estimated minimum flow based on observed mean winter flow, (2) the reservoir-filling phase, (3) a storage–area–depth relationship following the regular shape approximation of Liebe et al. (2005), and (4) a seasonally varying hydropower-production economic value that can be calibrated for hydropower production instead of a constant one; van Beek et al. (2011) further modified the retrospective inflow assumption to the prospective model by approximating the upcoming operational year inflow based on previous years' inflow (requires historical inflow observation) and then adjusting the release and demand every month using the actual inflow as estimated from a hydrologic model.

Solander et al. (2016) tested and compared six generic equations to represent reservoir release and storage simulations. The complexity of equations tested varies from the simplest case that assumes that reservoir outflow equals inflow (no-reservoir assumption) to a more complex representation using separate linear functions during reservoir-filling and release periods. While the reservoir-filling and release seasons were identified using long-term mean temperature, their respective release equations are configured as a function of reservoir inflow, storage, and optimized seasonal empirical parameters. Their results on California reservoirs showed that the equation dependent on inflow is best for the recharge season, while release during the drawdown season was better represented as a function of storage. Despite failing for highly regulated reservoirs, their study demonstrated the possibility of generalizing the seasonal empirical parameters as a function of the ratio between winter inflows and storage capacity. However, further testing is required to examine the usefulness of the Solander et al. (2016) method in different regions, such as cold regions, with different filling and release seasonality.

Although the inflow- and demand-based models provide improved results compared to the natural lake approach, these models do not accurately reproduce observed flows (Adam et al., 2007; Haddeland et al., 2006; Coerver et al., 2018). Overall, while the above methods have better flexibility for coupling with global hydrological and land-surface models, the methods have limitations in accounting for details of reservoir operation. For an adequate representation of reservoirs, particularly multipurpose reservoirs and/or those with multi-year carry-over capacity, it is important to consider reservoir zoning and adjust reservoir-release formulations for different storage levels. The absence of this consideration may limit the capability of this group of methods in representing complex reservoir operations.

2.3 Neural-network-based methods

Artificial NN models have been applied to establish data-driven rules that relate reservoir storage, inflow, and release data. This type of model includes (1) extensive data on reservoir release, storage and inflow, but minimal prior expert knowledge of the reservoir operation, and (2) extensive training of a model for each individual reservoir to deduce the reservoir operation rules. Neural-network techniques have been widely used beyond reservoir operation applications (e.g., flood forecasting, streamflow simulation, and water quality; Maier and Dandy, 2000; Razavi and Karamouz, 2007) and more recently have shown promise in reproducing historical reservoir operations (Coerver et al., 2018).

The study of Coerver et al. (2018) provides detailed background on NN applications for deduction of reservoir operation rules and also demonstrates the performance of NN-based fuzzy rules to describe the reservoir-release decisions. The analysis of Coerver et al. (2018) involves different levels of input complexity for the neural-network setups, such as the importance of accounting for inflow prediction and time of the season on the reservoir operation performance. Another similar application was shown by Ehsani et al. (2016), who demonstrated a general reservoir operation scheme that uses an NN technique to map the general input–output relationships to actual operating rules of 17 dams. Ehsani et al. (2016) demonstrated the possibility of aggregating multiple reservoirs that are closely located so that their integrated effect can be accounted for in large-scale hydrological modeling studies. In a subsequent study, Ehsani et al. (2017) integrated the reservoir model of Ehsani et al. (2016) into a global water security model to study reservoir operations under climate change.

While these studies demonstrated that the NN-based models can reproduce historical reservoir operation data and possibly outperform the widely used reservoir simulation models such as those of Hanasaki et al. (2006) and Wisser et al. (2010), the user of such models may have to deal with a fundamental limitation, i.e., their “black-box” nature. This limits their ability to provide insight into the underlying mechanisms of reservoir operation and might mask possible shortcomings in a derived NN model. Further, the credibility of their performance in extrapolation beyond the historical data can be in question, as they ignore the expert knowledge available on the actual physical and socio-economic processes that govern reservoir operations. Together, these limit the interpretation of results and their applicability in a changing environment. There have been some recent research efforts to reformulate neural networks such that they can overcome these limitations (e.g., see Razavi and Tolson, 2011).

2.4 Target storage-and-release-based methods

The target storage-and-release-based methods aim to emulate actual rule curves (i.e., reservoir target storage and release for different times of the year) that guide reservoir operators to decide on downstream releases (Burek et al., 2013; Yates et al., 2005; Neitsch et al., 2005). The target levels of storage divide the total reservoir storage capacity into multiple zones. For example, in the SWAT model (Arnold et al., 1998), a reservoir model is available in which the total storage of a reservoir is divided into sediment, principal, flood control, and emergency flood control zones, where each zone is either specified by the user or as a function of soil moisture wetness (Neitsch et al., 2005). Wu and Chen (2012) modified this approach by changing the reservoir zoning model and developed a reservoir-release simulation strategy that uses a decision-based parameterization to better fit both storage and release of multipurpose reservoirs. However, they reported only one application of this strategy to a local-scale reservoir, and its comprehensive evaluation needs to be performed on other reservoirs in other regions with different climates, levels of regulation, and allocation objectives.

Zhao et al. (2016) integrated a reservoir regulation module into a hydrology model, requiring user-specified (based on observed data) values to divide the reservoir into inactive, conservation, and flood control zones. In their module, the release from the conservation zone is determined using water demand, which includes multi-sectorial demand and environmental flow. The release from the flood storage zone is decided as a function of inflow (classified as flood inflow or non-flood inflow), downstream channel current discharge, and downstream maximum discharge. At the time of flood, if the downstream discharge is below the maximum limit, release from flood storage zone is estimated using available storage above the conservation zone, multiplied by a weight parameter which allows the release of more water. If the downstream discharge is at maximum capacity, there is no release from flood storage zone. Finally, any storage above the flood storage zone is automatically released. Additionally, Zhao et al. (2016) added the possibility to operate reservoirs conjunctively by giving release priorities to immediate downstream demands and by limiting the release if the downstream reservoir is within flood storage zone. The results of reservoir integration showed improved capability of the hydrological model to simulate storage and release for Lake Whitney and Aquilla Lake in Texas. The limitation of Zhao et al. (2016) for wider application is that there is no generic formulation of reservoir zoning (requires user specification), and evaluation was only performed on two reservoirs.

Similarly, Burek et al. (2013) divided the total reservoir storage into conservative, normal, and flood zones within the LISFLOOD model and defined releases in accordance with these storage zones using multiple linear regression. Zajac et al. (2017) showed the applicability of this method in capturing the effects of lakes and reservoirs globally using a parameterization that depends on naturalized inflow and maximum storage. Their results showed that the inclusion of reservoirs and lakes in a hydrologic model through this method helped improve streamflow simulation for many stations, but the performance in replicating observed storage dynamics was not reported.

Overall, the primary advantage of methods in this category is that they allow approximation of reservoir-release policy and have the potential of making use of detailed data on a reservoir when available. Their main limitation, however, is their relatively high data demands. When data are available, methods under this category have the potential to enhance the representation of dams and reservoirs in terms of both reservoir storage and release while adapting to the seasonality and change in operations on different timescales from daily to seasonal. These methods seem advantageous compared to NN-based models, as their functioning is transparent, accounting for the governing processes, while requiring similar data.

Given the advances in the field and the growing availability of data sources, the target storage-and-release methods seem to be the most promising, as they can better simulate the reservoir operation dynamics (the dynamics of both storage and release). The data requirement includes data on observed inflow, observed release, observed storage (level), and reservoir physical characteristics. Reservoir-level data are available for most lakes and reservoirs in the public domain, particularly in North America. These data can be converted to reservoir storage using reservoir elevation–area–volume relationships or by using area–volume relationships approximated by regular geometric shapes (Yigzaw et al., 2018; Liebe et al., 2005; Lehner et al., 2011). Inflows to and releases from a reservoir can be approximated by streamflow stations located upstream and downstream of the reservoir, respectively. Further, satellite missions such as MODIS (Savtchenko et al., 2004) and satellite radar altimetry provide information on lake and reservoir surface area dynamics and reservoir water elevation for some large reservoirs. The combination of MODIS and satellite radar altimetry allows deriving storage–area–depth relationships (Gao et al., 2012; Andreadis et al., 2007; Zhang et al., 2014; Yoon and Beighley, 2015). The planned SWOT (2021) mission (Garambois and Monnier, 2015; Biancamaria et al., 2016) will increase the availability of water-level data for smaller rivers (with widths going down to 100 m) that can be potentially converted to discharge to estimate reservoir inflows and downstream reservoir releases.

3 Material and methods

This study aimed to develop an improved reservoir model that better emulates reservoir operation for large-scale hydrologic modeling application in terms of both reservoir storage and release following the previous advances in target-storage-and-release-based methods reviewed in Sect. 2.4. In this section, we present the characteristics and formulation of our reservoir model. The reservoir water balance is maintained using the continuity equation, as shown in finite difference form in Eq. (1). The aim is to estimate unknown storage S_t and release Q_t at the current time step based on the storage at the previous time step S_t−1 and precipitation (P) over the reservoir, evaporation (E) from the reservoir, and inflows (I) during the current time step. When integrated within an H-LSM model, the inflow will be the modeled value of the upstream catchment that would account for delays in the precipitation-runoff generation and routing. This equation is solved in conjunction with the parameterization equations presented in the next section for reservoir releases to compute S_t and Q_t:

$$\begin{array}{}\text{(1)}& \begin{array}{rl}{\displaystyle \frac{S_t-S_{t-\mathrm{1}}}{\mathrm{\Delta }t}}& ={\displaystyle \frac{I_t+I_{t-\mathrm{1}}}{\mathrm{2}}}-{\displaystyle \frac{Q_t+Q_{t-\mathrm{1}}}{\mathrm{2}}}+{\displaystyle \frac{P_t+P_{t-\mathrm{1}}}{\mathrm{2}}}\\ & -{\displaystyle \frac{E_t+E_{t-\mathrm{1}}}{\mathrm{2}}}.\end{array}\end{array}$$

Figure 1 The schematic representation of reservoir zoning and the storage-release function: (a) four (active) reservoir zones with inflows and outflows, and (b) piecewise-linear reservoir-release function. m₁ and m₂ control the slope of the release curve, and they change monthly. The upward blue arrow is to indicate that inflow to the reservoir may also be considered in determining the release in Zone 3.

3.1 Proposed reservoir operation model

A detailed description of our proposed target storage-and-release model (or target-release model for brevity) is provided here. This model is formulated in the form of parametric piecewise-linear functions that approximate the reservoir-release rules that may be used by reservoir operators. This model can be set up on any timescale; in the case studies reported here, we define the target levels to dynamically change over time. We call the model the dynamically zoned target release (DZTR) model. Piecewise-linear-function-based reservoir operation models have already been used to solve complex reservoir operations and water resource management problems (e.g., Razavi et al., 2014; Asadzadeh et al., 2014). A systematic integration of such models into large-scale hydrological modeling has been reported in Burek et al. (2013), as implemented in the LISFLOOD hydrological model, and in Neitsch et al. (2005), as implemented in the SWAT model. Our DZTR model is a generalization of the method developed by Razavi et al. (2014), which may also be viewed as a modification to the model proposed by Burek et al. (2013) in terms of parameterization and reservoir zoning. Figure 1 shows the schematic representation of DZTR; Fig. 1a shows the reservoir zoning, and Fig. 1b shows the piecewise-linear functions to estimate the release for each zone based on DZTR.

The DZTR model divides reservoir storage into five zones in a similar fashion to Wu and Chen (2012) and Burek et al. (2013), namely dead storage, critical storage, normal storage, flood storage, and emergency storage. Whenever storage is below the emergency storage zone, release only occurs through the bottom outlet, but when the storage is within that zone, release happens through both a bottom outlet and the spillway. In the absence of data, the dead storage (Zone 0) is assumed to be 10 % of the maximum storage after Döll et al. (2009). To estimate the remaining storage zones in cases where no operational information on a reservoir is available, we propose two alternative strategies: (1) setting the zones based on suggested exceedance probabilities on historical reservoir storage time series and (2) optimizing these zones to reproduce the observed storage and release time series. Target releases for each zone can be obtained in a similar fashion. These target storages and releases are allowed to vary each month (or on any other arbitrarily selected time resolution) to allow a better representation of the seasonality of reservoir operation.

When reservoir storage is within the dead storage zone (Zone 0), the reservoir release is zero (Eq. 2). In Zone 1 (critical storage zone), the reservoir release is a function of storage at a given time step and the critical release target value (Eq. 3). In this zone, the reservoir operates to avoid storage depletion while trying to support environmental flow requirements defined as a critical (or minimum) release. In Zone 2 (normal storage zone), the reservoir release is purely governed by reservoir storage and varies between critical and normal release targets (Eq. 4). In this zone, the downstream release is greater for higher levels of storage. In Zone 3, the release decision considers both reservoir storage and inflow in that time step as well as the normal and maximum release targets (Eq. 5). When in this zone, two scenarios may occur: (A) the amount of inflow in a time step is equal to or less than the normal release rate or (B) the amount of inflow in this time step is greater than the normal release rate. As formulated in Eq. (5), in the case of scenario B, the inflow rate comes into play to augment the release in an attempt to keep the reservoir level within the normal storage zone. Scenario B is expected to occur more frequently in smaller reservoirs that only have “within-year” storage capacity, while scenario A should be more commonly seen with larger reservoirs that have “multi-year” carry-over capacity. Hanasaki et al. (2006) suggested that reservoirs that have a ratio of storage capacity to mean annual inflow (referred to as c) of less than 0.5 be assumed as within-year reservoirs and that the ones with a ratio of 0.5 and above be considered to be multi-year reservoirs. Other values for this threshold were also suggested in the literature; e.g., Wu and Chen (2012) used a c value of 0.3. In this study, scenario A is used for reservoirs that have multi-year capacity (c>0.5), and scenario B is used for reservoirs that have within-year capacity (c<0.5). Lastly, in Zone 4 (emergency storage zone) the reservoir algorithm operates to avoid reservoir overtopping by releasing the larger of the maximum release target or all excess storage above the maximum storage value (the flood storage) constrained to the downstream channel capacity Q_mc. If not specified, a rough estimate of the downstream channel capacity value could be the 99th percentile of non-exceedance probabilities of discharges from historical data.

$$\begin{array}{}\text{(2)}& {\displaystyle }\mathrm{Zone}\mathrm{0}{Q}_t=\mathrm{0}\left[S_t<\mathrm{0.1}{S}_{\max }\right],\text{(3)}& {\displaystyle }\begin{array}{rl}\mathrm{Zone}\mathrm{1}{Q}_t=& \min \left(Q_{{\mathrm{c}}_i},{\displaystyle \frac{S_t-\mathrm{0.1}{S}_{\max }}{\mathrm{\Delta }t}}\right)\\ & \left[\mathrm{0.1}{S}_{\max }<S_t\le S_{{\mathrm{c}}_i}\right],\end{array}\text{(4)}& {\displaystyle }\begin{array}{rl}\mathrm{Zone}\mathrm{2}{Q}_t=& Q_{{\mathrm{c}}_i}+\left(Q_{{\mathrm{n}}_i}-Q_{{\mathrm{c}}_i}\right){\displaystyle \frac{\left(S_t-S_{{\mathrm{c}}_i}\right)}{\left(S_{{\mathrm{n}}_i}-S_{{\mathrm{c}}_i}\right)}}\\ & \left[S_{{\mathrm{c}}_i}<S_t\le S_{{\mathrm{n}}_i}\right],\end{array}\end{array}$$

$$\begin{array}{}\text{(5a)}& {\displaystyle }\begin{array}{rl}\mathrm{Zone}\mathrm{3}\mathrm{A}{Q}_t=& Q_{{\mathrm{n}}_i}+\left(Q_{{\mathrm{m}}_i}-Q_{{\mathrm{n}}_i}\right){\displaystyle \frac{\left(S_t-S_{{\mathrm{n}}_i}\right)}{\left(S_{{\mathrm{m}}_i}-S_{{\mathrm{n}}_i}\right)}}\\ & \left[S_{{\mathrm{n}}_i}<S_t\le S_{{\mathrm{m}}_i}\right],\end{array}\text{(5b)}& {\displaystyle }\begin{array}{rl}\mathrm{Zone}\mathrm{3}\mathrm{B}{Q}_t=& Q_{{\mathrm{n}}_i}+\max \left(\left(I_t-Q_{{\mathrm{n}}_i}\right),\left(Q_{{\mathrm{m}}_i}-Q_{{\mathrm{n}}_i}\right)\right)\\ & {\displaystyle \frac{\left(S_t-S_{{\mathrm{n}}_i}\right)}{\left(S_{{\mathrm{m}}_i}-S_{{\mathrm{n}}_i}\right)}}\left[S_{{\mathrm{n}}_i}<S_t\le S_{{\mathrm{m}}_i}\right],\end{array}\end{array}$$

$$\begin{array}{}\text{(6)}& {\displaystyle }\begin{array}{rl}\mathrm{Zone}\mathrm{4}{Q}_t=& \min \left(\left[\max \left({\displaystyle \frac{\left(S_t-S_{{\mathrm{m}}_i}\right)}{\mathrm{\Delta }t,Q_{{\mathrm{m}}_i}}}\right)\right],Q_{\mathrm{mc}}\right)\\ & \left[S_{{\mathrm{m}}_i}<S_t\right].\end{array}\end{array}$$

I_t, Q_t, and S_t are inflow, release, and storage at time step t. $S_{{\mathrm{c}}_i}$, $S_{{\mathrm{n}}_i}$, and $S_{{\mathrm{m}}_i}$ are critical, normal, and maximum storage targets for month i. $Q_{{\mathrm{c}}_i}$, $Q_{{\mathrm{n}}_i}$, and $Q_{{\mathrm{m}}_i}$ are critical, normal, and maximum release targets for month i. Q_mc is the maximum channel capacity parameter. The unit for inflow, release, and target release parameters are in cubic meters per second. The unit for storage and target storage parameters are in cubic meters.

3.2 Evaluation criteria

We evaluated the performance of the proposed reservoir operation model in emulating the outflow and storage data collected for many reservoirs around the world. As this model was intended to be integrated into large-scale H-LSMs, we further evaluated it when embedded in the MESH (Modélisation Environmentale-Surface et Hydrologie; Pietroniro et al., 2007) model. For all of these evaluations, we used Nash–Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970) and Kling–Gupta efficiency (KGE; Gupta et al., 2009) as the metrics to assess the goodness of fit of the model to observed reservoir outflow and storage data.

3.3 Identification of reservoir operation model parameters

As demonstrated in Sect. 3.1, the proposed reservoir operation model has six parameters ($S_{{\mathrm{c}}_i}$, $S_{{\mathrm{n}}_i}$, $S_{{\mathrm{m}}_i}$, $Q_{{\mathrm{c}}_i}$, $Q_{{\mathrm{n}}_i}$, and $Q_{{\mathrm{m}}_i}$) that can vary for different times of the year. We recommend varying these parameters on a monthly basis, while other time resolutions are also possible. To normalize the parameters and their ranges across different types and sizes of reservoirs, for every reservoir, we use cumulative distribution functions (CDFs) of historical storage and release values; see Fig. 2 for an example of CDFs of the Lake Diefenbaker reservoir (Gardiner Dam) in the Saskatchewan River basin, Canada. Our preliminary analysis indicated that target storage and release values corresponding to 10 %, 45 %, and 85 % non-exceedance probabilities generally perform reasonably well. We call these our generalized parameterization.

Figure 2 Cumulative distribution function (CDF). (a) Storage CDF of Gardiner Dam. (b) Reservoir-release CDF of Gardiner Dam. Double arrows on y axis show parameterizations ranges for each generalized parameter.

However, optimal values of parameters for a given reservoir can be identified, when data are available, through optimization and parameter-identification techniques (Maier et al., 2019; Guillaume et al., 2019). For this purpose, we used a bi-objective optimization approach, as follows, that begins with the generalized parameter values as the starting point and optimizes the model fit to both storage and release data simultaneously:

$$\begin{array}{}\text{(7)}& F\left(\mathit{x}\right)=\left(f_{\mathrm{1}}\left(\mathit{x}\right),f_{\mathrm{2}}\left(\mathit{x}\right)\right),\end{array}$$

where x is a vector of decision variables (parameter values), Ω is decision space, f₁(x) is NSE (flow) measuring the goodness of fit in reproducing observed release, and f₂(x) is NSE (storage) measuring the goodness of fit in reproducing observed storage dynamics.

For parameter identification on a monthly basis, a total of 72 decision variables were used in the optimization. We chose rather arbitrarily the storage and release target intervals that correspond to 5 %–35 %, 35 %–75 %, and 75 %–95 % non-exceedance probabilities as the ranges of variation for critical, normal, and maximum (flood) storage and release, respectively.

The bi-objective optimization problem to calibrate 72 reservoir target release and storage parameters was conducted using the AMALGAM evolutionary multi-objective optimization algorithm (Vrugt and Robinson, 2007). AMALGAM was selected because it provides effective and reliable solutions for multi-objective optimization using multiple search operators (genetic algorithm, particle swarm optimization, adaptive metropolis search, and differential evolution) and self-adaptive offspring creation. Vrugt et al. (2009), Wöhling and Vrugt (2011), Zhang et al. (2011), Raad et al. (2009), Dane et al. (2011), and others showed that the performance of AMALGAM model parameter calibration was better than, or equivalent to, some other calibration algorithms across different complex response surfaces. AMALGAM was run using an initial population size of 100, resulting in a total of 15 000 model evaluations to estimate final Pareto solutions for every single reservoir.

3.4 Comparison of reservoir operation models

We compared the performance of our DZTR model against the performances of Hanasaki et al. (2006) and Wisser et al. (2010) using NSE and KGE performance metrics defined on both storage and release simulations. The comparisons were made only for selected non-irrigation reservoirs because their irrigation reservoir formulation requires additional data on water demands. For the method of Wisser et al. (2010), reservoir release was estimated under two conditions as shown in Eq. (8):

$$\begin{array}{}\text{(8)}& Q_t=\left\{\begin{array}{lc}\mathit{\kappa }{I}_t& I_t\ge I_{\mathrm{m}},\\ \mathit{\lambda }{I}_t+\left(I_{\mathrm{m}}-I_t\right)& I_t<I_{\mathrm{m}},\end{array}\right.\end{array}$$

where κ and λ are empirical constants set to 0.16 and 0.6, respectively, I_m is the mean annual inflow (m³ s⁻¹), and I_t is inflow to the reservoir (m³ s⁻¹) at time t.

In the method of Hanasaki et al. (2006), the release from non-irrigation reservoirs was estimated by multiplying the mean annual inflow by release-constraining coefficients (Eq. 9). The release-constraining coefficients for every given operational year were estimated by dividing the initial storage of that year by the maximum storage (Eq. 10). The start of the operational year was considered to be the month when the mean monthly inflow shifts from being greater to being lower than the mean annual inflow:

$$\begin{array}{}\text{(9)}& {\displaystyle }& {\displaystyle }{r}_{\mathrm{m},y}=\left\{\begin{array}{ll}{k}_{\mathrm{rls},y}\cdot r_{\mathrm{m},y}^{\prime }& (c\ge \mathrm{0.5}),\\ {\left(\frac{c}{\mathrm{0.5}}\right)}^{\mathrm{2}}\cdot k_{\mathrm{rls},y}\cdot r_{\mathrm{m},y}^{\prime }+\left(\mathrm{1}-{\left(\frac{c}{\mathrm{0.5}}\right)}^{\mathrm{2}}\right)\cdot i_{\mathrm{m},y}& (\mathrm{0}\le c<\mathrm{0.5}),\end{array}\right.\text{(10)}& {\displaystyle }& {\displaystyle }{k}_{\mathrm{rls},y}={\displaystyle \frac{S_{\mathrm{first},y}}{\mathit{\alpha }\cdot S_{\max }}},\end{array}$$

where c is the ratio of maximum reservoir storage to the mean total annual inflow, k_rls,y is the release coefficient, $r_{\mathrm{m},y}^{\prime }$ is the provisional monthly release (m³ s⁻¹) which is equal to mean annual inflow (m³ s⁻¹), and α is a dimensionless constant set to 0.85. Equation (12) differentiates between multi-year and single-year storage reservoirs based on a threshold value of 0.5 for c.

3.5 MESH modeling system

MESH is Environment and Climate Change Canada's land-surface–hydrology modeling system (Pietroniro et al., 2007) and has been widely used in different parts of Canada (Davison et al., 2016; Haghnegahdar et al., 2017; Yassin et al., 2017; Sapriza-Azuri et al., 2018; Berry et al., 2017). MESH is a grid-based modeling system composed of three components: (1) the Canadian Land Surface Scheme (CLASS; Verseghy, 1991; Verseghy et al., 1993), (2) lateral movement of surface (overland) runoff and subsurface water (interflow) to the channel system within a grid cell, and (3) hydrological routing using WATROUTE from the WATFLOOD hydrological model (Kouwen et al., 1993).

Currently, the reservoir representation in the MESH model is rudimentary. MESH offers two approaches to account for reservoir operations. In the first approach, the observed reservoir release rate at the reservoir location is provided as input to the model. In this approach, the flow from the catchment upstream of the reservoir is discarded as the release is replaced by observations, a process referred to as streamflow insertion, which limits the utility of the model in simulating future scenarios for which releases are not yet known. This approach violates the water conservation law in the model and also creates discontinuities within the model setup, especially if there are reservoir cascades. Nevertheless, streamflow insertion could be used when coupling water management models with MESH, and these coupled models could be used to formulate scenarios for reservoir operations. As mentioned in the objectives, however, model coupling is not the focus of this study, as we are looking to examine the internal representation of reservoir operations within CMs and LSMs. The second approach is a natural lake or uncontrolled reservoir representation model similar to that of Döll et al. (2003), which was shown to be unsuitable for highly managed reservoirs. To improve the reservoir representation in MESH, this study aims to incorporate the DZTR model for controlled reservoirs into the MESH framework and evaluate its performance.

Table 1 Summary of reservoirs.

^* Main purpose: WS – water supply, HP – hydropower, IR – irrigation, and FC – flood control. ¹ Monthly data and simulation. ² Multiple reservoir models that are compared on this reservoir. NA – not available.

3.6 Case studies and data

The dataset required to build and evaluate a reservoir operation model includes (1) reservoir physical characteristics such as the volume–level–area relationship and maximum capacity, which are static (in the absence of sedimentation or dam heightening), (2) time series of hydrologic variables such as inflow, release, and water level (or storage), and (3) environmental flows. In this study, we assembled such a dataset for 37 reservoirs located in several regions across the globe (Fig. 3) to test the model. These dams represent a wide range of storage sizes, from 0.132×10⁹ to 162×10⁹ m³, spanning multiple orders of magnitudes. Most of these are located in the western US and western Canada, while some are located in Vietnam, central Asian countries, and Egypt. Table 1 provides a summary of reservoir locations, construction years, main purposes, data periods, and other dam characteristics. Measured inflow, release, and storage time series were collected from different sources. For reservoirs located in Canada, the data were acquired from Water Survey Canada, Alberta Environment and Parks, and the Saskatchewan Water Security Agency. Data for the High Aswan Dam were acquired from the Nile Basin Encyclopedia via the Nile Basin Initiative. The data for other reservoirs were provided by the authors of previous studies (Hanasaki et al., 2006; Coerver et al., 2018). Additional information about the degree of regulation, dam height, and catchment area were obtained from the GRanD database (Lehner et al., 2011). Reservoir operation simulations were performed on a daily and monthly basis, with simulation periods varying from 8 to 62 years. The choice of simulation period and timescale was based on data availability (Table 1). The first year of the reservoir simulations was used for spin-up, while the first half of the remaining data periods was used for calibration and the second half for model validation.

Figure 3 Locations of dams used to evaluate the reservoir routing model.

We also evaluated the integration of our reservoir model into the MESH model on seven reservoirs in two major basins in western Canada. Six of the test reservoirs (Gardiner, St. Mary, Waterton, Oldman, Ghost, and Dickson dams) are located within the heavily regulated Saskatchewan River basin (SaskRB), and one reservoir (W. A. C. Bennett Dam) is located in the Mackenzie River basin (MRB). For both of the basins, the MESH model was set up on a grid resolution of 0.125^∘, and the data required to build the MESH model were obtained from different sources. The topographic data are based on the Canadian Digital Elevation Data (CDED) at a scale of 1:250 000 and were obtained from the GeoBase website (http://www.geobase.ca/, last access: February 2018). The data on seven climate forcing variables were obtained from a Global Environmental Multiscale (GEM) numerical weather prediction (NWP) model (Côté et al., 1998) and the Canadian Precipitation Analysis (CaPA; Mahfouf et al., 2007). The land-cover data used are based on a 2005 land-cover map from the Canada Centre for Remote Sensing (CCRS). Soil texture data were obtained from Soil Landscapes of Canada (SLC) data of Agriculture and Agri-Food Canada. The MESH parameter values were taken from previous studies for calibration to streamflow at major subbasins of the SaskRB and MRB.

4 Results and discussion

4.1 Evaluation of the dynamically zoned target release (DZTR) model with generalized parameters

Individual reservoir simulations were conducted using the DZTR model with generalized monthly storage and release parameter values set at non-exceedance probabilities recommended in Sect. 3.3 for representing the reservoir storage zones and their respective target releases. The evaluation of the DZTR model was based on the performance metrics and a comparison with the other reservoir operation approaches and a base case where the existence of a reservoir was ignored in a model, referred to as the “no-reservoir assumption”. Under the no-reservoir assumption, the release was considered equal to inflow, storage was considered constant, and, as such, the performance metrics were computed by directly comparing inflow with observed release.

Figure 4 shows performance metrics results of the DZTR model in terms of NSE and KGE for storage and release simulations compared to those of the base case. As shown in Fig. 4a, both NSE (flow) and NSE (storage) results are greater than 0.25 and 0.5 for 90 % and 50 % of reservoirs, respectively. Although NSE (flow) results are greater than zero for all reservoirs, 1 % of reservoirs resulted in a negative NSE (storage) values. The no-reservoir assumption resulted in NSE (base case) values of greater than 0.25 and 0.5 for 45 % and 30 % of reservoirs, respectively, which, in general, are much lower than those of the DZTR model. Under the no-reservoir assumption, 48 % of the reservoirs resulted in a negative NSE (base case). Almost all positive NSE (base case) results were observed in reservoirs with c<0.5, such as Dickson, E. B. Campbell, Kayrakkum, Oldman, and Tyuyamuyun (as explained in Sect. 3, c is the ratio of storage capacity to annual inflow volume). However, for reservoirs with c>0.5, such as Bhumibol, Flaming Gorge, Fort Peck, High Aswan, and W. A. C. Bennett, the NSE (base case) is negative, which indicates the significant influence of their regulations on the hydrograph shape. Similarly, Fig. 4b shows the evaluation of the different reservoir models based on the KGE metric (Gupta et al., 2009). The values of KGE (flow) and KGE (storage) are greater than 0.25 and 0.5 for 100 % and 86 % of the reservoirs, respectively. The KGE (base case) values of 21 % of reservoirs are less than zero, while those of 57 % and 49 % of the reservoirs are greater than 0.25 and 0.5, respectively. The NSE and KGE results show that the DZTR with the generalized parameter values is capable of simulating flow and storage simulation well.

Figure 4 Performance evaluation result of the DZTR model reservoir operation algorithm: (a) NSE performance metrics and (b) KGE performance metrics.

Figure 5 shows scatter plots between KGE, NSE, and the regulation level represented by c. These plot orientations of the scatter plot between NSE and KGE on flow and storage show a strong positive correlation between the evaluation metrics, which indicates that both metrics provide somewhat similar evaluation information. Figure 5a and b show that both the no-reservoir assumption and DZTR estimate the release more accurately for lower levels of regulation. As expected, the degradation of performance was pronounced for the no-reservoir assumption as the regulation level increased, while DZTR performance reduced by a much smaller extent (still positive values). Almost all low-regulation-level reservoirs (c<0.5) showed positive performance metrics, which means that the reservoir regulation does not strongly modify the flow regime, whereas the opposite case is true for highly regulated reservoirs (c>0.5) in which the reservoir regulation strongly changes the reservoir release. Coerver et al. (2018) also noted that low-regulation-level reservoirs are more dependent on the current time-step inflow knowledge because of their smaller influence on the flow regime. The method of Hanasaki et al. (2006) also recognizes the strong dependence of c<0.5 reservoirs on inflow to determine the release by configuring the release as a function of monthly mean inflow. Conversely, the relationship between the regulation level and the storage simulation performance – in terms of both KGE (storage) and NSE (storage) – did not show a strong correlation (Fig. 5c).

Figure 5 Scatter plot between KGE and NSE, with regulation scale represented in terms of c. (a) KGE and NSE on no-reservoir condition, (b) KGE and NSE on DZTR release, and (c) KGE and NSE on DZTR storage.

Figure 6 compares the reservoir simulation and observation time series for the whole simulation period, while Fig. 7 shows the long-term average of these simulations. Inflows are also included in Figs. 6 and 7 to show the regulation pattern and changes caused by reservoir operation. Both figures indicate that the DZTR model captures both release and storage dynamics well, reproducing the daily and monthly seasonality as well as the magnitude and timing of storage and releases for almost all reservoirs, especially for reservoirs with high regulation (multipurpose and multi-year reservoirs) such as American Falls, Bhumibol, High Aswan, Sirikit, Trinity, and W. A. C. Bennett dams. However, the simulations also show some systematic over- and underestimations; for example, the simulations of Bhumibol, Fort Peck, High Aswan, Int. Falcon, Navajo, Bennett, and Int. Amistad reservoirs show continuous underestimation and overestimation of reservoir storage. Some reservoirs, such as Trinity, Palisades, Kayrakkum, Flaming Gorge, and Garrison, show underestimation and overestimation of reservoir storage only for some seasons. A closer look at American Falls, Flaming Gorge, Fort Peck, Glen Canyon, and Navajo dams in Fig. 6 indicates that the DZTR model reliably captured storage and release seasonality, inter-annual trends, and release pattern shifts during the consecutive wet years 1982–1986 followed by consecutive dry years 1987–1993. Similar patterns can be observed for the Gardiner Dam, with good simulation results during both dry years (1984–1986, 1988–1989, and 1999–2004) and wet years (1993, 2005, and 2010–2011). Furthermore, as expected, Fig. 7 shows that lowly regulated reservoirs (c<0.5) have less of an impact on the flow regime, but with fairly significant storage seasonality (Oldman, E. B. Campbell, Palisades, and Andijan). In general, the DZTR model with the generalized parameterization of reservoir zones and releases showed an improved performance and can be applied to any hydrological model (CM or H-LSM) that involves reservoir simulation.

Figure 6 Daily and monthly reservoir simulations using the DZTR model with a generalized parameterization; x axes show month and year, the primary y axes show release (m³ s⁻¹), and the secondary y axis shows storage (m³).

Figure 7 Long-term average daily or monthly reservoir simulations with generalized parameterization; the x axes show days (1–365) or months (1–12), the primary y axes show release (m³ s⁻¹), and the secondary y axis shows storage (m³).

It is important to note that for the case of a cascade of reservoirs, the parameterization of the DZTR model implicitly accounts, to some extent, for the upstream regulation effects of the upstream cascade reservoirs. This is because the regulated inflow is used for parametrizing downstream reservoirs, which reflects the regulation information of upstream reservoirs in the cascade. In reality, the operations of some cascade reservoirs are highly interlinked, particularly during the flood season. The decision regarding the release from one reservoir accounts for the (forecasted) state of other reservoirs. Such dual- or multi-linked operation is, however, not accurately accounted for in the presented algorithm because it assumes that each reservoir operates using its own storage state, inflow, and target storage and releases. Such systems require detailed modeling of operations that is not usually attainable in large-scale hydrological models. Depending on the purpose of the model, the modeler may decide to lump those reservoirs together to improve simulations downstream, e.g., Ehsani et al. (2016).

4.2 Comparison with previously developed reservoir operation models

To further illustrate the reliability of the DZTR model in simulating reservoir operation, a comparison with the methods of Hanasaki et al. (2006) and of Wisser et al. (2010) was conducted, as shown in Fig. 8. The comparison shows that the DZTR model provides a considerable improvement according to all of the performance criteria, notably NSE (storage) and NSE (flow), except in the case of the E. B. Campbell dam, where the method of Hanasaki et al. (2006) showed similar performance to DZTR. The method of Hanasaki et al. (2006) outperformed that of Wisser et a. (2010). Out of the 13 reservoirs compared, the DZTR resulted in positive values for both NSE (storage) and NSE (flow) for all simulations except for that of E. B. Campbell storage. The method of Hanasaki et al. (2006) and Wisser et al. (2010) resulted in eight and five reservoirs with positive NSE (flow), respectively (Fig. 8a), while both produced negative values for NSE (storage) for all the reservoirs compared (Fig. 8c). A similar performance pattern was observed for KGE metrics for flow and storage. In addition, we compared the DZTR result shown in Figs. 7 and 8 with the results reported in Coerver et al. (2018), who applied a fuzzy-neural-network model to extract 11 operating rules. This comparison showed that the performance of our generalized parameterization is comparable to that of Coerver et al. (2018) in simulating reservoir release; note that performance on storage is not reported in Coerver et al. (2018). This indicates that the simple parameterization applied in the DZTR model can provide a solution that is at least as effective as that of a neural-network-based model. Equally importantly, the DZTR model is transparent, as opposed to neural-network methods that are often criticized as being a black box.

Figure 8 A comparison of our proposed reservoir operation model with generalized parameters with the models of Hanasaki et al. (2006) and Wisser et al. (2010): (a) NSE (flow), (b) KGE (flow), (c) NSE (storage), and (d) KGE (storage).

The above comparisons were conducted for non-irrigation reservoirs because water demand data are needed to use the Hanasaki et al. (2006) method for irrigation reservoirs. In the case of the DZTR approach, the idea is that the DZTR model operates in such a way that it infers existing operational rules that cater to those demands. Thus, the release from DZTR accounts implicitly for downstream demands as per the intended purpose of the reservoir, whether it is for flood control, irrigation, hydropower, etc., or any combination of these. The case study dams include reservoirs with different purposes, as shown in Table 1. The DZTR approach showed good performance for these reservoirs.

If the reservoir purpose is irrigation, the target releases from DZTR are to satisfy irrigation demands because the parameterization is optimized based on observed releases. The release from an irrigation dam will be available for abstraction at the predefined abstraction points downstream of the dam. The abstraction and distribution can be implemented as separate modules, as done within the MESH land-surface model (Yassin et al., 2019). In such an implementation, MESH takes care of (1) calculation of actual irrigation demand for a configured irrigation area, (2) water abstraction from the defined abstraction point along the river below the dam, and (3) distribution across the irrigation fields. Regarding the return flow, the excess water flows from the irrigation areas are assumed to join the nearest stream within the model grid cell.

The DZTR model can in principle handle multipurpose reservoirs, e.g., a reservoir that is used simultaneously for hydropower generation, irrigation water supply, and flood control (e.g., High Aswan Dam in Egypt, which is one of the studied reservoirs); the DZTR provides the release based on the inflow and storage conditions that will be available for irrigation downstream. Hydropower does not consume water but returns it back to the river (except in rare cases where it returns to a different channel). Flood control is directly accounted for in the scheme and becomes relevant when storage is within the flood storage zone. Further, the flexible formulation of DZTR allows for implicitly changing the priorities in operation for selected time periods (e.g., months or seasons) by changing the target storage values during flood periods (e.g., the storage target before the onset of snowmelt). During these flood months, lowering the target storage would increase the buffer for flood control. Conversely increasing the target storage during other months would be desirable to store water and release during irrigation months. When the scheme is optimized using inflow, release, and storage data, the parameterizations capture these priorities implicitly, as expressed in the data. When inflow data are lacking, the generalized parameterization will set the storage zones based on the suggested exceedance probabilities (that were deduced based on all reservoirs used in the study), and the priorities can be assumed as predefined.

4.3 Initial storage and inflow sensitivity test

The initial storage at the beginning of the simulation is an input that needs to be specified for the model. The initial values can be prescribed from the observations, if available. However, the simulation of a hydrological and land-surface model could start at any point in time when there is no observation is available (e.g., some time in far past, a future scenario simulation, or a hypothetical scenario). Additionally, in a long-term simulation, the initial storage may result from a previous model simulation, which may not be as close to observations as desired. The aim of the experiment is to examine and show the extent to which the initial storage value affects the simulation performance.

To test the effect of initial storage used in the reservoir simulation performance, two experiments were conducted on three reservoirs with different scales of regulations: (1) Charvak (c=0.28), (2) Gardiner (c=1.46), and (3) High Aswan (c=2.84). In the first experiment, the initial storage was allowed to vary between 10 % of maximum storage (0.1⋅S_max) and maximum storage (S_max). In the second experiment, the initial storage range was narrowed to starting simulation-month minimum and maximum historical observations. In both tests, 150 simulations were conducted by sampling the initial storage using uniform random sampling from the defined storage range.

Figure 9 Reservoir initial storage effect on storage and release simulation: (a) Charvak storage case 1, (b) Charvak release case 1, (c) Charvak storage case 2, (d) Charvak release case 2, (e) Gardiner storage case 1, (f) Gardiner release case 1, (g) Gardiner storage case 2, (h) Gardiner release case 2, (i) High Aswan storage case 1, (j) High Aswan release case 1, (k) High Aswan storage case 2, and (l) High Aswan release case 2.

Table 2 Reservoir initial storage effect on storage and release simulation.

Figure 9 and Table 2 show the results of these initial storage perturbation experiments. For both experiments the simulations on the Charvak dam showed a similar range for NSE (flow) [0.79, 0.83] and NSE (storage) [0.61, 0.74]. Using 1 year as a spin-up period on Charvak dam simulations stabilized the initial storage effects, resulting in NSE (flow) of 0.82 and NSE (storage) of 0.74. The simulations on Gardiner Dam in the first experiment showed a range of [0.35, 0.51] for NSE (flow) and a range of [−0.43, 0.88] NSE (storage), while in the second experiment the ranges were narrowed to [0.44, 0.49] for NSE (flow) and [0.87, 0.88] for NSE (storage). For a 1-year spin-up period on the Gardiner Dam, this simulation reduced the NSE (flow) range to [0.49, 0.51] and the NSE (storage) range to [0.76, 0.87] in the first experiment and to 0.49 NSE (flow) and 0.87 NSE (storage) for the second experiments. On the other hand, the simulation on the High Aswan Dam showed a range of [−0.28, 0.85] for NSE (flow) and [0.38, 0.91] for NSE (storage) for the first experiment and [0.52, 0.85] for NSE (flow) and [0.42, 0.91] for NSE (storage) for the second experiment. Excluding a 1-year spin-up period from the metric calculation on the High Aswan Dam simulation narrowed the NSE (flow) range to [0.62, 0.85] and the NSE (storage) range to [0.58, 0.91] for both experiments. Overall, as expected, the experiments suggest that the effect of initial storage on reservoir simulation performance depends on the regulation scale. Starting from observed storage values and using a 1-year warm-up period allows stabilization of the initial storage effect for low and medium regulated reservoirs. However, for highly regulated reservoirs, as in the case of High Aswan, longer spin-up periods are needed to stabilize the simulations. For example, a 5-year spin-up period was required to fully stabilize the performance for the High Aswan Dam simulations.

The existence of inflow bias is inevitable in any hydrological modeling practice. To understand the behavior of the DZTR model under biased inflow conditions, we conducted a sensitivity experiment on the Charvak, Gardiner, and High Aswan reservoirs. To do so, the DZTR model performance was tested using five simulations in which the entire inflow time series was changed by −50 %, −25 %, 0 %, +25 %, and +50 %. The sensitivity of simulations to bias in inflow was evaluated using the NSE (flow) and NSE (storage) performance metrics.

Figure 10 Inflow bias sensitivity test on storage and release simulation: (a) Charvak storage, (b) Gardiner storage, (c) High Aswan storage, (d) Charvak release, (e) Gardiner release, and (f) High Aswan release.

Table 3 Inflow bias sensitivity test on storage and release simulation.

Figure 10 and Table 3 show the results of the inflow bias test and that the reservoir simulation performance significantly changes as a result of this bias. Reducing the inflow by 50 % considerably reduced the reservoir storage and release and led to negative values of NSE (flow) and NSE (storage) for all reservoirs. For such a large negative inflow bias, the reservoir operation tries to recover the storage to the target (observed) level by releasing as low as possible. Conversely, the positive inflow bias increased simulated storage and releases for all reservoirs, which led to negative performance metrics for all reservoirs except for Gardiner NSE (storage). As shown in Fig. 10, with large positive inflow bias, storage quickly moves towards flood and maximum storage targets, resulting in insufficient storage left to attenuate flood peaks, and the operation model starts discharging large releases through the spillway to maintain the storage at the maximum storage target. An inflow bias of −25 % and +25 % showed similar behavior to a −50 % and +50 % bias for all reservoirs, but the simulation performance metrics during −25 % and +25 % provide significant positive NSE values for the Charvak and Gardener dams, except for the Gardiner NSE (flow) for +25 %, which resulted in a negative NSE value. However, on the highly regulated High Aswan Dam, the ±25 % inflow bias significantly reduced the performance to negative values.

Figure 11 Reservoir-release-parameter multi-objective calibration result. x axes show NSE (flow) multiplied by −1, and the y axes show NSE (storage) multiplied by −1.

4.4 Parameter calibration and validation of the DZTR model

We tried to improve upon the generalized parameterization by calibrating the DZTR parameters via bi-objective optimization for two objective functions, Nash–Sutcliffe on reservoir storage – NSE (storage) – and Nash–Sutcliffe on reservoir release – NSE (flow). This is an important step when the data and computational resources for optimization are available to enhance reservoir simulation and consequently hydrological modeling of the region of interest. Figure 11 shows the multi-criteria reservoir calibration (yellow circles) and validation (red circles) Pareto solutions for all reservoirs. The Pareto solutions show strong trade-offs between fitting observed reservoir storage versus downstream release, which also reflects the fact that the problem is multi-objective by nature and that it is required to consider both storage and release instead of fitting one at the cost of degrading the other. The generalized parameterization solution for the calibration (yellow square with blue border) and validation periods (red square with blue border) is also added in Fig. 11 for each reservoir to show the improvement gained through parameter calibration. Relative to the generalized solution for the calibration period, reservoir parameter calibration improved both NSE (flow) and NSE (storage) for all reservoirs, with a median improvement of 0.11 and 0.21, respectively. The NSE (flow) improvement ranged from 0.017 to 0.575, and NSE (storage) improvement ranged from 0.02 to 0.66. The parameter calibration has shown significant improvement in reservoirs that have lower performance with generalized parameterization. The best examples of this case are Fort Randall, Int. Amistad, Trinity, Int. Falcon, and E. B. Campbell, as shown in Fig. 11. Small improvements in performance have also been observed in reservoirs that have greater performance with generalized parameterization, such as American Falls, Andijan, Nurek, High Aswan, Waterton, and Charvak. The validation of calibrated solutions improved the NSE (flow) and NSE (storage) for 56 % of the reservoirs, with a median improvement of 0.035 and 0.092, respectively. The NSE (flow) improvement in the validation period ranged from 0.001 to 0.335, and NSE (storage) improvement ranged from 0.004 to 1.02. During validation, the remaining reservoirs (44 % of them) resulted in NSE (flow) and NSE (storage) reductions, with a median reduction of 0.032 and 0.089, respectively. The reductions of NSE (flow) ranged from 0.001 to 0.073, and those of NSE (storage) ranged from 0.001 to 0.257.

Overall, considerable improvement was achieved for both calibration and validation periods for several reservoirs, such as the Dickson, Gardiner, Ghost, Int. Amistad, Int. Falcon, Kayrakkum, Sirikit, Yellowtail, and Glenmore. However, as shown in Fig. 11, the improvements of DZTR model performance during calibration do not guarantee performance improvement in validation. This is because, as well as for any other type of model, the properties of the calibration and validation periods might differ significantly. In particular, the calibrated Pareto solution does not show the same trade-off or level of performance during validation when there is considerable change in inflow properties as a result of consecutive wet or dry years. Examples of this condition are shown for Glen Canyon (similarly Bhumibol, Fort Randall, and Fort Peck), where the calibration period had more wet and high-inflow years than the validation period. Such considerable changes of inflow, storage, and release result in performance degradation during the validation period. In general, a small change in inflow, storage, or release for the validation period can change the shape of the trade-off. However, the calibrated parameters in most cases were still capable of producing good performance during validation that was close to, or better than, that of the generalized parameterization for the same period.

To further test the role of the calibration period, we calibrated all reservoirs using the whole observational record. The result of this test is shown in Fig. 12, which demonstrates the strong role of the calibration period. All reservoirs showed trade-off between storage and release fitting. The solution resulted in a consistent Pareto pattern similar to the split-sample calibration results. The median NSE (flow) and NSE (storage) improvement when using the whole observational record for calibration is approximately 0.1 and 0.12 respectively, while the maximum improvement reached 0.45 and 0.55 for some reservoirs. High improvements in storage and flow simulations in the case of whole-period calibration are mostly observed in reservoirs that have considerable shift of observed storage and flow across the period of observation period. Figure 13 shows some example reservoirs that had considerable improvements, such as Bhumibol, Canyon Ferry, Int. Amistad, Int. Falcon, Navajo, and Trinity dams, compared to generalized parameters (Fig. 6). Similarly, for the remaining reservoirs, calibrating the whole period showed (Fig. 13) better agreement of daily and monthly simulations with the observations, even for years with extreme deviations that are most likely associated with extreme dry and wet conditions. Additionally, the long-term average simulations (Fig. 14) showed that calibrating using the whole period reduced the deviation between simulations and observations, and in most cases the Pareto simulation range encompasses the observation. Overall, the calibration period test indicates the benefit of using long-term observation for parameterization (even for generalized parameterization) to allow the parameterization to represent behavior in extreme periods. Thus, we recommend using as many data as available to parameterize the model for a specific reservoir so that all information on reservoir operation will be accounted for.

Figure 12 Reservoir-release-parameter multi-objective calibration using all available data for each reservoirs; x axes show NSE (flow) multiplied by −1, and the y axes show NSE (storage) multiplied by −1.

Figure 13 Daily and monthly reservoir simulations using DZTR model with a generalized parameterization; x axes show the month and year, the primary y axes show release (m³ s⁻¹), and the secondary y axis shows storage (m³).

Figure 14 Long-term average daily or monthly reservoir simulations with generalized parameterization; the x axes show days (1–365) or months (1–12), the primary y axes show release (m³ s⁻¹), and the secondary y axes show storage (m³).

The DZTR scheme introduces more parameters to the host land-surface model. However, its parameters are external to those of land-surface model and are determined a priori using storage and release data. The decision of the timescale to use for specifying the parameters is left to the modeler. The user has the ability to investigate the seasonal patterns in the storage and release data and decide whether a monthly or a coarser timescale (e.g., quarterly) would be sufficient. In fact, the configuration of DZTR is also flexible in using any user-specified zoning that is available from observation, reservoir information, or zoning values specified in other studies, such as that of Zhao et al. (2016).

4.5 DZTR model test within the MESH model

Finally, the generalized parameterization of the DZTR model was integrated into the MESH model and tested to simulate six reservoirs in the Saskatchewan River basin (Gardiner, St. Mary, Waterton, Oldman, Ghost, and Dickson dams) and one reservoir (W. A. C. Bennett Dam) in the Mackenzie River basin, both in Western Canada. The reservoir simulation was run using MESH-modeled inflows at a half-hourly time step, the usual MESH time step, and the performance metrics were calculated at a daily time step. The MESH-modeled inflows are considered to represent the base-case scenario, and the inflow can be assumed to be regulated or natural, depending on whether there are dams upstream or not.

Figure 15 illustrates that the generalized DZTR model generally improves upon having no representation of the reservoirs in the model. This improvement is apparent in the NSE (flow) values, which increase with the DZTR model. The only exception is Dickson Dam, with a small reduction in NSE. The importance of integration of the DZTR model was predominant for the Gardiner and W. A. C. Bennett dams, which are highly regulated reservoirs (c>0.5) when compared to the other reservoirs tested in MESH.

Figure 15 Reservoir simulation results within MESH model run for selected reservoirs; x axis shows time (d), the primary y axes show release (m³ s⁻¹), and the secondary y axes show storage (m³).

This general improvement of flow simulation when comparing a reservoir model to the no-reservoir assumption is, of course, not surprising. What is important to note, however, is that the improvement in NSE can be substantial without calibration of the DZTR parameters. This is important for many LSM applications where calibration is generally not performed. Hanasaki et al. (2006) illustrated that their method is superior to the natural lake (or unregulated reservoir) method applied in many CMs and H-LSMs, and this paper shows that the DZTR model improves upon the results of Hanasaki et al. (2006). Therefore, it is natural to assume that the DZTR model would also be an improvement in uncalibrated H-LSM applications.

However, calibration is very common in CM or H-LSM applications in which the DZTR model would likely be employed. A full comparison of calibrated results between a no-reservoir case, natural lake (or unregulated reservoir), and the DZTR model (and the other reservoir models) is beyond the scope of this paper. Again, given the improvements shown with the uncalibrated DZTR model when compared with other uncalibrated models, and the general improvements shown here when calibrating the DZTR model, it is assumed that calibrating the DZTR model within a CM or H-LSM would improve upon calibrating an unregulated reservoir model or the other reservoir models compared in this paper.

The storage simulation showed a low NSE (storage) value for the St. Mary and Waterton dams and a negative NSE (storage) for Oldman and Ghost dams. However, the simulation showed a reasonable representation of storage variability but with considerable underestimation. This underestimation in storage in Fig. 15 is attributable to the fact that the modeled inflow is underestimated. It is expected that calibration of the land-surface parameters in conjunction with the DZTR parameters in MESH would improve the modeled inflows and resulting modeled reservoir storage.

It is worth mentioning again that H-LSMs, such as MESH, can also be used for the original purpose of LSMs, which is to represent fluxes from the land-surface to the atmosphere. If the approach improves modeled flows where reservoirs operate, it could result in a better parameterization of the LSM, which should in turn improve land-surface fluxes and feedbacks to the atmosphere.

4.6 Uncertainties in reservoir operations and DZTR parameterization

Reservoir operation involves considerable uncertainties that several factors are attributed to. One major source of uncertainty in reservoir operations is future inflows (long-term and short-term inflow forecast). The forecast contains errors are rooted in the forecast method, the driving climate forecast, snowpack measurements, timing of snowmelt, and the statistical (stationarity) assumptions to generate inflows based on historical inflows. The inflow forecast uncertainty is more significant during flood seasons because it involves subjective decisions of operators to avoid the risk of dam overtopping and downstream flooding. Other sources of uncertainty in reservoir operations include changes in demand over time because of increases in demand for irrigation, power, water supply, etc. The purpose of the reservoir can also change from its initial intended purpose (e.g., adding a hydropower station to an irrigation dam). These changes are only implicitly captured by the DZTR scheme, as implied in the storage and release time series used for parameterizing it for a specific reservoir.

Given the above uncertainties, even the actual reservoir operation may deviate from the designed reservoir operation rule curve. Some of the decisions of reservoir operators are spontaneous, ad hoc, and depend on experiences that are not usually documented. Thus, there are difficulties in accurately representing the historical operation or establishing accurate relationships between reservoir storage, inflow, and release. These relationships typically contain considerable noise, e.g., different release values for the same storage level during the same season. As a result, these uncertainties influence considerably the parameterization of the model derived to represent the reservoir operation based on historical observations of each reservoir. This is particularly true for the algorithm presented because of two main factors. Firstly, the presented reservoir algorithm assumes that the relationship between reservoir storage and releases follows piecewise-linear functions. There is a chance that other functional forms represent such relationships better for some reservoirs. Secondly, in the case of the generalized parameterization, the bending points in the piecewise linear functions (zone classification points) are estimated based on fixed probabilities of exceedance extracted from historical data for all reservoirs. A different dataset (of reservoirs and/or time periods) could result in different quantiles. The assumption of having similar bending points of the piecewise-linear functions for all reservoirs cannot provide optimal zones for each reservoir. However, we showed that the generalized parameterization performs better compared to other widely used algorithms.

Optimizing storage and release parameters allows us to overcome the limitation of generalized bending points of the piecewise-linear function by adjusting the bending points so that the best fit can be identified. However, optimization usually does not provide a perfect storage-release relationship (i.e., in general, the trade-off between objectives never converges to a single point) because the perfect representation only occurs in the case of a perfect reservoir model and perfect data. The proposed model, like many other types of models, is not an exception because of the uncertainties highlighted in the previous point. Thus, the trade-off between storage and release objectives can be viewed as a measure of the limitation of the reservoir algorithm (piecewise linear functions, fixed number of zones, etc.) and observation errors. To examine the level of uncertainty of the trade-off, it is important to look at the shape and range of the trade-off on each objective function axis.

As shown in Figs. 11 and 12, apart from a few reservoirs, the range of Pareto solutions for each objective function is generally narrow, with good NSE values (Figs. 11 and 12). In such cases, the associated uncertainties are few, and the trade-off between improving simulated releases and improving simulated storage is minimal. Conversely, in some cases, an extended spread of the trade-off along one of the axes (objective function) was observed, indicating a higher uncertainty of the algorithm for the process that the axis represents, i.e., reservoir storage or release. This requires further investigation of the datasets and parameterization for those reservoirs and their history of operations. Shifts in operational management of reservoirs do occur, and these may obscure the parameterization. These may be detected by careful examination of the available records as well as metadata records of the reservoir history if accessible. The level of noise when determining the parameters could be an indicator of changes in operation.

4.7 Implementation strategies to overcome data limitation

The data requirement is the main limitation of the DZTR model for application at continental and global scales. One approach to overcome data limitations is to integrate our proposed method in land-surface and catchment models along with other reservoir operation methods (e.g., Hanasaki et al., 2006). Then, within the land-surface and catchment models, identifier flags can be used to indicate which method applies to which reservoirs. The DZTR approach can only be activated for reservoirs with data support, while the remaining reservoirs can use other approaches as dictated by data availability. We follow such an implementation within the MESH model.

As shown in our results, reservoir regulation has a large impact on downstream flows if the reservoir is highly regulated and/or is of the multi-year type (c>0.5). Thus, more emphasis can be placed on those reservoirs with c>0.5. At the moment, such methods will be more effective at the regional than global scale (for example for the Saskatchewan River basin in our case) because modelers at the regional scales have better access to inflow–storage–outflow data and have a better understanding of the system to acquire the necessary reservoir data. In a land-surface hydrologic model, important reservoirs are those causing large changes to the downstream flows, and those tend to be the larger ones with generally better data availability.

Data on reservoir storage, inflow, and release exist for most reservoirs, but sometimes they are not made publicly available. Storage data can be obtained from water-level data, which are generally available for major reservoirs and can be converted to storage. Release data can be deduced from the nearest downstream station. In addition, new initiatives are needed to gather and archive such reservoir datasets and move beyond the information on reservoir characteristics that is currently available in databases (e.g., GRanD database – Lehner et al., 2011). One of our recommendations is that the target release and storage data be archived for public use at least for highly regulated and multi-year type dams (c>0.5).

The possibility of estimating storage and release data from different satellite data products is promising; such new data sources will potentially improve the use of methods like the presented reservoir operation (optimized or generalized). More recently, Busker et al. (2019) showed an estimation of volume for 130 reservoirs using surface water dataset and satellite altimetry; this is an encouraging approach for reducing data limitation.

5 Summary and conclusions

Human interventions in hydrologic systems through dams and reservoirs significantly change the flow regime of many rivers. In this paper, we presented an improved reservoir operation model, called the dynamically zoned target release (DZTR) model, that can be integrated into any large-scale hydrological model; here we integrated it into the MESH land-surface–hydrology model. The DZTR model is based on parametric piecewise-linear functions that approximate reservoir-release rules used by reservoir operators. We proposed two strategies to identify the parameters of this model: one based on the distributions of historical storage and release to generate the so-called generalized parameters and the other one based on direct calibration to observed storage and release time series via multi-objective optimization. We first tested the DZTR model individually across a number of reservoirs around the globe and then tested its performance when plugged into the MESH model for a subset of those reservoirs. Our conclusions can be summarized as follows:

• The DZTR reservoir operation model performed well in reproducing observed storage and release time series in (almost) all reservoirs tested and outperformed the existing reservoir models proposed by Hanasaki et al. (2006) and Wisser et al. (2010). The model was capable of capturing inter- and intra-annual variability in both reservoir storage and release.

• As expected, calibration significantly improved the performance of the DZTR model compared with the performance of the generalized parameters. However, a significant trade-off exists between fitting reservoir storage versus release, signifying the importance of accounting for both storage and release in a multi-objective fashion.

• The integration of the DZTR reservoir model into the MESH land-surface–hydrology modeling system was straightforward and improved the overall model performance compared with the traditional methods of accounting for reservoirs in H-LSMs. This integration can be viewed as a successful example for improving the representation of reservoir operation in CMs, LSMs, and GWSMs.

Future research work may include (1) examining the applicability of the DZTR model for regions with severely limited data by examining the utility of other data sources such as those derived from satellite-based observations (Savtchenko et al., 2004; Garambois and Monnier, 2015; Gao et al., 2012) and using the area–volume relationship approximated by regular geometric shapes (e.g., Yigzaw et al., 2018) and (2) examining direct one-way and/or two-way coupling of WMMs with CMs and LSMs for developing a seamless coupled framework for the simulation of naturally engineered watershed systems.