Estimation of Reservoir Volumes at Drafts of 40% - 90%: Drought Magnitude Method Applied on Monthly River Flows from Canadian Prairies

Tribeni C. Sharma, Umed S. Panu^*

Department of Civil Engineering, Lakehead University, Thunder Bay, ON, Canada.

DOI: 10.4236/jwarp.2022.148030                

Abstract

The draft ratios for sizing the reservoirs can vary within a wide range (40% - 90% of the mean annual flow, MAF), depending upon the demands for water by various users, and environmental and ecological considerations. The reservoir volumes based on the drought magnitude (DM) method were assessed at aforesaid draft ratios using monthly-standardized hydrological index (SHI) sequences of 10 Canadian rivers located in the Canadian prairies and northwestern Ontario. These rivers are typified by a high level of persistence lag-1 autocorrelation, ρ_1m ≥ 0.50 and up to 0.94) and coefficient of variation (cv_o) in the range of 0.42 to 1.48. The moving average (MA) smoothing of monthly SHI sequences formed the basis of the DM method for estimating reservoir volumes. The truncation or cutoff level in the SHI sequences was found as SHI_x [=(α - 1)μ_o/σ_o], [(α - 1)μ_o/σ_max], or [(α - 1)μ_o/σ_av], where α (=0.40 to 0.90) is the draft ratio i.e. proportion of the MAF, μ_o and σ_o are the overall mean and standard deviation of the monthly flows, σ_max is the maximum value of standard deviations and σ_av the average of 12 monthly values. The failure probability levels (PF) were fixed at 5%, 2.5% and 0% (corresponding reliability of 95%, 97.5% and 100%). The study revealed that the coefficient of variation is the most important parameter that influences the reservoir size while the role of lag-1 autocorrelation (ρ_1m) appears more pronounced at high draft ratios, α such as 0.90, 0.80 and 0.70 in increasing the reservoir size. The DM based method can be regarded as an alternative to Behavior analysis for sizing reservoirs at the desired probability of failure or reliability level.

Keywords

Draft Ratio, Extreme Number Theorem, Markov Chain, Moving Average Smoothing, Reliability, Standardized Hydrological Index, Sequent Peak Algorithm

1. Introduction

The reservoir volumes can be estimated involving annual, monthly or weekly flow sequences. It has been shown by Sharma and Panu [1] that the monthly scale of flow sequences is adequate for estimating the reservoir volume (V_R) corresponding to a certain draft ratio (draft in proportion to the mean annual flow, MAF). The annual scale may mask the seasonal variability and persistence characteristics needed in the estimation process; whereas, the weekly scale only adds a little refinement over the results derived using the monthly scale and also renders the analysis rigorous and complex [1]. The draft ratio (α) is expressed as the proportion to the mean annual flow (MAF) such as 100% (α = 1 MAF), 90% (0.9 MAF), 75% (0.75 MAF), 60% (0.60 MAF) etc. Knowing that in terms of magnitude, the mean annual flow (MAF) and the overall mean monthly flow (µ_o) are equivalent to each other and thus, MAF and µ_o are used interchangeably at the monthly scale. The analysis can be done in terms of observed (recorded) flow sequences (without standardization) or in terms of statistically standardized (month-by-month) flow sequences. The standardized flow sequence is termed here as the standardized hydrologic index (SHI) sequence [1] [2]. In brief, SHI is an entity with mean = 0 and standard deviation = 1, while retaining the probabilistic structure/character of the flow sequence under consideration. In the ensuing text, the notations µ_o and σ_o stand for the mean and the standard deviation of the non-standardized (i.e., observed or recorded) monthly flow series and are also termed as an overall mean and standard deviation; σ_av is the average (arithmetic mean) of the 12 calendar monthly values of the standard deviation of the observed monthly flows; ρ_1m is the lag-1 autocorrelation in the SHI sequences, and other symbols are explained wherever they appear first. The term coefficient of variation, cv_o is defined as (σ_o/µ_o), which refers to the observed monthly flow series (non-standardized), and cv_av as (σ_av/µ_o), which also refers to the observed monthly flows. In the domain of hydrologic design of reservoirs, there are other methods such as the Sequent Peak Algorithm (SPA), Behavior Analysis, Dincer Procedure, Gould Gamma Procedure, and Gould Transitional Probability method, and details on these methods are well documented in hydrological and water resources books [3] - [8] and journal papers [9] - [14], among others. The SPA can be touted as the universal method for sizing reservoirs because of its popularity and wide familiarity and can be taken as a reference for comparing and complementing the results from other methods.

In terms of the hydrological drought analyses and more specifically for the DM method, two main parameters of interest in a drought episode are duration (length, L) and magnitude (M, also termed as severity), which have been the subject of intensive study [15] - [25]. The term M is a standardized entity (dimensionless) such that the drought deficiency volume D = σ × M [16]. For the estimation of reservoir volume, the term D can be regarded as a counterpart term to V_R. The modelling activity for drought magnitude is of paramount importance in terms of management of waters, and consequently for sizing and operation of reservoirs. Using the concept advanced by Dracup et al. [15], Sharma and Panu [23] have suggested a simple model for estimating drought magnitude by coupling it with the drought length through a third parameter namely, drought intensity, i.e. magnitude (M) = intensity (I_d) × duration (L). They have carried out analyses in terms of SHI sequences [1] [2]. The concept of SHI has performed well in the modelling of drought lengths and magnitudes on Canadian river flow sequences at annual, monthly and weekly scales [23] [24].

From an earlier study (Sharma and Panu [1]), the analysis using the draft at 100% of the mean annual flow (or α = 1) resulted in an excessively high volume of reservoirs. A median value of the draft for design purposes has been identified as 75% of MAF (α = 0.75) [3] [5] [8], which was analyzed for sizing reservoirs in Canadian rivers using the monthly SHI sequences by the DM method [2]. In these investigations, the DM method was compared to SPA at the probability of failure (PF) level of 0%. The term PF is defined as the ratio of the number of time units the reservoir became empty to the total number of time units used in the analysis [3]. Accordingly, reliability is defined as (100-PF). The crucial parameter in the DM method was a scaling parameter (Φ) for drought length, which emerged equal to 0 and 0.5. Also, the concept of moving average (MA) of SHI sequences was introduced in the DM method with smoothing steps MA1 and MA2. The MA2 smoothing means the moving average of two consecutive values of SHI and MA1 means the SHI sequence itself. Succinctly, the results of the DM method were based on the aforesaid fixed values of Φ and MA smoothing (MA1 and MA2 steps) and were found to be satisfactory and compatible with SPA based estimates.

In a worldwide study of reservoirs, McMahon et al. [14] noted that the draft ratio, α varies from 0.90 to 0.10 with the median value of 0.47 in Australia and 0.29 in South Africa. In the US, the α displays a large variability between the eastern and western regions because of wide variability in hydrologic conditions. However, α has been reported [14] to vary from 0.40 to 0.90 in the eastern region. In view of the proximity of Canada to the USA, and the hydrologic conditions of the Canadian rivers resembling the eastern USA, this paper attempts to demonstrate the DM method in assessing the size of reservoirs for a range of α from 0.90 down to 0.40. Three levels of the probability of failure (i.e., 5%, 2.5% and 0% respectively corresponding to the reliability of 95%, 97.5% and 100%) are considered. The motivation for the current study is to apply the DM method in the above range of draft ratios, which are used in the North American settings and evaluate its promise. In particular, the modus operandi of implementation of the DM method concerning the fixed or variable nature of Φ values in combination with MA smoothing steps was the main focus of the study.

A set of ten rivers from the Canadian prairies and northwestern Ontario exhibiting a high degree of persistence (ρ_1m ranging from 0.50 to 0.94) and a wide range (0.42 to 1.48) of cv_o (coefficient of variation) are the focus of analysis. In addition, such a wide range in autocorrelation and coefficient of variation offered an opportunity to identify the exclusive role of these parameters in influencing the reservoir volume by the DM method.

2. Preliminaries on Methods for Sizing Reservoirs

The popular method for sizing reservoirs is the sequent peak algorithm (SPA) which is well documented in [3] [6]. The SPA requires the historical or synthesized river flow data as inflows and demand levels as outflows, and data are analyzed using the original flows (or without standardization) in a sequential format. Computations are numerically conducted using the differences between cumulative inflows and demand values to assess the size of a reservoir. The procedure is fully amenable to computerized computations to arrive at the required reservoir volume (V_R) for a given situation. The SPA hinges on the reservoir water balance equation in which the reservoir is assumed to be initially full and thus does not allow the determination of the storage at a prefixed level of probability of failure (PF). However, the PF can be evaluated for the SPA based storage volume, which would turn out to be zero. The other variant of the SPA is the Behavior analysis, in which various storage levels in the reservoir (not initially full) are evaluated at the desired level of PF by a trial and error procedure.

In the DM method, SHI sequences are obtained after standardization (month-by-month) and in turn are chopped at a suitable truncation level corresponding to draft ratio, α. The truncation (or cutoff) level considered are equivalent to $SH{I}_o=\left(\alpha {\mu }_o-{\mu }_o\right)/{\sigma }_o=\left(\alpha -1\right){\mu }_o/{\sigma }_o$. The other variants of the truncation level SHI_x are considered as: $SH{I}_{\max }=\left(\alpha {\mu }_o-{\mu }_o\right)/{\sigma }_{\max }=\left(\alpha -1\right){\mu }_o/{\sigma }_{\max }$, $SH{I}_{av}=\left(\alpha -1\right){\mu }_o/{\sigma }_{av}$, $SH{I}_{gm}=\left(\alpha -1\right){\mu }_o/{\sigma }_{gm}$, and $SH{I}_{har}=\left(\alpha -1\right){\mu }_o/{\sigma }_{har}$ in which σ_max is the maximum value, σ_av is the mean (arithmetic), σ_gm is the geometric mean and σ_har is the harmonic mean of standard deviations of the 12 monthly flows [2]. The σ_gm is defined as = ${\left[{\sigma }_1\times {\sigma }_2\times {\sigma }_3\times \cdots \times {\sigma }_{12}\right]}^{1/12}$ ; σ_har as = ${\left[\left(1/12\right)\left\{{\sigma }_1^{-1}+{\sigma }_2^{-1}+{\sigma }_3^{-1}+\cdots +{\sigma }_{12}^{-1}\right\}\right]}^{-1}$.

The SHIs below the truncation level are referred to as the deficit (d), whereas above the level are referred to as the surplus (s). In a historical record of N (=T) months, several spells of deficit and surplus shall be encountered and the length of the longest spell (representing L_T) is recorded. Likewise, the corresponding SHIs are added to represent the largest magnitude (M_T). The foregoing mode of computation of L_T and M_T is hereafter referred to as the counting procedure. The largest deficit volume (D_T) during the longest drought period (L_T) is computed as D_T = σ_av × M_T. Since the unit of D_T is the same as that of σ_av because M_T is a dimensionless entity. It is worthy to emphasize that the computation of D_T in the DM method using the counting or the modelling methodology is equivalent to the computation of V_R in the SPA or the Behavior analysis.

While using models for estimating a value of M_T, the following probabilistic relationship through the use of the extreme number theorem [17] [23] provides the basic conceptual link.

$P\left(M_T\le Y\right)=\exp \left[-Tq\left(1-q_q\right)\left(1-P\left(M\le Y\right)\right)\right]$ (1)

In which, q represents the simple probability of drought and q_q represents the conditional probability that the present period is a drought given the past period was also a drought, and T is equivalent to the return period; M stands for the drought magnitude which takes on non-integer values represented by Y, and P(.) represents the notation of cumulative probability. The parameters, viz. simple probability, q, and conditional probability, q_q can be obtained through a counting procedure or a relationship as documented in [23]. Since Y’s (such as Y₁, Y_2, Y₃, Y₄, …) correspond to values of M, thus the largest of them will represent M_T. In the above expression, M is construed to follow a normal probability distribution function (pdf) with mean and variance related to the mean and variance of drought intensity (I_d) and a characteristic drought length (L_c). The term L_c is related to the mean drought length (L_m) and the extreme drought length (L_T) through the linkage relationship: $L_c=\Phi L_m+\left(1-\Phi \right)L_T$. The parameter Φ varies between and includes 0 and 1 (i.e. 0 ≤ Φ ≤ 1) and is very crucial in the estimation process. Further details on computing M_T-e (which is the expected value of M_T and hereafter is designated as the estimated value) are shown in Appendix A. The final expressions for computational purposes are expressed as follows.

1) For M_T-e based only on the mean of the drought intensity:

$M_{T\text{-}e}=\text{abs}\left[-\left\{\frac{\exp \left(-0.5z_0^2\right)}{q\sqrt{2\pi }}\right\}-z_0\right]L_c$ (2)

In which, z₀ is the standard normal variate at the drought probability, q.

2) For M_T-e based on both the mean and variance of the drought intensity:

$M_{T\text{-}e}={\displaystyle {\sum }_{j=0}^{n_1}\frac{Y_{j+1}+Y_j}{2}}\left[P\left(M_T\le Y_{j+1}\right)-P\left(M_T\le Y_j\right)\right]$ (3)

The reliability criteria generally adopted for the design of reservoirs are based on the probability of failure (PF) defined earlier in the text [3], i.e. PF is the ratio of the number of months the reservoir went dry to the total number of months used in the analysis. The crucial elements in the process of sizing reservoirs in the DM method are: 1) selection of the cutoff level (SHI_x) for truncating the monthly SHI sequences at the aforesaid levels of PF and 2) selection of a correct value of the parameter Φ.

3. Data Acquisition and Calculations of Reservoir Volumes

Seven rivers ({1} to {8}) from Canadian prairies (Table 1, Figure 1) and 2 rivers ({9} and {10}) from northwestern Ontario bordering the prairies were involved in the analysis. The rivers encompassed drainage areas ranging from 97 to

Figure 1. Location of the river flow gauging stations used in the analysis [source: Environment Canada].

Table 1. Summary of statistical properties of monthly flows of the selected ten rivers.

A value indicated by an asterisk (*) is the overall coefficient of variation (cv_o = σ_o/µ_o) based on non-standardized monthly flows. ^† indicates the average coefficient of variation (cv_av = σ_av/µ_o). Suffix “gm” stands for the geometric and “har” stands for the harmonic mean of 12 values of monthly standard deviations.

119,000 km² with the database spanning from 35 to 108 years. The chief criterion for the selection of these river-gauging stations is the high level of persistence (ρ_1m ≥ 0.50). The flow data for these 10 rivers were extracted from the Canadian Hydrological database, Environment Canada [26]. The statistical parameters viz. µ_o, σ_o, σ_max, σ_av, σ_gm, σ_har and ρ_1m were computed and are shown in Table 1. At the monthly scale, the valueσ_av was found to be the best estimator of standard deviation for the calculation of drought deficiency volumes, D = σ_av × M [23].

The analysis required the moving averages of the monthly SHI sequences for the evaluation of M, particularly at high draft ratios. A monthly SHI sequence, without any moving average, is designated as moving average 1 (MA1), which essentially is the monthly SHI sequence itself. At times, there arose a need to perform the moving averages of two consecutive values or three consecutive values of monthly SHIs, which are respectively designated as MA2 or MA3. The MA1 sequence of SHI was used for counting as well as estimating values of L_T and M_T. The M_T values so counted are essentially the observed values and thus are designated as M_T-o. Likewise, the estimated M_T values are labelled as M_T-e.

The aforesaid statistics were also evaluated for MA2, MA3 and MA4 sequences, whenever such a need arose. The MA2 sequence is characterized by mean = 0 with σ < 1 and lag-1 autocorrelation (designated as ρ_1m2) > than ρ_1m. The MA2 sequence was further standardized and all computations were accomplished in this new standardized (designated as MA2’) domain along with deriving of M values. A similar strategy would also apply to MA3 and MA4 sequences. For all MA smoothing (MA1 through MA4), the multiplier to M for obtaining D would be σ_av (D_T-o = σ_av × M_T-o; and likewise, D_T-e = σ_av × M_T-e). An illustration of computations for sizing reservoirs by SPA and the DM method for three typical rivers is provided viz. the Beaver River, the English River, and the Churchill River (Table 2).

It may be highlighted that the Beaver River flows tend to be highly variable with cv_o = 1.48 and ρ_1m = 0.78, whereas the Churchill River is characterized by

Table 2. Values of parameters used in relevant equations for calculating M_T-e in the DM method.

Asterisk (*) with SHI_m means SHI_max as the cutoff level.

cv_o = 0.42 but with high ρ_1m = 0.94. These two rivers thus exhibit so many contrasting features in terms of persistence and variability characteristics. The parameters viz. probabilities (q, q_q, q_p and z₀) needed for computations in Equations (1) through (3) were accomplished by the counting method and are shown in Table 2. In the strict sense, z₀ is the standard normal variate corresponding to probability q. In this paper, since q is derived by the counting method corresponding to the cutoff level equal to SHI_x, therefore, SHI_x can also be approximated by z₀.

Computational Strategy

In the DM method, SHI sequences were chopped at the level of SHI_x and the D_T-o values were counted (column 5, Tables 3-5). Thus, the evaluated values of D_T-o only yielded one value of PF (column 6, Table 3) and were away from the desired PF criteria (5% and 2.5%). To overcome this difficulty, the D_T-e values were estimated by plugging a suitably assumed value of Φ including and between 0 and 1 such that the criterion of PF = 5%, 2.5% or 0% was met.

Towards achieving the desired PF, one needs to consider three variables, i.e. 1) cutoff level (SHI_x = SHI_o, SHI_max, SHI_av or SHI_gm); 2) MA step (MA1, MA2, MA3 and MA4); 2) the use of an optimal value of Φ inclusive and between 0 and 1 in the relevant equations. A suitable combination of the three foregoing variables

Table 3. Estimation of reservoir volumes for the Beaver River (cv_o = 1.48, ρ_1m = 0.78) at PF = 5%, 2.5% and 0%.

Table 4. Estimation of reservoir volumes for the English River (cv_o = 0.74, ρ_1m = 0.76) at PF = 5%, 2.5% and 0%.

Table 5. Estimation of reservoir volumes for the Churchill River (cv_o = 0.49, ρ_1m = 0.94) at PF = 5.0%, 2.5% & 0%.

Asterisk (*) in Tables 3-5 means values in parentheses are the standardized values of V_R (=V_R/σ_av). Asterisk (**) in Tables 3-5 means the DM method did not perform well and hence, the Behavior analysis was used.

that fulfilled the criterion of meeting the desired PF is listed in Tables 3-5 for α ranging from 0.40 to 0.90. Three typical computational cases are presented viz. the Beaver River (Table 3) with the highest cv_o = 1.48, the English River (Table 4) with modest cv_o = 0.74, and the Churchill River (Table 5) with cv_o = 0.49.

In the case of the Beaver River, the SHI_max proved to be the best cutoff level in combination with step MA1 (except, α = 0.40). However, the computed PF values were not exactly 5% and 2.5% (column 9) but were within the close vicinity i.e., 5.07% and 2.58% etc. level with MA steps from MA1 to MA4 (Table 2).

The Churchill River performed well up to α = 0.50 but it encountered a problem at α = 0.4, where it could not reach the level of PF = 5%. Therefore, at this level of PF, the Behavior analysis was conducted to obtain a value of D_T-e and its corresponding value of M_T-e (Table 5). Equally, an intriguing spectacle arose with the Gods River (lowest cv_o of 0.42 among all rivers) in which at α = 0.90, 0.80 and with SHI_o performed well with MA1 step and for α varying from 0.70 to 0.60, the SHI_av proved better with MA1 step (for brevity, the table for the Gods River is not included). However, for this river, a difficult spectacle erupted with α ≤ 0.50 where the cutoff needed was SHI_gm. Even with such low levels of cutoff, the criteria of meeting the 5% and 2.50% of PF were not attainable. At this point, the cutoff with SHI_har (lowest value of cutoff) was considered, which also did not improve the results. Faced with this intrigue, the Behavior analysis was used which succeeded in providing a value of the D_T-e at the 0% PF level but was unable to provide D_T-e at the PF levels of 2.5% and 5%. The D_T-e so computed met the PF criterion up to 2.36% and 3.99% only. That means for highly correlated flows such as that of the Gods River, the DM method proved less useful and the Behavior analysis was able to provide approximate results within the vicinity of the PF levels of 2.5% and 5% at draft ratios (α) ≤ 0.50. However, at α > 0.50, the DM method worked well in the entire range of autocorrelations reported in this paper (ρ_1m from 0.50 to 0.94).

In short, for highly correlated flows with cv_o of 0.42 at α = 0.40, the storage requirement is greatly decimated to meet the criterion of PF = 2.5% and 5%. Similar behavior was observed for the Island Lake, Sturgeon, and Churchill rivers at α = 0.4. This value of storage volume was counterchecked with SPA (column 2, Tables 3-5), which tallied very well for all rivers (Table 1) with cv_o < 0.54. Under such a dilemma, it is desirable that the volume may be estimated with PF = 0% only, which can be used for sizing the reservoir at all PF levels, albeit such storage volumes will be superfluously larger than optimally desired.

4. Results and Discussion

The results of analyses indicated that rivers, in general, can be categorized into two groups: 1) exhibiting ρ_1m within the range of 0.50 to 0.78 viz. Bow, Athabasca, Beaver, South Saskatchewan, English and Pipestone rivers and 2) exhibiting ρ_1m ≥ 0.87 such as the Churchill, Sturgeon, Gods and Island Lake rivers (Table 1). The DM method can easily handle the estimation detail of sizing reservoirs for rivers in category (1) along with draft conditions considered. Out of all 10 rivers, the characteristics of the Beaver River (Alberta, Canada) differ significantly from the remaining nine rivers. For example, the Beaver River originates on the boreal plains rather than on the eastern slopes of the Rocky Mountains and as a result, the runoff is not subject to the stabilizing influence of mountain snowmelt and thus displays a considerable variability from year to year and within the year. The river rises to a peak flow in the latter part of April because of spring snowmelt and occurrences of rainstorms during the snowmelt period, and thereafter the flow, generally, recedes through the remainder of the year. Unlike, the prairie streams of southern Alberta, the Beaver River tends to respond to summer rainstorms with dramatic increases in flows. Indeed, in some years, the summer flows may be greater than the peak flow during the spring runoff (Beaver River Alliance, 2013) [27]. Such unique flow characteristics are ramified in cv_o = 1.48, requiring the cutoff level = SHI_max for all values of α. Besides, at α = 0.90, MA3 and MA4 steps were needed with the cutoff level of SHI_max, whereas for other values of α, the MA1 was found to be sufficient, though with the cutoff level of SHI_max.

At high draft ratios, in a majority of cases, the variance-based relationship (Equation (3)) was required whereas at low draft ratios only the mean-based relationship (Equation (2)) served the purpose. In the same vein, on several occasions at low draft ratios for the modestly correlated SHI sequences, the MC0 performed well with Equation (2). In the computations of reservoir volumes for rivers (with high values of ρ_1m) and belonging to category (2), the MA1 step with cutoff level as SHI_o was found to easily handle α from 0.60 to 0.90. However, the computational complexities erupted for α ≤ 0.50, where a recourse had to be taken with lower levels of truncation such as SHI_av, SHI_gm and in rare cases with SHI_har. In extreme situations, such as α = 0.40, SHI_har also could not salvage the situation, though, some relief was accorded by the Behavior analysis at PF = 0% and 2.5% (for example, the Gods River, which is not presented in the table for brevity). The Behavior analysis also could not yield satisfactory results because reservoir volumes were generally very low at a 5% level of PF. In general, SHI_o with a variance-based relationship (Equation (3)) and the MA1 step should be considered, as a first approximation, with suitable modifications exercised while meeting the desired PF criterion. Likewise, one should first consider the extreme values of Φ (0 or 1) and then a suitably modified value within the range (0 and 1) be used to meet the desired PF criterion. Briefly, the most important decisive factor to meet the desired PF criterion is the identification of the appropriate value of parameter Φ.

Since a majority of reservoirs on a global basis are designed for the draft conditions within the range of α = 0.70 and 0.80, where the DM method has performed satisfactorily with the truncation level equal to SHI_x (=SHI_o) or (=SHI_max) in combination with Equation (3) and by use of the easily adjustable value of Φ. The value of Φ works like a cursor, which should be moved up and down within the range from 0 to 1. Numerical computations can be easily performed via Excel-based programming with a button of Φ, which can easily be varied until the desired PF criterion is achieved.

The DM method is capable of assessing the effect of coefficient of variation (cv_o and cv_av) and lag-1 autocorrelation (ρ_1m) affecting the size (volume) of reservoirs. Towards achieving this objective, the ratio (R_s, named as storage ratio) of the reservoir volume (computed as D_T-e, m³) to MAF volume (m³) was evaluated (Figures 2-4). Also, the R_s values at PF levels of 0% and 2.5% are summarized in Table 6. It can be seen in Table 6 that as expected, R_s values increase

Table 6. Values of volume ratios, R_s at different draft ratios (α) at the PF levels of 2.5% and 0% for the rivers in Canadian prairies and Atlantic Canada.

The rivers affixed by an asterisk (*), (**) and sign (^†) refer to as pairs used for comparison. The bold italic numbers indicate the commonality of parameters in the pair.

(a) (b)

Figure 2. Effect of lag-1 autocorrelation (ρ_1m) in the monthly SHI sequences on reservoir volume at increasing draft ratios and PF levels: (a) 2.5% and (b) 0%.

(a) (b)

Figure 3. Effect of coefficient of variation (cv_o) in the monthly flow sequences on reservoir volumes at increasing draft ratios and PF levels: (a) PF = 2.5% and (b) PF = 0%.

(a) (b)

Figure 4. Effect of coefficient of variation (cv_av) in the monthly SHI sequences on the reservoir volumes at increasing draft ratios and PF levels: (a) 2.5% and (b) 0.0%.

consistently with an increase in the value of α for all rivers. The R_s values are primarily affected by 1) the autocorrelation, ρ_1m, and 2) the coefficient of variation cv_o (in case of monthly flow sequences (non-standardized), and cv_av (in case of the SHI sequences) and are discussed below.

4.1. Assessing the Effect of Autocorrelation and Coefficient of Variation

4.1.1. Effect of Autocorrelation on the Size of Reservoir Storage

To discern the effect of autocorrelation (ρ_1m) in river flows, a pairwise approach was adopted. The advantages of the pairwise approach can be cited as 1) the ease of comparison of the results based on the existing framework of contrasting data, 2) a quick visual display of the crucial information, and 3) empirically or experimentally confirming the behavior of the anticipated process outputs. Therefore, the pair of the Bow River [{1}, Table 6] and the South Saskatchewan River [{4}, Table 6] was chosen in which both rivers display nearly the same level of cv_o (≈1.05 to 1.06, Table 6) but ρ_1m differs significantly (0.50 to 0.64). A graph of R_s versus α shown in Figure 2(a) (Bow River) and Figure 2(b) (South Saskatchewan River) exhibit a stark effect of lag-1 autocorrelation on the reservoir volume. In particular, R_s rises sharply at the high draft ratios, implying that rivers with highly auto-correlated flows require large volumes, especially during extended droughts to meet the draft requirements. Likewise, during high flows, such rivers also need large reservoirs to store the floodwaters for flood control purposes. This observation is in agreement with the findings reported in the literature by Phatarfod [28], Otieno and Ndiritu [29]. Both the aforesaid authors unanimously agreed that the role of serial correlation becomes pronounced at high values of α, which compares very well with the graphical presentation shown in Figure 2(a) and Figure 2(b), computed by the DM method. It can be inferred that the pairwise approach is a simple and straightforward way of looking at the behavior of storage characteristics in relation to the autocorrelation in the inflows to the reservoir.

To further validate the effect of autocorrelation, two additional rivers from an earlier publication [2], the Upper Humber River (Newfoundland, Canada) and the North Margaree River (Nova Scotia, Canada) were analyzed. The statistical parameters of these rivers are summarized in Table 6. These rivers are characterized by low but significant values of ρ_1m (=0.13 and 0.17) [24]. These values of ρ_1m are far below compared to rivers from Canadian Prairies (Table 1 and Table 6). The R_s values of these rivers were computed by estimating D_T-e values by following the procedure used for other rivers. It is apparent from Table 6 that the R_s values of these rivers are very low and less than 0.50 even for α = 0.90, which is indicative that the influence of autocorrelation is very strong in increasing the reservoir volumes. In particular, with the near similar conditions of cv_o and cv_av between the North Margaree River and the English River, but with a large difference between values of ρ_1m (Table 6), the impact of autocorrelation can be appreciated because R_s values are much larger for the English River. Similar behavior is exhibited by a pair of the Upper Humber River and the Pipestone River with near comparable values of cv_o and cv_av but with a stark difference in ρ_1m, the R_s values are much higher for the Pipestone River at high draft ratios.

4.1.2. Effect of Coefficient of Variation on the Size of Reservoir Storage

The effect of the coefficient of variation on reservoir volumes was investigated by choosing a pair of rivers in which values of ρ_1m are close to each other while the values of cv_o and cv_av differ significantly. The identical conditions of autocorrelation structure of monthly flows (ρ_1m = 0.78) displayed by the Beaver River ({3}, Table 6) and the English River ({9}, Table 6), offered an opportunity to reveal the effect of the coefficient of variation on reservoir volumes. Knowing that monthly flows of the Beaver River are highly variable with cv_o (=1.48) compared to the English River with cv_o (=0.76), the reservoir volume for the Beaver River tends to be exceptionally high at high draft ratios. Reservoirs, by virtue, damp out the variable inflows to reasonably steady outflows but construction of such large reservoirs involve exorbitant cost.

The foregoing impact of variability in river inflows on the reservoir storage is further investigated by considering the coefficient of variation (cv_av) of monthly SHI sequences for two rivers as a pair viz. the Sturgeon River ({6}, Table 6), cv_o = 0.53, ρ_1m = 0.94,) and the Island Lake River ({7}, Table 6), cv_o = 0.54 and ρ_1m = 0.87). A review of the statistical parameters of these rivers indicates that the values of cv_o are almost identical but there exists some degree of dissimilarity in the autocorrelation structure of these rivers. Since the values (0.87 and 0.94) of ρ_1m lie within proximity of each other and for pragmatic reasons, these rivers can be perceived to belong to the same group of highly auto-correlated flows. The cv_av varies from 0.38 (Island Lake River) to 0.49 (Sturgeon River), and such a difference can be construed to be significant because its effect is reflected in Figure 4(a) and Figure 4(b) at the PF levels of 2.5% and 0% where an apparent rise in the values of R_s is exhibited.

There are two measures of variability between the monthly data sets considered in this paper, viz. cv_o and cv_av. Since the DM method requires stationarity in the data structure and therefore cv_av (standardized or the SHI sequence) plays a greater role. However, for the SPA and Behavior analysis, the data stationarity is not required and thus the role of cv_o becomes more apparent. In general, it can be perceived that high values of cv_av are associated with high values of cv_o, and the effect of cv_o is reflected in Figure 3, where R_s values rise steeply with increasing values of α. The role of cv_av seems to erupt in sync with cv_o as displayed in Figure 4, i.e. R_s increases with rising α in tandem with Figure 3.

Succinctly, the analysis amply demonstrates that the DM method is applicable over the entire range of draft ratios commonly encountered in the design of reservoirs at the desired PF levels. However, in the implementation of the DM method, the parameter Φ is variable (it is not confined to 0 and 0.5 as is the case with α = 1 and 0.75) and has to be calibrated for individual cases of a river and α being considered. Further, the MA smoothing from MA1 to MA4 steps are essential ingredients of the method with a greater number of MA steps needed at draft ratios equal to and greater than 0.7 associated with high cv_o. The results thus allude to the fact that though the role of coefficient of variation is well understood in the design of reservoir capacity, the impact of autocorrelation is also no less important. This finding underscores the significance of autocorrelation in the design of reservoirs and the DM method handles the role of autocorrelation lucidly through the use of the Markov chain, the extreme number theorem and MA smoothing based concepts.

In the paper, all the parameters (probabilities) were determined by the counting method after choosing a suitable cutoff level. No consideration was given to the gamma pdf of the flow and SHI sequences in the calculations. In fact, calculations are strictly applicable to the normal pdf of flows and SHI sequences. However, the procedural detail can be applied to the gamma pdf of flows and subsequently to SHI sequences by transforming SHI’s into the normal probability domain [24]. From the latter with rigorous computations involving the gamma pdf, it was observed that the outcome was only marginally improved from the one adopted in this paper. In other words for gamma probability distributed flows the calculations can be performed directly in terms of statistically standardized SHI sequences, as was conducted and demonstrated in this paper. In general, the DM method satisfactorily computes the reservoir volumes with α > 0.5, whereas at α ≤ 0.5, the DM method is less efficient, particularly with highly auto-correlated river flows. Under this situation, the Behavior analysis can be used to provide the needed results on the reservoir volumes while meeting the target PF criterion.

5. Conclusion

The DM method can be successfully utilized to estimate the reservoir volumes equivalent to drought deficiency volumes, D_T-e. The method utilizes the SHI sequences as a platform where the parameters viz. probabilities: q, q_q, q_p can be computed using the counting method by truncating the SHI sequences at the cutoff level equal to SHI_x. The SHI_x can take the form of SHI_o, SHI_max or SHI_av. The counterpart standard normal variate z₀ needed for computing the drought intensity at the probability, q can be treated equivalent to the aforesaid cutoff level equal to SHI_x. The DM method can be implemented with the proper choice of parameter Φ (within the bound of 0 and 1) and MA smoothing steps to meet the probability of failure criterion at 5%, 2.5% and 0%. Most importantly, the method can be used to assess the effect of cv (cv_o and cv_av) and lag-1 autocorrelation (ρ_1m) influencing the size of reservoirs. It was observed that the coefficient of variation (cv_o) is the most important parameter to influence the reservoir size, that is, a higher value of the coefficient of variation leads to a higher storage volume. Likewise, the autocorrelation in the river inflows increases the requirement for higher reservoir capacity, which rises steeply at higher draft ratios.

Acknowledgements

The partial financial support from the Natural Sciences and Engineering Research Council of Canada for this paper is gratefully acknowledged.

Appendix A

The following probability-based relationship from the first principles can be used to estimate the expected value of the largest drought-magnitude, E(M_T).

$E\left(M_T\right)={\displaystyle {\sum }_{j=0}^{\infty }\left(M_T=Y_j\right)\cdot p\left(M_T=Y_j\right)}$ (A1)

The notation P(.) stands for the cumulative probability and p(.) stands for the simple probability. Since M_T is a continuous random variable, therefore M_T has a continuous pdf so p (M_T = Y_j) can be evaluated as P (M_T ≤ Y_j+1) − P (M_T ≤ Y_j) with M_T = Y_j replaced by the mean value $\left(Y_{j+1}+Y_j\right)/2$. Equation (A1) can therefore be expressed as

$E\left(M_T\right)={\displaystyle {\sum }_{j=0}^{n_1}\frac{Y_{j+1}+Y_j}{2}\left[P\left(M_T\le Y_{j+1}\right)-P\left(M_T\le Y_j\right)\right]}$ (A2)

The upper limit of integration, ∞ in Equation (A1), is replaced by some finite number n₁. The general equation for evaluating P (M_T ≤ Y_j) based on the extreme number theorem [17] can be expressed as

$P\left(M_T\le Y_j\right)=\exp \left[-Tq\left(1-q_q\right)\left(1-P\left(M\le Y_j\right)\right)\right]$ (A3)

In which M takes on non-integer values represented by Y_j. Since Y’s (such as Y₁, Y₂, Y₃, Y₄, …) correspond to values of M, thus the largest of them will correspond to M_T. It is to be noted that M is the sum of several drought deficit spikes encountered in a drought spell. Each spike has a negative sign because they are derived by truncating SHI sequences and lay on the downside (negative side) of the cutoff level. The deficit spikes can be construed to obey a truncated normal distribution because a normal pdf encompasses values of a random variable from −∞ to +∞. Thus, the standard normal pdf is truncated at various levels of z₀ (standard normal number) corresponding to counterpart values of the probabilityq. For example, a standard normal pdf can be truncated at z₀ = −0.52 with the corresponding probability of 0.3, which can be obtained from the standard normal probability table or a polynomial equation documented in Chow et al. [30].

The truncated normal pdf version will have a mean and variance respectively different from 0 and 1. Applying the basic axioms for the evaluation of moments, expressions for the mean (denoted by μ_d) and variance (denoted by ${\sigma }_d^2$ ) of the truncated normal pdf version can be deduced [17] [24] as follows.

${\mu }_d=-\frac{\exp \left(-0.5z_0^2\right)}{q\sqrt{2\pi }}-z_0$ (A4)

${\sigma }_d^2=1-\frac{z_0\exp \left(-0.5z_0^2\right)}{q\sqrt{2\pi }}-\frac{\exp \left(-z_0^2\right)}{q^{2}2\pi }$ (A5)

Because drought episodes lay below the desired truncation level, therefore an absolute value of the term μ_d is an estimator of the mean value of drought intensity (I_d), as it represents the mean of the several deficit spikes. Likewise, the term ${\sigma }_d^2$ is an estimator of the variance of drought intensity, whose value is unaffected by the sign. For example, at the cutoff level of z₀ = −0.52, q = 0.3; and substitution of these values in Equations (A4) and (A5) yields μ_dequal to −0.64 and ${\sigma }_d^2$ equal to 0.25. As explained earlier, the absolute value of μ_dwill be taken as the mean value of drought intensity (I_d), i.e. = 0.64. Likewise at the cutoff level, z₀ = −0.10, q = 0.16; μ_d and ${\sigma }_d^2$ will work out as −0.51 and 0.23.

Given the central limit theorem and since M is the sum of deficit spikes, its probability structure can be approximated by a normal pdf with meanμ_M and variance ${\sigma }_M^2$ [17]. Such a consideration reduces the expression for the term P (M ≤ Y_j) in Equation. (A3) as follows.

$P\left(M\le Y_j\right)=\frac{1}{{\sigma }_M\sqrt{2\pi }}{\displaystyle {\int }_0^{Y_j}\exp \left[-0.5{\left(\frac{M-{\mu }_M}{{\sigma }_M}\right)}^2\right]\text{d}M}$ (A6)

It was noted that the parameters μ_M and ${\sigma }_M^2$ are related to the extreme drought length, L_T and the mean drought length, L_M [24]. Such a representative drought length is named herein as characteristic drought length, L_C, which can be expressed as follows.

$L_C=\Phi L_M+\left(1-\Phi \right)L_T$ (A7)

The parameter Φ can be designated as a weighing parameter as it weighs the mean drought length L_M and the longest drought length L_T. The value of Φ varies from 0 to 1 and is determined through optimization or trial and error procedure. For first-order dependence or Markov Chain order 1 (MC1) situation, the mean length, L_M, can be expressed as follows

$L_M=\frac{1}{1-q_q}$ (A8)

The expression for the expected value of L_T at MC1 situation in drought periods can be obtained as follows [23]

$L_T=1-\frac{\log \left[FT\left(1-q\right)q_p\right]}{\log \left(q_q\right)}$ (A9)

where,F is the factor to account for the plotting position in the empirical estimation of the exceedance probability. That is, in the Hazen plotting position formula, the exceedance probability = 0.5/T (T = sample size), so the return period is equal to T/0.5 = 2T or F = 2. Likewise in the Weibull plotting position formula F = 1. In this analysis, the plotting position formula [31] as developed for Canadian rivers has been used. The formula evaluates the exceedance probability = 0.75/(T + 0.25), so F = 1.33(1 + 0.25/T) ≈ 1.33 as T is generally large. The term q_q stands for the conditional probability of the present period being drought given the previous period was also a drought and likewise, q_p stands for the present period being drought given the previous period was wet. The conditional probabilities q_q and q_p can be computed from the counting method.

The expressions for μ_M and ${\sigma }_M^2$ are related as follows [17].

${\mu }_M=L_C{\mu }_d$ (A10)

${\sigma }_M^2=L_C{\sigma }_d^2\left[\frac{1+\rho }{1-\rho }-\frac{2\rho \left(1-{\rho }^{L_C}\right)}{L_C{\left(1-\rho \right)}^2}\right]$ (A11)

Once a proper value of L_C has been determined, Equation (A6) is integrated numerically to evaluate P (M ≤ Y_j) and then the value of integrand is inserted in Equation (A3) to yield an estimate of P (M_T ≤ Y_j). Letting these values of Y_j as Y₀ (j = 0) = 0, Y₁ (j = 1) = 0.05, Y₂ (j = 2) = 0.10, Y₃ (j = 3) = 0.15, Y_n1 = 150 (for the monthly scale) with an increment, ∆ = 0.05, E(M_T) can be computed using Equation (A3). It should be noted that in Equation (A11), ρ should be replaced by ρ_1m in the relevant calculations.

A particular version for the estimation of E(M_T) can be taken as (A10) itself and Equation (A10) takes the following form

$E\left(M_T\right)={\mu }_d{L}_c=\text{abs}\left\{-\left[\frac{\exp \left(-0.5z_0^2\right)}{q\sqrt{2\pi }}\right]-z_0\right\}{L}_c$ (A12)

Note Equation (A2) involves both mean and variance of drought intensity (I_d) for the estimation of E(M_T) and therefore is more comprehensive.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Sharma, T.C. and Panu, U.S. (2021) A Drought Magnitude Based Method for Reservoir Sizing: A Case of Annual and Monthly Flows from Canadian Rivers. Journal of Hydrology: Regional Studies, 36, Article ID: 100829. https://doi.org/10.1016/j.ejrh.2021.100829
[2]	Sharma, T.C. and Panu, U.S. (2021) Reservoir Sizing at the Draft Level of 75% of Mean Annual Flow Using Drought Magnitude-Based Method on Canadian Rivers. Hydrology, 8, 79. https://doi.org/10.3390/hydrology8020079
[3]	McMahon, T.A. and Mein, R.G. (1978) Reservoir Capacity and Yield, Development in Water Science #9. Elsevier, Amsterdam, 16.
[4]	Loucks, D.P., Stedinger, J.R. and Haith, D.A. (1981) Water Resources Systems Planning and Analysis. Prentice-Hall, Englewood Cliffs.
[5]	McMahon, T.A. and Mein. R.G. (1986) River and Reservoir Yield. Water Resources Publications, Littleton, 102.
[6]	Linsley, R.K., Franzini, J.B., Freyburg, D.L. and Tchobanoglous, G. (1992) Water Resources Engineering. 4th Edition, Irwin McGraw-Hill, New York, 192.
[7]	Nagy, I.V., Asante-Duah, K. and Zsuffa, I. (2002) Hydrological Dimensioning and Operation of Reservoirs: Practical Design Concepts and Principles. Kluwer Academic Publishers, Boston, 127. https://doi.org/10.1007/978-94-015-9894-1
[8]	McMahon, T.A. and Adeloye, A.J. (2005) Water Resources Yield. Water Resources Publications, Littleton, 67.
[9]	Parks, Y.P. and Gustard, A. (1982) A Reservoir Storage Yield Analysis for Arid and Semiarid Climate. IAHS. Optimal Allocation of Water Resources, 135, 49-57.
[10]	Lele, S.M. (1987) Improved Algorithms for Reservoir Capacity Calculation Incorporating Storage-Dependent and Reliability Norms. Water Resources Research, 23, 1819-1823. https://doi.org/10.1029/WR023i010p01819
[11]	Vogel, R.M. and Stedinger, J.R. (1987) Generalised Storage-Reliability-Yield Relationships. Journal of Hydrology, 89, 303-327. https://doi.org/10.1016/0022-1694(87)90184-3
[12]	Montaseri, M. and Adeloye, A.J. (1999) Critical Period of Reservoir Systems for Planning Purposes. Journal of Hydrology, 224, 115-136. https://doi.org/10.1016/S0022-1694(99)00126-2
[13]	Adeloye, A.J., Lallemand, F. and McMahon, T.A. (2003) Regression Models for Within-Year Capacity Adjustment in Reservoir Planning. Hydrological Sciences, 48, 539-552. https://doi.org/10.1623/hysj.48.4.539.51409
[14]	McMahon, T.A., Pegram, G.G.S., Vogel, R.M. and Peel, M.C. (2007) Revisiting Reservoir Storage-Yield Relationships Using a Global Streamflow Database. Advances in Water Resources, 30, 1858-1872. https://doi.org/10.1016/j.advwatres.2007.02.003
[15]	Dracup, J.A., Lee, K.S. and Paulson Jr., E.G. (1980) On the Statistical Characteristics of Drought Events. Water Resources Research, 16, 289-286. https://doi.org/10.1029/WR016i002p00289
[16]	Yevjevich, V. (1983) Methods for Determining Statistical Properties of Droughts. In: Yevjevich, V., da Cunha, L. and Vlachos, E., Eds., Coping with Droughts, Water Resources Publications, Littleton, 22-43.
[17]	Sen, Z. (1980) Statistical Analysis of Hydrological Critical Droughts. ASCE Journal of Hydraulic Engineering, 106, 99-115. https://doi.org/10.1061/JYCEAJ.0005362
[18]	Sen, Z. (2015) Applied Drought Modelling, Prediction and Mitigation. Elsevier Inc., Amsterdam, 154.
[19]	Guven, O. (1983) A Simplified Semi-Empirical Approach to Probabilities of Extreme Hydrological Droughts. Water Resources Research, 19, 441-453. https://doi.org/10.1029/WR019i002p00441
[20]	Tallaksen, L.M., Madsen, H. and Clausen, B. (1997) On the Definition and Modelling of Streamflow Drought Duration and Deficit Volume. Hydrological Sciences Journal, 42, 15-33. https://doi.org/10.1080/02626669709492003
[21]	Salas, J., Fu, C., Cancelliere, A., Dustin, D., Bode, D., Pineda, A. and Vincent, E. (2005) Characterizing the Severity and Risk of Droughts of the Poudre River, Colorado. ASCE Journal of Water Resources Planning and Management, 131, 383-393. https://doi.org/10.1061/(ASCE)0733-9496(2005)131:5(383)
[22]	Nalbantis, I. and Tsakiris (2009) Assessment of Hydrological Drought Revisited. Water Resources Management, 23, 881-897. https://doi.org/10.1007/s11269-008-9305-1
[23]	Sharma, T.C. and Panu, U.S. (2014) A Simplified Model for Predicting Drought Magnitudes: A Case of Streamflow Droughts in Canadian Prairies. Water Resources Management, 28, 1597-1611. https://doi.org/10.1007/s11269-014-0568-4
[24]	Sharma, T.C. and Panu, U.S. (2014) Modelling of Hydrological Drought Durations and Magnitudes: Experiences on Canadian Streamflows: A Case of Streamflow Droughts in Canadian Prairies. Journal of Hydrology: Regional Studies, 1, 92-106. https://doi.org/10.1016/j.ejrh.2014.06.006
[25]	Akyuz, D.E., Bayazit, M. and Onoz, B. (2012) Markov Chain Models for Hydrological Drought Characteristics. Journal of Hydrometeorology, 13, 298-309. https://doi.org/10.1175/JHM-D-11-019.1
[26]	Environment Canada (2018) HYDAT CD-ROM Version 96-1.04 and HYDAT CD-ROM User’s Manual. Surface Water and Sediment Data. Water Survey of Canada.
[27]	Beaver River Watershed Alliance (2013) The State of the Beaver River Watershed—A Summary, Alberta, Canada.
[28]	Phatarfod, R.M. (1986) The Effect of Serial Correlation on Reservoir Size. Water Resources Research, 22, 927-934. https://doi.org/10.1029/WR022i006p00927
[29]	Otieno, F.A.O. and Ndiritu, J.G. (1997) The Effect of Serial Correlation on Reservoir Capacity Using the Modified Gould Probability Matrix Method. Water South Africa, 23, 63-70.
[30]	Chow, V.T., Maidment, D.R. and Mays, L.W. (1988) Applied Hydrology. McGraw-Hill, New York, 367.
[31]	Adamowski, K. (1981) Plotting Formula for Flood Frequency. Water Resources. Bulletin, 17, 197-202. https://doi.org/10.1111/j.1752-1688.1981.tb03922.x