Uncertainties in Flow-duration-frequency Relationships of High and Low Flow Extremes in Lake Victoria Basin

This paper focuses on uncertainty analysis to aid decision making in applications of statistically modeled flow-duration-frequency (FDF) relationships of both daily high and low flows. The analysis is based on 24 selected catchments in the Lake Victoria basin in Eastern Africa. The FDF relationships were derived for aggregation levels in the range 1–90 days for high flows and 1–365 days for low flows. The validity of the projected FDF quantiles for high return periods T was checked using growth factor curves. Monte Carlo simulations were used to construct confidence intervals CI on both the estimated Ts for given flows and the estimated FDF quantiles for given T. The average bias of the modeled T of high and low flows are for all catchments and Ts up to 25 years lower than 8%. Despite this relatively small average bias in the modeled T, the limits of the CI on the modeled 25-year flows go up to more than 100% for high flows and more than 150% for low flows. The assessed FDF relationships and accompanied uncertainties are useful for various types of risk based water engineering and water management applications related to floods and droughts.


Introduction
In support of water related risk analysis, there is a great need for frequency analysis for both high and low flow extremes for proper management applications related to floods and droughts.Examples of applications include reservoir operations, irrigation, hydropower scheduling, industrial planning, flow control for ecological purposes e.g., compensation flows and dilution flows for improving the quality of water for treatment plants or power generating plants.Proper management of water resources under global climate change and/or anthropogenic influence is an important key to development.It requires an accurate descriptive study of hydrological extremes and their recurrence rates, at the relevant scales, based on long-term time series of observations of rainfall intensities, discharges or water levels.One approach to obtain substantially compressed frequency information on such extremes from a hydrological time series is through extreme value analysis for a range of aggregation levels to constitute Amplitude-Duration-Frequency (ADF) relationships (FDF or IDF, for discharges or rainfall respectively).Aggregation levels are simply durational intervals over which the hydrological values are averaged.Premised on such durations, the conditional relationships are essentially cumulative functions of the amplitude values in the time series [1].
ADF relationships have been presented in a number of studies (see e.g., [2][3][4][5][6][7][8][9][10]).These relationships are very important in water engineering.According to Nhat [11], ADF relationships are among the most commonly used tools in water resources engineering, either for planning, designing and operating of water resource projects, or for various engineering projects against floods.They are used to construct design storms for hydrological modeling applications [1,12].Another application is the calibration and validation of stochastic rainfall generators [13].Several studies have used ADF relationships to assess the impact of climate change and/or variability on hydrological extremes (see e.g., [8,[14][15][16][17][18]).
These studies did, however, not include uncertainty analysis on the calibrations or communication of confidence intervals (CIs) on the ADF relationships.Although there are some studies including [19][20] that take into account uncertainty analysis on flow duration curves, uncertainties in the FDF relationships can, however, be large in data scarce regions, as is shown in this paper for the basin of Lake Victoria in Eastern Africa.It would in such cases be important to take these uncertainties into account in the water resources management or decision making.
According to Ayyub [21], the need to model and analyze uncertainties stems from the awareness that data abundance does not necessarily give us certainty, and sometimes can lead to overwhelmingly confusing situations, and/or a sense of over-confidence leading to an improper information use.The former can be an outcome of the limited capacity of our human mind in some situations to deal with complexity and data abundance whereas the latter can be attributed to a higher order of ignorance, called the ignorance of self-ignorance.Predictive uncertainty, its quantification and its reduction, is a key issue in statistical modeling of a hydrological variable that allows judgment of the degree of confidence in estimated results.It then can be taken into account in decision making in water resources management.Several approaches exist to assess errors in statistical models, in which the Monte Carlo technique is commonly used.Examples of Monte Carlo simulations applied in statistical modeling can be found in many researches (see e.g., [8,[22][23][24][25][26]).Other approaches are the Jackknife method [27] and the typical-value principle [28], which both construct CIs by employing subsample values of a general statistic as the building block.Bootstrapping was introduced by Efron [29] to further strengthen the jackknife method of estimating bias or standard error.Other than probabilistic approaches, the generalized likelihood uncertainty estimation (GLUE) technique of Beven and Binley [30] and Bayesian methods are commonly employed.The aforementioned methods deal with a number of uncertainties as elaborated in [31].However, this study considers sampling uncertainty of quantile estimates by applying parametric bootstrapping.Consequentially, this study is aimed at not only statistically modeling the FDF relationships but also carrying out an in-depth analysis of errors and uncertainty for communication of the research findings to water resources managers and environmental policy makers.The parametric bootstrapping Monte Carlo method is used to establish CIs on the FDF based flow quantiles.Root mean square error and average bias on these quantiles and on the parameters describing the FDF relationships are also estimated to understand uncertainty.

Study Area and Data Series
Lake Victoria is the world's second largest freshwater lake but has a relatively small drainage basin of about 184,000 km 2 , being slightly less than three times the Lake's surface in area.The Lake's basin is situated at an altitude of 1134 m above mean sea level and stretches 355 km from west to east between 31°37' E to 34°53' E and 412 km from north to south between longitudes 00°30' N and 3°12' S. Low-lying parts close to Lake Victoria are characterized by episodes of floods, for instance, downstream of River Nzoia, around Budalang'i [32][33][34].Despite the episodes of flooding, low flows tend to subsequently punctuate the basin's hydrology often for long periods of time.Severe droughts in the Lake Victoria basin can be expected on average about every seven to eight years during the hot dry season from December to February [35].According to Otieno and Awange [36], these droughts affected power production with resulting economic losses.According to Awange et al. [35], Lake Victoria and its environment is more recently under threat from declining water levels, which has had a number of social and economic effects.
This study is based on 24 selected catchments in the Lake Victoria basin.Long discharge time series, preferably above 25 years were used.Figure 1 shows the locations of the discharge measurement stations used in this study and Table 1 shows for these stations the flow record length, coordinates, area of the catchment upstream of each station, and percentage of missing records.

FDF Modeling
The extreme value analysis (EVA) and FDF modeling are based on nearly independent high and low flow extremes extracted from the full time series.Independent high flows are selected using independence criteria based on threshold values for the time between two successive independent flow (F) peaks, the ratio of the minimum flow between the two peaks over the peak value, and the peak height; see Willems [37] for details on the method.Extraction of low flows was carried out by applying the same method but on the inverted flow (1/F) series.Prior to the extraction of the extreme values from the full time series for each of the selected stations, n-day moving averaging window was passed through the series.The aggregation levels considered for high flows were 1, 3, 5, 7, 10, 30, 60, 90 days while for low flows 1, 10, 30, 90, 150, 180, 240 and 365 days were taken.This is the range covered by the relevant water engineering or management applications as agriculture, irrigation, hydropower, domestic water supply, pollution control, etc.The highest aggregation levels of three months considered for high flows and one year for low flows are based on the differences in time scale of the high/low flow related hydrological processes.Peak flows are sudden and result in immediate effects due to the excess of water, whereas low flows are due to progressive low rainfall periods with long term effects such as shortage in water availability.
To come up with the FDF relationships, for the selected range of aggregation levels, EVA was carried out and the suitable extreme value distribution (EVD) selected.To enable an adequate selection of the most optimal threshold level and to avoid systematic over-/underestimation in the tail of the EVD, quantile plots or Q-Q plots were considered.As seen in the principle suggested by Csörgö et al. [38] and Beirlant et al. [39], and used by Taye and Willems [8], Onyutha [9], Willems et al. [40], Willems [41], and Willems [18], the extreme value index γ (or k = −γ) enables identification of the shape of the EVD.This extreme value index describes the tail heaviness of the EVD.It is a parameter in the Generalized Extreme Value (GEV) distribution of Jenkinson [42] or Generalized Pareto Distribution (GPD) of Pickands [43] commonly applied as EVDs.
The class of the GEV distribution or GPD is identified as heavy tail (when γ >0 or k <0), normal tail (when γ = k = 0) and light tail (when γ <0 or k >0).The principle of calibrating the GPD by a weighted linear regression in Q-Q plots, with the primary forms of testing the tail shape of the GPD was adopted in this study.It has been shown in some asymptotic sense that for independent peak flow extremes as used in this study, the conditional distribution of these extremes follows the GPD [43].
Considering x t , α and k as threshold, scale and shape parameters respectively; the cumulative distribution function G(x) of the GPD is given by: ( ) ( ) This distribution is valid for values of x above the threshold x t .
The relationship between the T-year flow F T and the return period T is given by: [ ] When based on the calibrated GPD, or by: When based on the empirical data.In Equations ( 2) and (3), n is the data record length in years; t the number of observed flows above the threshold x t that is considered in the EVD; j the rank of the events (j = 1 for the highest).The relationships between F T and T for the GPD can be given by Equations ( 5) and (6).For the exponential distribution of Equation (1), Equation (3) transfers to a linear relationship between F T and log(T) as in Equation ( 6): where F T0 is the flow at return periodT0, and T0 is equal or higher than the return period n/t of the threshold x t .The quantiles F T are hereafter called the growth factors (Gf T ).
Figure 2 shows examples of such calibrated normal tailed GPDs as linear regression lines in exponential Q-Q plots for station 16.The slopes of the linear calibrated EVD were quantified by the weighting factors proposed by Hill [44].The slope is computed as the difference between two ordered flows divided by the difference between their corresponding ranks.By considering j to be the rank of events, the MSE of the weighted linear regression in the exponential Q-Q plot is: and the slope estimate provided by: Because of high fluctuations which occur in the slope of the Q-Q plots for high thresholds (see Figure 3) stemming from randomness of the available dataset, the slope estimates for these high thresholds have high statistical uncertainty.Instead for very low thresholds the slope estimates might result in pronounced bias because according to Pickands [43] the slope asymptotically converges to a constant one (in the exponential Q-Q plot for normal tails) for higher thresholds.The selection of optimal threshold values x t above which the EVD is calibrated was ensured to be at a point above which the mean squared error (MSE) of the linear regression is minimal, i.e., within nearly horizontal sections in the plot of the slope versus the number of observations above the threshold.In cases where it is possible for MSE to reach local minima at different threshold orders within the range where the  What followed next after carefully selecting, in a consistent way, the optimal thresholds for the different aggregation levels , was the calibration of the parameters of the EVD and analysis of the relationship between the model parameters and the aggregation levels using the formula presented by Willems [45] and used by Taye and Willems [8], and Onyutha [9] as expressed below: ( ) In this formula, A is the area of the catchment upstream of the discharge measuring station considered.The formula is based on scaling properties for the rainfall intensities and consequently of the river discharges.The scaling property indicates that the same EVD is valid for different aggregation levels after application of a scaling factor to the rainfall or discharges values.The scaling factor is different for different aggregation levels.The formula has five parameters: c, w, H, z and a.The last three parameters are called "scaling exponents" in the scaling theory and a specific interpretation can be given to these parameters.The parameter H is called "Hurst-exponent"; while z represents the dynamic scaling exponent; and a, the scaling exponent applied for the aggregation level.
For parameters α and t, threshold discharge q (mean discharge value calculated based on the complete time series) equals 0. After calibrating of parameters of the GPD using Equation ( 9), FDF relationships are derived using Equations ( 10) and ( 11): The parameter-aggregation level relationships, together with the analytical description of the EVD, finally constituted the FDF relationships that are used to estimate high or low flow quantiles as a simultaneous function of different Ts and aggregation levels.

Uncertainty and Error Analysis
The difference between the observed (Q f ) and the FDF-based (Q d ) flow quantiles were used as a measure of bias.Considering i to be the rank of POT events (i = 1 for the highest); and R the number of POT events above optimal threshold event; the mean of values obtained from expression (Q f,I − Q d,i ) as percentage of Q d,i for i = 1 to R is considered the average percentage bias as in Equation (12).The overall differences between Q f,i and Q d,i values for each of the selected catchments were also evaluated in terms of the root mean squared error [RMSE; Equation ( 13)].
, , %] 100 ( ) Since there are no empirical values for T higher than the length of the available flow series, the validity of the flow quantiles higher than that T obtained from the FDF relationships was checked using the at-site EVD based quantiles.
To check for the bias in the theoretical T in comparison with the empirical T, the distribution of the residuals of FDF based versus empirical or Gf T based T values, and the CIs of these residuals were assessed.On assumptions that a small sample size of observed flow extremes was obtained and that the residuals on the modeled T are random and follow a normal distribution, a t-mean test was conducted on the null hypothesis of an unbiased model hence on the bias (mean residual) in FDF based T or flow quantile estimates of each selected catchment.
Probability distributions and/or CIs on the FDF based T or flow quantile can be applied to provide the FDF based estimates with uncertainty measures.The biases in fitting of theoretical T to empirical quantiles, and that of calibrating parameter-aggregation level relationships were considered.The CIs were computed applying the parametric bootstrapping method, based on samples of river flows randomly generated from the EVD using the Monte Carlo method.Random samples of equal size as the dataset of observed flow extremes were randomly generated.This was repeated 1000 times.The CIs were computed by the percentile method after ranking the flow quantiles in the generated samples and picking the 25th and 975th quantiles as the upper and lower limits of the 95% CI respectively.
Figure 4 shows examples of the FDF relationships obtained after compiling the exponential EVD calibration results for river flows at various aggregation levels for station 14.Up to Ts equal to the length of the available time series, empirical quantiles were derived as well.Because the lengths of the available river flow series were all longer than 25 years but less than 100 years, empirical T-year flow quantiles are only shown for curves up to 25 years in Figure 4.For higher T values, due to the randomness involved in the empirical data, the empirical quantiles can be far more inaccurate in comparison with the theoretical quantiles.Differences between the empirical and theoretical quantiles can, for the higher T values, also be explained by the influence of river flooding and the higher observation errors for higher flows.One of the reasons of the latter increase in observation errors for higher flows are due to bias in rating curve extrapolation or the difference between the river discharge and the catchment rainfall-runoff discharge.Another reason is the increasing statistical model uncertainty with increasing T as seen later in this paper.

Evaluation of the FDF Relationships
Figure 5 shows the graphical goodness-of-fit of the flow quantiles after calibration of the EVDs for the daily aggregation level and for T of 5, 10 and 25 years.Considering the full range of aggregation levels, it can be visualized that the theoretical quantiles fit well the empirical ones.This shows that the FDF calibrations are highly acceptable.
Figure 6 shows a representative example of comparison between the at-site EVD (or Gf T ) based curves and the corresponding curve derived from the FDF relationships of station 1.The figure shows increasing deviations between flow quantiles derived from the two types of curves.Given that no bias could be observed in Figure 5, the systematic deviation between the two curves must be attributed to the EVD versus FDF extrapolations for larger T. Results also show that there is increasing uncertainty in flow quantiles and T estimates for T higher than the length of the observed data series.Extrapolation of data measured over relatively short record length introduces large uncertainties [46].This clearly indicates that caution must be exercised in using data of short record length in estimating very high T.
It has been recommended that T should be extrapolated only for values higher than about two or three times the series length due to the uncertainties introduced by the finite sample size [47].

Uncertainty in Return Periods
Figure 7 shows comparison between the theoretical (FDF modeled) and empirical T values.The figure shows increasing statistical modeling uncertainty with increasing T.This again can be explained by the higher uncertainty in EVA for higher T. The alignment of the data points along straight lines parallel to the bisector shows that the deviations between the empirical and theoretical T values are related to the parameter α.The shifts of the data points from the bisector are proportional to the parameter α.The bias as shown in Figure 7 can be explained by the uncertainty in EVD parameter α, which controls the vertical differences between two successive T-year curves on the FDF relationships.The labels of the stations from 1 to 24 are in the order arranged in Table 1.
Figure 8 shows the Monte Carlo simulation results on the FDF modeled T for daily high and low hydrological extremes.It can be seen that the CIs increase in width with increase in T, which obviously is the result of the high uncertainty associated with EVA for very high T.The widths of the CIs were also noted to vary significantly from one catchment to the next.This reflects the difference in the degree of temporal variability in river flows for the different catchments of the study area.
Figure 9 shows the bias and RMSE of the deviation between FDF modeled and empirical T averaged over the entire EVD fitted to daily flows.It can be seen in Figure 9 that for both high and low flows, the average biases for all the selected catchments of the study area are less than 8%.This explains the acceptability of the quantiles from the statistically modeled FDF relationships of the hydrological extremes in the study area.However, as seen from Figure 9, catchments 2, 6 and 7 border each other in the North Western quadrant of the study area (see Figure 1).This indicates the reduced variability of the river flows in this portion of the study area compared to other areas, e.g., the western or eastern quadrant (see stations 5, 15, and 24 in Figure 1).
The differences in the bias and uncertainty in the FDF modeled flow and T values indicate spatial differences across the study area.The catchments with wider CIs in general present higher temporal variability in the flows due to higher differences between low and high flows or stronger differences between short-duration values and longer duration values.This could be explained by the strong variation in the rainfall extreme intensities in the Nile basin as seen from the regional extreme value analysis by Nyeko-Ogiramoi et al. [48].
As can be seen from the probability distributions of the residuals on the FDF modeled T for daily flows in Figure 10, the residuals on the T estimates (FDF based versus empirical) are wider for high flows than for low flows.This was also seen in Figures 8 and 9.It is noted that the zero value is within the CIs on the T residuals for all the catchments indicating that the FDF modeled Ts are unbiased.

Uncertainty in Flow Quantiles
Figure 11 shows the validity of the projected FDF quantiles for very high T checked using Gf T curves.It can be seen that the bias increases for higher T of both high and low flows.Whereas general underestimations are obtained for high flows, the statistical FDF models for low flows are characterized by overestimations.The overall variation of the biases across the catchments (min, max) for high flows are (−24.49%,10.29%), (−29.98%,13.03%), and (−35.14%,22.02%) for T = 25, 100 and 500 years respectively; while the corresponding figures for low flows are (−32.14%,67.41%), (−47.05%,81.23%), and (−59.95%,110.71%) for T = 25, 100 and 500 years respectively.On average, higher biases (for T = 25, 100 and 500) were obtained in FDF models of low flows (11.61%, 19.01%, and 28.08%) than high flows (5.37%, 11.28%, and 16.85%).This suggests that the FDF models for low flows are less capable of capturing the flow variability than those of high flows; this is somewhat contradictory to the conclusion in previous section that the FDF based T estimates are less biased for low flows than high flows.The reason for this is the lower low flow values, which leads to higher relative errors for the same absolute errors.The average difference between CI (lower, upper) limits and empirical quantiles (EQ) as percentages of EQ on daily flows of FDF relationships of high flows are (−51.9%,60.5%), (−61.0%,82.5%) and (−70.7%,116.7%) for T of 5, 10 and 25 years respectively.Correspondingly, for low flows, CIs of (−56.1%,91.6%), (−65.5%,116.3%) and (−77.7%,151.2%) are obtained.The highest differences are found at stations 20 and 16 with 155.41% and 179.08% for high and low flows respectively; and the minimum % are respectively 37.74% and 91.47% at stations 23 and 11 for high and low flows respectively.These differences are expected to be due to the variation in the influence of local climate bringing about uneven wet and dry periods across the study area.
The differences are noted to become narrower and wider with increase in aggregation levels of high flows and low flows respectively (see Figure 12).This explains that the CIs on the FDF based extreme flow quantiles depend on the magnitude of the aggregated river flows.With increase in the aggregation level, the magnitudes of the river flows reduce (for high flows) and increase (for low flows).This can be seen in Figure 12, which shows CIs constructed using Monte Carlo simulations on the FDF curve of T = 5 years.Importantly, as shown in Figure 12, the CI for any selected T-year quantile on FDF relationship can be estimated.

Conclusions
This paper has developed statistically modeled FDF relationships for high and low flow extremes and Monte Carlo based uncertainty estimates on both T-year flow quantiles and return period estimates.It was possible to adequately fit the exponential case of the GPD as the EVD for all the selected study catchments.This was based on the analysis of the tail's shape of the EVD in Q-Q plots, and calibration of the EVD by weighted regression in the exponential Q-Q plot.These results suggest that the GPD family of distributions represented in this case by the exponential distribution is most suited for hydrological extremes in most of the catchments of the Lake Victoria basin.However, no definite generalization can be made for the entire Lake Victoria basin.More studies need to be carried out in other catchments of the region with longer and up-to-date time series.It can also be recommended that instead of MSE in Q-Q plots, other goodness-of-fit measures e.g., Anderson-Darling test, Kolmogorov-Smirnov test, Z-statistics, L-Moment ratio diagrams, etc., be applied.
The average bias on the modeled return periods of high flows and low flows are all less than 8% (averaged over the entire EVD fitted to daily flows).This confirmed the acceptability of the established FDF relationships.Despite this relatively small value for the average bias in the modeled return period, the bias for individual locations and the limits of the 95% CIs on the modeled T-year flows can differ much more from the observed values.These limits go up to 117% for high flows and return periods up to 25 years, and up to 152% for low flows up to the same return period.The 95% CI of the average bias in the T-year flow ranges from −24.49% to +10.29% for high flows and return period of 25 years, and from to 67.41% for low flows.
In addition, the validity of the projected FDF quantiles for return periods higher than the length of the available flow series was checked.When the Gf T curves were taken as the reference, the 95% CI of the average bias widens to (−29.98%, 13.03%) and (−35.14%,22.02%) for high flows and return periods of 100 and 500 years respectively.They widen to (−47.05%, 81.23%) and (−59.95%,110.71%) for low flows.For individual stations, the CIs become even wider.This result shows that if FDF relationships calibrated for data scarce regions are used for estimating projected T-year flows for high return periods the uncertainty is large, hence should not be ignored.Quantification of this uncertainty then becomes important.
The assessed uncertainty will be useful for decision making to use the constructed FDFs for various applications to estimate cumulative volumes of water during drought or flood periods at various aggregation levels or return periods, as clarified before in the introduction section.

Figure 1 .
Figure 1.Locations of the discharge stations; seeTable 1 for details on these stations.

Figure 2 .
Figure 2. (a) Daily high flows; (b) Daily low flows.Symbol (○) shows observations in exponential Q-Q plots; (□) denotes the selected optimal threshold.The regression lines are the calibrated EVDs.
(m³/s)-1   Return Period (years) b) optimal threshold is situated, the thresholds at the different local minima are calibrated to the EVD separately for visual aid of selecting the most suitable one.The examples in Figure2are for the daily flows at station 16.The optimal thresholds are determined as the flow values with threshold ranks t = 100 and 69 [the 100th and 69th highest flow values for high and low (1/F) flow values respectively] as shown in Figure3.A linear tail behavior in the exponential Q-Q plot can be observed towards the higher F or (1/F) values.

Figure 3 .
Figure 3. (a) Daily high flows; (b) Daily low flows.The symbol (♦) shows Hill-type estimation of the slope in the exponential Q-Q plot; (◊) is for mean squared error (MSE) of the Hill-type regression in the exponential Q-Q plot; and (□) represents selected optimal threshold.

Figure 4 .
Figure 4. (a) FDF of high flows; (b) FDF of low flows.The horizontal axis of (a) is logarithmic.

Figure 6 .
Figure 6.(a) Daily high flow quantiles; (b) Daily low flow quantiles.Solid lines are for Gf T curves and dashed lines are for FDF relationships.

Figure 7 .
Figure 7. (a) Daily high flow quantiles; (b) Daily low flow quantiles.The labels of the stations from 1 to 24 are in the order arranged in Table1.

Figure 8 .
Figure 8.(a) and (c) are for daily high flows; (b) and (d) are for daily low flows.(a) and (b) are for station 7; (c) and (d) are for station 6. Round dotted line is the upper limit of the simulated 95% CI.Thin solid line is the lower limit of the CI.Thick solid line is the modeled T. Vertical axis of each graph is logarithmic.

Figure 10 .
Figure 10.(a) and (c) are for daily high flows; (b) and (d) are for daily low flows.(a) and (b) are for station 20; (c) and (d) are for station 22.The blue line (bell-shaped) is the assumed Gaussian probability density function (pdf).The dark line (sigmoid) is the cumulative probability distribution (cdf).(a), (b) and (d) all have the same legend as (c).

Figure 11 .
Figure 11.(a) High flows; (b) Low flows.The labels of the catchments from 1 to 24 are in the order arranged in Table1.
Figure 11.(a) High flows; (b) Low flows.The labels of the catchments from 1 to 24 are in the order arranged in Table1.

Figure 12 .
Figure 12.(a) and (c) are high flow quantiles; (b) and (d) are low flow quantiles.(a) and (b) are for station 16; (c) and (d) are for station 24.(b), (c) and (d) all have the same legend as (a)."T5 Emp." stands for empirical FDF quantiles for T = 5 years.

Table 1
for details on these stations.

Table 1 .
Overview of selected stations and some characteristics.