Flood-Frequency Analysis for Dams in Romania

Abstract

Accurately determining the maximum designed water discharges of dams is extremely important, considering the economic costs of carrying out these types of hydrotechnical works and the possible disastrous consequences resulting from their incorrect design. This article describes and applies probability distributions used in hydrology, with some recommended by Romanian legislation standard NP 129-2011. The methods for estimating the parameters presented in this article, as well as the establishment of directions for correlating the normative with international regulations, resulting from the research on many rivers with different characteristics, conducted within the Faculty of Hydrotechnics, were completed with specialized computer applications for applying the normative. In this article, two case studies reflecting this research are presented. The verification of the proposed recommendations, on rivers with hydrographic basins with different physiographic characteristics, confirmed the opportunity to implement rigorous and simple criteria. The presentation of the quantile form of some distributions (especially Pearson III) and of the expressions of moments (central and raw) of high order, as well as the presentation of the frequency factors of each analyzed distribution necessary to calculate the confidence interval, constitute novelties, thus facilitating the ease of use of these distributions.

1. Introduction

According to the Romanian legislation, the water dams are designed and verified at the values of the maximum water discharges with different probabilities of being exceeded [1], depending on the importance class of the dams [2].

Rigorous determination of these maximum flows is of particular importance given the costs and risks associated with these types of hydrotechnical works, especially in the context of climate change. The availability of hydrological data from the last period allows the recalculation of maximum flows for the imposition of constructive/non-constructive measures regarding the increase of the capacity of large water discharges.

For a better prediction of floods and the adoption of appropriate mitigation and protection measures, it is essential to know the hydrological and hydraulic characteristics of watercourses [3].

The current legislation [4] requires a correlation with international recommendations and modern practices, as well as the implementation of computer applications for the normative, to be used by engineers without advanced knowledge of mathematics and computer science.

In Romania, the method of ordinary moments (MOM) is established, but taking into account the fact that hydrometric monitoring is deficient, in most cases, the data series are not older than 25 years, and it is necessary to achieve a regionalization based on the L-moments method [5,6], which is generally a method less influenced by the length of the data series and the existence of extreme values. Another solution is the use of partial series (POT), which is much more laborious to achieve, requiring additional checks of data independence and criteria for establishing threshold values [7,8,9].

The comparative presentation, in this article, of the two methods for estimating the parameters, MOM and L-moments, for certain distributions used in Romania and correlated with international regulations [10,11,12] represents a starting point in achieving a regionalization based on the L-moments method. It should be emphasized that the L-moments method has not been presented, analyzed and applied to specific hydrography in Romania. Switching to the L-moment can be done by linear regressions between the coefficient of variation based on the MOM (C_V) and the coefficient of variation based on L-moments (L-C_V), HC_S and L-C_S, where HC_S is Hazen’s unbiased skewness or C_s adopted depending on the genesis of the flows, and L-C_S is the L-moment’s skewness. The genesis of the flows represents the generating mechanism of the maximum flows given by the spatio-temporal conditions of the precipitations and the physiographic characteristics of the hydrographic basin [13,14,15].

This article analyses Pearson III, the most used three parameters of statistical distribution in Romania (parent distribution), as well as LogPearson and GEV distributions, which are recommended by the Romanian legislation [4]. In addition to these distributions, two distributions commonly used internationally are analyzed and applied, the LogNormal and the Wakeby distributions [5,10,11,12,16,17,18,19]. This article describes the latest methods for estimating distribution parameters using engineering mathematics software, by presenting, for the first time, some expressions of the functions, moments and frequency factors, facilitating the ease of use of these distributions.

Thus, novelty elements, such as the Pearson III expressions of the probability density function using the dgamma built-in function, the cumulative probability function using pgamma and pchisq built-in functions, the quantile function expressed with qgamma and qchisq functions and the form of the L-skewness coefficient expressed with ibeta and pbeta built-in functions; the LogPearson III expressions of the probability density function using the dgamma function, the cumulative probability function using the pgamma function, the quantile function expressed with the qgamma function, the form of the three L-moments and the expression of the frequency factor; the GEV expression for the third L-moment; the LogNormal expressions of the probability density function using the dnorm built-in function and dlnorm function, the cumulative probability function (using pnorm, plnorm and cnorm functions) and the expressions of the three conditions for L-moments using the quantile function and the Wakeby expressions for high order moments and the frequency factor.

All of these novelty elements help hydrology researchers understand and apply the processes behind dedicated softwares, in which they select some options without knowing the mathematics behind them. We consider that the use of these kinds of softwares (not knowing the mathematics behind them) is not beneficial in the long run.

2. Flood Frequency Analysis—Determination of Maximum Flows

In Romania, the series of maximum flows used to estimate the parameters of statistical distributions are those of the annual maximum series (AM). A disadvantage of this series is that the values resulting from the theoretical distributions are conservative and, in some cases, exaggerated, especially in the field of low probabilities, and there are underestimated values for exceeding probabilities higher than 80%. The probability of exceeding 80% has a special importance because, according to Romanian legislation Order no. 326/2007, it represents a criterion for determining the bankfull channel. Currently, according to Romanian legislation Order no. 2.115/2021, the annual probability of exceeding 50% is used, but it is too high for mountain river areas.

The annual maximum series consists of the values of the maximum flows characterizing each hydrological year. The advantage of this analysis is the certainty of data independence. Large datasets with length n ≥ 20 [12,13,14,15] are required for these analyses. Based on these considerations, it is necessary to impose an analysis based on the L-moments method.

The determination of the maximum flows was carried out in stages, according to Figure 1. The verification of the character of outliers, normality and homogeneity were carried out in the data-curation phase.

In the next section, the theoretical distributions used in this article for the calculation of the maximum flows are presented [5,16,17,19]. Only the notation of MOM and L-moments will be used to represent the two methods of estimating the parameters of statistical distributions.

2.1. Pearson III Distribution (PE3)

Pearson III distribution is a particular case of the four-parameter exponential gamma distribution (FPEGD) under the condition of b = 1. The probability density function of FPEGD is [20]:

$$\mathrm{f}\left(\mathrm{x}\right)=\frac{{\mathsf{\beta }}^{\mathsf{\alpha }}{\left(\mathrm{x}-\mathsf{\gamma }\right)}^{\frac{{\alpha }}{\mathrm{b}}-1}}{\mathsf{\beta }\cdot \mathsf{\Gamma }\left(\mathsf{\alpha }\right)}\cdot \exp {\left(-\frac{\mathrm{x}-\mathsf{\gamma }}{\mathsf{\beta }}\right)}^{\frac{1}{\mathrm{b}}}$$

(1)

The probability density function, $\mathrm{f}\left(\mathrm{x}\right)$; the complementary cumulative distribution function, $\mathrm{F}\left(\mathrm{x}\right)$, and quantile function, $\mathrm{x}\left(\mathrm{p}\right)$, of the Pearson III distribution are:

$$\mathrm{f}\left(\mathrm{x}\right)=\frac{{\left(\mathrm{x}-\mathsf{\gamma }\right)}^{\mathsf{\alpha }-1}}{{\mathsf{\beta }}^{\mathsf{\alpha }}\cdot \mathsf{\Gamma }\left(\mathsf{\alpha }\right)}\cdot \exp \left(-\frac{\mathrm{x}-\mathsf{\gamma }}{\mathsf{\beta }}\right)=\frac{1}{\mathsf{\beta }}\cdot \mathrm{d}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(\frac{\mathrm{x}-\mathsf{\gamma }}{\mathsf{\beta }},\mathsf{\alpha }\right)$$

(2)

$$\mathrm{F}\left(\mathrm{x}\right)=1-\frac{1}{\mathsf{\beta }\cdot \mathsf{\Gamma }\left(\mathsf{\alpha }\right)}\cdot {{\displaystyle \underset{\mathsf{\gamma }}{\overset{\mathrm{x}}{\int }}\left(\frac{\mathrm{x}-\mathsf{\gamma }}{\mathsf{\beta }}\right)}}^{\mathsf{\alpha }-1}\cdot \exp \left(-\frac{\mathrm{x}-\mathsf{\gamma }}{\mathsf{\beta }}\right)\mathrm{d}\mathrm{x}=\frac{\mathsf{\Gamma }\left(\mathsf{\alpha },\frac{\mathrm{x}-\mathsf{\gamma }}{\mathsf{\beta }}\right)}{\mathsf{\Gamma }\left(\mathsf{\alpha }\right)}$$

(3)

$$\mathrm{F}\left(\mathrm{x}\right)=1-\mathrm{p}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(\frac{\mathrm{x}-\mathsf{\gamma }}{\mathsf{\beta }},\mathsf{\alpha }\right)=1-\mathrm{p}\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{s}\mathrm{q}\left(2\cdot \frac{\mathrm{x}-\mathsf{\gamma }}{\mathsf{\beta }},2\cdot \mathsf{\alpha }\right)$$

(4)

$$\begin{array}{l}\mathrm{x}\left(\mathrm{p}\right)=\mathsf{\gamma }+\mathsf{\beta }\cdot \mathrm{q}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(1-\mathrm{p},\mathsf{\alpha }\right)=\mathsf{\mu }+\mathsf{\sigma }\cdot \frac{\mathrm{q}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(1-\mathrm{p},\mathsf{\alpha }\right)-\mathsf{\alpha }}{\sqrt{\mathsf{\alpha }}}\\ =\mathsf{\gamma }+\mathsf{\beta }\cdot \mathrm{q}\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{s}\mathrm{q}\left(1-\mathrm{p},2\cdot \mathsf{\alpha }\right)=\mathsf{\mu }+\mathsf{\sigma }\cdot \left(\frac{\mathrm{q}\mathrm{c}\mathrm{h}\mathrm{i}\mathrm{s}\mathrm{q}\left(1-\mathrm{p},2\cdot \mathsf{\alpha }\right)}{2\cdot \sqrt{\mathsf{\alpha }}}-\sqrt{\mathsf{\alpha }}\right)\end{array}$$

(5)

where $\mathsf{\alpha },\mathsf{\beta },\mathsf{\gamma }$ are the shape, the scale and the position parameters, with conditions $\mathsf{\gamma }<\mathrm{x}<\infty$ if $\mathsf{\beta }>0$ or $-\infty <\mathrm{x}<\mathsf{\gamma }$ if $\mathsf{\beta }<0$ and $\mathsf{\alpha }>0$; $\mathsf{\mu },\mathsf{\sigma }$ represent the mean and standard deviation.

Appendix C presents the built-in mathematical functions in Mathcad.

2.2. LogPearson Distribution (LP3)

The probability density function, $\mathrm{f}\left(\mathrm{x}\right)$; the complementary cumulative distribution function, $\mathrm{F}\left(\mathrm{x}\right)$, and the quantile function, $\mathrm{x}\left(\mathrm{p}\right)$, of the LogPearson distribution are:

$$\mathrm{f}\left(\mathrm{x}\right)=\frac{{\left(\ln \left(\mathrm{x}\right)-\mathsf{\gamma }\right)}^{\mathsf{\alpha }-1}}{\mathrm{x}\cdot {\mathsf{\beta }}^{\mathsf{\alpha }}\cdot \mathsf{\Gamma }\left(\mathsf{\alpha }\right)}\cdot \exp \left(-\frac{\ln \left(\mathrm{x}\right)-\mathsf{\gamma }}{\mathsf{\beta }}\right)=\frac{1}{\mathsf{\beta }\cdot \mathrm{x}}\cdot \mathrm{d}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(\frac{\ln \left(\mathrm{x}\right)-\mathsf{\gamma }}{\mathsf{\beta }},\mathsf{\alpha }\right)$$

(6)

$$\mathrm{F}\left(\mathrm{x}\right)=\frac{\mathsf{\Gamma }\left(\mathsf{\alpha },\frac{\ln \left(\mathrm{x}\right)-\mathsf{\gamma }}{\mathsf{\beta }}\right)}{\mathsf{\Gamma }\left(\mathsf{\alpha }\right)}=1-\mathrm{p}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(\frac{\ln \left(\mathrm{x}\right)-\mathsf{\gamma }}{\mathsf{\beta }},\mathsf{\alpha }\right)$$

(7)

$$\mathrm{x}\left(\mathrm{p}\right)=\exp \left(\mathsf{\gamma }+\mathsf{\beta }\cdot \mathrm{q}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(1-\mathrm{p},\mathsf{\alpha }\right)\right)$$

(8)

where $\mathsf{\alpha },\mathsf{\beta },\mathsf{\gamma }$ are the shape, the scale and the position parameters. For case $\mathsf{\beta }>0$, the quantile takes the form of $\mathrm{x}\left(\mathrm{p}\right)=\exp \left(\mathsf{\gamma }+\mathsf{\beta }\cdot \mathrm{q}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(\mathrm{p},\mathsf{\alpha }\right)\right)$.

2.3. Generalized Extreme Value Distribution (GEV)

The probability density function, $\mathrm{f}\left(\mathrm{x}\right)$; the complementary cumulative distribution function, $\mathrm{F}\left(\mathrm{x}\right)$, and the quantile function, $\mathrm{x}\left(\mathrm{p}\right)$, of the GEV distribution are [17]:

$$\mathrm{f}\left(\mathrm{x}\right)={\left(1-\frac{\mathsf{\alpha }}{\mathsf{\beta }}\cdot \left(\mathrm{x}-\mathsf{\gamma }\right)\right)}^{\frac{1}{\mathsf{\alpha }}-1}\cdot \frac{1}{\mathsf{\beta }}\cdot \exp \left(-{\left(1-\frac{\mathsf{\alpha }}{\mathsf{\beta }}\cdot \left(\mathrm{x}-\mathsf{\gamma }\right)\right)}^{\frac{1}{\mathsf{\alpha }}}\right)$$

(9)

$$\mathrm{F}\left(\mathrm{x}\right)=1-\exp \left(-{\left(1-\frac{\mathsf{\alpha }}{\mathsf{\beta }}\cdot \left(\mathrm{x}-\mathsf{\gamma }\right)\right)}^{\frac{1}{\mathsf{\alpha }}}\right)$$

(10)

$$\mathrm{x}\left(\mathrm{p}\right)=\mathsf{\gamma }+\frac{\mathsf{\beta }}{\mathsf{\alpha }}\cdot \left(1-\left(-\ln \left(1-\mathrm{p}\right)\right)\right)$$

(11)

where $\mathsf{\alpha },\mathsf{\beta },\mathsf{\gamma }$ are the shape, the scale and the position parameters; $\mathrm{x}<\left(\mathsf{\gamma }+\frac{\mathsf{\beta }}{\mathsf{\alpha }}\right)$ if $\mathsf{\alpha }>0$, and $\mathrm{x}>\left(\mathsf{\gamma }+\frac{\mathsf{\beta }}{\mathsf{\alpha }}\right)$ if $\mathsf{\alpha }<0$.

2.4. LogNormal Distribution (LN3)

The probability density function, $\mathrm{f}\left(\mathrm{x}\right)$; the complementary cumulative distribution function, $\mathrm{F}\left(\mathrm{x}\right)$, and the quantile function, $\mathrm{x}\left(\mathrm{p}\right)$, of the LN3 distribution are:

$$\begin{array}{l}\mathrm{f}\left(\mathrm{x}\right)=\frac{\exp \left(-\frac{{\left(\ln \left(\mathrm{x}-\mathsf{\gamma }\right)-\mathsf{\alpha }\right)}^2}{2\cdot {\mathsf{\beta }}^2}\right)}{\left(\mathrm{x}-\mathsf{\gamma }\right)\cdot \mathsf{\beta }\cdot \sqrt{2\cdot \mathsf{\pi }}}=\frac{1}{\mathrm{x}-\mathsf{\gamma }}\cdot \mathrm{d}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\left(\ln \left(\mathrm{x}-\mathsf{\gamma }\right),\mathsf{\alpha },\mathsf{\beta }\right)\\ =\mathrm{d}\mathrm{l}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\left(\mathrm{x}-\mathsf{\gamma },\mathsf{\alpha },\mathsf{\beta }\right)\end{array}$$

(12)

$$\begin{array}{l}\mathrm{F}\left(\mathrm{x}\right)=1-\frac{1}{2}\cdot \left(\mathrm{e}\mathrm{r}\mathrm{f}\left(\frac{1}{\sqrt{2}\cdot \mathsf{\beta }}\cdot \left(\ln \left(\mathrm{x}-\mathsf{\gamma }\right)-\mathsf{\alpha }\right)\right)+1\right)=1-\mathrm{p}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\left(\ln \left(\mathrm{x}-\mathsf{\gamma }\right),\mathsf{\alpha },\mathsf{\beta }\right)\\ =1-\mathrm{c}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\left(\frac{\ln \left(\mathrm{x}-\mathsf{\gamma }\right)-\mathsf{\alpha }}{\mathsf{\beta }}\right)=1-\mathrm{p}\mathrm{l}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\left(\mathrm{x}-\mathsf{\gamma },\mathsf{\alpha },\mathsf{\beta }\right)\end{array}$$

(13)

$$\mathrm{x}\left(\mathrm{p}\right)=\mathsf{\gamma }+\exp \left(\mathsf{\alpha }+\mathsf{\beta }\cdot \mathrm{q}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\left(1-\mathrm{p},0,1\right)\right)=\mathsf{\gamma }+\mathrm{q}\mathrm{l}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\left(1-\mathrm{p},\mathsf{\alpha },\mathsf{\beta }\right)$$

(14)

where $\mathsf{\alpha },\mathsf{\beta },\mathsf{\gamma }$ are the shape, the scale and the position parameters; $\mathsf{\mu },\mathsf{\sigma }$ represent the mean and standard deviation. There are also other forms of expression of the functions [16,17,19].

2.5. Four-Parameters Wakeby Distribution (WK4)

The four-parameters Wakeby distribution represents an alternative to the LogNormal distribution. The four-parameters Wakeby distribution has no form for density and cumulative function, being classified as a quantile function.

The quantile function of WK4 distributions is [5,17]:

$$\mathrm{x}\left(\mathrm{p}\right)=\frac{\mathsf{\alpha }}{\mathsf{\beta }}\cdot \left(1-{\mathrm{p}}^{\mathsf{\beta }}\right)-\frac{\mathsf{\gamma }}{\mathsf{\delta }}\cdot \left(1-{\mathrm{p}}^{-\mathsf{\delta }}\right)$$

(15)

where $\mathsf{\alpha },\mathsf{\gamma }$ are the scale parameters, and $\mathsf{\beta },\mathsf{\delta }$ are the shape parameters.

3. Method of Parameters Estimation

In this article, two methods of parameters estimation are studied: the MOM and L-moments methods.

3.1. Method of Ordinary Moments (MOM)

In the MOM, the function parameters are obtained from the formulas expressed for the expected value, variance and skewness. In the case of the four parameters Wakeby distribution, it is necessary to express kurtosis as an additional condition.

The formulas for these indicators can be obtained according to the central moments or depending on the raw moments, as follows [17,21,22]:

$$\mathrm{expected} \mathrm{value}: \mathsf{\mu }={\mathrm{m}}_1={\mathrm{m}}_1^{\prime }$$

(16)

$$\mathrm{variance}: {\mathsf{\sigma }}^2={\mathrm{m}}_2={\mathrm{m}}_2^{\prime }-{\left({\mathrm{m}}_1^{\prime }\right)}^2$$

(17)

$$\mathrm{skewness}: {\mathsf{\gamma }}_1=\frac{{\mathrm{m}}_3}{{\mathrm{m}}_2{}^{3/2}}=\frac{{\mathrm{m}}_3^{\prime }-3\cdot {\mathrm{m}}_2^{\prime }\cdot {\mathrm{m}}_1^{\prime }+2\cdot {\left({\mathrm{m}}_1^{\prime }\right)}^3}{{\left[{\mathrm{m}}_2^{\prime }-{\left({\mathrm{m}}_1^{\prime }\right)}^2\right]}^{3/2}}$$

(18)

$$\mathrm{kurtosis}: {\mathsf{\gamma }}_2=\frac{{\mathrm{m}}_4}{{\mathrm{m}}_2{}^2}=\frac{{\mathrm{m}}_4^{\prime }-4\cdot {\mathrm{m}}_3^{\prime }\cdot {\mathrm{m}}_1^{\prime }+6\cdot {\mathrm{m}}_2^{\prime }\cdot {\left({\mathrm{m}}_1^{\prime }\right)}^2-3\cdot {\left({\mathrm{m}}_1^{\prime }\right)}^4}{{\left[{\mathrm{m}}_2^{\prime }-{\left({\mathrm{m}}_1^{\prime }\right)}^2\right]}^2}$$

(19)

$\mathrm{E}\left[{\mathrm{x}}^{\mathrm{r}}\right]={\mathrm{m}}_{\mathrm{r}}^{\prime }={\displaystyle \underset{-\infty }{\overset{\infty }{\int }}{\mathrm{x}}^{\mathrm{r}}\cdot \mathrm{f}\left(\mathrm{x}\right)\mathrm{d}\mathrm{x}}$, which represents the raw moment of the r order, and ${\mathsf{\gamma }}_1=\mathsf{\xi }\cdot {\mathrm{C}}_{\mathrm{v}}$, in which $\mathsf{\xi }$ is chosen depending on the origin of the flows, namely, 2 if the flows come exclusively from the melting of snow, 3 if the flows have a mixed origin of melting snow and rain and 4 if the flows come exclusively from rain [13,14,15].

3.1.1. Pearson III Distribution

The formulas for expected value, variance and skewness are as follows [10,16,17]:

$$\mathrm{expected} \mathrm{value}: \mathsf{\mu }=\mathsf{\gamma }+\mathsf{\alpha }\cdot \mathsf{\beta }$$

(20)

$$\mathrm{variance}: {\mathsf{\sigma }}^2=\mathsf{\alpha }\cdot {\mathsf{\beta }}^2$$

(21)

$$\mathrm{skewness}: {\mathsf{\gamma }}_1=\frac{2}{\sqrt{\mathsf{\alpha }}}$$

(22)

3.1.2. LogPearson Distribution (LP3)

The formulas for expected value, variance and skewness are [16]:

$$\mathrm{expected} \mathrm{value}: \mathsf{\mu }=\frac{\exp \left(\mathsf{\gamma }\right)}{{\left(1-\mathsf{\beta }\right)}^{\mathsf{\alpha }}}$$

(23)

$$\mathrm{variance}: {\mathsf{\sigma }}^2=\exp \left(2\cdot \mathsf{\gamma }\right)\cdot \left(\frac{1}{{\left(1-2\cdot \mathsf{\beta }\right)}^{\mathsf{\alpha }}}-\frac{1}{{\left(1-\mathsf{\beta }\right)}^{2\cdot \mathsf{\alpha }}}\right)$$

(24)

$$\mathrm{skewness}: {\mathsf{\gamma }}_1=\frac{\frac{1}{{\left(1-3\cdot \mathsf{\beta }\right)}^{\mathsf{\alpha }}}-\frac{3}{{\left(1-2\cdot \mathsf{\beta }\right)}^{\mathsf{\alpha }}\cdot {\left(1-\mathsf{\beta }\right)}^{\mathsf{\alpha }}}+\frac{2}{{\left(1-\mathsf{\beta }\right)}^{3\cdot \mathsf{\alpha }}}}{{\left(\frac{1}{{\left(1-2\cdot \mathsf{\beta }\right)}^{\mathsf{\alpha }}}-\frac{1}{{\left(1-\mathsf{\beta }\right)}^{2\cdot \mathsf{\alpha }}}\right)}^{\frac{3}{2}}}$$

(25)

3.1.3. Generalized Extreme Value Distribution (GEV)

The formulas for expected value, variance and skewness have the following expressions [6,17]:

$$\mathrm{expected} \mathrm{value}: \mathsf{\mu }=\mathsf{\gamma }+\frac{\mathsf{\beta }}{\mathsf{\alpha }}\cdot \left(1-\mathsf{\Gamma }\left(\mathsf{\alpha }+1\right)\right)$$

(26)

$$\mathrm{variance}: {\mathsf{\sigma }}^2=\frac{{\mathsf{\beta }}^2}{{\mathsf{\alpha }}^2}\cdot \left(\mathsf{\Gamma }\left(2\cdot \mathsf{\alpha }+1\right)-\mathsf{\Gamma }{\left(\mathsf{\alpha }+1\right)}^2\right)$$

(27)

skewness:

$${\mathsf{\gamma }}_1=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\mathsf{\alpha }\right)\cdot \frac{3\cdot \mathsf{\Gamma }\left(2\cdot \mathsf{\alpha }+1\right)\cdot \mathsf{\Gamma }\left(\mathsf{\alpha }+1\right)-\mathsf{\Gamma }\left(3\cdot \mathsf{\alpha }+1\right)-2\cdot \mathsf{\Gamma }{\left(\mathsf{\alpha }+1\right)}^3}{\sqrt{{\left(\mathsf{\Gamma }\left(2\cdot \mathsf{\alpha }+1\right)-\mathsf{\Gamma }{\left(\mathsf{\alpha }+1\right)}^2\right)}^3}}$$

(28)

3.1.4. LogNormal Distribution (LN3)

The formulas for expected value, variance and skewness are [6,16,17]:

$$\mathrm{expected} \mathrm{value}: \mathsf{\mu }=\mathsf{\gamma }+\exp \left(\mathsf{\alpha }+\frac{1}{2}\cdot {\mathsf{\beta }}^2\right)$$

(29)

$$\mathrm{variance}: {\mathsf{\sigma }}^2=\exp \left({\mathsf{\beta }}^2+2\cdot \mathsf{\alpha }\right)\cdot \left(\exp \left({\mathsf{\beta }}^2\right)-1\right)$$

(30)

$$\mathrm{skewness}:{\mathsf{\gamma }}_1=3\cdot \sqrt{\exp \left({\mathsf{\beta }}^2\right)-1}+{\left(\sqrt{\exp \left({\mathsf{\beta }}^2\right)-1}\right)}^3$$

(31)

3.1.5. Wakeby Distribution (WK4)

Due to complicated expressions, only the formulas for expected value and variance are presented in this section [17].

The expression for skewness is presented in Appendix B.

The expression for kurtosis is obtained based on the raw moments presented in Appendix B.

The WK4 equations are:

$$\mathrm{expected} \mathrm{value}: \mathsf{\mu }=\frac{\mathsf{\alpha }}{\mathsf{\beta }+1}-\frac{\mathsf{\gamma }}{\mathsf{\delta }-1}$$

(32)

variance:

$$\begin{array}{l}{\mathsf{\sigma }}^2=\frac{2\cdot {\mathsf{\alpha }}^2}{\left(\mathsf{\beta }+1\right)\cdot \left(2\cdot \mathsf{\beta }+1\right)}-{\left(\frac{\mathsf{\alpha }}{\mathsf{\beta }+1}-\frac{\mathsf{\gamma }}{\mathsf{\delta }-1}\right)}^2\\ -\frac{2\cdot \mathsf{\gamma }\cdot \left(2\cdot \mathsf{\alpha }\cdot {\mathsf{\delta }}^2-\left(\left(\mathsf{\beta }+1\right)\cdot \mathsf{\gamma }+2\cdot \mathsf{\alpha }\cdot \mathsf{\beta }+5\cdot \mathsf{\alpha }\right)\cdot \mathsf{\delta }+\left({\mathsf{\beta }}^2+2\cdot \mathsf{\beta }+1\right)\cdot \mathsf{\gamma }+\mathsf{\alpha }\cdot \mathsf{\beta }+2\cdot \mathsf{\alpha }\right)}{\left(\mathsf{\beta }+1\right)\cdot \left(\mathsf{\delta }-\left(\mathsf{\beta }+1\right)\right)\cdot \left(\mathsf{\delta }-1\right)\cdot \left(2\cdot \mathsf{\delta }-1\right)}\end{array}$$

(33)

3.2. Method of L-Moments

Method of L-Moments is based on linear combinations of probability-weighted moments [5,17]. The flow data must be in ascending order. The formulas for sample L-moments are defined in Appendix A.

3.2.1. Pearson III Distribution

The equations necessary to obtain the exact solution for the parameters of the PE3 function are the following:

$${\mathrm{L}}_1=\mathsf{\gamma }+\mathsf{\alpha }\cdot \mathsf{\beta }$$

(34)

$${\mathrm{L}}_2=\frac{\mathsf{\beta }\cdot \mathsf{\Gamma }\left(\mathsf{\alpha }+\frac{1}{2}\right)}{\sqrt{\mathsf{\pi }}\cdot \mathsf{\Gamma }\left(\mathsf{\alpha }\right)}$$

(35)

$${\mathsf{\tau }}_3=3\cdot \left(2\cdot \mathrm{i}\mathrm{b}\mathrm{e}\mathrm{t}\mathrm{a}\left(\frac{1}{3},\mathsf{\alpha },2\cdot \mathsf{\alpha }\right)-1\right) \mathrm{or} {\mathsf{\tau }}_3=3\cdot \left(2\cdot \mathrm{p}\mathrm{b}\mathrm{e}\mathrm{t}\mathrm{a}\left(\frac{1}{3},\mathsf{\alpha },2\cdot \mathsf{\alpha }\right)-1\right)$$

(36)

where $\mathrm{i}\mathrm{b}\mathrm{e}\mathrm{t}\mathrm{a}\left(\mathrm{a},\mathrm{x},\mathrm{y}\right)$ and $\mathrm{p}\mathrm{b}\mathrm{e}\mathrm{t}\mathrm{a}\left(\mathrm{x},{\mathrm{s}}_1,{\mathrm{s}}_2\right)$ are explained in Apendix C.

An approximate solution of the parameters depending on the parameter $\mathsf{\alpha }$ can also be adopted, obtained based on the solution presented by Hosking [5] and RAO [17].

3.2.2. LogPearson Distribution (LP3)

The exact estimation of the parameters is done based on the following system of equations, solvable based on the Gauss–Quadrature method:

$${\mathrm{L}}_1={\displaystyle \underset{0}{\overset{1}{\int }}\exp \left(\mathsf{\gamma }+\mathsf{\beta }\cdot \mathrm{q}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(1-\mathrm{p},\mathsf{\alpha }\right)\right)}\mathrm{d}\mathrm{p}$$

(37)

$$\frac{{\mathrm{L}}_1-{\mathrm{L}}_2}{2}={\displaystyle \underset{0}{\overset{1}{\int }}\exp \left(\mathsf{\gamma }+\mathsf{\beta }\cdot \mathrm{q}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(1-\mathrm{p},\mathsf{\alpha }\right)\right)}\cdot \mathrm{p}\mathrm{d}\mathrm{p}$$

(38)

$$\frac{2\cdot {\mathrm{L}}_1-3\cdot {\mathrm{L}}_2+{\mathrm{L}}_3}{6}={\displaystyle \underset{0}{\overset{1}{\int }}\exp \left(\mathsf{\gamma }+\mathsf{\beta }\cdot \mathrm{q}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(1-\mathrm{p},\mathsf{\alpha }\right)\right)\cdot {\mathrm{p}}^2}\mathrm{d}\mathrm{p}$$

(39)

For $\mathsf{\beta }<0$, the inverse form will be $\mathrm{q}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(\mathrm{p},\mathsf{\alpha }\right)$.

3.2.3. Generalized Extreme Value Distribution (GEV)

The three equations necessary to obtain the exact solution for the parameters of the GEV function are the following:

$${\mathrm{L}}_1=\mathsf{\gamma }-\frac{\mathsf{\beta }}{\mathsf{\alpha }}\cdot \left(\mathsf{\Gamma }\left(1+\mathsf{\alpha }\right)-1\right)$$

(40)

$${\mathrm{L}}_2=\mathsf{\Gamma }\left(\mathsf{\alpha }\right)\cdot \left(1-2^{-\mathsf{\alpha }}\right)\cdot \mathsf{\beta }$$

(41)

$${\mathrm{L}}_3=\mathsf{\Gamma }\left(\mathsf{\alpha }\right)\cdot \left(1-2^{-\mathsf{\alpha }}\right)\cdot \mathsf{\beta }\cdot \left(\frac{\left(1-3^{-\mathsf{\alpha }}\right)\cdot 2}{1-2^{-\mathsf{\alpha }}}-3\right)$$

(42)

An approximate solution of the parameters depending on the parameter $\mathsf{\alpha }$ can also be adopted, obtained based on the solution presented by Hosking [5] and RAO [17].

3.2.4. LogNormal Distribution (LN3)

The exact estimation of the parameters is performed based on the following system of equations:

$${\mathrm{L}}_1={\displaystyle \underset{0}{\overset{1}{\int }}\mathsf{\gamma }+\mathrm{q}\ln \mathrm{o}\mathrm{r}\mathrm{m}\left(1-\mathrm{p},\mathsf{\alpha },\mathsf{\beta }\right)}\mathrm{d}\mathrm{p}$$

(43)

$$\frac{{\mathrm{L}}_1-{\mathrm{L}}_2}{2}={\displaystyle \underset{0}{\overset{1}{\int }}\left(\mathsf{\gamma }+\mathrm{q}\ln \mathrm{o}\mathrm{r}\mathrm{m}\left(1-\mathrm{p},\mathsf{\alpha },\mathsf{\beta }\right)\right)}\cdot \mathrm{p}\mathrm{d}\mathrm{p}$$

(44)

$$\frac{2\cdot {\mathrm{L}}_1-3\cdot {\mathrm{L}}_2+{\mathrm{L}}_3}{6}={\displaystyle \underset{0}{\overset{1}{\int }}\left(\mathsf{\gamma }+\mathrm{q}\ln \mathrm{o}\mathrm{r}\mathrm{m}\left(1-\mathrm{p},\mathsf{\alpha },\mathsf{\beta }\right)\right)\cdot {\mathrm{p}}^2}\mathrm{d}\mathrm{p}$$

(45)

An approximate solution of the parameters depending on the parameter $\mathsf{\beta }$ can also be adopted, obtained based on the solution presented by RAO [17]:

$$\begin{array}{l}\mathsf{\beta }=0.999281\cdot \mathrm{z}-0.006118\cdot {\mathrm{z}}^3+0.000127\cdot {\mathrm{z}}^5\\ \mathrm{u}\mathrm{n}\mathrm{d}\mathrm{e}\hspace{1em}\mathrm{z}=\sqrt{\frac{8}{3}}\cdot \mathrm{q}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}\left(\frac{1+{\mathsf{\tau }}_3}{2},0,1\right)\end{array}$$

(46)

$$\begin{array}{l}\mathsf{\alpha }=\ln \left(\frac{{\mathrm{L}}_2}{\mathrm{e}\mathrm{r}\mathrm{f}\left(0.5\cdot \mathsf{\beta }\right)}\right)-\frac{{\mathsf{\beta }}^2}{2}\\ \end{array}$$

(47)

$$\begin{array}{l}\mathsf{\gamma }={\mathrm{L}}_1-\exp \left(\mathsf{\alpha }+0.5\cdot {\mathsf{\beta }}^2\right)\\ \end{array}$$

(48)

3.2.5. Wakeby Distribution (WK4)

The WK4 parameters are obtained based on the following formulas [5]:

$${\mathrm{L}}_1=\frac{\mathsf{\alpha }}{\mathsf{\beta }+1}-\frac{\mathsf{\gamma }}{\mathsf{\delta }-1}$$

(49)

$${\mathrm{L}}_2=\frac{\mathsf{\alpha }}{\left(\mathsf{\beta }+1\right)\cdot \left(\mathsf{\beta }+2\right)}+\frac{\mathsf{\gamma }}{\left(\mathsf{\delta }-2\right)\cdot \left(\mathsf{\delta }-1\right)}$$

(50)

$${\mathrm{L}}_3=-\frac{\mathsf{\gamma }\cdot \left(\mathsf{\delta }+1\right)}{\left(\mathsf{\delta }-3\right)\cdot \left(\mathsf{\delta }-2\right)\cdot \left(\mathsf{\delta }-1\right)}-\frac{\mathsf{\alpha }\cdot \left(\mathsf{\beta }-1\right)}{\left(\mathsf{\beta }+1\right)\cdot \left(\mathsf{\beta }+2\right)\cdot \left(\mathsf{\beta }+3\right)}$$

(51)

$${\mathrm{L}}_4=-\frac{\mathsf{\gamma }\cdot \left(\mathsf{\delta }+1\right)\cdot \left(\mathsf{\delta }+2\right)}{\left(\mathsf{\delta }-4\right)\cdot \left(\mathsf{\delta }-3\right)\cdot \left(\mathsf{\delta }-2\right)\cdot \left(\mathsf{\delta }-1\right)}+\frac{\mathsf{\alpha }\cdot \left(\mathsf{\beta }-2\right)\cdot \left(\mathsf{\beta }-1\right)}{\left(\mathsf{\beta }+1\right)\cdot \left(\mathsf{\beta }+2\right)\cdot \left(\mathsf{\beta }+3\right)\cdot \left(\mathsf{\beta }+4\right)}$$

(52)

4. Confidence Intervals

In the previous Romanian normatives, the statistical distributions were selected that for probabilities of exceeding less than 0.1%, the values did not exceed ±20% of the value determined by genetic methods [1]. This approach is difficult to apply, because there are no well-founded hydrological syntheses and regionalizations that are valid for this domain of small probabilities. Thus, it is recommended to use a reference distribution, which is scientifically confirmed over time, and to report on it using the confidence interval. The confidence interval for the MOM is defined, in the World Meteorological Organization (WMO 718) [10], for a statistical distribution, as the 90% confidence level. This assumes that the confidence interval is variable depending on the exceedance probability and standard error specific to each statistical distribution. The confidence limits represent the upper and lower bounds of the interval. The confidence interval can be expressed in three ways: with the frequency factor (${\mathrm{K}}_{\mathrm{T}}$), with the standard error of the theoretical distribution or based on the Kite equation [17].

In this article, only the confidence interval based on the frequency factor is analyzed, due to the ease of application in a normative.

$$\mathrm{x}\left(\mathrm{p},\mathsf{\alpha },{\mathrm{K}}_{\mathrm{T}}\right)=\mathsf{\mu }+\left(\frac{{\mathrm{K}}_{\mathrm{T}}\pm \sqrt{{\mathrm{K}}_{\mathrm{T}}{}^2-\left(1-\frac{\mathrm{q}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}{\left(\frac{1-\mathsf{\alpha }}{2},0,1\right)}^2}{2\cdot \left(\mathrm{n}-1\right)}\right)\cdot \left({\mathrm{K}}_{\mathrm{T}}{}^2-\frac{\mathrm{q}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}{\left(\frac{1-\mathsf{\alpha }}{2},0,1\right)}^2}{\mathrm{n}}\right)}}{1-\frac{\mathrm{q}\mathrm{n}\mathrm{o}\mathrm{r}\mathrm{m}{\left(\frac{1-\mathsf{\alpha }}{2},0,1\right)}^2}{2\cdot \left(\mathrm{n}-1\right)}}\right)\cdot \mathsf{\sigma }$$

(53)

where ${\mathrm{K}}_{\mathrm{T}}$ is the frequency factor of the theoretical distribution, and $\mathsf{\alpha }$ is the confidence level; $\mathsf{\mu }$ represents the arithmetic mean, and $\mathsf{\sigma }$ is the standard deviation.

The frequency factor of theoretical distributions analyzed for MOM is:

-

Pearson III:

$${\mathrm{K}}_{\mathrm{T}}=\frac{\mathrm{q}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(1-\mathrm{p},\mathsf{\alpha }\right)-\mathsf{\alpha }}{\sqrt{\mathsf{\alpha }}}$$

(54)

The frequency factor can also be expressed with the qchisq function or, for approximate solutions, using the Kite, Wilson Hilferty of Cornish–Fisher frequency factor [17].

-

LogPearson:

$${\mathrm{K}}_{\mathrm{T}}=\frac{\exp \left(\mathsf{\gamma }+\mathsf{\beta }\cdot \mathrm{q}\mathrm{g}\mathrm{a}\mathrm{m}\mathrm{m}\mathrm{a}\left(1-\mathrm{p},\mathsf{\alpha }\right)\right)-\exp \left(\mathsf{\gamma }\right)\cdot \left(\frac{1}{{\left(1-\mathsf{\beta }\right)}^{\mathsf{\alpha }}}\right)}{\sqrt{\exp \left(2\cdot \mathsf{\gamma }\right)\cdot \left(\frac{1}{{\left(1-2\cdot \mathsf{\beta }\right)}^{\mathsf{\alpha }}}-\frac{1}{{\left(1-\mathsf{\beta }\right)}^{2\cdot \mathsf{\alpha }}}\right)}}$$

(55)

-

Generalized extreme value distribution (GEV):

$${\mathrm{K}}_{\mathrm{T}}=\mathrm{s}\mathrm{i}\mathrm{g}\mathrm{n}\left(\mathsf{\alpha }\right)\cdot \frac{\left(\mathsf{\Gamma }\left(1+\mathsf{\alpha }\right)-{\left(-\ln \left(1-\mathrm{p}\right)\right)}^{\mathsf{\alpha }}\right)}{{\left(\mathsf{\Gamma }\left(1+2\cdot \mathsf{\alpha }\right)-\mathsf{\Gamma }\left(1+\mathsf{\alpha }\right)\right)}^{0.5}}$$

(56)

-

LogNormal distribution (LN3):

$${\mathrm{K}}_{\mathrm{T}}=\frac{\mathrm{q}\ln \mathrm{o}\mathrm{r}\mathrm{m}\left(1-\mathrm{p},\mathsf{\alpha },\mathsf{\beta }\right)-\exp \left(\mathsf{\alpha }+\frac{1}{2}\cdot {\mathsf{\beta }}^2\right)}{\sqrt{\left(\exp \left({\mathsf{\beta }}^2\right)-1\right)\cdot \exp \left({\mathsf{\beta }}^2+2\cdot \mathsf{\alpha }\right)}}$$

(57)

-

Wakeby distribution (WK4):

$${\mathrm{K}}_{\mathrm{T}}=\frac{\mathsf{\alpha }\cdot \left(\frac{1-{\mathrm{p}}^{\mathsf{\beta }}}{\mathsf{\beta }}-\frac{1}{\mathsf{\beta }+1}\right)-\mathsf{\gamma }\cdot \left(\frac{1-{\mathrm{p}}^{-\mathsf{\delta }}}{\mathsf{\delta }}-\frac{1}{\mathsf{\delta }-1}\right)}{\sqrt{\frac{\mathsf{\gamma }\cdot \left\{\begin{array}{l}4\cdot \mathsf{\alpha }\cdot {\mathsf{\delta }}^2-\left[\mathsf{\gamma }\cdot \left(\mathsf{\beta }+1\right)+6\cdot \mathsf{\alpha }\right]\cdot \mathsf{\delta }+\\ \mathsf{\gamma }\cdot \left({\mathsf{\beta }}^2+2\cdot \mathsf{\beta }+1\right)+2\cdot \mathsf{\alpha }\end{array}\right\}}{\left(\mathsf{\beta }+1\right)\cdot \left(\mathsf{\delta }-\mathsf{\beta }-1\right)\cdot {\left(\mathsf{\delta }-1\right)}^2\cdot \left(2\cdot \mathsf{\delta }-1\right)}+\frac{{\mathsf{\alpha }}^2}{{\left(\mathsf{\beta }+1\right)}^2\cdot \left(2\cdot \mathsf{\beta }+1\right)}}}$$

(58)

For situations where the quantile is derived by only two parameters, which generally refer to the mean and the dispersion, a standard error is obtained depending on the C_V. The confidence interval obtained is a simplification, being easy to determine and narrow (small). The confidence interval with the derivation of skewness is very wide for low probabilities. Because, in Romania, the skewness is chosen according to the genesis of the maximum flows, it is recommended to apply the confidence interval that does not take into account the errors given by the skewness. In Romania, the confidence interval (depending on the standard error) can be defined with the Pearson III distribution by choosing the skewness coefficient from regionalization studies regarding the origin of the maximum flows.

5. Case Studies

The case studies presented (Figure 2) consist of determining the maximum flows in two rivers with different geneses of maximum flows. The study on the Ialomita river uses data from the Romanian legislation NP 129/2011 [4], because the influence of the proposed improvements on official data can be observed.

The verification of the proposed solutions was performed on other data series, especially on mountain rivers that have a natural regime of high flows, represented by case study number 2 on the Prigor river.

The distribution functions presented in the NP 129/2011 were used and added to the distributions frequently used in international recommendations [10,11,12,23,24,25,26,27]. All the calculations were performed in Mathcad Prime 7 for the transparency of algorithms and methods.

5.1. Tandarei Section, Ialomita River, Romania

The Ialomita river, code XI, is part of the Danube river basin, being its left tributary. Located in the southern part of Romania, it has a total length of 417 km, a river basin area of 10,350 km² and an average slope of 15%, an average altitude of 327 m and a sinuosity coefficient of 1.88 [28]. The Tandarei section is situated about 20-km upstream of the confluence with the Danube river.

In

Table 1. Observed data from the Tandarei station.

Data	[yr]	1976	1977	1978	1979	1980	1981	1982	1983	1984	1985	1986
Flow	[m3/s]	220	136	152	381	401	273	159	405	346	161	89
Data	[yr]	1987	1988	1989	1990	1991	1992	1993	1994	1995	1996	1997
Flow	[m3/s]	106	308	47.5	65.3	306	85	228	72	224	317	424
Data	[yr]	1998	1999	2000	2001	2002	2003	2004	2005	2006	2007	2008
Flow	[m3/s]	341	180	94.5	270	192	103	251	468	237	249	104

, the data series (observed data) is presented. The data series has a length of 33 records. The mean ($\mathsf{\mu }$), the standard deviation ($\mathsf{\sigma }$), the coefficient of variation (C_v), the skewness (${\mathsf{\gamma }}_1$) and the kurtosis (${\mathsf{\gamma }}_2$) for observed data were 224 m³/s, 118 m³/s, 0.527, 0.327 and −0.926. The sample L-moments ratio L-C_v, L-skew and L-kurt were 0.306, 0.089 and 0.025.

Taking into account the origin of the flows, a coefficient of 2 (${\mathsf{\gamma }}_1=2\cdot {\mathrm{C}}_{\mathrm{v}}=1.054$) was used to establish the skewness for estimating the parameters with the MOM.

Table 2. Parameter estimates by the MOM and L-moments.

Frequency Distribution	MOM				L-Moments
	$\mathsf{\alpha }$	$\mathsf{\beta }$	$\mathsf{\gamma }$	$\mathsf{\delta }$	$\mathsf{\alpha }$	$\mathsf{\beta }$	$\mathsf{\gamma }$	$\mathsf{\delta }$
Pearson III	3.602	62.2	0	-	13.36	33.6	−224.2	-
LogPearson	18.243	−0.131	7.65	-	2.852	−0.412	−6.396	-
GEV	0.015	93.8	171	-	0.129	110	173.2	-
LogNormal	5.801	0.329	−124	-	6.48	0.183	−440.7	-
Wakeby	117	0.063	1187	−9.49	3.06 × 102	0.672	3.8 × 1011	−8.5 × 109

shows the values of the distributions parameters.

Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7 show the frequency curves using the annual maximum flow for the Ialomita river. For plotting positions, the Nguyen formula was used [29].

In

Table 3. Comparison of the estimated values for the usual probability for the Ialomita river.

P	Frequency Distributions
	Pearson III		LogPearson		GEV		LogNormal		Wakeby
	MOM	L	MOM	L	MOM	L	MOM	L	MOM	L
[%]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]
0.01	941	828	909	582	979	765	999	851	948	495
0.1	768	699	757	560	787	675	789	711	785	492
1	585	557	585	511	588	554	586	560	596	475
5	446	443	448	441	443	444	443	443	446	435
40	233	244	232	247	234	244	234	244	230	250
80	123	118	123	111	126	119	126	120	136	104
95	70.1	43.1	72.5	49.1	67.5	43.5	67.4	42.8	54.3	56.3

, the values for the frequency distribution are estimated for the two methods analyzed.

5.2. Prigor Section, Prigor River, Romania

The Prigor river is a tributary of the Nera river, code VI, which is part of the Danube river basin, being its left tributary. Located in the southwestern part of Romania, the Nera river has a total length of 143 km, a river basin area of 1380 km² and an average slope of 0.9%. The Prigor river has a length of 33 km, with an average slope of 22 ‰, a sinuosity coefficient of 1.83, an average altitude of 713 m and a hydrographic basin area of 153 km² [28]. In

Table 4. Observed data from the Prigor station.

Data	[yr]	2000	2001	2002	2003	2004	2005	2006	2007	2008	2009	2010
Flow	[m3/s]	35	19.9	22.1	11.6	80.3	88	51.6	72.2	16.2	42.6	28.5
Data	[yr]	2011	2012	2013	2014	2015	2016	2017	2018	2019	2020
Flow	[m3/s]	12.8	31.2	24.1	52.2	20.2	18.9	6.49	24.9	15	36.6

, the observed data for the Prigor river are presented. The data series has a length of 21 records. The mean ($\mathsf{\mu }$), the standard deviation ($\mathsf{\sigma }$), the coefficient of variation (C_v), the skewness (${\mathsf{\gamma }}_1$) and the kurtosis (${\mathsf{\gamma }}_2$), are 33.9 m³/s, 23.0 m³/s, 0.679, 1.19 and 0.592. The sample L-moments ratio L-C_v, L-skew and L-kurt are 0.370, 0.312 and 0.153.

In

Table 5. Parameter estimates by the MOM and L-moments.

Frequency Distribution	MOM				L-Moments
	$\mathsf{\alpha }$	$\mathsf{\beta }$	$\mathsf{\gamma }$	$\mathsf{\delta }$	$\mathsf{\alpha }$	$\mathsf{\beta }$	$\mathsf{\gamma }$	$\mathsf{\delta }$
Pearson III	0.964	23.4	11.3	-	1.139	23.1	7.5	-
LogPearson	187	−0.047	12.1	-	3461	0.012	−36.6	-
GEV	−0.112	15.1	23.3	-	−0.210	14.3	21.9	-
LogNormal	3.482	0.559	−4.13	-	3.35	0.654	−1.27	-
Wakeby	299	23.7	21.1	0.03	3776	453	26.7	−0.044

, the values of the parameters for each frequency distribution for the two methods of estimation are displayed. Taking into account the mixed origin of the flows, a coefficient of 3 (${\mathsf{\gamma }}_1=3\cdot {\mathrm{C}}_{\mathrm{v}}=2.04$) was used to establish the skewness.

Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12 show the frequency curves using the annual maximum flow for the Prigor river. For plotting positions, the Weibull formula was used.

In

Table 6. Comparison of the estimated values for usual probability for the Prigor river.

P	Frequency Distributions
	Pearson III		LogPearson		GEV		LogNormal		Wakeby
	MOM	L	MOM	L	MOM	L	MOM	L	MOM	L
[%]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[m3/s]
0.01	224	229	249	352	266	424	255	322	226	210
0.1	170	176	177	226	180	243	178	213	171	167
1	117	121	116	132	114	132	115	129	117	119
5	79.8	82.9	77.7	82.6	76.5	80.8	77.4	81.9	79.6	83.2
40	31.8	32.3	33.0	31.8	33.8	32.2	33.3	32.2	32.9	32.3
80	16.1	14.4	16.2	15.2	16.3	15.5	16.2	15.1	16.2	14.3
95	12.3	9.34	9.35	8.85	7.72	7.94	8.84	8.40	12.8	9.70

, the estimated values for the frequency distribution are displayed, for the two methods analyzed.

6. Discussion

For the current situation of hydrology in Romania, it is recommended to estimate the parameters by the MOM, because there are studies, including regionalization, based on the coefficient of variation and skewness. The transition to L-moments or LH-moments [30] can only be done after regionalization studies are conducted based on the L-moment ratio diagrams, which requires considerable but necessary efforts.

In this article, the maximum flows were determined for two rivers that have different origins of their maximum flows, through the MOM and L-moments, using statistical distributions suitable for Romania. The data series of maximum annual flows, according to Bulletin 17C [12], must have a minimum of 20 records. In the case of the Prigor river, the data series has a length of 21 records, and in the case of the Ialomita river, the data series has a length of 33 records.

Pearson III was chosen as the reference statistical distribution for the case studies. This was chosen for two reasons: the long-term use in Romania and the values close to the L-moments method with the MOM. The relative mean error (RME) and the relative absolute error (RAE) criteria [31] were used to compare the results, as well as the framing of the quantile values in the confidence interval of the reference distribution.

$$\mathrm{R}\mathrm{M}\mathrm{E}=\frac{1}{\mathrm{n}}\cdot \sqrt{{{\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}\left(\frac{{\mathrm{x}}_{\mathrm{i}}-\mathrm{x}\left(\mathrm{p}\right)}{{\mathrm{x}}_{\mathrm{i}}}\right)}}^2}$$

(59)

$$\mathrm{R}\mathrm{A}\mathrm{E}=\frac{1}{\mathrm{n}}\cdot {\displaystyle \sum _{\mathrm{i}=1}^{\mathrm{n}}\left|\frac{{\mathrm{x}}_{\mathrm{i}}-\mathrm{x}\left(\mathrm{p}\right)}{{\mathrm{x}}_{\mathrm{i}}}\right|}$$

(60)

where $\mathrm{n},{\mathrm{x}}_{\mathrm{i}},\mathrm{x}\left(\mathrm{p}\right)$ represent the sample size, the observed value and the estimated value for a given probability.

The results for the RME and RAE are presented in

Table 7. RME and RAE score values for the Ialomita river.

Frequency Distribution	RME		RAE
	MOM	L-Moments	MOM	L-Moments
Pearson III	0.0205	0.0338	0.0934	0.1038
LogPearson	0.0208	0.0191	0.0969	0.0643
GEV	0.0227	0.0340	0.0233	0.0183
LogNormal	0.0219	0.0349	0.0990	0.1067
Wakeby	0.0332	0.0103	0.1426	0.0418

and

Table 8. RME and RAE score values for the Prigor river.

Frequency Distribution	RME		RAE
	MOM	L-Moments	MOM	L-Moments
Pearson III	0.0473	0.0262	0.1039	0.0759
LogPearson	0.0280	0.0213	0.0893	0.0673
GEV	0.0237	0.0204	0.0149	0.0175
LogNormal	0.0258	0.0211	0.0900	0.0670
Wakeby	0.0284	0.0318	0.0770	0.1047

.

The validation of some distributions using these two criteria is not recommended, because they are determined by the differences only in the probability domain of the observed values, so it is recommended to approach validation with the criterion of the confidence interval for the MOM. The results of the performance indicators for L-moments were presented to note that their values were very small, but the significance was not valid.

For the Prigor river, where there is certainty of the natural registration of hydrometric data, the results were similar in the case of estimating the parameters with the MOM, with all distributions falling within the confidence interval of the reference distribution. For the L-moments method, all distributions fell within the confidence interval up to the probability of exceeding 0.1%, except for the GEV, which fell within 0.5%. Poor results for low probabilities are too influenced by small annual maximum flows in the case of a short data series [30].

Figure 13 shows the statistical distributions analyzed compared to the chosen reference distribution for the Prigor case study.

For the Ialomita river (Tandarei hydrometric station), the observed data were in the influenced regime, and the reconstruction of the data in the natural regime was with relatively large errors. In general, for large river basins and long river lengths, especially in low-altitude areas, there is also a high degree of natural attenuation of maximum flows. Furthermore, in this case, the results with the MOM were in accordance with NP 129/2011, all falling within the confidence interval of the reference distribution. For the L-moments method, the Wakeby, Log-Pearson and GEV distributions gave different shapes, possibly due to the high degree of attenuation of the maximum flows, with there being a threshold of the maximum values at low probabilities.

Figure 14 shows the statistical distributions analyzed compared to the chosen reference distribution for the Ialomita case study.

In the case of the Prigor river, the L-skewness (τ₃) and L-kurtosis (τ₄) values of the observed data did not differ much from those of the characteristic values of the τ₄–τ₃ variation of the analyzed theoretical distributions (Figure 15). Thus, all analyzed distributions had a similar graphic appearance.

The resulting values for τ₃ and τ₄ for the two rivers are presented in

Table 9. The resulting values for τ₃ and τ₄.

Distributions	Prigor		Ialomita
	τ3	τ4	τ3	τ4
PE3	0.312	0.160	0.089	0.125
LPE3	0.312	0.196	0.089	0.070
GEV	0.312	0.222	0.089	0.124
LN3	0.312	0.199	0.089	0.129
WK4	0.312	0.153	0.089	0.025

.

In the case of the Ialomita river, the τ₃ and τ₄ values of the observed data differed greatly from those of the characteristic values of the τ₄–τ₃ variation of the analyzed theoretical distributions. Thus, for the three-parameter distributions, with the calibration being done only as a function of τ₃, the L-kurtosis τ₄ took values consistent with the variation of τ₄–τ₃ of the theoretical distribution, disregarding the τ₄ of the observed data, because for larger moments, they became unstable, making the resulting values unrealistic [32]. The Wakeby distribution is a distribution that was introduced in the flood frequency analysis in order to fulfil the “separation effect”, described by Matalas et al. 1975 [33], as much as possible, namely, to carry out an analysis so that the maximum flows were not excessively influenced by much of the small flows. The Wakeby distribution separates the right-hand side from the left-hand side of the distribution [33]. The Wakeby distribution has the property of a very thick left-hand tail (high flows) and a right-hand tail that is thick (small flow) enough to decrease the average skew [33], which makes the middle part of the distribution steeper than traditional skewed curves.

In the case of the four-parameter Wakeby distribution, the resulting values of the parameters matched those of a particular case of the distribution, namely, the generalized Pareto [17]. Thus, in the particular case of the Ialomita river dataset, the Wakeby distribution turned into the generalized Pareto distribution, which has the expression of the quantile function, as follows:

$$\mathrm{x}(\mathrm{p})=\mathsf{\gamma }+\frac{\mathsf{\beta }}{\mathsf{\alpha }}\cdot \left(1-{\mathrm{p}}^{\mathsf{\alpha }}\right)$$

(61)

The term on the right of the Wakeby quantile has a constant value of −40.9 m³/s, which, due to the “–“ sign in front of the term, acquires positive values. The value of this term represents the value of the position parameter γ from the Pareto quantile, with the rest of the terms remaining unchanged. It can also be seen on the graph that the value of τ₄ corresponding to τ₃ is that of the particular Wakeby distribution, namely, the generalized Pareto, respecting the τ₄ of the observed data.

In the case of the Log-Pearson distribution, calculated based on the moments obtained from the density function, it had a similar appearance to the Wakeby distribution, but the latter had, in some cases, better results [32]. The τ₄–τ₃ variation in the case of the LogPearson (Figure 16) distribution was varied depending on the asymmetry values in the log space [34].

The graphs between τ₃ and τ₄ also represent a criterion for choosing the best distributions to use in achieving a regionalization based on L-moments [5,17].

7. Conclusions

This article briefly described the most used statistical distributions for the maximum discharge and the latest methods for estimating the distribution parameters using engineering mathematical software.

The use of the Pearson III distribution in Romania must be maintained for the MOM, because the estimation of the parameters was simple and gave good results. A statistical distribution that can be a good alternative to Pearson III is LogNormal (LN3), being the most widely used distribution in Europe for maximum flows [10,24].

The introduction of other statistical distributions should be done after the regionalization studies are conducted, and it should be a simpler alternative, if this proves to be the case.

The method of estimating the parameters must remain the method of ordinary moments, because it has proven to be effective for Romania. The adoption of other estimation methods (L-moments) is useful but requires a long transition period, in which the two methods must be used in parallel.

In the regulations for calculating the maximum flows, it is recommended that statistical distributions should be presented in full, i.e., the functions of density, distribution and quantile, where they exist, as well as the inclusion of computer applications for estimating the parameters of statistical distributions and calculation examples.

The calibration of the three-parameter distributions is easy with the method of ordinary moments. Since asymmetry does not characterize short data strings, it is generally correct, highlighting two methods. A method of correcting the asymmetry according to the genesis of the data was presented in this article, which is based on choosing the asymmetry as a multiple of the coefficient of variation. Another method is to use skewness correction coefficients depending on the relatively short length of the observed data string. A disadvantage of this method is that it is not possible to calibrate the kurtosis function (non-linear function), because it has large differences for a sample compared to that of a population.

The L-moments for the samples were very close to those of the population; thus, no correction of τ₃ and τ₄ was necessary, being thus more robust and less influenced by the effects of sampling variability [5] due to the fact that the L-moments are linear functions. The use of a four-parameter distribution allows calibration as a function of τ₃ and τ₄, which is an advantage in achieving a regionalization based on L-moments.