An unbiased estimator of coefficient of variation of streamflow

Highlights

• Daily streamflow exhibits periodic, skewed, multimodal and intermittent behavior.

• Those factors lead to severe bias and variability in coefficient of variation C, estimates.

• A zero-inflated lognormal monthly mixture, $\mathrm{\Delta }$LN3MM, model mimics observations.

• A $\mathrm{\Delta }$LN3MM estimator leads to reliable and unbiased estimates of C.

• Comparisons reveal an entirely new level of streamflow variability never before witnessed.

Abstract

Given increasing demand for high frequency streamflow series (HFSS) at daily and subdaily time scales there is increasing need for reliable metrics of relative variability for such series. HFSS can exhibit enormous relative variability especially in comparison with low frequency streamflow series formed by aggregation of HFSS. The product moment estimator of the coefficient of variation C, defined as the ratio of sample standard deviation to sample mean, as well as ten other common estimators of C, are shown to provide severely downward biased and highly variable estimates of C for very long records of highly skewed and periodic HFSS particularly for rivers which exhibit zeros. Resorting to the theory of compound distributions, we introduce an estimator of C corresponding to a mixture of monthly zero-inflated lognormal distributions denoted as a delta lognormal monthly mixture $\mathrm{\Delta }$LN3MM model. Through monthly stratification, our $\mathrm{\Delta }$LN3MM model accounts for the seasonality, skewness, multimodality, and the possible intermittency of HFSS. In comparisons among estimators, our $\mathrm{\Delta }$LN3MM based C estimator is shown to yield much more reliable and approximately unbiased estimates of C not only for small samples but also for very large samples (tens of thousands of observations). We document values of C in the range of [0.18, 42,000] with a median of 1.9 and an interquartile range of [1.34, 3.75] for 6807 daily streamflow series across the U.S. from GAGES-II dataset, with the highest values of C occurring in arid and semiarid regions. A multivariate analysis and national contour map reveal that extremely large values of C, never previously documented, tend to occur in arid watersheds with low runoff ratios, which tend to also exhibit a considerable number of zero streamflows.

Introduction

Streamflow variability has a profound impact on nearly every aspect of water resource design, planning and management, and therefore its quantification plays a key role. Reliable and unbiased metrics for streamflow variability are needed in a range of activities including, but not limited to: classification of regional hydrologic homogeneity; goodness-of-fit assessments of hydrologic models; evaluation of impacts of climatic variability and change on hydrologic systems; and in a wide range of research activities which seek to understand the hydroclimatic mechanisms which give rise to hydrologic variability. Recent research has shown that high levels of hydrologic variability play a dominant role and create considerable obstacles relating to our ability to estimate very common hydrologic statistics such as the Pearson correlation coefficient (Barber et al., 2019) and the Nash-Sutcliffe efficiency goodness-of-fit metric (Lamontagne et al., 2020).

Among the existing measures of relative variability, the coefficient of variation, C, introduced by Pearson (1896, pp. 276–277), is now perhaps the most widely used metric. C is defined as $C=\sigma /\mu$ where μ and σ denote the population mean and standard deviation respectively, of the random variable of interest. The index C has been applied in many other research areas including business, engineering, science, medicine, economics, psychology and other social sciences as well as many other fields (Nairy and Rao, 2003, Kelley, 2007, Soliman et al., 2012). Alternatively, the standard deviation σ is often used for comparing the variability of samples, because it has the same units as the random variable of interest. However, when one’s interest is in comparing the relative variability of several samples, each with different mean values, C is preferred to σ because it is nondimensional and thus accounts for differences in the mean of the samples. In finance, the inverse of C is commonly used as a measure of the performance of an investment portfolio (Knight and Satchell, 2005) and is often termed the risk to reward ratio.

We note at the outset, that C should only be computed for data measured on a ratio scale and would have no meaning for summarizing data on an interval scale (Velleman and Wilkinson, 1993). Ratio scale data are equally spaced data which exhibit a zero, whereas interval scale data are equally spaced data without a predefined zero point. Examples of ratio scale variables include streamflow, precipitation, and temperature measured in Kelvin. In contrast, most common temperature scales (e.g., Celsius, Fahrenheit etc.) are interval scales with arbitrary zeros, so C would be different depending on which scale is used. Most statistics including C, are meaningful for ratio data because their interpretation is unchanged when linear transformations are applied to the data. Other statistics such as means, standard deviations, and product moment correlations are meaningful for summarizing data on both ratio and interval scales.

Streamflow variability and skewness are linked. Vargo et al. (2010, Fig. 3) document the theoretical relationship between C and the coefficient of skewness γ for 36 probability distribution functions (pdfs). For two-parameter pdfs a unique relationship usually exists between C and γ (see Fig. 4 and discussion in Vogel and Fennessey, 1993). For example, for positively skewed Gamma and LN2 variables, γ is related to C via the relations $\gamma =2C$, and $\gamma =3C+C^3$, respectively. For more complex pdfs, with more than two parameters, there is usually no unique relationship between C and γ, however there still remains a linkage between the two, given by a two-dimensional region within the plot of C versus γ, as depicted in Fig. 3 of Vargo et al. (2010).

It is widely understood that streamflow observations exhibit both variability and skewness, yet due to numerous factors discussed in Section 1.5, the most common and widely used product moment estimators of the coefficients of variation C and skewness γ, exhibit both severe downward bias and variability, and are generally not to be trusted, even for sample sizes in the tens of thousands (see Vogel and Fennessey, 1993).

Obtaining reliable estimates of C for daily flow series constitutes the central challenge of this study. Initial efforts to obtain reliable estimates of C for daily streamflow series were made by Limbrunner et al., 2000, Vogel et al., 2003 who applied an L-moment estimator of C based on a three-parameter lognormal distribution to flow series at 1571 watersheds across the U.S. Vogel et al. (2003) reported values of C for daily flow series ranging from approximately 0.5 to 10,000 with a median value of 10, and an interquartile range from 3 to 33. We were unable to find examples of such high values of C reported for any other variables, across multiple disciplines, which in part explains the need for new methods introduced here.

In the past, many water resources design, planning and management problems relied on low frequency streamflow series (LFSS) such as annual and monthly series resulting from the temporal aggregation (average) of daily, hourly, or subhourly high frequency streamflow series (HFSS). All such LFSS exhibit much less variability than the HFSS from which they were created. Regardless of the random variable of interest, aggregation leads to a reduction in variability. Given our focus on streamflow variability it is instructive to first consider the impact of aggregation of independent and identically distributed series. The aggregation (i.e. taking the average) of any independent and identically distributed (iid) random variable X over $\mathrm{n}$ intervals leads to a drop in its standard deviation $\sigma$ to $\sigma /\sqrt{\mathrm{n}}$, so that the coefficient of variation of the aggregated variable denoted as $C_n$ is reduced to $C_n=C/\sqrt{n}$. Assume a daily time scale as a reference and denote the corresponding C as $C_1$. Since average annual streamflows are known to exhibit values of $C_{365}$ in the range of [0.2, 1.5] across the conterminous U.S. (see Vogel et al., 1998), one would expect C for daily streamflows to be in the range [3.8, 28.5] resulting from the relationship $C_1=\sqrt{365}{C}_{365}\approx 19C_{365}\to \left[3.8,28.5\right]$. These initial results are only very crude approximations because daily streamflows are neither identically distributed nor independent, and both of these factors affect the estimation of C, a central focus of this study.

Most streamflow statistics attempt to provide a summary of the statistical behavior of streamflow, as distinguished from its physical or deterministic characteristics. Daily streamflow is subject to numerous deterministic characteristics including a periodic component which is dominated by seasonal climatic conditions which can lead to intermittent and ephemeral streamflows and the occurrence of observations equal to zero. HFSS, such as hourly streamflow, may be subject to other forms of periodic behavior such as diurnal variations.

In a recent study closely related to this study, the need to account for periodicity when estimating sample statistics of daily streamflow series was documented by Lamontagne et al. (2020) in their Monte-Carlo experiments which evaluated the sampling variability of estimates of the goodness-of-fit metric termed efficiency E. Fig. 4 of Lamontagne et al. (2020) documents that accounting for the periodic behavior of daily streamflow led to marked reductions in the variability of estimates of E when compared with estimators of E which did not account for streamflow periodicity.

Zero streamflows are defined as streamflow below the measurement threshold, which is approximately 0.01 cfs in the U.S. (Granato et al., 2017). Of the 20,438 U.S. Geological Survey (USGS) river gages evaluated by Granato et al. (2017), 36% of those gages had at least one occurrence of zero streamflow and 2.6% of those gages had more than 297 days per year (or 81.3%) of zero streamflow. According to Levick et al. (2008), ephemeral and intermittent streams make up approximately 59% of all streams in the United States (excluding Alaska), and over 81% in the arid and semi-arid Southwest according to the USGS National Hydrography Dataset. Such streams usually reside in the headwaters or major tributaries of perennial streams in the Southwest.

Since the occurrence of zero daily streamflow is so common, we introduce a model which accommodates their occurrence; such a model that allows for frequent zero-valued observations is known as a zero-inflated model. We document later using both a zero-inflated model and streamflow observations, that the occurrence of zero streamflows leads to considerable increases in C, requiring estimators of C that accommodate their occurrence.

All statistics have both a theoretical and empirical interpretation. For example, the sample mean $\stackrel{-}{x}$ computed from a single sample of length n, is a sample estimate of the true or population mean $\mu$. The theoretical sampling properties of $\stackrel{-}{x}$ are known for any iid variable and can be summarized by the mean and variance of $\stackrel{-}{x}$. Since $E\left[\stackrel{-}{x}\right]=\mu$, $\stackrel{-}{x}$ is said to be an unbiased estimator of $\mu$. Barber et al., 2019, Lamontagne et al., 2020 provided examples of the sampling properties of estimates of the correlation coefficient and the Nash-Sutcliffe efficiency, respectively.

The difference between the population C and the mean of estimates of C is often referred to as ‘sampling bias’ because it results from estimating the true population value by a finite sample. Wallis et al. (1974) first exposed the importance of sampling bias in their seminal paper “Just a moment”. Bias results from the combination of numerous phenomena acting together. The occurrence of zeros, skewness, persistence, and seasonality characterizing HFSS lead to considerable challenges associated with estimation of C. Kirby (1974) also derived an upper bound on the product moment estimator of C which, until this study, was only considered to be important for small samples.

Values of population C exceeding unity correspond to samples with extremely high skewness, leading to enormous downward bias in all ratio estimators commonly used in hydrology. For example, for the case of population C = 10, Vogel and Fennessey (1993) reported downward bias associated with conventional product moment estimator of C of about 40% and 80% for samples of length 10,000 from synthetic iid samples from lognormal and generalized Pareto distributions, respectively. Note that Limbrunner et al., 2000, Vogel et al., 2003 and this study all document that a value of C = 10 for daily streamflows is not uncommon.

It is well known that autocorrelation inflates the sampling variance of most statistics. Vogel et al., 1998, Lombardo et al., 2014 analyzed the effect of the autocorrelation on the sampling properties of various moments and moment ratios, expanding the results of Wallis et al. (1974). One can think of the impact of autocorrelation as decreasing the effective sample size, so that in the limit, the sample size approaches unity as the autocorrelation approaches unity (which is the case for HFSS). Moreover, all ratio estimators are known to exhibit bias, which is induced by the fact that the numerator and denominator are often correlated random variables each with different sampling properties of their own, resulting in bias in estimation of the overall ratio.

All product moment ratio estimators are known to exhibit bias due to outliers, because very large/small observations, which are far away from the sample mean, exert much more influence than the other observations, due to the exponentiation involved in higher-order moments. In such instances, which occur frequently in daily streamflow series, observations do not exert the same weight, and single or few observations can dictate the value of the sample estimates. Vogel and Fennessey (1993, Fig. 3) documented the enormous impact of the largest observation even for very large samples of daily streamflows in the tens of thousands; they also discussed remarkably large downward bias associated with product moment estimates of skewness and by analogy the same would be true for kurtosis and all other higher order moment ratios. For highly skewed bivariate lognormal samples, Lai et al. (1999) concluded that significant upward bias in estimates of the Pearson correlation coefficient exist, and only begins to disappear for sample sizes in the range of 3–4 million observations. Barber et al. (2019) extend the results of Lai et al. (1999) to highly skewed and periodic hydrologic series, showing clearly that the ordinary product moment estimator of the Pearson correlation coefficient should generally be avoided for use with daily and sub daily streamflow series. Thus, there is ample evidence in the literature that new methods are needed to better estimate C for highly skewed hydrologic data, even for very large samples.

The study of the sampling bias associated with commonly used product moment ratio estimators of HFSS has received very little attention, which is surprising when one considers the increasing attention being given to the application of HFSS in water resource management activities. It is now commonplace for hydrologists to model streamflow at sub daily scales (including hourly and even sub hourly scales) for use in flood forecasting as well as for real-time stormwater and water quality management activities and in hillslope hydrology applications. For example, the National Water Model (Office of Water Prediction, 2017) provides hourly streamflow forecasts for any location within the U.S. With increasing focus on challenges relating to big data, combined with increasing access to graphical processing units, supercomputer resources and the internet of things, one can expect to see continuing development of new HFSS modeling and data acquisition approaches which will result in profound challenges associated with estimation of various summary statistics corresponding to zero-inflated, periodic and highly skewed HFSS. Examples of research which address such challenges include this study, as well as the two recent studies by Barber et al., 2019, Lamontagne et al., 2020.

The primary goals of this paper are (1) to develop and compare approximately unbiased estimators of C for use with highly skewed, periodic, and possibly intermittent daily flow series and (2) to apply those estimators to large samples of daily streamflow data to enable a better understanding of the influence of skewness, zeros, periodicity and other physical factors, on the relative variability of streamflow. We begin with a literature review of estimators of C, followed by introduction of zero-inflated and monthly mixture models to deal with periodicity, skewness and non-perennial rivers. Three classes of estimators of C are introduced. The first class of estimators account for the high degree of variability and skewness associated with daily streamflows. The second class of estimators account for both skewness and zero-inflation, while the third class of estimators account for skewness, zero-inflation and the periodicity of daily streamflows. Then, we perform Monte Carlo experiments which compare the behavior of these estimators of C, and finally we apply those estimators to observed daily flow series across the U.S. A multivariate analysis and contour map of estimates of C enable us to summarize some of the hydroclimatic mechanisms which drive daily streamflow variability. Summary and recommendations conclude this study.