Hydrological Drought Frequency Analysis in Water Management Using Univariate Distributions

Abstract

The study of extreme phenomena in hydrology generally involves frequency analysis and a time series analysis. In this article we provide enough mathematics to enable hydrology researchers to apply a wide range of probability distributions in frequency analyses of hydrological drought. The article presents a hydrological drought frequency analysis methodology for the determination of minimum annual flows, annual drought durations and annual deficit volumes for exceedance probabilities common in water management. Eight statistical distributions from different families and with different numbers of parameters are used for the frequency analysis. The method of ordinary moments and the method of linear moments are used to estimate the parameters of the distributions. All the mathematical characteristics necessary for the application of the eight analyzed distributions, for the method of ordinary moments and the method of linear moments, are presented. The performance of the analyzed distributions is evaluated using relative mean error and relative absolute error. For the frequency analysis of the annual minimum flows, only distributions that have a lower bound close to the annual minimum value should be used, a defining characteristic having the asymptotic distributions at this value. A case study of hydrological drought frequency analysis is presented for the Prigor River. We believe that the use of software without the knowledge of the mathematics behind it is not beneficial for researchers in the field of technical hydrology; thus, the dissemination of mathematical methods and models is necessary. All the research was conducted within the Faculty of Hydrotechnics.

1. Introduction

Drought is a complex phenomenon and can have various approaches. Based on Wilhite and Glantz [1], Tallaksen and van Lanen [2], usually different types of droughts are distinguished, such as meteorological drought, soil moisture drought, hydrological drought and socio-economical drought. In general, drought analysis is based on drought indices (Standardized Precipitation Index, Effective Drought Index, Streamflow Drought Index, Surface Water Supply Index, etc.) or drought frequency analysis. Matalas defines drought as an extended period of dryness, and drought analysis should be performed according to magnitude and duration [3]. Mehmet defines drought “as a natural disaster that occurs when precipitation is significantly less than its normal time” [4].

The objective of this article is the frequency analysis of hydrological drought, important in the management of water resources, in the context of economic development, but also of climate change. According to A. F. Van Loon [5], hydrological drought is defined as “a lack of water in the hydrological system, manifesting itself in abnormally low streamflow in rivers and abnormally low levels in lakes, reservoirs, and groundwater”.

In this article we analyze hydrological drought represented by water deficit which is characterized by minimum flows, duration of drought and volume of water deficit. They are defined according to various thresholds generally chosen as fractions of monthly average flows or values from the duration curve of daily average flows [6,7,8,9].

The Romanian methodology regarding the analysis of minimum flows is usually carried out by extracting some data from the duration curve of average daily flows or minimum monthly flows. Another methodology is statistical processing to obtain insurance curves for the minimum annual flows, volumes and durations of the annual deficit [3,10], but it is not yet implemented in Romania, and there is no standard regarding their statistical determination.

International methodology recommends performing statistical analyses using the minimum data of the 1-day, 7-day and 30-day moving-average chronological series of minimum annual flows [3,7,9]. For proper analysis, droughts are also defined by other variables, namely durations and deficit volumes. Thus, a statistical analysis of these variables is required. In water management, durations and annual cumulative volumes are important. The time series of these two variables, resulting from the imposition of thresholds, may contain many zero values. This aspect leads to the need to use the conditional probability model [7,11,12].

In water management, the statistical analysis of the water deficit is important for the management of water resources, including the improvement of the ecological state of rivers. In Romania, until the year 2020, flow requirements of aquatic ecosystems represented the flow with a duration of 95% of the duration curve of average daily flows. The new legislation, H.G. 148/2020 [13], chooses to define periods of drought with the maximum value between the minimum monthly mean discharge with a 95% probability of occurrence and the lowest value among the monthly ecological flows; this leads in many cases to flows much lower than those before 2020. For example, for a multitude of rivers in protected areas (Porumbacu, Sambata, Vistisoara, Uz, …), the values based on the new methodology affect the hydromorphological indicators [14] that characterize the good ecological condition of rivers for most of the year, constituting a real threat to the environment.

The Faculty of Hydrotechnics, following hydrological research, started a procedure of statistical analysis of the minimum flows using a multitude of statistical distributions but also mathematical methods for determining periods of hydrological drought, being part of more extensive research regarding the correction of the ecological flow legislation.

Water deficit analysis was performed for calendar-year periods, as the use of a hydrological-year period is relevant to the water balance period but interrupts drought periods. Seasonal analysis of drought is not appropriate, because the freezing period leads to hydrological droughts smaller in magnitude and less important in water management, because there is not much water consumption, mainly only hydropower being critical. Through the annual summation of drought periods, water deficit analysis proves useful in water management decisions.

Probability distributions with three and five parameters are presented in this article, for the first time, in the statistical analysis of hydrological drought. In other materials [7,8,9,10], the Pearson III, Pearson V, log-normal, GEV, Frechet, Weibull, LogLogistic, etc., distributions are recommended for the frequency analysis of hydrological drought, especially of the minimum annual flows, all of which are three-parameter distributions. Pearson III, Pearson V, log-normal, and Weibull distributions are used in this article. In addition, the three-parameter distributions, Wilson–Hilferty and fatigue lifetime, and the five-parameter Wakeby and lambda distributions are used. Based on [7,8,9,10], in the case of perennial rivers, only statistical distributions that have a positive lower bound should be used for frequency analysis, all being three-parameter distributions and applied with the MOM or maximum likelihood parameter estimation method. Thus, it is required that the frequency analysis be performed using the L-moments parameter estimation method, which is much more stable, and that distributions of at least four parameters be used. It is recommended to use some distributions that still have a lower bound that is not only positive, but also very close to the minimum value recorded.

The methods for estimating the parameters of the presented distributions are MOM and L-moments [15]. The use of some software in hydrology limits the user’s freedom of choice, but, above all, it distances the user from understanding the phenomenon. The elimination of these inconveniences requires the detailed presentation of the distributions used and the equations for estimating the parameters of the distributions for the two methods. Other novel elements are also presented, such as the presentation of some distributions which have received very little attention in hydrology, the expressions of the quantiles of some distributions, parameter approximation relations, and frequency factors for estimating the quantile with MOM and L-moments (Appendix A). The estimation of the parameters leads to the solution of some systems of nonlinear equations, often having the integral operator, so it was necessary to obtain approximation relations of the parameters to some probability distributions. For this, approximation tests were carried out, based on [16], with polynomial, exponential or rational functions. In some cases, the variable function by which the parameter is approximated is a logarithm. The relative approximation errors of the parameters are below 10⁻³. For some distributions, no approximation relations were calculated because the skewness depends on two parameters, resulting in a function of two variables with poor results. For these, numerical methods were used to solve systems of equations.

Thus, all these novelty elements for these distributions presented in

Table 1. Novelty elements.

A. Distribution	New Elements
Fatigue Lifetime	inverse function; approximate relation for the shape parameter estimation with the L-moments method
Wakeby	the first four central moments for parameter estimation with the method of ordinary moments
Lambda	the first four central moments for parameter estimation with the method of ordinary moments
B. Frequency factors	the frequency factors for expressing the inverse functions of the distributions with MOM and L-moments The international materials do not contain references regarding the expression of the inverse function of statistical distributions with the L-moments method.
C. Distribution equivalence relations LCs-Cs	The equivalence relations between the two methods are presented for an easy transition from MOM to L-moments.
D.${\tau }_4-{\tau }_3$relations	The variation is presented for a number of 15 distributions. In the case of the existing relations, presented in other materials, for PE3, GEV, W3, LN3, the relations have been improved, and in the case of the Pearson V, WH3, Chi, Inverse Chi, Pseudo-Weibull, Fréchet, Paralogistic, Inverse Paralogistic distributions, the relations are new.
E. Duration curve	Approximation of duration curves with polynomials with logarithmic argument

will help hydrology researchers to use these distributions easily.

The fatigue lifetime and Wilson–Hilferty distributions are used for the first time in the frequency analysis of hydrological drought.

In the research carried out in the Faculty of Hydrotechnics, other distributions were also analyzed (four-parameter Johnson SU, Kappa, Burr and Kumaraswamy distributions) but the results were unsatisfactory. For transparency and to verify the results obtained, the equations necessary to apply these distributions as well as the results obtained after the frequency analysis based on the case study are presented in Tables S1–S4, Fgure S1 in the supplementary material.

The objective of the article is the presentation of rigorous criteria for the methodology for performing a frequency analysis of hydrological drought, by detailing the data selection criteria, the presentation of all the elements of the distributions, and the methods of estimating their parameters recommended in the frequency analysis of hydrological drought. Since in Romania there is no norm regarding hydrological drought, this article aims to provide the methodological basis aligned with current concerns. Modern water management principles are presented for the hydrological characteristics and particularities of Romania.

The article is organized as follows. The description of methodology, the statistical distributions by presenting the density function, the complementary cumulative function and the quantile function are found in Section 2.1. The presentation of the relations for exact calculation and the approximate relations for determining the parameters of the distributions are found in Section 2.2. Case study of hydrological drought frequency analysis for the Prigor river is in Section 3. Results, discussions and conclusions are in Section 4 and Section 5.

2. Methodology

The data analysis period must be at least the World Meteorological Organization (WMO) reference period [17].

In general, low flow frequency analysis is performed for values of the minimum annual flows (AM), the values being the minimum flows recorded in one day of a year [3,10]. In water management, data on the water deficit is of more interest, being a series of minimum 7-day low flows or minimum 30-day low flows, both of which are moving-average daily annual discharges. The 7-day time period eliminates the effect of daily fluctuations [3], being useful for run-of-river facilities, and the 30-day time period is useful for dams with water reservoirs. The analysis was performed for the three variables, namely minimum 1-day low flows, minimum 7-day low flows, and minimum 30-day low flows. These values are obtained from the chronological series of daily flows observed over the entire analysis period, by selecting minimum 1-day low flows, minimum 7-day-moving-average low flows, and minimum 30-day-moving-average low flows. Figure 1 shows the duration curves used in establishing the thresholds that define hydrological drought.

These values represent the magnitude of the drought, but a drought is defined as an extended period of dryness, so in addition to the magnitude presented, the duration is also necessary [3,6]. The duration is obtained using thresholds that define the drought period. These thresholds are obtained from the duration curve of average daily flows. Usual values of thresolds are 80%, 90% and 95% of the duration curve [6]. With the help of these thresholds, the cumulative annual drought durations and annual volumes can be determined (Figure 2).

The duration curve of average daily flows can be approximated with a polynomial having logarithmic terms of the duration (d) and flow (k = Q/Q_m) as follows:

$$k\left(d\right)={\displaystyle \sum _{i=1}^{n_p+1}{C}_i\cdot \ln {\left(d\right)}^{i-1}}$$

(1)

$$D\left(k\right)=\exp \left({\displaystyle \sum _{j=1}^{n_k+1}{C}_j\cdot k^{j-1}}\right)$$

(2)

If the duration curve is tabulated, cubic spline interpolation can also be used.

The determination of the volumes for the thresholds (80%, 90% and 95%) is achieved by the annual summation of the threshold differences and time series (of daily average flows, 7-day moving average, 30-day moving average). Figure 2 shows how to determine the drought volume for the time series of average daily flows corresponding to the Q_80% threshold.

The duration is determined by the annual summation of the days that have values of the chronological series lower than the established thresholds.

A problem with extreme drought events is the much longer duration of low flow compared to the duration of other extreme values, such as floods. Moreover, if in floods there is a form of the flood that can be simplified (flood hydrograph), in droughts, a valid hydrograph cannot be obtained for these types of events; only recession curves are obtained (lack of precipitation). The complexity of the drought phenomenon compared to other hydrological phenomena is noticeable; it requires various types of analyses.

In this article, only frequency analysis will be analyzed for minimum annual flows, durations and annual deficit volumes. The analysis is carried out according to Figure 3.

For the study of hydrological drought, the time series of average daily flows and annual minimum flows are used. The chronological series of moving averages (7-day moving average, 30-day moving average) are obtained from the time series of daily average flows through moving averaging.

From the time series of average daily flows, the duration curve of average daily flows is obtained by descending ordering. The thresholds for hydrological drought analysis (80%, 90%, 95%) are extracted from the duration curve of average daily flows.

The annual deficit volumes corresponding to the three thresholds are obtained from the time series by daily summation of the volumes lower than the chosen threshold.

The annual drought durations corresponding to the three thresholds are obtained from the time series by summing the days with flows lower than the chosen threshold. These annual data represent the statistical data for frequency analysis.

In the frequency analysis of low flow, when the observed values are lower than 100 values, a correction of the skewness coefficient ($C_s$) is necessary to estimate the parameters with MOM [11,18,19,20,21].

In Romania, for the low flow frequency analysis, the skewness coefficient is established according to the perennial/ephemeral character of the river. The skewness coefficient of the flows is established based on the coefficient of variation of the minimum annual flows by multiplying with a coefficient $C_v$, depending on the perennial/ephemeral character of the river, chosen as follows [22,23,24]:

-

3 for rivers with strong supply from the aquifer;

-

2 for rivers that never dry up and with a significant supply from the aquifer in periods of low water;

-

1.5 for rivers with a deficient supply from the aquifer and with a strong trend of decreasing flows;

-

0–1 in the case of rivers that dry up (small rivers, without water supply).

This method for the calculation of the corrected skewness is an outdated method, needing an update by aligning with modern norms and methodologies. A solution is the use of the L-moments method for estimating the parameters of statistical distributions, knowing that this estimation method is more stable and less influenced by the small lengths of the observed data [15,16,17,18,19,25,26,27].

To estimate the parameters with MOM, for some three-parameter distributions, the L-moments method can be used to obtain an equivalent skewness. This involves the calculation of the shape parameter for each distribution using the L-moments method, and this will substitute the shape parameter from the MOM, obtaining the skewness, for the $C_s/C_v$ ratio. This method preserves the mean and coefficient of variation of the observed data. This skewness estimation method is better than the arbitrary choice, many times, of the flow genesis function.

The use of MOM with skewness determined from L-moments is useful for a period of transition from MOM to L-moments in Romania, because the interpretation of the $C_s/C_v$ ratio is well known. For the application of this method, the equivalence relations of the $C_s/C_v$ ratio according to ${\tau }_3$ are presented in

Table A2. Distribution equivalence relations.

Distribution	Skewness
Generalized Extreme Value	$C_s=sign\left(\alpha \right)\cdot \frac{3\cdot \Gamma \left(2\cdot \alpha +1\right)\cdot \Gamma \left(\alpha +1\right)-\Gamma \left(3\cdot \alpha +1\right)-2\cdot \Gamma {\left(\alpha +1\right)}^3}{\sqrt{{\left(\Gamma \left(2\cdot \alpha +1\right)-\Gamma {\left(\alpha +1\right)}^2\right)}^3}}$ $\alpha =\frac{0.283759107-1.669931462\cdot \left|{\tau }_3\right|}{1+0.441588375\cdot \left|{\tau }_3\right|-0.071007671\cdot {\tau }_3^2+0.015634368\cdot {\left|{\tau }_3\right|}^3}$
Log-Normal	$C_s=3\cdot \sqrt{\exp \left({\beta }^2\right)-1}+{\left(\sqrt{\exp \left({\beta }^2\right)-1}\right)}^3$ if $0<\left|{\tau }_3\right|\le \frac{1}{3}$ $\beta =\frac{0.283759107-1.669931462\cdot \left|{\tau }_3\right|}{1+0.441588375\cdot \left|{\tau }_3\right|-0.071007671\cdot {\tau }_3^2+0.015634368\cdot {\left|{\tau }_3\right|}^3}$ $\frac{1}{3}<\left|{\tau }_3\right|\le \frac{2}{3}$ $\begin{array}{l}\beta =-5.3943247\cdot {10}^{-2}+2.6732827\cdot \left|{\tau }_3\right|-2.9211411\cdot {\tau }_3^2+\\ 7.3388138\cdot {\left|{\tau }_3\right|}^3-8.1997299\cdot {\tau }_3^4+4.3299046\cdot {\left|{\tau }_3\right|}^5\end{array}$ $\frac{2}{3}<\left|{\tau }_3\right|<1$ $\beta =\frac{0.81822527+7.818908\cdot {\tau }_3^4-8.1529453\cdot {\tau }_3^8+0.1554959\cdot {\tau }_3^{12}}{1+2.5188862\cdot {\tau }_3^4-4.9742123\cdot {\tau }_3^8+1.592552\cdot {\tau }_3^{12}}$
Frechet	$C_s=\frac{\Gamma \left(1-\frac{3}{\alpha }\right)+2\cdot \Gamma {\left(1-\frac{1}{\alpha }\right)}^3-3\cdot \Gamma \left(1-\frac{2}{\alpha }+\right)\cdot \Gamma \left(1-\frac{1}{\alpha }\right)}{\sqrt{{\left(\Gamma \left(1-\frac{2}{\alpha }\right)-\Gamma {\left(1-\frac{1}{\alpha }\right)}^2\right)}^3}}$ $\alpha =\frac{1+0.441588375\cdot \left|{\tau }_3\right|-0.071007671\cdot {\tau }_3^2+0.015634368\cdot {\left|{\tau }_3\right|}^3}{-0.283759107+1.669931462\cdot \left|{\tau }_3\right|}$
Weibull	$C_s=\frac{\Gamma \left(1+\frac{3}{\alpha }\right)-3\cdot \Gamma \left(1+\frac{2}{\alpha }\right)\cdot \Gamma \left(1+\frac{1}{\alpha }\right)+2\cdot \Gamma {\left(1+\frac{1}{\alpha }\right)}^3}{{\left(\Gamma \left(1+\frac{2}{\alpha }\right)-\Gamma {\left(1+\frac{1}{\alpha }\right)}^2\right)}^{1.5}}$ $\alpha =\frac{3.528107902-6.294082546\cdot \left|{\tau }_3\right|+2.767652838\cdot {\tau }_3^2}{1+4.599024923\cdot \left|{\tau }_3\right|-7.993601572\cdot {\tau }_3^2+2.423742593\cdot {\left|{\tau }_3\right|}^3}$
Log-Logistic	$C_s==\frac{\Gamma \left(1+\frac{3}{\alpha }\right)\cdot \Gamma \left(1-\frac{3}{\alpha }\right)+2\cdot \Gamma {\left(1+\frac{1}{\alpha }\right)}^3\cdot \Gamma {\left(1-\frac{1}{\alpha }\right)}^3-3\cdot \Gamma \left(1+\frac{2}{\alpha }\right)\cdot \Gamma \left(1-\frac{2}{\alpha }\right)\cdot \Gamma \left(1+\frac{1}{\alpha }\right)\cdot \Gamma \left(1-\frac{1}{\alpha }\right)}{{\left[\Gamma \left(1+\frac{2}{\alpha }\right)\cdot \Gamma \left(1-\frac{2}{\alpha }\right)-\Gamma {\left(1+\frac{1}{\alpha }\right)}^2\cdot \Gamma {\left(1-\frac{1}{\alpha }\right)}^2\right]}^{1.5}}$ $\alpha =\frac{1}{{\tau }_3}$
Pareto	$C_s=\frac{\alpha }{\left|\alpha \right|}\cdot \frac{2\cdot \left(\alpha -1\right)\cdot {\left(2\cdot \alpha +1\right)}^{0.5}}{3\cdot \alpha +1}$ $\alpha =\frac{1-3\cdot {\tau }_3}{{\tau }_3+1}$
Pearson V	$C_s=4\cdot \frac{{\left(\alpha -3\right)}^{0.5}}{\alpha -4}$ $\alpha =\exp \left(\begin{array}{l}0.692740135-0.466055586\cdot \ln \left(\left|{\tau }_3\right|\right)+0.460508385\cdot \ln {\left(\left|{\tau }_3\right|\right)}^2+\\ 0.470172127\cdot \ln {\left(\left|{\tau }_3\right|\right)}^3+0.938027763\cdot \ln {\left(\left|{\tau }_3\right|\right)}^4+\\ 0.805363294\cdot \ln {\left(\left|{\tau }_3\right|\right)}^5+0.359846034\cdot \ln {\left(\left|{\tau }_3\right|\right)}^6+\\ 0.08945247\cdot \ln {\left(\left|{\tau }_3\right|\right)}^7+0.011775528\cdot \ln {\left(\left|{\tau }_3\right|\right)}^8+0.000641664\cdot \ln {\left(\left|{\tau }_3\right|\right)}^9\end{array}\right)$
Wilson–Hilferty	$C_s=\frac{\alpha -\frac{3}{\Gamma {\left(\alpha \right)}^2}\cdot \Gamma \left(\alpha +\frac{2}{3}\right)\cdot \Gamma \left(\alpha +\frac{1}{3}\right)+\frac{2}{\Gamma {\left(\alpha \right)}^3}\cdot \Gamma {\left(\alpha +\frac{1}{3}\right)}^3}{\sqrt{\frac{\Gamma \left(\alpha +\frac{2}{3}\right)}{\Gamma \left(\alpha \right)}-\frac{\Gamma {\left(\alpha +\frac{1}{3}\right)}^2}{\Gamma {\left(\alpha \right)}^2}}\cdot \left(\Gamma \left(\alpha +\frac{2}{3}\right)-\frac{\Gamma {\left(\alpha +\frac{1}{3}\right)}^2}{\Gamma \left(\alpha \right)}\right)}$ if $0<\left|{\tau }_3\right|\le \frac{1}{3}$ $\alpha =\exp \left(\begin{array}{l}-4.16202506-3.261604018\cdot \ln \left({\tau }_3\right)-1.783702334\cdot \ln {\left({\tau }_3\right)}^2-0.770946644\cdot \ln {\left({\tau }_3\right)}^3-\\ 0.221698815\cdot \ln {\left({\tau }_3\right)}^4-0.041191426\cdot \ln {\left({\tau }_3\right)}^5-0.004665295\cdot \ln {\left({\tau }_3\right)}^6-\\ 0.000287712\cdot \ln {\left({\tau }_3\right)}^7-0.000007207\cdot \ln {\left({\tau }_3\right)}^8\end{array}\right)$ $\frac{1}{3}<\left|{\tau }_3\right|\le \frac{2}{3}$ $\alpha =\exp \left(\begin{array}{l}-5.264027693-9.134993593\cdot \ln \left({\tau }_3\right)-15.845477811\cdot \ln {\left({\tau }_3\right)}^2-19.874003352\cdot \ln {\left({\tau }_3\right)}^3-\\ 15.495267042\cdot \ln {\left({\tau }_3\right)}^4-6.752325319\cdot \ln {\left({\tau }_3\right)}^5-1.255615645\cdot \ln {\left({\tau }_3\right)}^6\end{array}\right)$ $\frac{2}{3}<\left|{\tau }_3\right|<1$ $\alpha =\exp \left(\begin{array}{l}-7.712526023-93.06660109\cdot \ln \left({\tau }_3\right)-1.519099681\cdot {10}^3\cdot \ln {\left({\tau }_3\right)}^2-1.602587644\cdot {10}^3\cdot \ln {\left({\tau }_3\right)}^4-\\ 1.041173897\cdot {10}^5\cdot \ln {\left({\tau }_3\right)}^4-4.152482998\cdot {10}^5\cdot \ln {\left({\tau }_3\right)}^5-9.887282\cdot {10}^5\cdot \ln {\left({\tau }_3\right)}^6-\\ 1.287749298\cdot {10}^6\cdot \ln {\left({\tau }_3\right)}^7-7.050767642\cdot {10}^5\cdot \ln {\left({\tau }_3\right)}^8\end{array}\right)$
Pseudo-Weibull	$C_s=\frac{2\cdot \Gamma {\left(\frac{1}{\alpha }\right)}^2\cdot \Gamma \left(\frac{4}{\alpha }\right)+8\cdot \Gamma {\left(\frac{2}{\alpha }\right)}^3-9\cdot \Gamma \left(\frac{1}{\alpha }\right)\cdot \Gamma \left(\frac{2}{\alpha }\right)\cdot \Gamma \left(\frac{3}{\alpha }\right)}{4\cdot {\left(\frac{3}{4}\Gamma \left(\frac{1}{\alpha }\right)\cdot \Gamma \left(\frac{3}{\alpha }\right)-\Gamma {\left(\frac{2}{\alpha }\right)}^2\right)}^{1.5}}$ $\begin{array}{l}\alpha =3.3784105-24.3298763\cdot {\tau }_3+1.3320721\cdot {10}^2\cdot {\tau }_3^2-6.002157\cdot {10}^2\cdot {\tau }_3^3+2.0552631\cdot {10}^3\cdot {\tau }_3^4-\\ 5.037878\cdot {10}^3\cdot {\tau }_3^5+8.5230901\cdot {10}^3\cdot {\tau }_3^6-9.6216631\cdot {10}^3\cdot {\tau }_3^7+6.8802849\cdot {10}^3\cdot {\tau }_3^8-\\ 2.8078183\cdot {10}^3\cdot {\tau }_3^9+4.9671387\cdot {10}^2\cdot {\tau }_3^{10}\end{array}$
K3— Generalized Gumbel	$C_s={\displaystyle \underset{0}{\overset{1}{\int }}\frac{{\left[\ln \left(\frac{\alpha }{1-{\left(1-p\right)}^{\alpha }}\right)-\psi \left(\frac{1}{\alpha }+1\right)-\ln \left(\alpha \right)-{\gamma }_e\right]}^3}{{\left[{\displaystyle \underset{0}{\overset{1}{\int }}ln{\left(\frac{\alpha }{1-{\left(1-p\right)}^{\alpha }}\right)}^{2}dp-{\left(\psi \left(\frac{1}{\alpha }+1\right)+\ln \left(\alpha \right)+{\gamma }_e\right)}^2}\right]}^{1.5}}dp}$ if $0<{\tau }_3\le 1/3$ $\alpha =-971.10886182\cdot {10}^{-3}+6.3606944\cdot {\tau }_3-6.35405524\cdot {\tau }_3^2+15.03572719\cdot {\tau }_3^3$ $1/3<{\tau }_3<1$ $\alpha =\frac{-661.76303922\cdot {10}^{-3}+4.30832382\cdot {\tau }_3-1.09504916\cdot {\tau }_3^2+373.12384033\cdot {10}^{-3}\cdot {\tau }_3^3}{1-999.98952576\cdot {10}^{-3}\cdot {\tau }_3}$
K3—Park	$C_s=\frac{\begin{array}{l}2\cdot \Gamma {\left(\frac{2}{\alpha }\right)}^3\cdot {\alpha }^{\frac{3}{\alpha }}\cdot \Gamma {\left(1-\frac{1}{\alpha }\right)}^3+\Gamma {\left(\frac{1}{\alpha }\right)}^2\cdot \Gamma \left(\frac{4}{\alpha }\right)\cdot {\alpha }^{\frac{3}{\alpha }}\cdot \Gamma \left(1-\frac{3}{\alpha }\right)-\\ 3\cdot \Gamma \left(\frac{1}{\alpha }\right)\cdot \Gamma \left(\frac{2}{\alpha }\right)\cdot \Gamma \left(\frac{3}{\alpha }\right)\cdot \Gamma \left(1-\frac{2}{\alpha }\right)\cdot \Gamma \left(1-\frac{1}{\alpha }\right)\cdot {\alpha }^{\frac{3}{\alpha }}\end{array}}{\Gamma {\left(\frac{1}{\alpha }\right)}^3\cdot {\left[\frac{\Gamma \left(\frac{3}{\alpha }\right)\cdot {\alpha }^{\frac{2}{\alpha }}\cdot \Gamma \left(1-\frac{2}{\alpha }\right)}{\Gamma \left(\frac{1}{\alpha }\right)}-\frac{\Gamma {\left(\frac{2}{\alpha }\right)}^2\cdot {\alpha }^{\frac{2}{\alpha }}\cdot \Gamma {\left(1-\frac{1}{\alpha }\right)}^2}{\Gamma {\left(\frac{1}{\alpha }\right)}^2}\right]}^{1.5}}$ $\alpha =\exp \left(\begin{array}{l}0.000753505-0.983447042\cdot \ln \left(\left|{\tau }_3\right|\right)+0.162877024\cdot \ln {\left(\left|{\tau }_3\right|\right)}^2+0.247492727\cdot \ln {\left(\left|{\tau }_3\right|\right)}^3+\\ 0.116264397\cdot \ln {\left(\left|{\tau }_3\right|\right)}^4+0.029183504\cdot \ln {\left(\left|{\tau }_3\right|\right)}^5+0.004031622\cdot \ln {\left(\left|{\tau }_3\right|\right)}^6+\\ 0.000270699\cdot \ln {\left(\left|{\tau }_3\right|\right)}^7+0.000005628\cdot \ln {\left(\left|{\tau }_3\right|\right)}^8\end{array}\right)$

from Appendix C.

For frequency analysis of the annual volume deficit and duration, due to the fact that the data can contain zero values, it is necessary to use a conditional probability model [3,11,28]. The same model is recommended for application in the statistical analysis of annual low flows when the observed data contains zero values (rivers that dry up). This model assumes a correction of the cumulative function and the quantile of the analyzed distribution, as follows:

$$F{\left(x\right)}^{*}=\left(1-p_0\right)\cdot F\left(x\right)$$

(3)

$$x{\left(p\right)}^{*}=\left\{\begin{array}{ll}0,& \hfill p>F{\left(0\right)}^{*}\\ x\left(\frac{p}{1-p_0}\right)& \hfill \end{array}\right.$$

(4)

where $p_0$ is the empirical probability from which the values are zero.

The next section presents the theoretical distributions analyzed in the research of the Faculty of Hydrotechnics regarding the regionalization studies of the low flows in Romania.

The inverse forms of the distributions (quantile functions) can be expressed both for MOM and for L-moments with the frequency factor.

The frequency factors of the analyzed distributions are presented in Appendix A.

2.1. Probability Distributions

The probability density function, $f\left(x\right)$, the complementary cumulative distribution function, $F\left(x\right)$, and quantile function, $x\left(p\right)$, for analyzed distributions are:

2.1.1. Log-Normal (LN3)

It represents one of the most used statistical distributions in hydrology. The three-parameter log-normal distribution (LN3) is a generalized form of the two-parameter log-normal distribution (LN2), with x shifted by a position parameter that characterizes the lower bound [18,20,21,29].

$$f\left(x\right)=\frac{\exp \left(-\frac{{\left(\ln \left(x-\gamma \right)-\alpha \right)}^2}{2\cdot {\beta }^2}\right)}{\left(x-\gamma \right)\cdot \beta \cdot \sqrt{2\cdot \pi }}$$

(5)

$$F\left(x\right)=1-\frac{1}{2}\cdot \left(erf\left(\frac{1}{\sqrt{2}\cdot \beta }\cdot \left(\ln \left(x-\gamma \right)-\alpha \right)\right)+1\right)$$

(6)

$$x\left(p\right)=F^{-1}\left(x\right)=\gamma +\exp \left(\alpha +\beta \cdot qnorm\left(1-p,0,1\right)\right)=\gamma +qlnorm\left(1-p,\alpha ,\beta \right)$$

(7)

where $\alpha ,\beta ,\gamma$ are the shape, the scale, and the position parameters, respectively; $x$ can take any values of range $\gamma <x<\infty$.

2.1.2. Pearson III (PE3)

The Pearson III represents a generalized form of the two-parameter gamma distribution and a particular case of the four-parameter gamma distribution [18,20]. It represents the parent distribution in Romania in statistical analysis in hydrology [30,31].

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

(8)

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

(9)

$$x\left(p\right)=\gamma +\beta \cdot qgamma\left(1-p,\alpha \right)$$

(10)

where $\alpha ,\beta ,\gamma$ are the shape, the scale, and the position parameters, respectively, and $x$ can take any values of range $\gamma <x<\infty$ if $\beta >0$, or $-\infty <x<\gamma$ if $\beta <0$ and $\alpha >0$; $\mu ,\sigma$ represent the mean (expected value) and standard deviation. If $\beta <0$ (negative skewness), then the first argument of the inverse of the distribution function gamma, ${\Gamma }^{-1}\left(1-p;\alpha \right)$ becomes ${\Gamma }^{-1}\left(p;\alpha \right)$.

The built-in function from Mathcad, $qgamma\left(1-p,\alpha \right)={\gamma }^{-1}\left(\left(1-p\right)\cdot \Gamma \left(\alpha \right),\alpha \right)$, returns the inverse cumulative probability distribution for probability p, for the gamma distribution, where ${\gamma }^{-1}$ is the inverse of the lower incomplete gamma function [32].

2.1.3. Pearson V (IPV)

The Pearson V distribution represents the inverse of the Pearson III distribution [3,20].

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

(11)

$$F\left(x\right)=\frac{\Gamma \left(\alpha -1,\frac{\beta }{x-\gamma }\right)}{\Gamma \left(\alpha -1\right)}$$

(12)

$$x\left(p\right)=\gamma +\frac{\beta }{{\Gamma }^{-1}\left(p;\alpha -1\right)}=\gamma +\frac{\beta }{qgamma\left(p,\alpha -1\right)}$$

(13)

where $\alpha ,\beta ,\gamma$ are the shape, the scale, and the position parameters, respectively; $\alpha ,\beta >0;\gamma <x_{\min }$, and $x$ can take any values of range $\gamma <x<\infty$.

2.1.4. Wilson–Hilferty (WH3)

The three-parameter Wilson–Hilferty distribution is a generalized form of the two-parameter Wilson–Hilferty distribution. Both are cases of Amoroso distribution [20,29].

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

(14)

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

(15)

$$x\left(p\right)=\gamma +\beta \cdot \sqrt[3]{qgamma\left(1-p,\alpha \right)}$$

(16)

where $\alpha ,\beta ,\gamma$ are the shape, the scale, and the position parameters, respectively; $\alpha ,\beta >0$; $x$ can take any values in the range $\gamma <x<\infty$.

2.1.5. Fatigue Lifetime (FL3)

This distribution is also known as the Birnbaum–Saunders distribution [33].

$$f\left(x\right)=\frac{1+\left(x-\gamma \right)\cdot \beta }{2\cdot \alpha \cdot \sqrt{{\left(x-\gamma \right)}^3\cdot \beta }}\cdot \frac{\exp \left(-\frac{{\left(\beta \cdot \left(x-\gamma \right)-1\right)}^2}{2\cdot \left(x-\gamma \right)\cdot {\alpha }^2\cdot \beta }\right)}{\sqrt{2\cdot \pi }}$$

(17)

$$F\left(x\right)=1-{\displaystyle \underset{\gamma }{\overset{x}{\int }}\frac{1+\left(x-\gamma \right)\cdot \beta }{2\cdot \alpha \cdot \sqrt{{\left(x-\gamma \right)}^3\cdot \beta }}\cdot \frac{\exp \left(-\frac{{\left(\beta \cdot \left(x-\gamma \right)-1\right)}^2}{2\cdot \left(x-\gamma \right)\cdot {\alpha }^2\cdot \beta }\right)}{\sqrt{2\cdot \pi }}dx}$$

(18)

$$F\left(x\right)=1-\frac{1}{2}\cdot \left(1+erf\left(\frac{\beta \cdot \left(x-\gamma \right)-1}{\alpha \cdot \sqrt{2\cdot \left(x-\gamma \right)\cdot \beta }}\right)\right)$$

(19)

$$x\left(p\right)=\gamma +\frac{\alpha \cdot \sqrt{qnorm{\left(p,0,1\right)}^2}\cdot \sqrt{{\alpha }^2\cdot qnorm{\left(p,0,1\right)}^2+4}+{\alpha }^2\cdot qnorm{\left(p,0,1\right)}^2+2}{2\cdot \beta }, \mathrm{if} p\le 0.5$$

(20)

$$x\left(p\right)=\gamma +\frac{-\alpha \cdot \sqrt{qnorm{\left(p,0,1\right)}^2}\cdot \sqrt{{\alpha }^2\cdot qnorm{\left(p,0,1\right)}^2+4}+{\alpha }^2\cdot qnorm{\left(p,0,1\right)}^2+2}{2\cdot \beta }, \mathrm{if} p>0.5$$

(21)

where $\alpha ,\beta ,\gamma$ are the shape, the scale and the position parameters; $\alpha ,\beta >0$; $x$ can take any values in the range $\gamma <x<\infty$.

2.1.6. Weibull (W3)

The distribution represents a particular case of the generalized extreme value distribution when the shape parameter is positive. It is also known as the Type III extreme value distribution [7,11,12,16,18,20,34].

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

(22)

$$F\left(x\right)=e^{-{\left(\frac{x-\gamma }{\beta }\right)}^{\alpha }}$$

(23)

$$x\left(p\right)=\gamma +\beta \cdot {\left(-\ln \left(p\right)\right)}^{1/\alpha }$$

(24)

where $\alpha ,\beta ,\gamma$ are the shape, the scale, and the position parameters, respectively; $\alpha ,\beta >0$; $x$ can take any values in the range $\gamma <x<\infty$.

2.1.7. The Five-Parameter Wakeby Distribution (WK5)

The five-parameter Wakeby distribution has no form for density and cumulative function, being classified as a quantile function. The Wakeby distribution is a distribution that was introduced in the frequency analysis of extremes by Houghton [18,19,35] in order to fulfil, as much as possible, the “separation effect” described by Matalas et al. 1975 [35], namely to carry out an analysis so that the maximum flows are not influenced too much by the low flows, and vice versa. The Wakeby distribution separates the right-hand side from the left-hand side of the distribution. It has the property of a very thick left-hand tail (high flows) and a right-hand tail (low flow) thick enough to decrease the average skew, which makes the middle part of the distribution steeper than traditional skewed curves. A particular case of the five-parameter Wakeby distribution represents the Wakeby distribution of four parameters when $\xi =0$ [18] and the Pareto Type II distribution when $\gamma =0$, or in particular cases when $\frac{\alpha }{\beta }\cdot \left(1-p^{\beta }\right)$ or $\frac{\gamma }{\delta }\cdot \left(1-p^{-\delta }\right)$ has a constant value [19].

The quantile function of WK5 distributions is [18,35]

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

(25)

where $\alpha ,\gamma$ are the scale parameters, $\beta ,\delta$ are the shape parameters, and $\xi$ is the position parameter.

2.1.8. The Five-Parameter Lambda Distribution (L5)

The five-parameter lambda distribution seems part of the same family as the Wakeby distribution, being also classified as a quantile function. It has no form for density and cumulative function. It represents a distribution of a generalization of the Tukey lambda quantile function [36] and an alternative to the Wakeby distribution to fulfil the “separation effect”. A particular case of the five-parameter lambda distribution represents the lambda distribution of four parameters when $\gamma =0$ [36], and the Tukey distribution when $\lambda =0$, or in particular cases when $\alpha \cdot p^{\beta }$ or $\lambda \cdot {\left(1-p\right)}^{\delta }$ has a constant value.

The quantile function of L5 distributions is [36]:

$$x\left(p\right)=\gamma +\alpha \cdot p^{\beta }-\lambda \cdot {\left(1-p\right)}^{\delta }$$

(26)

where $\alpha ,\lambda$ are the scale parameters, $\beta ,\delta$ are the shape parameters, and $\gamma$ is the position parameter.

2.2. Parameter Estimation

The parameter estimation of the analyzed statistical distributions is presented for MOM and L-moments, two of the most used methods in hydrology for parameter estimation.

2.2.1. Log-Normal (LN3)

For estimation with MOM, the distribution parameters have the following expressions [37]:

$$\alpha =\ln \left(\sigma \cdot \frac{\sqrt[3]{2\cdot \left(\sqrt{C_s^2+4}-C_s\right)}}{\sqrt[3]{4}-\sqrt[3]{{\left(\sqrt{C_s^2+4}-C_s\right)}^2}}\right)-\frac{1}{2}\cdot \ln \left({\left(\frac{\sqrt[3]{4}-\sqrt[3]{{\left(\sqrt{C_s^2+4}-C_s\right)}^2}}{\sqrt[3]{2\cdot \left(\sqrt{C_s^2+4}-C_s\right)}}\right)}^2+1\right)$$

(27)

$$\beta =\ln {\left({\left(\frac{\sqrt[3]{4}-\sqrt[3]{{\left(\sqrt{C_s^2+4}-C_s\right)}^2}}{\sqrt[3]{2\cdot \left(\sqrt{C_s^2+4}-C_s\right)}}\right)}^2+1\right)}^{0.5}$$

(28)

$$\gamma =\mu -\sigma \cdot \frac{\sqrt[3]{2\cdot \left(\sqrt{C_s^2+4}-C_s\right)}}{\sqrt[3]{4}-\sqrt[3]{{\left(\sqrt{C_s^2+4}-C_s\right)}^2}}$$

(29)

The parameter estimation with the L-moments method is performed numerically (definite integrals) based on the equations using the quantile of the function [19].

An approximate form of parameter estimation can be adopted.

Thus, for the estimation with the L-moments, the shape parameter $\beta$ can be evaluated numerically with the following approximate forms, depending on L-skewness (${\tau }_3$) [20]:

if $0<\left|{\tau }_3\right|\le \frac{1}{3}$:

$$\begin{array}{l}\beta =4.379498\cdot {10}^{-4}+2.0295824\cdot \left|{\tau }_3\right|+0.23041762\cdot {\tau }_3^2-\\ 0.92166328\cdot {\left|{\tau }_3\right|}^3+3.8546644\cdot {\tau }_3^4-3.6560389\cdot {\left|{\tau }_3\right|}^5\end{array}$$

(30)

if $\frac{1}{3}<\left|{\tau }_3\right|\le \frac{2}{3}$:

$$\begin{array}{l}\beta =-5.3943247\cdot {10}^{-2}+2.6732827\cdot \left|{\tau }_3\right|-2.9211411\cdot {\tau }_3^2+\\ 7.3388138\cdot {\left|{\tau }_3\right|}^3-8.1997299\cdot {\tau }_3^4+4.3299046\cdot {\left|{\tau }_3\right|}^5\end{array}$$

(31)

if $\frac{2}{3}<\left|{\tau }_3\right|<1$:

$$\beta =\frac{0.81822527+7.818908\cdot {\tau }_3^4-8.1529453\cdot {\tau }_3^8+0.1554959\cdot {\tau }_3^{12}}{1+2.5188862\cdot {\tau }_3^4-4.9742123\cdot {\tau }_3^8+1.592552\cdot {\tau }_3^{12}}$$

(32)

With parameter $\beta$ known, the parameter $\alpha$ and the parameter $\gamma$ are determined with the following expressions:

$$\alpha =\ln \left(\frac{L_2}{erf\left(0.5\cdot \beta \right)}\right)-\frac{{\beta }^2}{2}$$

(33)

$$\gamma =L_1-\exp \left(\alpha +0.5\cdot {\beta }^2\right)$$

(34)

2.2.2. Pearson III (PE3)

For estimation with MOM, the distribution parameters have the following expressions [17,18,19,20]:

$$\alpha ={\left(\frac{2}{C_s}\right)}^2$$

(35)

$$\beta =\frac{\sigma }{2}\cdot C_s$$

(36)

$$\gamma =\mu -\alpha \cdot \beta$$

(37)

where $C_s$ represents the skewness coefficient.

The parameter estimation with the L-moments method is performed numerically (definite integrals) based on the equations using the quantile of the function. An approximate form of parameter estimation can be adopted. Thus, for the estimation with the L-moments, the shape parameter $\alpha$ can be evaluated numerically with the following approximate forms, depending on L-skewness (${\tau }_3$) [20]:

if $0<\left|{\tau }_3\right|\le \frac{1}{3}$:

$$\alpha =\exp \left(\begin{array}{l}-3.164791927-5.108735285\cdot \ln \left(\left|{\tau }_3\right|\right)-4.116014079\cdot \ln {\left(\left|{\tau }_3\right|\right)}^2-\\ 2.985250105\cdot \ln {\left(\left|{\tau }_3\right|\right)}^3-1.327399577\cdot \ln {\left(\left|{\tau }_3\right|\right)}^4-0.373944875\cdot \ln {\left(\left|{\tau }_3\right|\right)}^5-\\ 0.065421611\cdot \ln {\left(\left|{\tau }_3\right|\right)}^6-0.006508037\cdot \ln {\left(\left|{\tau }_3\right|\right)}^7-0.000281969\cdot \ln {\left(\left|{\tau }_3\right|\right)}^8\end{array}\right)$$

(38)

if $\frac{1}{3}<\left|{\tau }_3\right|\le \frac{2}{3}$:

$$\alpha =\exp \left(\begin{array}{l}-3.9918551-10.781466\cdot \ln \left(\left|{\tau }_3\right|\right)-21.557807\cdot \ln {\left(\left|{\tau }_3\right|\right)}^2-\\ 33.8752604\cdot \ln {\left(\left|{\tau }_3\right|\right)}^3-35.0641585\cdot \ln {\left(\left|{\tau }_3\right|\right)}^4-22.921163\cdot \ln {\left(\left|{\tau }_3\right|\right)}^5-\\ 8.5491823\cdot \ln {\left(\left|{\tau }_3\right|\right)}^6-1.3855653\cdot \ln {\left(\left|{\tau }_3\right|\right)}^7\end{array}\right)$$

(39)

if $\frac{2}{3}<\left|{\tau }_3\right|<1$:

$$\alpha =\frac{\begin{array}{l}5.17817436-26.209448756\cdot \left|{\tau }_3\right|+62.12494027\cdot {\tau }_3^2-84.39423264\cdot {\left|{\tau }_3\right|}^3+\\ 67.08589624\cdot {\tau }_3^4-29.150288079\cdot {\left|{\tau }_3\right|}^5+5.364968945\cdot {\tau }_3^6\end{array}}{1+0.0005134\cdot \left|{\tau }_3\right|+0.00063644\cdot {\tau }_3^2}$$

(40)

The scale parameter $\beta$ and the position parameter $\gamma$ are determined with the following expressions:

$$\beta =L_2\cdot \sqrt{\pi }\cdot \frac{\Gamma \left(\alpha \right)}{\Gamma \left(\alpha +\frac{1}{2}\right)}$$

(41)

$$\gamma =L_1-\alpha \cdot \beta$$

(42)

2.2.3. Pearson V (IPV)

The equations needed to estimate the parameters with MOM have the following expressions [20]:

$$\mu =\gamma +\frac{\beta }{\alpha -2}$$

(43)

$${\sigma }^2=\frac{{\beta }^2}{\alpha -2}\cdot \left(\frac{1}{\alpha -3}-\frac{1}{\alpha -2}\right)$$

(44)

$$C_s=4\cdot \frac{{\left(\alpha -3\right)}^{0.5}}{\alpha -4}$$

(45)

For estimation with MOM, parameter $\alpha$ can be evaluated numerically with the following approximate forms, depending on $C_s\ge 0.2$:

$$\alpha =\exp \left(\begin{array}{l}3.04183273-1.532335271\cdot \ln \left(C_s\right)+0.350524466\cdot \ln {\left(C_s\right)}^2+\\ 0.114266418\cdot \ln {\left(C_s\right)}^3-0.018372078\cdot \ln {\left(C_s\right)}^4-0.026493586\cdot \ln {\left(C_s\right)}^5-\\ 0.000919818\cdot \ln {\left(C_s\right)}^6+0.004375959\cdot \ln {\left(C_s\right)}^7+\\ 0.000310213\cdot \ln {\left(C_s\right)}^8-0.000357638\cdot \ln {\left(C_s\right)}^9\end{array}\right)$$

(46)

The scale parameter $\beta$ and the position parameter $\gamma$ are determined with the following expressions:

$$\beta =\sigma \cdot \left(\alpha -2\right)\cdot \sqrt{\alpha -3}$$

(47)

$$\gamma =\mu -\frac{\beta }{\alpha -2}$$

(48)

The parameter estimation with the L-moments method is performed numerically (definite integrals) based on the equations using the quantile of the function.

An approximate form of parameter estimation can be adopted. Thus, for the estimation with the L-moments, the shape parameter $\alpha$ can be evaluated numerically with the following approximate forms, depending on L-skewness (${\tau }_3$) [20]:

$$\alpha =\exp \left(\begin{array}{l}0.692740135-0.466055586\cdot \ln \left(\left|{\tau }_3\right|\right)+0.460508385\cdot \ln {\left(\left|{\tau }_3\right|\right)}^2+\\ 0.470172127\cdot \ln {\left(\left|{\tau }_3\right|\right)}^3+0.938027763\cdot \ln {\left(\left|{\tau }_3\right|\right)}^4+\\ 0.805363294\cdot \ln {\left(\left|{\tau }_3\right|\right)}^5+0.359846034\cdot \ln {\left(\left|{\tau }_3\right|\right)}^6+\\ 0.08945247\cdot \ln {\left(\left|{\tau }_3\right|\right)}^7+0.011775528\cdot \ln {\left(\left|{\tau }_3\right|\right)}^8+0.000641664\cdot \ln {\left(\left|{\tau }_3\right|\right)}^9\end{array}\right)$$

(49)

The scale parameter $\beta$ and the position parameter $\gamma$ are determined with the following expressions [16]:

$$\beta =\frac{L_2}{\frac{1}{\alpha -2}-2\cdot z}$$

(50)

$$\gamma =L_1-\frac{\beta }{\alpha -2}$$

(51)

where $z$ has the following expression:

$$z=\exp \left(\begin{array}{l}4.083852739-11.442833473\cdot \ln \left(\alpha \right)+10.701793287\cdot \ln {\left(\alpha \right)}^2-\\ 6.513639901\cdot \ln {\left(\alpha \right)}^3+2.527142191\cdot \ln {\left(\alpha \right)}^4-\\ 0.640389591\cdot \ln {\left(\alpha \right)}^5+0.105546499\cdot \ln {\left(\alpha \right)}^6-\\ 0.010903032\cdot \ln {\left(\alpha \right)}^7+0.000641094\cdot \ln {\left(\alpha \right)}^8-0.000016371\cdot \ln {\left(\alpha \right)}^9\end{array}\right)$$

(52)

2.2.4. Wilson–Hilferty (WH3)

The equations needed to estimate the parameters with MOM have the following expressions [20]:

$$\mu =\gamma +\frac{\beta }{\Gamma \left(\alpha \right)}\cdot \Gamma \left(\alpha +\frac{1}{3}\right)$$

(53)

$${\sigma }^2=\frac{{\beta }^2}{\Gamma \left(\alpha \right)}\cdot \left[\Gamma \left(\alpha +\frac{2}{3}\right)-\frac{1}{\Gamma \left(\alpha \right)}\cdot \Gamma {\left(\alpha +\frac{1}{3}\right)}^2\right]$$

(54)

$$C_s=\frac{\alpha -\frac{3}{\Gamma {\left(\alpha \right)}^2}\cdot \Gamma \left(\alpha +\frac{2}{3}\right)\cdot \Gamma \left(\alpha +\frac{1}{3}\right)+\frac{2}{\Gamma {\left(\alpha \right)}^3}\cdot \Gamma {\left(\alpha +\frac{1}{3}\right)}^3}{\sqrt{\frac{\Gamma \left(\alpha +\frac{2}{3}\right)}{\Gamma \left(\alpha \right)}-\frac{\Gamma {\left(\alpha +\frac{1}{3}\right)}^2}{\Gamma {\left(\alpha \right)}^2}}\cdot \left(\Gamma \left(\alpha +\frac{2}{3}\right)-\frac{\Gamma {\left(\alpha +\frac{1}{3}\right)}^2}{\Gamma \left(\alpha \right)}\right)}$$

(55)

The shape parameter can be obtained approximately, depending on the skewness coefficient, using the following exponential function [20]:

$$\alpha =\exp \left(\begin{array}{l}-1.6047146-1.2117058\cdot \ln \left(C_s\right)-2.4627986\cdot {10}^{-1}\cdot \ln {\left(C_s\right)}^2-\\ 3.0754515\cdot {10}^{-2}\cdot \ln {\left(C_s\right)}^3+1.3529125\cdot {10}^{-2}\cdot \ln {\left(C_s\right)}^4+\\ 5.4495596\cdot {10}^{-3}\cdot \ln {\left(C_s\right)}^5+6.0310303\cdot {10}^{-7}\cdot \ln {\left(C_s\right)}^6-\\ 3.5860178\cdot {10}^{-4}\cdot \ln {\left(C_s\right)}^7-7.3564689\cdot {10}^{-5}\cdot \ln {\left(C_s\right)}^8-4.7318329\cdot {10}^{-6}\cdot \ln {\left(C_s\right)}^9\end{array}\right)$$

(56)

$$\beta =\frac{\sigma }{\sqrt{\frac{\Gamma \left(\alpha +\frac{2}{3}\right)}{\Gamma \left(\alpha \right)}-\frac{\Gamma {\left(\alpha +\frac{1}{3}\right)}^2}{\Gamma {\left(\alpha \right)}^2}}}$$

(57)

$$\gamma =\mu -\frac{\beta \cdot \Gamma \left(\alpha +\frac{1}{3}\right)}{\Gamma \left(\alpha \right)}$$

(58)

The parameter estimation with the L-moments method is performed numerically (definite integrals) based on the equations using the quantile of the function.

An approximate form can be adopted based on the parameter estimation, depending on L-skewness (${\tau }_3$), as follows [20]:

if $0<{\tau }_3\le 1/3$:

$$\alpha =\exp \left(\begin{array}{l}-4.16202506-3.261604018\cdot \ln \left({\tau }_3\right)-1.783702334\cdot \ln {\left({\tau }_3\right)}^2-\\ 0.770946644\cdot \ln {\left({\tau }_3\right)}^3-0.221698815\cdot \ln {\left({\tau }_3\right)}^4-0.041191426\cdot \ln {\left({\tau }_3\right)}^5-\\ 0.004665295\cdot \ln {\left({\tau }_3\right)}^6-0.000287712\cdot \ln {\left({\tau }_3\right)}^7-0.000007207\cdot \ln {\left({\tau }_3\right)}^8\end{array}\right)$$

(59)

if $1/3<{\tau }_3\le 2/3$:

$$\alpha =\exp \left(\begin{array}{l}-5.264027693-9.134993593\cdot \ln \left({\tau }_3\right)-\\ 15.845477811\cdot \ln {\left({\tau }_3\right)}^2-19.874003352\cdot \ln {\left({\tau }_3\right)}^3-\\ 15.495267042\cdot \ln {\left({\tau }_3\right)}^4-6.752325319\cdot \ln {\left({\tau }_3\right)}^5-1.255615645\cdot \ln {\left({\tau }_3\right)}^6\end{array}\right)$$

(60)

if $2/3<{\tau }_3<1$:

$$\alpha =\exp \left(\begin{array}{l}-7.712526023-93.06660109\cdot \ln \left({\tau }_3\right)-1.519099681\cdot {10}^3\cdot \ln {\left({\tau }_3\right)}^2-\\ 1.602587644\cdot {10}^3\cdot \ln {\left({\tau }_3\right)}^4-1.041173897\cdot {10}^5\cdot \ln {\left({\tau }_3\right)}^4-\\ 4.152482998\cdot {10}^5\cdot \ln {\left({\tau }_3\right)}^5-9.887282\cdot {10}^5\cdot \ln {\left({\tau }_3\right)}^6-\\ 1.287749298\cdot {10}^6\cdot \ln {\left({\tau }_3\right)}^7-7.050767642\cdot {10}^5\cdot \ln {\left({\tau }_3\right)}^8\end{array}\right)$$

(61)

$$\beta =\frac{L_2}{\frac{\Gamma \left(\alpha +\frac{1}{3}\right)}{\Gamma \left(\alpha \right)}-2\cdot z}$$

(62)

$$\gamma =L_1-\frac{\beta \cdot \Gamma \left(\alpha +\frac{1}{3}\right)}{\Gamma \left(\alpha \right)}$$

(63)

where $z={\displaystyle \underset{0}{\overset{1}{\int }}qgamma{\left(1-p,\alpha \right)}^{1/3}\cdot p\cdot dp}$, which can be approximated with the following equation:

$$z=\exp \left(\begin{array}{l}-1.037385169+0.592727202\cdot \ln \left(\alpha \right)-0.107494558\cdot \ln {\left(\alpha \right)}^2+\\ 0.027616773\cdot \ln {\left(\alpha \right)}^3-0.002977204\cdot \ln {\left(\alpha \right)}^4-\\ 0.000546413\cdot \ln {\left(\alpha \right)}^5+0.00123125\cdot \ln {\left(\alpha \right)}^6+\\ 0.000420922\cdot \ln {\left(\alpha \right)}^7+0.000052295\cdot \ln {\left(\alpha \right)}^8+0.000002315\cdot \ln {\left(\alpha \right)}^9\end{array}\right)$$

(64)

2.2.5. Fatigue Lifetime (FL3)

The equations needed to estimate the parameters with MOM have the following expressions [33]:

$$\mu =\gamma +\frac{2+{\alpha }^2}{2\cdot \beta }$$

(65)

$${\sigma }^2=\frac{{\alpha }^2\cdot \left(4+5\cdot {\alpha }^2\right)}{4\cdot {\beta }^2}$$

(66)

$$C_s=\frac{4\cdot {\alpha }^2\cdot \left(11\cdot {\alpha }^2+6\right)}{\alpha \cdot \sqrt{{\left(4+5\cdot {\alpha }^2\right)}^3}}$$

(67)

The parameter estimation with the L-moments method is performed numerically (definite integrals), based on the equations using the quantile of the function.

An approximate form can be adopted based on the parameter estimation, depending on L-skewness (${\tau }_3$), as follows:

if $0<{\tau }_3\le 0.3$:

$$\alpha =\frac{0.000329932+2.035420236\cdot \left|{\tau }_3\right|}{1-0.060182808\cdot \left|{\tau }_3\right|-1.000075667\cdot {\tau }_3^2-0.630055585\cdot {\left|{\tau }_3\right|}^3}$$

(68)

if $0.3<{\tau }_3\le 0.7$:

$$\alpha =\frac{0.413986909-0.600184554\cdot \left|{\tau }_3\right|}{1-3.179942891\cdot \left|{\tau }_3\right|+3.026068085\cdot {\tau }_3^2-0.754825029\cdot {\left|{\tau }_3\right|}^3}$$

(69)

$$\beta =\frac{1}{2\cdot L_2}\cdot \frac{0.000467677+1.13508777\cdot \alpha +0.514463902\cdot {\alpha }^2+0.532996863\cdot {\alpha }^3}{\begin{array}{l}1+0.503786644\cdot \alpha +0.036559683\cdot {\alpha }^2-0.005584778\cdot {\alpha }^3+\\ 0.000507637\cdot {\alpha }^4-0.00002505\cdot {\alpha }^5+0.000000515\cdot {\alpha }^6\end{array}}$$

(70)

$$\gamma =L_1-\frac{2+{\alpha }^2}{2\cdot \beta }$$

(71)

2.2.6. Weibull (W3)

The equations needed to estimate the parameters with MOM have the following expressions [18,20]:

$$\mu =\gamma +\beta \cdot \Gamma \left(1+\frac{1}{\alpha }\right)$$

(72)

$${\sigma }^2={\beta }^2\cdot \left[\Gamma \left(1+\frac{2}{\alpha }\right)-\Gamma {\left(1+\frac{1}{\alpha }\right)}^2\right]$$

(73)

$$C_s=\frac{\Gamma \left(1+\frac{3}{\alpha }\right)-3\cdot \Gamma \left(1+\frac{2}{\alpha }\right)\cdot \Gamma \left(1+\frac{1}{\alpha }\right)+2\cdot \Gamma {\left(1+\frac{1}{\alpha }\right)}^3}{{\left(\Gamma \left(1+\frac{2}{\alpha }\right)-\Gamma {\left(1+\frac{1}{\alpha }\right)}^2\right)}^{1.5}}$$

(74)

The shape parameter can be obtained approximately, depending on the skewness coefficient, using the following rational function [20]:

$$\alpha =\frac{\begin{array}{l}3.5947875-3.7374529\cdot C_s+2.7023151\cdot C_s^2-1.1835102\cdot C_s^3+\\ 0.4318107\cdot C_s^4-0.0042306\cdot C_s^5+0.0045545\cdot C_s^6-0.0046984\cdot C_s^7\end{array}}{\begin{array}{l}1+0.0931316\cdot C_s-0.12783\cdot C_s^2+0.0671082\cdot C_s^3+0.1113667\cdot C_s^4-\\ 0.0425475\cdot C_s^5+0.0673629\cdot C_s^6-0.0157628\cdot C_s^7\end{array}}$$

(75)

$$\beta =\frac{\sigma }{\sqrt{\Gamma \left(1+\frac{2}{\alpha }\right)-\Gamma {\left(1+\frac{1}{\alpha }\right)}^2}}$$

(76)

$$\gamma =\mu -\beta \cdot \Gamma \left(1+\frac{1}{\alpha }\right)$$

(77)

The equations needed to estimate the parameters with L-moments have the following expressions [18,20]:

$$L_1=\gamma +\beta \cdot \Gamma \left(1+\frac{1}{\alpha }\right)$$

(78)

$$L_2=\beta \cdot \Gamma \left(1+\frac{1}{\alpha }\right)\cdot \left(1-2^{-\frac{1}{\alpha }}\right)$$

(79)

$$L_3=\beta \cdot \Gamma \left(1+\frac{1}{\alpha }\right)\cdot \left(1-2^{-\frac{1}{\alpha }}\right)\cdot \left(2\cdot \frac{2^{-\frac{1}{\alpha }}-3^{-\frac{1}{\alpha }}}{2^{-\frac{1}{\alpha }}-1}+1\right)$$

(80)

An approximate form of the shape parameters can be adopted, depending on L-skewness (${\tau }_3$), as follows [20]:

$$\alpha =\frac{3.528107902-6.294082546\cdot \left|{\tau }_3\right|+2.767652838\cdot {\tau }_3^2}{1+4.599024923\cdot \left|{\tau }_3\right|-7.993601572\cdot {\tau }_3^2+2.423742593\cdot {\left|{\tau }_3\right|}^3}$$

(81)

With parameter $\alpha$ known, the parameters $\beta$ and $\gamma$ are determined with the following expressions:

$$\beta =\frac{L_2}{\Gamma \left(1+\frac{1}{\alpha }\right)\cdot \left(1-2^{-\frac{1}{\alpha }}\right)}$$

(82)

$$\gamma =L_1-\frac{\beta }{\alpha }\cdot \Gamma \left(\frac{1}{\alpha }\right)$$

(83)

2.2.7. The Five-Parameter Wakeby Distribution (WK5)

The equations needed to estimate the parameters with MOM have the following expressions:

$$\mu =\xi +\frac{\alpha }{\beta +1}-\frac{\gamma }{\delta -1}$$

(84)

$${\sigma }^2=\frac{{\alpha }^2}{{\left(\beta +1\right)}^2\cdot \left(2\cdot \beta +1\right)}-\frac{2\cdot \alpha \cdot \gamma }{\left(\beta +1\right)\cdot \left(\beta +1-\delta \right)\cdot \left(\delta -1\right)}-\frac{{\gamma }^2}{{\left(\delta -1\right)}^2\cdot \left(2\cdot \delta -1\right)}$$

(85)

The equations for skewness and kurtosis are presented in Appendix D.

The equations needed to estimate the parameters with L-moments have the following expressions [16,18,35,36]:

$$L_1=\xi +\frac{\alpha }{\beta +1}+\frac{\gamma }{1-\delta }$$

(86)

$$L_2=\frac{\alpha }{\left(\beta +1\right)\cdot \left(\beta +2\right)}+\frac{\gamma }{\left(2-\delta \right)\cdot \left(1-\delta \right)}$$

(87)

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

(88)

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

(89)

2.2.8. The Five-Parameter Lambda Distribution (L5)

The equations needed to estimate the parameters with MOM have the following expressions:

$$\mu =\frac{\gamma \cdot \left(\beta +1\right)+\alpha }{\beta +1}-\frac{\lambda }{\delta +1}$$

(90)

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

(91)

The equations for skewness and kurtosis are presented in Appendix D.

The equations needed to estimate the parameters with L-moments have the following expressions [36]:

$$L_1=\frac{\gamma \cdot \left(\beta +1\right)+\alpha }{\beta +1}-\frac{\lambda }{\delta +1}$$

(92)

$$L_2=-\frac{\alpha \cdot \beta }{\left(\beta +1\right)\cdot \left(\beta +2\right)}-\frac{\lambda \cdot \delta }{\left(\delta +1\right)\cdot \left(\delta +2\right)}$$

(93)

$$L_3=\frac{\alpha \cdot \beta \cdot \left(\beta -1\right)}{\left(\beta +1\right)\cdot \left(\beta +2\right)\cdot \left(\beta +3\right)}+\frac{\lambda \cdot \delta \cdot \left(1-\delta \right)}{\left(\delta +1\right)\cdot \left(\delta +2\right)\cdot \left(\delta +3\right)}$$

(94)

$$L_4=-\frac{\alpha \cdot \beta \cdot \left(\beta -2\right)\cdot \left(\beta -1\right)}{\left(\beta +1\right)\cdot \left(\beta +2\right)\cdot \left(\beta +3\right)\cdot \left(\beta +4\right)}+\frac{\lambda \cdot \delta \cdot \left(3\cdot \delta -{\delta }^2-2\right)}{\left(\delta +1\right)\cdot \left(\delta +2\right)\cdot \left(\delta +3\right)\cdot \left(\delta +4\right)}$$

(95)

3. Case Study

The presented case study consists in the determination of minimum annual flows, drought durations and annual volume deficits on the Prigor River, Romania, using the proposed methodology.

The Prigor River is the left tributary of the Nera River, and it is located in the south-western part of Romania, as shown in Figure 4.

The main morphometric characteristics of the river are presented in

Table 2. The morphometric characteristics.

Length [km]	Average Stream Slope [‰]	Sinuosity Coefficient [-]	Average Altitude, [m]	Drainage Area, [km2]
33	22	1.83	713	153

[38].

In the section of the hydrometric station, the drainage area is 141 km² and the average altitude is 729 m. The minimum 1-day data are presented in

Table 3. The observed data from Prigor hydrometric station.

Minimum 1-Day Low Flows
		1990	1991	1992	1993	1994	1995	1996	1997	1998	1999	2000
Flow	[m3/s]	0.257	0.335	0.250	0.249	0.170	0.310	0.235	0.796	0.400	0.500	0.258
		2001	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011
Flow	[m3/s]	0.268	0.330	0.294	0.300	0.450	0.610	0.450	0.241	0.350	0.520	0.193
		2012	2013	2014	2015	2016	2017	2018	2019	2020
Flow	[m3/s]	0.216	0.173	0.687	0.180	0.388	0.207	0.182	0.185	0.390

.

There are 31 annual observed data for low flows, with the values of the main statistical indicators presented in

Table 4. The statistical indicators of the data series.

Low Flows	Statistical Indicators
	$\mathit{\mu }$	$\mathit{\sigma }$	${\mathit{C}}_{\mathit{v}}$	${\mathit{C}}_{\mathit{s}}$	${\mathit{C}}_{\mathit{k}}$	${\mathit{L}}_1$	${\mathit{L}}_2$	${\mathit{L}}_3$	${\mathit{L}}_4$	${\mathit{\tau }}_2$	${\mathit{\tau }}_3$	${\mathit{\tau }}_4$
	[m3/s]	[m3/s]	[-]	[-]	[-]	[m3/s]	[m3/s]	[m3/s]	[m3/s]	[-]	[-]	[-]
minimum 1-day	0.335	0.157	0.469	1.359	4.593	0.335	0.084	0.026	0.013	0.251	0.307	0.152

Where $\mu ,\sigma ,C_v,C_s,C_k,L_1,L_2,L_3,L_4,{\tau }_2,{\tau }_3,{\tau }_4$ represent the mean, the standard deviation, the coefficient of variation, the skewness, the kurtosis, the four L-moments, the L-coefficient of variation, the L-skewness, and the L-kurtosis, respectively P.

.

For parameter estimation with L-moments, the data series must be in ascending order for the calculation of natural estimators, the respective L-moments.

In this case study, for low-flow frequency analysis with MOM, a coefficient of 3 was chosen for the choice skewness, to reflect the perennial/ephemeral character of the river. It can be seen from Figure 5 that the predominant droughts are in the summer, and in the winter period, the important droughts appear only in 2012, 2013 being the driest year as well. The droughts of 2012 and 2013 are predominantly in the summer period, having the highest magnitude in those years.

The deficit volumes calculated for the thresholds of 95% are presented in

Table 5. The annual volume deficit.

95% Threshold
		1990	1991	1992	1993	1994	1995	1996	1997	1998	1999	2000
W	[1000 m3]	28	0.0	49.6	0.4	164.4	0.0	4.4	0.0	0.0	0.0	1.4
		2001	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011
W	[1000 m3]	1.5	0.0	0.0	0.0	0.0	0.0	0.0	131.7	0.0	0.0	100.8
		2012	2013	2014	2015	2016	2017	2018	2019	2020
W	[1000 m3]	82	133.8	0.0	179.7	0.0	233	141.1	122.8	0.0
W	[1000 m3]	3511.8	1907.7	0.0	1276.2	359.1	2402.3	1642.2	1338.4	366.9

. The analysis was also conducted for other thresholds, such as 90% and 80%.

The drought durations determined for the thresholds of 95% are presented in

Table 6. The annual drought durations.

95% Threshold
		1990	1991	1992	1993	1994	1995	1996	1997	1998	1999	2000
D	[days/year]	23	0.0	23	2.0	54	0.0	1.0	0.0	0.0	0.0	2.0
		2001	2002	2003	2004	2005	2006	2007	2008	2009	2010	2011
D	[days/year]	1.0	0.0	0.0	0.0	0.0	0.0	0.0	59	0.0	0.0	34
		2012	2013	2014	2015	2016	2017	2018	2019	2020
D	[days/year]	59	39	0.0	35	0.0	67	33	27	0.0

. The analysis was also conducted for other thresholds, such as 90% and 80%.

4. Results

The frequency analysis of hydrological drought consisted in determining the minimum flows, annual volumes and durations for exceedance probabilities used in water management [39]. Eight probability distributions, with three and five parameters, were used. The distribution parameters were estimated for MOM and L-moments for the analysis of minimum annual flows, and for L-moments for the analysis of annual volume deficits and annual drought durations.

For the analysis of minimum flows with MOM, the skewness coefficient was chosen, depending on the ephemeral/perennial character of the river. In this case, the multiplication coefficient 3, applied to the coefficient of variation of the data, was used, resulting in a skewness of 1.351, different from 1.053 of the observed values.

Figure 6, Figure 7 and Figure 8 show the frequency curves using the annual minimum flow (minimum 1-day flow), minimum 7-day flow, and minimum 30-day flow, respectively. For plotting positions, the Landwehr formula was used [18].

If the presented graphs are for three time series, the data are presented in the tables only for annual minimums, being specific to the run-of-river facilities. A frequency analysis was performed with distributions of three, four and five parameters for the two estimation methods, and the distributions that have a lower bound close to the minimum recorded value were selected. In general, it is considered that with a larger number of parameters, anything can be fit. In the case of parameter estimation with L-moments (the method recommended in this article, being much less influenced by the length of the observed data), it is necessary to estimate the first four linear moments, an aspect fulfilled properly only by distributions of at least four parameters. In the case of three-parameter distributions, they properly calibrate only L-skewness (L₃/L₂), and L-kurtosis (L₄/L₂) is the natural one of the distributions, not being able to calibrate according to the fourth order moment.

It should be mentioned, however, that not every distribution of at least four parameters gives good results; they must be used carefully, depending on the nature of extreme events (low flow or flood). In the analysis of extreme events, the most used are the Wakeby and lambda distributions, these having the advantage of the fact that the left hand is less influenced by the right hand, fulfilling the so-called ‘’separation effect’’ described by Matalas. Their advantage in the frequency analysis of low flows is represented by their asymptotic character in the field of high probabilities, with a lower bound close to the minimum value of the data sets.

The supplementary file graphically presents the results obtained for other four-parameter distributions (four-parameter kappa, four-parameter Kumaraswamy, four-parameter Burr, four-parameter Johnson SU). However, these cannot be applied for the analysis of the minimum extreme values, because for high probabilities of exceeding (right hand), they have a decreasing trend, being characterized by a small lower bound (even zero or negative); this is an unrealistic aspect for perennial rivers, generating unrealistic values in this field of probabilities.

From the case study, it is observed that among all the distributions analyzed, the ones that meet the conditions mentioned above are Wilson–Hilferty, Weibull, Wakeby and lambda. The Pearson III distribution was presented comparatively, because this is the parent distribution in Romania in the frequency analysis of extreme events. For the L-moments method, this gives good results.

Table 7. Quantile results with MOM.

P [%]	Probability Distributions
	LN3	PE3	IPV	WH3	FL3	W3	WK5	L5
	[m3/s]
50	0.304	0.299	0.306	0.272	0.301	0.296	0.290	0.297
70	0.240	0.235	0.241	0.218	0.238	0.232	0.228	0.237
90	0.169	0.172	0.167	0.199	0.171	0.174	0.179	0.179
95	0.142	0.152	0.138	0.198	0.147	0.158	0.169	0.159
97	0.126	0.142	0.121	0.198	0.134	0.151	0.165	0.147

and

Table 8. Quantile results with L-moments.

P [%]	Probability Distributions
	LN3	PE3	IPV	WH3	FL3	W3	WK5	L5
	[m3/s]
50	0.290	0.287	0.291	0.279	0.284	0.287	0.287	0.286
70	0.234	0.228	0.236	0.215	0.230	0.227	0.227	0.227
90	0.180	0.180	0.180	0.186	0.185	0.180	0.180	0.181
95	0.163	0.169	0.160	0.184	0.172	0.169	0.170	0.170
97	0.153	0.164	0.149	0.183	0.165	0.165	0.166	0.166

present the result values of quantiles for some of the most common exceedance probabilities for 1-day low flow analysis.

The resulting values of the estimated parameters are presented in

Table 9. The estimated values of the parameters.

Parameter	Distributions
	LN3	PE3	IPV	WH3	FL3	W3	WK5	L5	LN3	PE3	IPV	WH3	FL3	W3	WK5	L5
	MOM								L-moments
$\alpha$	−1.121	2.024	13.0	0.129	0.482	1.264	2.7 × 10−4	−3.487	−1.64	1.176	6.142	0.153	0.801	1.07	0.027	−4.1
$\beta$	0.421	0.110	5.451	0.519	3.490	0.212	3.984	0.045	0.644	0.152	1.296	0.507	6.312	0.181	3.984	0.043
$\gamma$	−0.022	0.112	−0.161	0.198	14.77	0.138	0.199	2.338	0.095	0.155	0.022	0.183	125.6	0.158	0.177	2.77
$\lambda$	-	-	-	-	-	-	-	−1.353	-	-	-	-	-	-	-	−1.49
$\delta$	-	-	-	-	-	-	−0.130	0.013	-	-	-	-	-	-	−0.041	0.001
$\xi$	-	-	-	-	-	-	0.159	-	-	-	-	-	-	-	0.160	-

, for both estimation methods, for 1-day low flow analysis.

For the analysis of deficit volumes and durations of annual droughts, only three distributions of three parameters are presented. The calculations were made only for the estimation of parameters with L-moments, because MOM requires the choice of skewness, there being no criteria for this. The same considerations apply to the statistical analysis of drought durations. The distributions PE3 and W3 were chosen because they best estimate the data set, PE3 being also the parent distribution in Romania. The LN3 distribution is presented comparatively, as it is generally considered an alternative to PE3 [20].

Figure 9 shows the frequency curves used for the annual deficit volumes corresponding to the thresholds of 95%, 90% and 80%.

Table 10. Quantile results for annual deficit volume.

Distributions	Probability [%]
	5	10	20	30	40	50	60	70	80	90
95% threshold—[1000 m3]
LN3	193.01	154.51	106.01	66.71	23.51	0	0	0	0	0
PE3	196.16	154.69	106.06	66.63	23.35	0	0	0	0	0
W3	194.24	156.61	107.05	65.83	21.01	0	0	0	0	0

presents the result values of quantiles for annual deficit volume.

Figure 10 shows the frequency curves used for the annual drought durations corresponding to the thresholds of 95%, 90% and 80%. The calculations were made only to estimate the parameters with L-moments.

Table 11. Quantile results for annual drought duration.

Distributions	Probability [%]
	5	10	20	30	40	50	60	70	80	90
95% threshold—[days/year]
LN3	60.84	49.81	35.53	23.62	10.14	0	0	0	0	0
PE3	60.85	49.83	35.54	23.61	10.13	0	0	0	0	0
W3	61.11	50.33	35.83	23.42	9.45	0	0	0	0	0

presents the result values of quantiles for annual drought duration.

The performance of the analyzed distribution was evaluated using the relative mean error (RME) and relative absolute error (RAE) criteria [40].

$$RME=\frac{1}{n}\cdot \sqrt{{{\displaystyle \sum _{i=1}^n\left(\frac{x_i-x\left(p\right)}{x_i}\right)}}^2}$$

(96)

$$RAE=\frac{1}{n}\cdot {\displaystyle \sum _{i=1}^n\left|\frac{x_i-x\left(p\right)}{x_i}\right|}$$

(97)

where $n,x_i,x\left(p\right)$ represent sample size, observed value, and estimated value for a given probability, respectively.

The distribution performance values are presented in

Table 12. Distribution performance values for minimum annual flows.

Distribution	Methods of Parameter Estimation						Observed Data
	MOM		L-moments
	RME	RAE	RME	RAE	${\mathit{\tau }}_3$	${\mathit{\tau }}_4$	${\mathit{\tau }}_3$	${\mathit{\tau }}_4$
					[-]	[-]	[-]	[-]
LN3	0.0169	0.0577	0.0087	0.0342	0.307	0.197	0.307	0.152
PE3	0.0114	0.0431	0.0057	0.0257		0.159
IPV	0.0187	0.0632	0.0102	0.0392		0.212
WH3	0.0124	0.0573	0.0096	0.0444		0.083
FL3	0.0141	0.0495	0.0066	0.0306		0.186
W3	0.0087	0.0367	0.0056	0.0254		0.155
WK5	0.0058	0.0273	0.0055	0.0253		0.152
L5	0.0110	0.0401	0.0056	0.0256		0.152

,

Table 13. Distribution performance values for annual deficit volumes.

Distribution	Methods of Parameter Estimation
	L-moments
	RME	RAE
95% threshold
LN3	0.2118	0.7498
PE3	0.2119	0.7498
W3	0.2121	0.7497

and

Table 14. Distribution performance values for annual drought durations.

Distribution	Methods of Parameter Estimation
	L-moments
	RME	RAE
95% threshold
LN3	0.2087	0.7343
PE3	0.2087	0.7343
W3	0.2093	0.7357

, only for the analysis of minimum annual flows (1-day), annual deficit volumes and annual drought durations.

We observed good performances for all the distributions used, for the probabilities defined on the range of observed values. For MOM estimation, the best approximation is the Wakeby distribution. For estimation with L-moments, a good approximation of L-skewness and L-kurtosis is observed for three distributions, Wakeby, lambda and Weibull.

In the case of the analysis for annual deficit volumes and annual drought durations, the performances are good for the distributions LN3 and PE3.

5. Conclusions

This article presents hydrological drought analysis using frequency analysis. Considering the complexity of this hydrological phenomenon, it is necessary to establish well-founded criteria regarding the analysis of minimum flows, deficit volumes and drought durations. Thus, this article presents the necessary criteria regarding the frequency analysis of minimum annual flows, the establishment of drought thresholds based on the duration curve of the average daily flows, as well as the deficit volumes and drought durations corresponding to the thresholds. A preferred analysis in hydrology is the use of average monthly flows, which often do not reflect the phenomenon of drought. In this article, the analysis with time series moving averages for 7 and 30 days is proposed [7]. The methodology provides the necessary mathematical support for frequency analysis of hydrological drought for the elaboration of a norm in Romania. Mathematical support in statistical analysis is useful, because the use of software (EasyFit, HEC-SSP, …) without knowledge of mathematical foundations often leads to superficial analyses with negative consequences in water management.

The case study consisted of hydrological drought analysis for the Prigor River. The frequency analysis was performed using eight probability distributions, their parameters being estimated with two of the most used parameter estimation methods in the analysis of extreme phenomena in hydrology. All the mathematical characteristics necessary for the application of the eight analyzed distributions, for the method of ordinary moments and the method of linear moments are presented.

For ease of application of these distributions, parameter approximations were made for MOM and L-moments.

For confidence interval calculation using Chow’s relation [37], Appendix A presents the frequency factor of distributions analyzed for both parameter estimation methods. A criterion for choosing distributions for L-moments analysis is the variance of ${\tau }_3-{\tau }_4$. Figure A1 from Appendix B presents approximation relations of ${\tau }_4$ with regards to ${\tau }_3$ for a wide range of distributions used in hydrology.

A peculiarity in the analysis of the frequency of annual drought volumes and durations is the choice of thresholds. Many researchers avoid choosing minimum flow thresholds (Q90%,Q95%) that generate many zero values of volumes and durations, thus ignoring these values corresponding to these thresholds needed in water management. This impediment of threshold selection is solved by using conditional probabilities, which apply a correction to the cumulative function and implicitly to the inverse function.

All research was carried out within the Faculty of Hydrotechnics for elaboration of a norm regarding hydrological drought and for improving the legislation of ecological flow, based on mathematical methods for determining the water deficit.