Compatibility of Drought Magnitude Based Method with SPA for Assessing Reservoir Volume: Analysis Using Canadian River Flows

Tribeni C. Sharma, Umed S. Panu^*

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

DOI: 10.4236/jwarp.2022.141001                

Abstract

The traditional sequent peak algorithm (SPA) was used to assess the reservoir volume (V_R) for comparison with deficit volume, D_T, (subscript T representing the return period) obtained from the drought magnitude (DM) based method with draft level set at the mean annual flow on 15 rivers across Canada. At the annual scale, the SPA based estimates are larger, on an average of nearly 70%, compared to the DM based estimates. To ramp up the DM based estimates to be in parity with SPA based values, the analysis was conducted through the counting and the analytical procedures involving only the annual SHI (standardized hydrological index, i.e. standardized values of annual flows) sequences. It was found that MA2 or MA3 (moving average of 2 or 3 consecutive values) of SHI sequences was required to match the counted values of D_T to V_R. Further, the inclusion of mean, as well as the variance of the drought intensity in the analytical procedure, with the aforesaid smoothing led D_T comparable to V_R. The distinctive point in the DM based method is that no assumption is necessary such as the reservoir being full at the beginning of the analysis—as it is the case with the SPA.

Keywords

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

1. Introduction

A considerable amount of research can be traced to the hydrologic drought models utilizing the river flow data that focus on the estimation of drought duration and magnitude (previously termed as severity). Two major elements of the hydrologic drought studies have been the truncation level approach and the analysis by simulation and/or analytical methods. The analytical methods are pursued by the use of the frequency analyses of drought events in terms of duration and deficit volumes. The noteworthy contributions in this area of frequency analyses are that of [1] [2] [3] [4] , among others. The other route in the domain of analytical methods is the use of the theory of runs, which is well documented in [5] - [13] . Several hydrologic drought indices have been suggested such as standardized runoff index (SSI) [14] , streamflow drought index (SDI) [15] , and standardized hydrologic index (SHI) [12] [13] . These indices are essentially standardized (statistically) values of historical stream flows or in some transformed version (normalization in a probabilistic sense) at the desired time scale. The standardized hydrological index (SHI) is the standardized value (statistical) of river flows with the mean 0.0 and the standard deviation equal to 1.0, unlike the standardized precipitation index (SPI), which is normalized after standardization [16] . On the monthly time scale, it is the month-by-month standardization and so on at the weekly time scale.

The major application of the SPI refers to drought monitoring which is an essential element in the process of drought early warning and preparedness. Applications of SPI are amenable because of the widespread availability of precipitation data. Though some attempts have been made to classify the hydrological drought [15] on the lines of SPI, yet such uses of hydrological drought indices are limited. However, there have been investigations to use the SPI to relate the propagation of meteorological droughts to hydrological droughts in Spanish catchments [17] and for the U.K. catchments [18] , among others. Despite such limitations, hydrological drought indices have potential in the estimation of drought magnitude that plays an important role in the assessment of shortage of water in rivers and consequently in reservoirs. Even with the aforesaid studies, a few investigations other than Sharma and Panu [19] [20] have been made to link the deficit volume to reservoir volume, and also how and at what time scale of analysis would be aptly meaningful in this regard.

The term drought magnitude has been variously defined in the earlier literature such as the drought severity [5] [6] and the deficit volume [2] . In this paper, the term deficit volume (denoted by D) represents the deficiency or the shortage of water below the truncation level in a river flow sequence, and the drought magnitude (M) refers to the deficiency in terms of the SHI (standardized flow) sequences. The deficit volume and drought magnitude are related by the linkage relationship: D = σ × M [5] , in which σ is the standard deviation of the flow sequence. The analyses are usually conducted in the standardized domain to assess deficit volume, D through the above linkage relationship.

To the best of the authors’ knowledge, no research investigations other than those of authors [19] [20] have been reported in the literature on the application of drought indices and magnitude-based analyses, and models for sizing of reservoirs. This paper represents one of the pioneering attempts to address and bridge the above gap and demonstrate the utility of such analyses of drought magnitude in assessing the size of reservoirs. The standardized hydrological index (SHI) has been used in this analysis utilizing streamflow data from Canadian rivers. The data on annual, monthly and weekly flow sequences were analyzed using the draft at the mean annual flow for sizing of reservoirs. However, the authors’ preliminary investigations indicate that the detailed analysis related to the sizing of reservoirs be conducted at an annual scale in view of ease and simplicity in handling annual streamflow data.

2. Preliminaries on Methods for Sizing the Reservoirs

The two textbook-based methods [21] [22] for sizing the reservoirs are the Rippl graphical procedure and the sequent peak algorithm (SPA). In the Rippl method, the graphical plot of cumulative inflows as well as outflows is used to derive an estimate of the reservoir size. In the SPA, the calculations are conducted numerically using the cumulative or residual mass curve methods to obtain the estimate of reservoir volume, V_R. In the drought magnitude-based method, SHI sequences are obtained after standardization. It corresponds to truncating the annual river flow sequences at the mean level or SHI value = 0. Since the drought lengths and corresponding drought magnitudes yield conservative values for the design of reservoirs, therefore SHI = 0 as a truncation level is preferred and has been used in the analysis. The SHIs below the truncation level are referred to as the deficit (dubbed as d), whereas above the level are referred to as the surplus (dubbed as s). In a historical record of N (=T) years, there shall emerge several spells of deficit and surplus, and the longest spell length of deficits (representing L_T) is recorded. Likewise, the corresponding deficits are added to represent the largest magnitude (M_T). These deficits are being referred to as drought intensities and represent truncated values of SHIs below the truncation level. The foregoing approach of calculation of L_T and M_T is dubbed as the counting procedure in the ensuing sections. The largest deficit volume (D_T) during the drought period is computed as D_T = σ × M_T [5] . It is noted that the unit of D_T is the same as that of σ because M_T is a dimensionless entity. It is stated that either the quantity D_T obtained using the DM based counting or analytical procedure is perceived equivalent to V_R calculated by the SPA method. In the counting procedure, the entities M_T and D_T are respectively obtained from the historical or observed data and hence are denoted as M_T-o and D_T-o, where subscript “_o” stands for the observed.

Estimation of Deficit Volumes by DM Model

A majority of models for the estimation of drought magnitude implicitly involves the use of the frequency distribution of drought events [2] [4] . It is in this regard that the moving average (MA) and sequent peak algorithm (SPA) form the important tools for analysis [2] [9] . In other approaches, the probability-based relationships are hypothesized for estimating drought magnitude (M) using the relationship: drought magnitude = drought intensity × drought length [23] . As mentioned above, drought intensities essentially are deficit spikes and are derived by truncating a SHI sequence. The deficit spikes have a negative sign because each spike lay on the downside (negative side) of the truncation level, with the lower bound as −∞ and an upper bound as truncation level such as z₀, which is also a negative number with a maximum value of 0. It is tacitly assumed that the SHI sequences obey standard normal probability density function (pdf) which after truncating at the desired level (z₀) shall result in a truncated normal pdf, whose mean and variance would be different from 0 and 1. One can develop a probabilistic relationship for M_T, using the extreme number theorem [7] [24] that implicitly involves drought intensity and drought length (L_T) as expressed below.

$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 drought and T is equivalent to 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. Since Y’s (such as Y₁, Y₂, 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 pdf with mean and variance related to the mean and variance of drought intensity and a characteristic drought length. The characteristic drought length is related to the mean drought length and the extreme drought length, L_T.

At the annual level, the flow sequences in Canadian rivers have been found to follow the normal pdf [13] , leading SHI sequences to obey standard normal pdf. Therefore, the assumption of deficit spikes to obey truncated normal distribution is reasonably justified. Based on the above premises, a detailed derivation has been tracked by Sharma and Panu [19] [20] and is not reproduced here for brevity. The concluding expressions for the present paper are described as follows.

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

To compute $\left[P\left(M_T\le Y_{j+}{}_1\right)–P\left(M_T\le Y_j\right)\right]$ in Equation (2), the integration of the normal probability function is numerically performed as described in Sharma and Panu [12] . Theoretically, the upper limit of summation (n₁) in Equation (2) is ∞, but for numerical integration purposes, a finite value is chosen. For drought magnitude analysis based on annual flows, a value of n₁ = 30 (with an increment in j = 0.05 was found to be large enough to ensure sufficient accuracy in the process of numerical integration. For brevity, henceforth E(M_T) shall be written as $M_{T\text{-}e}$ , i.e., an estimated value of M_T. It may be noted that Equation (2) involves both the mean and variance of drought intensity to arrive at a value of M_T-e. Likewise, the estimated value of D_T is designated as D_T-e (=σ × M_T-e).

A particular version of M_T-e involving the mean of drought intensity (denoted as µ_d) only can be written as follows [12] .

$M_{T\text{-}e}={\mu }_d\cdot L_T=abs\left[-\frac{\exp \left(-0.5z_0^2\right)}{q\sqrt{2\pi }}-z_0\right]\cdot L_T$ (abs means absolute) (3)

where, L_T is the largest drought length obtained using Markov chain based algorithm, q is drought probability at the truncation level z₀. For example, a standard normal pdf respectively can be truncated at z₀ = 0.0 and z₀ = −0.52 and the corresponding drought probability q can be found from the standard normal probability table to be 0.5 and 0.3.

3. Data Acquisition and Calculations of Reservoir Volumes

Fifteen rivers from prairies to Atlantic Canada (Table 1, Figure 1) were involved

Table 1. Summary of statistical properties of annual flows of the rivers under consideration.

Note: The cv, γ, ρ respectively represent the coefficient of variation, skewness, lag-1 autocorrelation. Small values of skewness indicate the normal pdf of the annual flow sequences.

Figure 1. Location of the river gauging stations used in the analysis across Canada (not to the scale) [Source: Environment Canada].

in the analysis. The rivers encompassed drainage areas ranging from 97 to 56,369 km² with the data bank spanning from 38 to 108 years. The flow data for these 15 rivers were extracted from the Canadian hydrological database [25] . To increase the number of samples, some of the rivers with large data sizes such as the Bow, English, Lepreau, Bevearbank, and North Margaree were also analyzed by forming 2 - 4 subsamples with the data size of 40 years or more. This type of analysis created around 30 samples from 15 rivers to obtain a robust and reliable estimate of the performance statistics. Based on the above premises, the results of various analyses are described in the sections to follow.

The first step in the analysis was to discern the role of time scale in influencing the reservoir size. Therefore, reservoir volumes (V_R) were assessed using the SPA at the demand level equivalent to the mean flows at the annual, monthly, and weekly scales. The procedure advanced in Linsley et al. [21] was used to calculate the V_R. In turn, the V_R values were compared with the deficit volumes, D_T-o. The calculations were done by writing Macros in Visual Basic and coupling them with associated data in the Microsoft Excel framework. Therefore, flows were standardized at the above three time scales to obtain SHI sequences. In all three time scales, the values of the drought probability, q, were obtained by the counting procedure in which an SHI sequence was chopped at level 0 (mean level). In general, the annual flows tend to follow a normal pdf in the Canadian settings and thus, at the mean level, q values cluster around 0.50 (Table 2, column 2). In view of the gamma pdf of the flow sequences, at the monthly and weekly scales [13] , the q values are significantly larger than 0.50 (Table 2, columns 3 and 4).

Table 2. Calculations of storage volumes at the mean level of flows for varying time scales.

Note: The italicized values in parentheses are calculated at the mean levels (variable means) of the respective months and weeks without standardization of the flow sequences.

For each time scale, the drought magnitudes, M_T-o were computed and accordingly converted to D_T-o. On the annual scale, there is only one set of µ and σ, whereas there are respectively 12 and 52 such sets of µ and σ at the monthly and weekly scales. Therefore, in the calculations of D_T-o at the monthly and weekly scales, an averaged out (arithmetic average) value, denoted by σ_av was used in the analysis. Various versions of σ such as σ_min (subscript min for minimum), σ_max (max for maximum) and the geometric mean of 12 monthly and 52 weekly values were tested, and the arithmetic mean turned out to be the best estimator [12] . The D_T-o values were also estimated without standardization by truncating the flow series at the variable mean levels corresponding to the respective monthly and weekly time scales. In the standardized domain, the variable means are homogenized with a common mean value equal to 0.0.

The MA sequences can be formed from flows or the SHI sequence, alike. However, it is convenient to apply flow sequences to compute the V_R using SPA, whereas the DM based method explicitly requires SHI sequences. When the annual SHI (or flow) sequence is used without involving any moving average operations then such a sequence is designated as moving average 1 (MA1) sequence. In other words, a non-averaged value of SHI (or flow) is essentially the annual SHI (or flow). When consecutive 2 or 3 or annual SHIs (or flows) are averaged out then such a running sequence is termed as MA2 or MA3 sequence. Figure 2 displays MA1, MA2 and MA3 annual SHI sequences with the drought parameters for the South Saskatchewan River. The MA1 sequence (flows) was subjected to analysis to compute the mean (µ), standard deviation (σ) and lag-1 autocorrelation (ρ). The aforesaid statistics were also evaluated for the MA2 and MA3 (flows) sequences and are shown in Table 3. Using the above values of mean and

Figure 2. Redistribution of drought lengths and magnitudes with varying MA smoothing.

Table 3. Summarized V_R (SPA) and D_T (DM method-counting) on the MA smoothed annual SHI Sequences.

Note: *SPA denotes sequent peak algorithm, **V_R reservoir volume (bold letter) closest to SPA based value.

standard deviation, the MA1, MA2 and MA3 flow sequences were converted to respective SHI sequences. In the process of analysis, the number of drought spells (Ns) dropped from the MA1 through MA3 sequences and are presented in column 5 of Table 3. After a few MA smoothing, Ns attained nearly an equilibrium state and thus suggesting no further MA smoothing were warranted. For example, in Table 3, Ns values for MA3 smoothing marginally deviate from MA2 but significantly drop from MA1.

For a comparative analysis on V_R, the counting procedure was applied to the MA1 sequences. The V_R (Table 3, column 11 and in the subsequent text) were computed using SPA for comparison with D_T-o. The counting for D_T-o was done in terms of M_T-o (SHI sequences, Table 3), which were truncated at the level of 0, then converted to D_T-o (=σ × M_T-o). In the MA1 smoothing, there is only one standard deviation, so after calculating the value of M_T-o1, the value of D_T-o1 was obtained using the above relationship by replacing σ with σ₁, i.e. the standard deviation of the MA1 sequence.

When the D_T-o based on the MA1 analysis did not match or were far less than the V_R, and then the analysis was extended to MA2 and at times to MA3 sequences. For the MA2 smoothing, the value of D_T-o (denoted as D_T-o2) was obtained using the corresponding value of the standard deviation (denoted by σ₂) and M_T-o (say M_T-o2), i.e. D_T-o2 = σ₂ × M_T-o2. For further comparison, another value of D_T-o2 was also computed using σ₁ (based on MA1 or the original flow sequence), i.e. D_T-o2 = σ₁ × M_T-o2 (Table 3). Likewise, for the MA3 smoothing, two values of D_T-o3 were obtained; one based on σ₃ (D_T-o3 = σ₃ × M_T-o3; σ₃ being the standard deviation obtained after MA3 smoothing) and another value as D_T-o3 = σ₁ × M_T-o3. The aim is to choose the MA smoothing that will provide the best equivalence of D_T-o to V_R. The above sequence of computations is portrayed in a flow diagram (Figure 3). In the flow diagram, the symbols D_T-o and M_T-o are denoted by D_T and M_T for the sake of brevity and ease of writing and they can also represent D_T-e and M_T-e with the common multiplier σ₁. In the flow diagram ${V^{\prime }}_R$ stands for the standardized value of V_R, i.e. = V_R/σ₁, which corresponds to M_T.

Parallel to the counting procedure, the M_T values (denoted as M_T-e) were obtained by using the analytical procedure (Equations (2) and (3)). The corresponding values of D_T-e were computed (D_T-e = σ₁ × M_T-e). The V_R values were compared with the values of D_T-e. Two values of M_T-e were computed for each situation, i.e. one value using the mean as well as the variance in Equation (2) and the other value by simply using the mean based on Equation (3), thus yielding two values of D_T-e. Both the values were compared with V_R for arriving at an appropriate value of D_T-e for further analysis and use.

4. Results

4.1. Role of the Time Scale on the Estimates of Reservoir Size

To discern the role of the annual, monthly and weekly time scales, V_R (SPA) and

Figure 3. Flow diagram for computing deficit volumes under various options in the DM method.

D_T-o (the counting procedure in the DM method) from the respective flows and SHI sequences were computed and are summarized for 6 rivers in Table 2. It was found that the values of V_R at the aforesaid time scales were fairly close to each other with a tendency to slightly increase (0% to 20% with an average value of 10%) at the monthly and weekly scales compared to the annual scale. The monthly and weekly D_T-o values tended to decrease compared to annual values of D_T-o with an average reduction of 25% (range from 21% to 59%). The other point that emerged from the calculations was that the V_R was greater than D_T-o at all-time scales. For example, at the annual scale, V_R values were found to be larger (ranging from 0% to 200%), with an average of nearly 70%. The reduction in the storage requirement in terms of deficit volumes makes sense because at monthly and weekly scales the actual drought periods are estimated more accurately due to time scale effects and are usually shortened and thus requiring less amount of water to meet the demand. On the contrary, in the SPA based calculations for $V_R$ , the fluctuations at a shorter time scale would be larger requiring greater reservoir volume to damp out such fluctuations to meet the constant demand.

The above numbers displaying the large discrepancies between the SPA and the DM based (MA1) estimates highlight that either the SPA yields excessive values of reservoir volume or the DM method yields too small estimates at the draft level of mean annual flow (MAF, 1µ). At the annual scale, however, when the draft was lowered to 0.90µ or less, the estimates by the SPA and the DM method converged to the same value [20] . In other words, a region with a draft level between 0.90µ and 1µ requires special consideration for the estimation of reservoir volumes by the DM method. The SPA based estimates can be construed as fixed with a little scope to lower them in view of the inherent algorithm imbued in it. But the DM based estimates can be boosted by utilizing the MA procedure to attain parity with SPA based estimates. The SPA has been in vogue since the 1960s [26] and is universally accepted to design the reservoir capacity, therefore, the focus in this study is to arrive at a suitable MA smoothing that should yield D_T comparable to V_R at the draft level of 1µ.

The DM based estimates (italicized, Table 2) at the monthly and weekly scales without standardization were found to be slightly different (mostly smaller) in comparison to the standardization based values (SHI sequences). On average, the standardization based estimates were found about 12% larger than those based on the non-standardized values. This discrepancy can be perceived to arise because of ${\sigma }_{av}$ , which has been taken as the representative value of the standard deviation to convert the magnitude in deficit volume (D_T = σ_av × M_T). Needless to mention that σ_av is an estimator representing 12 values of monthly σ’s and is unlikely to be the best for all situations, but is construed to be a better option compared to other options mentioned earlier. Since standardization is purely statistical in this operation, the pdf of monthly sequences is less likely to play the role in explaining the aforesaid discrepancy. On an annual scale, however, it should be noted that the standard deviation has only one value, thus estimates by both routes turned out to be identical. Therefore, only one value of D_T-o is reported in Table 2. Since these estimates (i.e. without standardization) do not require the use of standard deviation in the calculations and thus can be deemed more accurate. However, the standardization procedure is better amenable to statistical analysis, and estimates based on this approach are more conservative (i.e. higher compared to the non-standardization, Table 2). Based on the foregoing reasoning, the route involving standardization (i.e., the use of SHI sequences) for the estimation of D_T-o was preferred in subsequent analyses and the annual scale is considered as a first choice.

4.2. Comparison of Reservoir Sizes Using the SPA and the DM Based Counting Procedure

In computing the D_T-o, the first step is to choose the right value of σ at each MA smoothing. For example, in the MA2 smoothing, there are two D_T-o: one based on σ₂ (i.e. 'D_T-02 = σ₂ × M_T-o2) and another based on σ₁ (i.e. D_T-02 = σ₁ × M_T-o2). For the MA1 smoothing, there is only one standard deviation and thus D_T-o1 = 'D_T-o1 (Table 3). It is apparent from Table 3 that 'D_T-o values either inconsistently decrease or increase in MA2 and MA3 smoothing; whereas D_T-o values are consistently increasing and hence σ₁ is the crucial parameter to be used for matching to the V_R to arrive at an appropriate MA smoothing. In other words, σ₁ must be used as a multiplier with M_T-o in every MA smoothing for estimation of D_T-o and the role of σ₂ and σ₃ is confined to the standardization of the smoothed MA2 and MA3 flow sequences. It should be noted that for consistency and ease, only 1.0 σ₁ is being used as a multiplier. Other proportion of σ₁ (such as 1.2, 1.1, 0.90 or 0.80) were neither considered nor tested for their efficacy in this study.

In assessing the efficacy of various smoothing, the values of D_T-o and V_R were compared on a 1:1 basis and the performance statistics, viz. the Nash-Sutcliffe efficiency (NSE), and the mean error (MER) were used [27] . To arrive at the above estimates of performance statistics, values of V_R and D_T-o were standardized dividing them by σ₁ (MA1). In other words, a 1:1 comparison was made between D_T-o/σ₁ = M_T-o and V_R/σ₁ (denoted by ${V^{\prime }}_R$ ). In doing so, the wild variation in these entities (i.e. D_T-o and V_R) from small to large rivers were homogenized while rendering them non-dimensional, and thus resulting in sensible estimates of NSE and MER. The efficacy is being tested using NSE and MER [27] as these statistics have been extensively used during the past 50 years and are time tested measures in hydrologic investigations.

Based on aforesaid calculations, it was found that M_T-o for MA1 sequences turned out to be significantly less than ${V^{\prime }}_R$ with a caveat that in a few cases, the values of M_T-o were found to be equal to ${V^{\prime }}_R$ . In other words, the values of the M_T-o compared poorly with ${V^{\prime }}_R$ which is also apparent from an utterly low value of NSE ≈ 24% and MER ≈ −39% (Table 4). In brief, the ${V^{\prime }}_R$ tended to be very conservative (meaning larger), whereas DM based estimates, i.e. M_T-o appeared to

Table 4. Performance statistics for comparison of SPA based V_R with DM based D_T.

be significantly smaller. Since the discrepancies in values of M_T-o and the ${V^{\prime }}_R$ were excessively large, therefore, the MA2 and MA3 smoothing were considered.

Firstly, the MA2 based M_T-o2 values were compared to the ${V^{\prime }}_R$ on a 1:1 basis. It was discovered that with the MA2 smoothing, the matching to ${V^{\prime }}_R$ significantly improved resulting in NSE ≈ 72% and MER ≈ −13% (Table 4). In other words, the underestimation was ameliorated significantly as the value of MER ascended from −39% to −13%. Although the NSE values have improved remarkably, there still existed a scope for improvement in the estimates of the D_T-o because underestimation was endemic as revealed by MER of −13%.

Thus, MA3 smoothing was undertaken (flow chart—Figure 3) and values of M_T-o3 were obtained. The MA3 sequences resulted in the over-estimation of D_T-o values with MER = 17.50% although the value of NSE dropped marginally to 69.36% (Table 4). In short, the MA2 smoothing led to the under-counting whereas the MA3 smoothing led to the over-counting of the D_T-o values with NSE being nearly the same. For comparison with values of V_R, therefore, it was considered reasonable to average out the D_T-o values based on the MA2 and the MA3 smoothing. Such a comparison was made by plotting the average values of the M_T-o2 and M_T-o3 against the values of ${V^{\prime }}_R$ and resulted in a remarkably improved match with NSE ≈ 87% and MER ≈ 2% (Figure 4(A), Table 4). The important point to be noted is that at every smoothing, a new value of M_T-o will emerge, which is multiplied by the MA1 smoothing based value of σ (=σ₁) to arrive at the new estimate of D_T-o.

In the process of moving from the MA1 smoothing to the MA2 smoothing, there has been a considerable reduction in the number of drought spells (column 5, Table 3). Such a reduction suggests that there is a significant increase in the drought length (Figure 2) and in turn, there is also a significant increase in the drought magnitude. In other words, the smoothing procedure led to the amalgamation of smaller drought episodes with the larger ones which resulted in enhanced values of the D_T-o (or M_T-o). Such enhanced values have been found to compare well with V_R (or ${V^{\prime }}_R$ ).

Figure 4. Comparison of SPA based V_R with (A) M_T-o by counting procedure (B) M_T by a hybrid procedure i.e. M_T = bigger between M_T-o and M_T-e.

4.3. Comparison of Reservoir Sizes Using the SPA and DM Based Model

The drought magnitudes (M_T-e) in the standardized domain were estimated using Equations (2) and (3). At the annual scale, the characteristic drought length was found equivalent to extreme drought length, L_T [12] , obtained from the Markov chain based relationship. In the first version, the M_T-e was estimated by involving only the mean of the drought intensity i.e. Equation (3) while in the second version, both the mean and variance of drought intensity were considered i.e. Equation (2) to arrive at estimates of the M_T-e hence estimates of the D_T-e. The calculation for D_T-e was done using σ₁, i.e. D_T-e = σ₁ × M_T-e, which is similar to the case of D_T-o. Because of similarity, the best multiplier was σ₁ for all MA smoothing in the estimation of D_T-e. For example, there are two estimates of D_T-e (viz. D_T-e2 = σ₂ × M_T-e2 and D_T-e2 = σ₁ × M_T-e2) if the MA2 smoothing was conducted, then the appropriate value will be D_T-e2 (= σ₁ × M_T-e2), which, in turn, should be comparable to V_R or the counting based D_T-o2. One can advance similar arguments to the MA3 smoothing. Succinctly, σ₁ is the multiplier for all MA smoothing chosen for estimating D_T-e in the analytical approach as was the case for D_T-o. It was found that the estimation by a simple version involving only the mean of the drought intensity proved too inadequate both in terms of NSE and MER. The calculations showed that values of M_T-e are nearly 50% of V_R with the MA1 smoothing and such an underestimation persisted even with the MA3 smoothing leading to the value of MER = −27% (Table 4). The NSE for the MA3 smoothing was also low with the highest value nearly equal to 57%. Similar was the case for the MA2 smoothing with NSE = 49% and MER = −36% (Table 4).

In view of the abysmal values of the performance statistics by Equation (3), Equation (2) was used to estimate M_T-e and its corresponding D_T-e. The performance statistics turned out to be encouraging. Although, there was a significant underestimation (≈−24%) for the MA1 smoothing, however, the underestimation improved remarkably (MER = 0.75%) with the corresponding NSE = 72% for the MA2 smoothing. A consideration of MA3 smoothing resulted in a significant overestimation of nearly 14% and a slight reduction of NSE to 65%, which suggested that the MA3 smoothing is less meaningful. However, the estimates of M_T-e, based on the MA2 smoothing and the MA3 smoothing were averaged out and the resultant performance statistics improved compared to those of the MA2 smoothing with an acceptable overestimation (7.30%). Likewise, the NSE of 71% was almost equal to 72% that was obtained for the MA2 smoothing. In nutshell, the analytical (model) approach also yielded estimates of M_T-e (or D_T-e) which are in agreement with those of the counting method. However, the estimation procedure proved a bit rigorous as it involved the numerical integration of relevant equations as is reported by Sharma and Panu [12] and the resultant output from which, in turn, became input into Equation (2).

It was observed that the MA2 smoothing resulted in similar values of NSE for both the counted D_T-o as well as the estimated D_T-e (Equation (2)). The counted values of D_T-o were ameliorated by averaging the values obtained from the MA2 and the MA3 smoothing. Such an averaging by the analytical estimates involving the MA2 and MA3 smoothing resulted in little improvement over the counted values. At this point, it was mooted that the MA2 smoothing be preserved and the larger value between D_T-o (M_T-o) and D_T-e (M_T-e) (Table 5) be used as the final estimate of the reservoir volume. For an evaluation of the performance statistics,

Table 5. Summary of the ${V^{\prime }}_R$ , M_T-o and M_T-e based on the hybrid procedure on MA2 smoothed annual SHI sequences for the selected rivers.

Note: *(asterisk) means the value is based on comparing M_T-o and M_T-e (Equation (3)) because of ρ being < 0.42.

viz. NSE and MER, the ${V^{\prime }}_R$ and M_T (larger between M_T-e and M_T-o) were compared on a 1:1 basis (Table 5, Figure 4(B)).

The foregoing selection criteria of the larger estimate between the M_T-e and the M_T-o (named as hybrid procedure) resulted in the value of NSE about 85% with an acceptable level of overestimation (MER = 4.7%). The criteria performed well in a majority of rivers except in cases where ρ was found to be less than 0.42 for the MA2 sequences. In such cases, the value of the M_T-e based on Equation (3) was compared with the estimate of the M_T-o and the bigger value was chosen to represent the reservoir volume. In nutshell, the MA2 smoothing is satisfactory for the evaluation of reservoir volumes at the annual scale and conversely, a little gain is achieved by invoking higher MA smoothing, with the caveat that the bigger one between D_T-e (=σ₁ × M_T-e) and D_T-o (=σ₁ × M_T-o) should be chosen as the estimator of deficit volume to correspond with reservoir volume.

5. Discussion

From the foregoing analysis, it was observed, that the discrepancy between the V_R (SPA) and D_T (DM) is significant at the draft level of mean annual flow when only MA1 smoothing is applied on SHI sequences. The large discrepancies between the SPA and DM based (MA1) estimates are an eye-opener in that either the SPA yields too large values of reservoir volume or the DM method yields too small estimates. No such estimates can be considered to be absolute as each method has its own logistics and limitations. In the case of SPA, the difference between the full reservoir level (reservoir is assumed to be full at the beginning) and the lowest level reached during the sampling period, is taken as the reservoir volume. During this intervening period, several droughts including the longest one may occur and the recovery in the water levels in the reservoir may succeed. Conversely in the DM based method, no such assumption is invoked and thus the total shortfall of water below the long term mean flow in a river during the longest drought period is regarded as the required reservoir volume. The reservoir should be designed to store the above volume of water during the period of excess flow in the river.

It turned out that at the demand level equivalent to the long term mean of the river flow, the V_R values are fairly constant no matter what time scale is chosen. In contrast, the drought magnitude-based methodology resulted in significant discrepancy among these estimates with the annual time scale yielding much higher values compared to those at the monthly and weekly time scales. In such a scenario, one can even be tempted to limit the analysis with the annual flow sequences only as it is trivial and the annual flow data can easily be synthesized or generated. The MA1 based estimates were found to be significantly smaller than the SPA based values but adequately catered for deficiencies arising in the wake of severe droughts over a return period of T years. One may infer that under such a scenario, the SPA based values of V_R are excessive and are not warranted to fulfil the water shortages. Such a reservoir would require a large investment for construction and maintenance with less tangible benefits otherwise called for. Despite the wide divergence between these two estimates of V_R, no value can be said to be absolutely perfect. At this point, there is also a need to examine other methods of estimating the reservoir capacity and compare them with the DM based estimates. These values can be treated as different estimates, which can further be fine-tuned by conducting the analyses at the monthly time scales. For example, in the case of South Saskatchewan River (Table 3), there are 4 estimates of V_R, i.e. 708.07 (SPA); 718.55 (average of D_T-o using MA2 and MA3 with DM method), 630.68 (larger between D_T-o and D_T-e for MA2 using DM method), 533.31(average of D_T-e using MA2 and MA3 with the DM method), and one may consider taking an average value of these 4 estimates resulting into 647.60 as a value for design considerations. Here only 4 estimates are considered and a final design value can be fine-tuned by evaluating other methods such as Behavior analysis, the Hardison Gamma method, the Alexander method, the Dincer method, the Gould Gamma method, and Gould probability method as documented in McMahon and Adeloye [22] . Finally, all the estimates may be averaged out to arrive at the final design value of the reservoir capacity.

In a bid to attain the same D_T values as V_R using the DM based method, the drought lengths and in turn, the magnitudes were amplified by a moving average procedure that resulted in the MA2 and MA3 sequences. The D_T values based on the MA2 sequences tend to undercount whereas those based on the MA3 sequences tend to over count compared to the V_R. On an annual basis, either MA2 or MA3 smoothing can only be conducted because there is no smoothing operation between these two, i.e. there is no integer number between MA2 and MA3 smoothing. Therefore, any further refinement of results stresses that the analysis should be conducted at a shorter time scale (i.e. the monthly scale). There exists an opportunity for a suitable match between the V_R and D_T values, provided MA smoothing such as 3-, 4-, 5-, 6- or higher monthly SHI sequences are utilized [28] . This is an area for further research justifying the use of monthly based analysis in the design of reservoirs.

In the present analysis, a draft at the level of mean annual flow was chosen, with the sole objective of demonstrating the application of the drought magnitude-based methodology for estimating the storage capacity of reservoirs both under independent and dependent (Markovian) river flow conditions. In practice, the majority of rivers worldwide are designed, based on the draft of 75% of the mean annual flow [22] , under such a situation the above-described methodology can be applied. Using the above methodology, Sharma and Panu [19] noted that at such a draft level, the DM method with no moving averaging (MA1) of SHI sequences turn out to be equal to SPA based estimates meaning that no higher-order averaging i.e. MA2 or MA3 is required for evolving storage estimates. This is an important observation, which speaks the worth of the DM based methodology as a viable method in tandem with SPA. The method can be extended to monthly SHI sequences, which yields more accurate estimates of reservoir capacity [19] while considering the higher level of dependence (autocorrelation) and skewness (gamma probability distribution) in monthly flows.

6. Conclusions

The analysis was carried out to compare the V_R and D_T respectively, using the SPA and the DM based method (the counting and the estimation procedures) on the flow data from 15 rivers across Canada. The estimates of V_R at the annual, monthly and weekly scales with the draft set at the mean annual flow level were observed close to each other with a slight tendency to increase at the monthly and weekly scales. On the contrary, estimates of D_T at monthly and weekly scales tended to decrease compared to the annual scale. At all three time scales, the D_T estimates turned out to be smaller than V_R. To ameliorate the D_T to the level of V_R, the counting and estimation (analytical) procedures (in the DM based Method) were applied to flow and SHI sequences on an annual scale. In the estimation procedure, the relationships were built on the extreme number theorem, the truncated normal probability distribution of the drought intensity, the normal distribution of the drought magnitude and a Markov chain based value of extreme drought length. In the counting procedure, the MA2 and MA3 smoothing of SHI sequences found the best parity with the averaged out values of D_T-o that were obtained. Likewise, the analytical procedure yielded similar results (in terms of D_T-e) when applied on SHI sequences resulting from the MA2 and MA3 smoothing.

The estimation of D_T-e was found inadequate when only the mean of the drought intensity was used. The consideration of the mean and the variance of the drought intensity in the estimation procedure turned out to be satisfactory and corroborated the results obtained from the counting procedure. Another finding of the study was that the MA2 smoothing of SHI sequences sufficiently provided the larger value between D_T-o and D_T-e is taken as a counterpart value of the reservoir volume for design purposes. The novel feature of the DM method lies in its ability to assess the reservoir volume without assuming the reservoir being full at the beginning of the analysis as it is the case with SPA. Further, the DM based method is capable of considering the return period and associated risk in the design process of reservoirs. However, the DM method requires flow sequences to be stationary unlike the SPA, which applies to stationary and nonstationary flow sequences alike. It is recommended that the study be extended to the annual and monthly scales at varying draft levels such as 80%, 75%, 70%, 60%, 50% etc. of the mean annual flow, which are largely used where environmental concerns are the overriding factors in the design of reservoirs across the globe.

Acknowledgements

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

NOTES

*Corresponding author.

Conflicts of Interest

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

References

[1]	Stedinger, J.R., Vogel, R.M. and Foufoula-Georgia, E. (1993) Frequency Analysis of Extreme Events. Chapter 18. In Maidment, D.R., Ed., Handbook of Hydrology, McGraw Hill, New York, 18.1-18.66.
[2]	Tallaksen, L.M., Madsen, H. and Clausen, B. (1997) On the Definition and Modeling of Streamflow Drought Duration and Deficit Volume. Hydrological Sciences Journal, 42, 15-33. https://doi.org/10.1080/02626669709492003
[3]	Tallaksen, L.M. and Van Lanen, H.A. (2004) Hydrological Drought: Processes and Estimation Methods for Streamflow and Groundwater, Developments in Water Science Vol. 48. Elsevier, Amsterdam.
[4]	Fleig, A.K., Tallaksen, L.M., Hisdal, H. and Demuth, S. (2006) A Global Evaluation of Streamflow Drought Characteristics. Hydrology and Earth System Science, 10, 535-552. https://doi.org/10.5194/hess-10-535-2006
[5]	Yevjevich, V. (1983) Methods for Determining Statistical Properties of Droughts, In: Yevjevich, V., da Cunha, L. and Vlachos, Eds., Coping With Droughts, Water Resources Publications, Littleton, 22-43.
[6]	Horn, D.R. (1989) Characteristics and Spatial Variability of Droughts in Idaho. ASCE Journal of Irrigation and Drainage Engineering, 115, 111-124. https://doi.org/10.1061/(ASCE)0733-9437(1989)115:1(111)
[7]	Sen, Z. (1980) Statistical Analysis of Hydrological Critical Droughts. ASCE Journal of Hydraulic Engineering, 106, 99-115. https://doi.org/10.1061/JYCEAJ.0005362
[8]	Sen, Z. (2015) Applied Drought Modelling, Prediction and Mitigation. Elsevier Inc., Amsterdam.
[9]	Zelenhasic, E. and Salvai, A. (1987) A Method of Streamflow Drought Analysis. Water Resources Research, 23, 156-158. https://doi.org/10.1029/WR023i001p00156
[10]	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 Management, 131, 383-393. https://doi.org/10.1061/(ASCE)0733-9496(2005)131:5(383)
[11]	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
[12]	Sharma, T.C. and Panu, U.S. (2013) A Semi-Empirical Method for Predicting Hydrological Drought magnitudes in the Canadian Prairies. Hydrological Sciences Journal, 58, 549-569. https://doi.org/10.1080/02626667.2013.772688
[13]	Sharma, T.C. and Panu, U.S. (2014) Modelling of Hydrological Drought Durations and Magnitudes: Experiences on Canadian streamflows. Journal of Hydrology: Regional Studies, 1, 92-106. https://doi.org/10.1016/j.ejrh.2014.06.006
[14]	Shukla, S. and Wood, A.W. (2008) Use of a Standardized Runoff Index for Characterizing Hydrologic Drought. Geophysical Research Letters, 35, L02405. https://doi.org/10.1029/2007GL032487
[15]	Nalbantis, I. and Tsakiris, G. (2009) Assessment of Hydrological Drought Revisited. Water Resources Management, 23, 881-897. https://doi.org/10.1007/s11269-008-9305-1
[16]	McKee, T.B., Doesen, N.J. and Kleist, J. (1993) The Relationship of Drought Frequency and Duration to Time Scales. 8th Conference on Applied Climatology, Anaheim, 17-22 January 1993, 179-184.
[17]	Vicente-Serrano, S.M. and Moreno, L. (2005) Hydrological Response to Different Time Scales of Climatological Drought: An Evaluation of the Standardized Precipitation Index in a Mountainous Mediterranean Basin. Hydrology and Earth System Science, 9, 523-533. https://doi.org/10.5194/hess-9-523-2005
[18]	Barker, K.J., Hannaford, J., Chiverton, A. and Stevenson, C. (2016) From Meteorological to Hydrological Droughts Using Standardized Indicators. Hydrology and Earth System Science, 20, 2483-2505. https://doi.org/10.5194/hess-20-2483-2016
[19]	Sharma, T.C. and Panu, U.S. (2021) Reservoir Sizing at Draft Level of 75% of Mean Annual Flow Using Drought Magnitude Based Method on Canadian Rivers. American Institute of Hydrology Journal Hydrology, 8, 79. https://doi.org/10.3390/hydrology8020079
[20]	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
[21]	Linsley, R.K., Franzini, J.B., Freyburg, D.L. and Tchobanoglous, G. (1992) Water Resources Engineering. 4th Edition, Irwin McGraw-Hill, New York, 192.
[22]	McMahon, T.A. and Adeloye A.J. (2005) Water Resources Yield. Water Resources Publications, Littleton.
[23]	Dracup, J.A., Lee, K.S. and Paulson Jr., E.G. (1980) On the Statistical Characteristics of Drought Events. Water Resources Research, 16, 289-296. https://doi.org/10.1029/WR016i002p00289
[24]	Todorovic, P. and Woolhiser, D.A. (1975) A Stochastic Model of n Day Precipitation. Journal of Applied Meteorology, 14, 125-127. https://doi.org/10.1175/1520-0450(1975)014<0017:ASMODP>2.0.CO;2
[25]	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.
[26]	Thomas, H.A. and Burden, R.P. (1963) Operation Research in Water Quality Management. Division of Engineering and Applied Physics, Harvard University, Cambridge.
[27]	Nash, J.E. and Sutcliffe, J.V. (1970) River Flow Forecasting through Conceptual Models: Part 1. A Discussion of Principles. Journal of Hydrology, 10, 282-289. https://doi.org/10.1016/0022-1694(70)90255-6
[28]	Vicente-Serrano, S.M. (2006) Differences in Spatial Patterns of Drought on Different Time Scales: An Analysis of the Iberian Peninsula. Water Resources Management, 20, 37-60. https://doi.org/10.1007/s11269-006-2974-8