*Article* **A Probabilistic Model for Maximum Rainfall Frequency Analysis**

**Maurycy Ciupak 1, Bogdan Ozga-Zieli ´nski 1, Tamara Tokarczyk 1,\* and Jan Adamowski <sup>2</sup>**


**Abstract:** As determining the probability of the exceedance of maximum precipitation over a specified duration is critical to hydrotechnical design, particularly in the context of climate change, a model was developed to perform a frequency analysis of maximum precipitation of a specified duration. The PMAXTP model (Precipitation MAXimum Time (duration) Probability) harbors a pair of computational modules fulfilling different roles: (i) statistical analysis of precipitation series, and (ii) estimation of maximum precipitation for a specified duration and its probability of exceedance. The input data consist of homogeneous 30-element series of precipitation values for 16 different durations: 5, 10, 15, 30, 45, 60, 90, 120, 180, 360, 720, 1080, 1440, 2160, 2880, and 4320 min, obtained through Annual Maximum Precipitation (AMP) and Peaks-Over-Threshold (POT) approaches. The statistical analysis of the precipitation series includes: (i) detecting outliers using the Grubbs-Beck test; (ii) checking for the random variable's independence using the Wald-Wolfowitz test and the Anderson serial correlation coefficient test; (iii) checking the random variable's stationarity using nonparametric tests, e.g., the Kruskal-Wallis test and Spearman rank correlation coefficient test for trends of mean and variance; (iv) identifying the trend of the random variables using correlation and regression analysis, including an evaluation of the form of the trend function; and (v) checking for the internal correlation of the random variables using the Anderson autocorrelation coefficient test. To estimate maximum precipitations of a specified duration and with a specified probability of exceedance, three-parameter theoretical probability distributions were used: a shifted gamma distribution (Pearson type III), a log-normal distribution, a Weibull distribution (Fisher-Tippett type III), a log-gamma distribution, as well as a two-parameter Gumbel distribution. The best distribution was selected by: (i) maximum likelihood estimation of parameters; (ii) tests of the hypothesis of goodness of fit of the theoretical probability distribution function with the empirical distribution using Pearson's *χ*<sup>2</sup> test; (iii) selection of the best-fitting function within each type according to the criterion of minimum Kolmogorov distance; (iv) selection of the most credible probability distribution function from the set of various types of best-fitting functions according to the Akaike information criterion; and (v) verification of the most credible function using single-dimensional tests of goodness of fit: the Kolmogorov-Smirnov test, the Anderson-Darling test, the Liao-Shimokawa test, and Kuiper's test. The PMAXTP model was tested on data from two meteorological stations in northern Poland (Chojnice and Bialystok) drawn from a digital database of high-resolution precipitation records for the period of 1986 to 2015, available for 100 stations in Poland (i.e., the Polish Atlas of Rainfall Intensities (PANDa)). Values of maximum precipitation with a specified probability of exceedance obtained from the PMAXTP model were compared with values obtained from the probabilistic Bogdanowicz-Stachý model. The comparative analysis was based on the standard error of fit, graphs of the density function for the probability of exceedance, and estimated quantile errors. The errors of fit were lower for the PMAXTP compared to the Bogdanowicz-Stachý model. For both stations, the smallest errors were obtained for the quantiles determined on the basis of maximum precipitation POT using PMAXTP.

**Citation:** Ciupak, M.; Ozga-Zieli ´nski, B.; Tokarczyk, T.; Adamowski, J. A Probabilistic Model for Maximum Rainfall Frequency Analysis. *Water* **2021**, *13*, 2688. https://doi.org/ 10.3390/w13192688

Academic Editor: Renato Morbidelli

Received: 1 September 2021 Accepted: 22 September 2021 Published: 28 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Keywords:** annual maximum precipitation; peaks-over-threshold methods; statistical analysis; maximum precipitation frequency analysis; gamma; Weibull; log-gamma; log-normal; Gumbel distributions; nonparametric tests

#### **1. Introduction**

A frequency analysis of values of maximum precipitation of a specified duration and probability of exceedance is an essential part of engineering [1]. Given the significant impact of maximum precipitation on various spheres of human activity (e.g., the economy, agriculture, industry, and the environment), such an analysis is widely applied, particularly in the context of observed climate change [2,3].

A widely used tool in the statistical description of rare meteorological (climatic) events is the extreme value theorem (EVT). Two probability distributions are used when employing the EVT: the generalized extreme value distribution (GEV) and the generalized Pareto distribution (GPD) [4,5]. Encompassing three families of distributions (Gumbel (G), Fréchet (F), and Weibull (WE)), the GEV distribution offers the advantage of high accuracy of fit to observed precipitation data [6]. Commonly used methods for the estimation of the unknown parameters of theoretical probability distributions include: maximum likelihood, L-moments, and the Bayesian method [7–9]. Ragulina and Reitan [10] proposed a Bayesian hierarchical model approach to the selection of a GEV distribution, where Bayesian inference was applied both to parameter estimation and model selection. For most locations in Japan investigated by Yuan et al. [11], a log-Pearson type 3 distribution (LGA) proved to be the best-fitting theoretical probability distribution for annual maximum hourly precipitation data. Mły ´nski et al. [12] found that among the G, GA, WE, log-normal (LN), and GEV distributions, the latter best described annual maximum daily precipitation in Poland's upper Vistula basin.

An assumption of the EVT is that the random variables subjected to analysis show stationarity, i.e., the statistical properties of the mechanism generating these variables remain unchanged over time. Such conditions are rarely encountered in nature, and extreme events are increasingly of a nonstationary nature. In the case of maximum precipitation, its natural variation is overlaid by changes in climate and human intervention in land use (e.g., reduction in soil drainage). In this situation, time series of maximum precipitation values exhibit non-stationarity in the form of long-term trends and/or periodic fluctuations. In recent years, it has become increasingly common to analyze the frequency of nonstationary phenomena using the theory of nonstationary extreme value (NSEV). Katz et al. [13] extended the traditional approach to a frequency analysis to deal with nonstationary cases, where it is assumed that there is a constant probability of the occurrence of an extreme event with values that vary with time. Likewise, Adlouni et al. [14] developed a method for estimating a GEV distribution under nonstationary conditions. Parameters of the distribution were estimated by the maximum likelihood method (MLM), and the covariance of the observed variables was included in the parameters of the probability distribution.

Another approach, used in engineering practice for estimating values of maximum precipitation with a specified duration and probability of exceedance, is regionalization. In Poland, Bogdanowicz and Stachý [15,16] used a clustering procedure for a series of annual maximum precipitation values to distinguish three precipitation regions. In these regions, annual maximum values were described using a WE extreme value distribution. Satisfying the assumptions of independence, stationarity, and identity of probability distribution, Shahzadi et al. [17] used a regional analysis of flooding frequency and a Monte Carlo method to divide the territory of Pakistan into three homogeneous subregions. The estimation of parameters followed the L-moments method, while quantile estimation was carried out using GA, GEV, GPA, generalized normal (GNO), and generalized logistic (GLO) distributions.

Quantiles of an extreme value distribution are usually estimated directly from a random sample of annual maximum precipitation (AMP) values. In view of the shortness of the time series, alternative solutions were used, thereby enabling statistical inference to be carried out based on a broader set of information than the annual maxima. Examples include analyses of seasonal maxima and models of annual maxima with different seasonal variances. In these models, the probabilistic description is usually based on mixed distributions. Earlier research on mixed distributions assumed the same probability density function for the distinguished seasons (homogeneous mixed distributions). An example of this approach is the two-population general extreme value distribution (TPGEV), based on the assumption of GEV-GEV distributions [18], gamma-gamma distributions (GA-GA), and log-normal-log-normal distributions (LN-LN) [19,20]. However, hydrometeorological variables are composed of different types of probability density functions.

Numerous studies on non-homogeneous mixed distributions have led to an improvement of the characteristics of the analyzed variables through the use of two-component models, such as the mixed gamma-Gumbel distribution (GA-G) [21] or the two-component generalized extreme value distribution (TCGEV) composed of a GEV and a Gumbel (G) distribution [22]. A GA-GP mixed distribution, incorporating a gamma distribution [23] and generalized Pareto distribution (GP), is commonly used. It serves mainly to model meteorological situations featuring both dry and wet periods. Another approach to the frequency analysis of maximum precipitation is the determination of the relationship between the intensity of precipitation and duration, and between duration and frequency of occurrence. For the modeling of two-dimensional dependences, the use of copula functions is recommended as a method of estimation of a two-dimensional distribution function [24,25]. In recent years, analyses have been made of a multidimensional dependence structure of extreme precipitation event variables using vine copula functions. The method involves the step-by-step mixing of two-dimensional copulas, which leads to a simplification of the estimation of multidimensional distribution functions [26].

Although there have been many attempts at using models for nonstationary series of extreme events [27–34], engineering practice shows that the assumption of the stationarity of time series is still widely adopted.

The purpose of this paper is to present the PMAXTP model for a frequency analysis of maximum precipitation with a specified duration and probability of exceedance, together with the results of testing the model against data from two meteorological stations located in northern Poland: Chojnice and Białystok. Values of maximum precipitation with a specified duration and probability of exceedance were estimated for two time series: (i) a 30-year series of annual maximum precipitation (AMP) values from the period 1986–2015 and (ii) a 30-element series of maximum precipitation values from the period 1986–2015 obtained by means of peaks-over-threshold (POT) analysis. The 30 highest values from the obtained set were used for further analyses. Computations were performed for 16 different durations: 5, 10, 15, 30, 45, 60, 90, 120, 180, 360, 720, 1080, 1440, 2160, 2880, and 4320 min. The results given by the PMAXTP model were compared with those obtained with the probabilistic Bogdanowicz-Stachý model of maximum precipitation [15,16], which is in common use in Polish engineering practice.

#### **2. Problem Formulation and Methodology**

The PMAXTP model for a frequency analysis of maximum precipitation with a specified duration and probability of exceedance was developed with the use of the method of alternative events (MAE), which serves to compute annual maximum flows with a specified probability of exceedance [35]. The overall scheme of the PMAXTP model is shown in Figure 1. The model contains two computational modules, one that performs a statistical analysis of series of precipitation data, and another that estimates maximum precipitation with a given duration and probability of exceedance. The latter includes an estimation of parameters of the distributions by the maximum likelihood method, verification of goodness of fit by Pearson's *χ<sup>2</sup>* test, selection of the best-fitting probability distribution function

within each distribution type according to the criterion of minimum Kolmogorov distance, selection of the most credible function according to the Akaike information criterion (AIC), and determination of the quantile confidence interval with regard to the randomness of the series of observations. The results returned by the PMAXTP model are values of maximum precipitation with a specified duration τ (min) ∈ {5, 10, 15, 30, 45, 60, 90, 120, 180, 360, 720, 1080, 1440, 2160, 2880, 4320} and a given probability of exceedance *p* (%) ∈ {99.9, 99.5, 99, 98.5, 98, 95, 90, 80, 70, 60, 50, 40, 30, 20, 10, 5, 3, 2, 1, 0.5, 0.3, 0.2, 0.1, 0.05, 0.03, 0.02, 0.01}.

**Figure 1.** Overall scheme of the PMAXTP model.

An analysis of the homogeneity of the random variables of series of maximum precipitation with different durations was performed by genetic (physical) methods and by statistical methods [35,36]. The identification of the trend of the analyzed random variables and evaluation of the form of the trend function were carried out by correlation and regression analysis, where the dependent variable is the maximum precipitation selected by the AMP or POT method, and the independent variable is the time (*τ*). The correlation was analyzed using the nonparametric Spearman rank correlation test [37] and the parametric Pearson linear correlation coefficient test [38]. In regression analysis, the global Fisher-Snedecor F-test [39] tests three equivalent null hypotheses: the significance of the slope, the significance of the coefficient of determination, and the significance of the linear relationship between the analyzed variables. Verification is performed for the null hypothesis that the independent variable (time *τ*) has no effect on the analyzed dependent variable, which here is the maximum precipitation (*P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* ). An evaluation of the form of the trend function is performed using scatter plots of the analyzed random variables with respect to time (*τ*). These provide a visual assessment and an evaluation of the form of the trend function: linear, power, exponential, etc.

The internal correlation of the analyzed random variable was checked using the Anderson autocorrelation coefficient test [40]. This analysis identifies the occurrence of periodic fluctuations and their effect on the variation of the analyzed variables. The results are presented numerically and graphically for a specified lag, with an indication of the autocorrelation coefficients and an evaluation of white noise (standard error) for the confidence level assumed (*α*).

The computation of the maximum precipitation with a specified probability of exceedance is performed using probabilistic models of the properties of the random variables *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* . An analysis of the properties of random maximum precipitations served as the basis for the acceptance of potential probability distribution models: e.g., G, GA, LN, log-gamma (LGA), and WE. The first four models are three-parameter distributions with the following parameters: *α* (*α* > 0), *λ* (*λ* > 0) or *μ* (*μ* > 0), and *ε* (*ε* ≤ *x* ≤ + ∞), representing, respectively, the parameters of scale, shape, and position, i.e., the lower (left-hand) limit of the probability distribution (see details in Appendix A).

The PMAXTP model assumes that each type of distribution is represented by a family of functions *fi*(*x*), shifted with respect to each other, each of which has a certain fixed lower limit (*εi*) satisfying 0 ≤ *ε<sup>i</sup>* < min 1≤*j*≤*n xj* , where *n* is the size of the random sample. The value of *ε<sup>i</sup>* may take values ranging from 0 up to the minimum value of the variable (*X*) in the random sample (*x*1, *x*2, ··· , *xn*). Hence, the lower limit (*εi*) of the *i*th specific function in the family of a selected type of distributions is the discriminant of that function within the family, and is not subject to estimation. In the G distribution, described by Equations (A9) and (A10) in Appendix A, only two parameters appear: the scale *α* and the shape *μ*.

The parameters of probability density functions were estimated by the MLM using dedicated software [41]. The procedure was as follows:


rejected. Sets of noncontradictory functions are formed separately for each selected probability distribution function type: GA, WE, LGA, and LN.


#### **3. Study Area and Data**

The PMAXTP model was tested on data from two meteorological stations located in Poland: Chojnice and Bialystok (Figure 2, black hexagons). The choice of stations was based on the availability of long series of historical data and current meteorological observations.

Data were drawn from the Rain-Brain database, created under the Development and Implementation of a Polish Atlas of Rainfall Intensities (PANDa) project [51] carried out in 2016 and 2017 by Poland's Institute of Meteorology and Water Management—National Research Institute (IMGW—PIB). Under the PANDa project, a series of depths of precipitation having specific durations were subjected to qualitative assessment, including a comparison of digital records with analog data (from Hellmann rain gauges), and information was drawn from a system of ground-based radars operating in the measurement and observation network of the IMGW—PIB. The observations were verified with respect to the occurrence of meteorological configurations which might cause rainfall of a given quantity in specified pressure conditions, characteristic of the analyzed region.

The study was based on the 30 highest precipitation depth values for 16 specified durations, *τ* = {5, 10, 15, 30, 45, 60, 90, 120, 180, 360, 720, 1080, 1440, 2160, 2880, 4320} (minutes) for the two precipitation stations mentioned above.

**Figure 2.** Location of the Chojnice and Bialystok meteorological stations in Poland.

Two methods were used to select maximum precipitation values: AMP [1,2,52] and POT [53]. Under the AMP method, a single maximum precipitation value was selected for the year, independent of its duration. A defect of the AMP method is that it fails to take into account all the high precipitation depth values occurring in a given year. In the POT method, it is possible to take into account all high precipitation depth values in a given year, i.e., the method selects these values that exceed a threshold determined a priori. The analyses were based on events with values not less than *P*POT min,*<sup>τ</sup>* = 3.5*τ*0.275 [51]. Thus, threshold values *P*MAX *<sup>τ</sup>* (mm) were set for precipitation with specified durations (*τ*), as given in Table 1 [51]. The subsequent analyses used 30-element series of maximum precipitation data, selected by both methods.

**Table 1.** Minimum quantity of precipitation *P*POT min,*<sup>τ</sup>* (mm) taken as a threshold in the POT method.


#### **4. Results and Discussion**

*4.1. Results of Analysis of Homogeneity for the PMAXTP Model*

An analysis was made of the genetic, time, and measurement homogeneity of the precipitation series from the stations in Chojnice and Bialystok. Based on a visual assessment of the measurement series and information contained in IMGW—PIB reports (Meteorological Yearbooks and Precipitation Yearbooks Report [51]), no significant factors were

found that might have an impact on the genetic homogeneity of the series of maximum precipitation values observed in the years 1986–2015.

An analysis was made of the statistical properties of the series of precipitation measurements from Chojnice and Bialystok using nonparametric significance tests [35,36]. The results are presented in Tables 2–6. Tables 2 and 3 contain the results of outlier detection using the Grubbs-Beck test [54,55], checking for the independence of the analyzed random variable using the Wald-Wolfowitz test (Test of Series) and Anderson serial autocorrelation coefficient test [40,55,56], and checking the stationarity of the analyzed random variable using the Kruskal-Wallis test and Spearman rank correlation coefficient test for the trends of mean and variance [57,58]. The final column of Tables 2 and 3 indicates genetically and statistically homogeneous series of maximum precipitation data selected by the AMP and POT methods.

In the case of *P*AMP *<sup>τ</sup>* , the Grubbs-Beck test detected outliers for precipitation with the duration *τ* = 360 and *τ* = 720 min, at both the Chojnice station (Table 2) and the Bialystok station (Table 3). In Tables 2 and 3, for a positive test result (+), the number of the outlier in the chronological sequence and the quantity of precipitation are also given. For the *P*POT *<sup>τ</sup>* series at Chojnice (Table 2), outliers were detected for *τ* ∈ {15, 30} and *τ* ∈ {120, . . . , 4320} min, while at Bialystok (Table 3), outliers were detected for *τ* ∈ {5, . . . , 15}, *τ* ∈ {60, ... , 360} and *τ* ∈ {2160, ... , 4320} min. Based on the theorem developed by Neyman and Scott [59] stating that the families of LN, G, and WE distributions—these being the distributions assumed as potential models describing the maximum precipitation values—are entirely susceptible to the occurrence of outliers in a random sample, it was concluded that the occurrence of the detected outliers should be considered entirely natural, and such elements were not removed from the measurement series.




**Table 3.** Results of nonhomogeneity analysis of AMP and POT precipitation series from Bialystok meteorological station; (−)/(+) denotes, respectively, negative and positive test results; <sup>√</sup>—denotes homogenous series.

> For all observed values of maximum precipitation *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* (Tables 2 and 3), the Wald-Wolfowitz test (Test of Series) and the Anderson serial correlation coefficient test showed that the analyzed measurement series were random and formed a simple sample, i.e., the random variables were independent variables. The significance level *α* = 0.05 used in the test took account of the size of the random sample, *n* = 30. For series of length greater than 30, a lower value may be taken as the test significance level (e.g., *α* = 0.01). For the detection of outliers with the Grubbs-Beck test, the higher value *α* = 0.10 was used, on the assumption that series of measurements of meteorological phenomena may be characterized by greater anthropogenic impact.

> The stationarity of the measurement series was checked using the Kruskal-Wallis test and Spearman rank correlation test for the trends of the mean and variance. According to the Kruskal-Wallis test, in the *P*AMP *<sup>τ</sup>* series from both Chojnice and Bialystok, jumps in the mean were detected, with the exception of the observations for *τ* = 5 and *<sup>τ</sup>* ∈ {720, . . . , 4320} min at Chojnice. In the case of the *<sup>P</sup>*POT *<sup>τ</sup>* precipitation values, most of the observations were stationary, with the exception of *τ* = 5 at Chojnice and *τ* ∈ {15, 30} and *τ* = 1440 min at Bialystok.

> The Spearman's rank correlation test for the trends of mean and variance revealed nonstationarity mainly for the *P*AMP *<sup>τ</sup>* precipitation values. In the case of *P*POT *<sup>τ</sup>* , nonstationary observations were the exception. For example, in the observations from Chojnice for *τ* = 10 min and *τ* ∈ {45, 60} min, a trend was detected in the mean and variance, respectively, while for the Bialystok data, such trends were detected, respectively, for *τ* ∈ {360, 1440} and *τ* = 120 min.

> The results of correlation testing and the identification of the trend of maximum precipitation for the AMP and POT series are given in Tables 4–6. The identification of the trend of the analyzed random variables was performed using the nonparametric Spearman rank correlation test [37] and the parametric Pearson linear correlation coefficient test [38]. An analysis was made of the correlation between the studied random variables (*P*AMP <sup>τ</sup> and *P*POT <sup>τ</sup> ) and the time variable *τ* (Table 4). Positive and negative values indicate upward and

downward trends, respectively. Spearman's coefficient also indicates the strength of the trend. The closer the values are to 1.0, the stronger is the relationship between the analyzed random variable and the time variable *τ*. Pearson's coefficient indicates proportionality, that is, linear dependence between variables, while Spearman's coefficient indicates any monotonic relationship, even if nonlinear. Figures shown in bold type in Table 4 indicate significant correlations, with the probability *p* ≤ 0.05. Strong dependences between the observed maximum precipitation values and the independent variable *τ* were recorded in the case of *P*AMP <sup>τ</sup> at both Chojnice and Bialystok.

**Table 4.** Correlations between the maximum precipitation variables and time *τ* for the Chojnice and Bialystok stations. Bold values of Spearman's rank correlation and Pearson's linear correlation coefficients are significant at *p* < 0.05 for *n* = 30, where *n* is the size of the sample.


The form of the trend function was assessed using regression analysis (Tables 5 and 6), where the dependent variable is the maximum precipitation and the independent variable is the time *τ*. Tables 5 and 6 give the results of the regression analysis, including the following indicators: Pearson's correlation coefficient *r*, the coefficient of determination *r*2, the Fisher-Snedecor global *F*-test [60], the test probability *p* resulting from the latter test, the size of the random sample *n*, and the standard error of estimation S(*E*). Statistically significant regression coefficients for the analyzed variables are identified according to the criterion for statistical significance adopted in the model, with *α* = 0.05. This means that the regression coefficients are significant for a test probability *p* ≤ 0.05.

The global *F*-test tests three equivalent null hypotheses: H0: *β*<sup>1</sup> = 0 (significance of the slope); H0: *r*<sup>2</sup> = 0 (significance of the coefficient of determination); and H0: *y* = *β*1*x*+ *β*<sup>0</sup> (significance of the linear relationship between the analyzed variables), where *β*<sup>1</sup> is the slope; *β*<sup>0</sup> is a free term; and *x* and *y* denote the independent and dependent variables, respectively. Verification is made of the null hypothesis that the independent variable *x* (in Tables 5 and 6, the independent variable is time, *τ*) does not influence the analyzed dependent variable *y* (in Tables 5 and 6, the dependent variables are *P*AMP <sup>5</sup> , ... , *<sup>P</sup>*AMP 4320 and *P*POT <sup>5</sup> , ... , *<sup>P</sup>*POT <sup>4320</sup> ). If, in the course of verification, the null hypothesis is rejected, the regression coefficient is assessed as significant, meaning that *τ* has a significant influence on the analyzed dependent variable. Examples of random variables with no trend and showing a trend are given in Tables 5 and 6, respectively, for observations from Chojnice and Bialystok.

**Table 5.** Results of simple regression analysis for the Chojnice station, where the dependent variables are *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* , and the independent variable is time (*τ*), for *n* = 30, where *r* is Pearson's correlation coefficient; *r*<sup>2</sup> is the coefficient of determination; *F*(1,*n*) is the Fisher-Snedecor test; *S*(*E*) is the standard error of estimation; and *p* (*p*-value) is the value of the test probability. Bold type indicates significance of regression parameters, namely the existence (for *p* ≤ 0.05) of a significant linear trend coefficient.


**Table 6.** Results of simple regression analysis for the Bialystok station, where the dependent variables are *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* , and the independent variable is time (*τ*), for *n* = 30, where *r* is Pearson's correlation coefficient; *r*<sup>2</sup> is the coefficient of determination; *F*(1,*n*) is the Fisher-Snedecor test; *S*(*E*) is the standard error of estimation; and *p* (*p*-value) is the value of the test probability. Bold type indicates significance of regression parameters, namely the existence (for *p* ≤ 0.05) of a significant influence of the variable *τ* on the analyzed dependent variable.


Values shown in bold type in Tables 5 and 6 indicate the presence of a significant influence of time *τ* on the analyzed random variable. In these cases, the estimated regression slope coefficients *β*<sup>1</sup> are significantly different from zero. At Chojnice, the observations of maximum precipitation showed a trend only in the case of *P*AMP <sup>τ</sup> for the durations *τ* ∈ {90, . . . , 360} min. At Bialystok, however, in all of the analyzed observations of maximum precipitation *P*AMP <sup>τ</sup> and in three cases of *P*POT <sup>τ</sup> (*τ* ∈ {720, ... , 1440} min), an upward

trend was detected. The test probability *p* determined for the computed regression coefficients was below the assumed significance level *α* = 0.05.

An assessment of the form of the trend function (linear, power, exponential, etc.) was made using scatter plots of the analyzed random variables with respect to time *τ* (Figure 3). The scatter plots of *P*AMP <sup>10</sup> , *<sup>P</sup>*AMP <sup>30</sup> , and *<sup>P</sup>*AMO <sup>60</sup> showed a clear linear upward trend, while those for the variables *P*POT <sup>10</sup> , *<sup>P</sup>*POT <sup>30</sup> , and *<sup>P</sup>*POT <sup>60</sup> showed, respectively, small upward and downward trends. In this case, the slope *β*<sup>1</sup> was close to 0, and the test probabilities (*P*POT <sup>10</sup> : *p* = 0.660; *P*POT <sup>30</sup> : *p* = 0.749; *<sup>P</sup>*POT <sup>60</sup> : *p* = 0.928) were substantially higher than the significance level *α* = 0.05 used in the analysis. In the annual data, seasonal (monthly or daily) fluctuations were not analyzed. If the analyzed series of values of *P*AMP *<sup>τ</sup>* or *P*POT *<sup>τ</sup>* contain a trend or periodic fluctuations, they cannot be used as an input in the computational procedures of the PMAXTP method.

**Figure 3.** Scatter plots of dependent random variables observed at the Bialystok station: *P*POT <sup>10</sup> , *<sup>P</sup>*POT *ao* , *P*POT <sup>60</sup> and *<sup>P</sup>*AMP <sup>10</sup> , *<sup>P</sup>*AMP <sup>30</sup> , *P*AMP <sup>60</sup> with respect to the independent variable time (*τ*), with indication of the simple regression equation, coefficient of determination (*r*2), linear correlation coefficient (*r*), and test probability (*p*) compared with the assumed significance level *α* < 0.05.

An analysis was made of the internal correlation of the series of random variables *P*AMP <sup>τ</sup> and *P*POT <sup>τ</sup> using Anderson's test [40]. An autocorrelation analysis was performed for lags up to 25 (Figure 4). The greatest autocorrelation coefficients were detected for *P*AMP <sup>1080</sup> with *lag* =1(*<sup>ρ</sup>* = 0.358) and for *<sup>P</sup>*POT <sup>90</sup> with *lag* =4(*ρ* = 0.417). Other autocorrelation values were not large and lay within the confidence interval for the assumed significance level *α* = 0.05. This is a sufficient condition to conclude a lack of correlation; that is, that the analyzed random variables are independent. An analysis of the autocorrelation plots (Figure 4) also showed an absence of periodic fluctuations.

Nonhomogeneity analysis, performed using genetic and statistical methods, showed that most of the observations of maximum precipitation selected by the POT method satisfied the homogeneity requirements, except for the observations for duration *τ* = {10, 45, 60} min at Chojnice and *τ* = {15, 30, 120, 1440} min at Bialystok (Tables 2 and 3). Most of the maximum precipitation observations selected by the AMP method are nonhomogeneous; exceptions are the *P*AMP <sup>τ</sup> observations from Chojnice with duration *τ* = 5 and *τ* = {720, . . . , 4320} min.

**Figure 4.** Autocorrelation function of random variables observed at Bialystok: *P*POT <sup>10</sup> , *<sup>P</sup>*POT <sup>30</sup> , and *<sup>P</sup>*POT <sup>60</sup> for lags of up to 25 elements in a series, with indication of autocorrelation coefficients, calculated white noise (standard error), and confidence level *α*.

#### *4.2. Computation of Maximum Precipitation with Specified Probability of Exceedance Using the PMAXTP Method*

Parameters of the probability distributions of the analyzed random variables were estimated for the two adopted methods of selection of maximum precipitations, *P*AMP *τ* and *P*POT *<sup>τ</sup>* (for details, see Section 2). The most credible distribution was selected for the analyzed random variable by minimizing the value of the Akaike information criterion (AIC) [45]. Calculations were performed for three-parameter (*α*, *λ* or *μ*, *ε*; Equations (A1), (A3), (A5) and (A7) in Appendix A) probability distributions GA, WE, LGA, and LN, and for the two-parameter (*α*, *μ*; Equation (A9) in Appendix A) G distribution. Sample results obtained at each stage of the procedure are given in Table 7. The most credible theoretical probability distribution for precipitation *P*AMP <sup>5</sup> at the Chojnice station was found to be GA, while for *P*POT <sup>5</sup> , it was found to be WE. At the Bialystok station, the most credible theoretical distribution for *P*POT <sup>5</sup> was determined to be LGA.

Verification of the distributions of maximum precipitation identified as most credible at the meteorological stations in Chojnice and Bialystok was performed by means of nonparametric tests of goodness of fit: DK−S, DA−D, DL−S, and DK (defined by Equations (A11)–(A14) in Appendix B). For purposes of inference, a significance level of *α* = 0.05 was arbitrarily selected. This is a consequence of the fact that the value of the significance level of a test is closely related to the size (length) of the random sample on whose basis the parameters of the theoretical distributions are estimated. In the present analysis, the series contained *n* = 30 elements, which means that the significance level can be taken to be at most *α* = 0.05. Verification was performed for the most credible theoretical probability distributions, which are shown in Table 8 for maximum precipitation with specified duration *τ*, together with the results obtained in single-dimensional statistical tests and the critical values, respectively for *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* at the Chojnice station and *P*POT *<sup>τ</sup>* at Bialystok. All of the tests failed to reject the null hypothesis on the goodness of fit of the theoretical distribution with the empirical distribution, for the analyzed variables *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* , with the exception of the DA−<sup>D</sup> test in relation to the maximum precipitation *<sup>P</sup>*POT <sup>90</sup> at Chojnice (value shown in bold type in Table 8). The least of the maximum distances between values of the theoretical and empirical cumulative probability distributions, particularly in

the tail part, was situated decidedly below the critical value of the DA−<sup>D</sup> test defined at a significance level of *α* = 0.05, which signifies rejection of the hypothesis of the goodness of fit of the theoretical and empirical distributions.

**Table 7.** Sample results of the procedure to select probability distributions for maximum precipitation values *P*AMP *<sup>τ</sup>* and *P*POT *τ* for *τ* = 5 min. GA—gamma distribution; WE—Weibull; LN—log-normal; LGA—log-gamma; G—Gumbel; *χ2*—Pearson's *χ<sup>2</sup>* goodness-of-fit test; min*(D*max)—Kolmogorov's minimum distance criterion. Bold values represent the most credible distributions according to the Akaike information criterion (AIC).


**Table 8.** Results of tests of fit of the theoretical probability distributions for *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* , where *τ* = {5, 10, 15, 30, 45, 60, 90, 120, 180, 360, 720, 1080, 1440, 2160, 2880, 4320} (min). DK−S—Kolmogorov-Smirnov test, DA−D— Anderson-Darling test, DL−S—Liao-Shimokawa test, DK—Kuiper's test, significance level *α* = 0.05. The value in bold type indicates rejection of the hypothesis of goodness of fit to the empirical distribution according to the statistic DA−<sup>D</sup> at *α* = 0.05.


*α*cr. = 0.05 for: DK−Scr. = 0.242; DA−Dcr. = 0.795; DL−Scr. = 1.505; DKcr. = 0.317.

The results obtained from the PMAXTP model for the values of maximum precipitation with a specified probability of exceedance were compared with the results from the Bogdanowicz-Stachý model [1,2]. In the latter model, the procedure for computing the values of maximum precipitation with a specified probability of exceedance *p* consisted of:


The procedure of the Bogdanowicz-Stachý model conforms to the recommendations of the World Meteorological Organization [61]. The input data originated from 20 meteorological stations situated in latitudinal strips running along the coast, lake districts, lowland parts, and southern upland parts of Poland. Mountain areas were omitted, due to the absence of stations monitoring precipitation at all altitudes. The maximum quantity of precipitation with a specified duration and specified probability of exceedance was determined using the formula (A16) in Appendix C, taking account of the regionalization of the meteorological stations in Chojnice and Bialystok.

Quantile values determined using the PMAXTP and Bogdanowicz-Stachý models were compared using statistical and graphical measures. According to the regionalization carried out by Bogdanowicz and Stachý, the Chojnice meteorological station belongs to the north-west region for precipitation with durations in the range <5, >60 min, to the central region for durations in the range <60, >720 min, and to the southern/coastal region for durations in the range <720, >4320 min. The Bialystok station, located in the northeast of Poland, belongs to the central region irrespective of the duration of precipitation being considered.

For a comparison of the results given by the two models, i.e., PMAXTP and Bogdanowicz-Stachý, various statistical measures can be used [62]. In our study, we used the standard error of fit S(*E*), which is shown in Table 9. The error is given by the following formula [63]:

$$\mathbf{S}(E) = \sqrt{\sqrt{\sum\_{i=1}^{i=m} \left(P\_{\overline{\tau\_i}}^{\text{MAX}} - P\_{\overline{\tau\_i}}^{\text{MAX}}\right)}} \tag{1}$$

where *P*MAX *<sup>τ</sup>*,*<sup>i</sup>* is the observed maximum precipitation selected by the AMP or POT method for a specified duration (*τ*); *P*ˆMAX *<sup>τ</sup>*,*<sup>i</sup>* is the estimated maximum precipitation from the PMAXTP or Bogdanowicz-Stachý model; *m* = 30 is the size of the random sample formed from empirical quantiles for *m* = 30 selected probabilities *p* ∈ {96.8, 93.6, 90.3, 87.1, 83.9, 80.7, 77.4, 74.2, 70.9, 67.7, 64.5, 61.3, 58.1, 54.8, 51.6, 48.4, 45.2, 41.9, 38.7, 35.5, 32.3, 29.0, 25.8, 22.6, 19.4, 16.1, 12.9, 9.7, 6.5, 3.2} % and the corresponding theoretical distributions computed using the PMAXTP and Bogdanowicz-Stachý methods. Finally, *l* is the number of parameters of the theoretical probability distribution according to the density function (Equations (A1), (A3), (A5), (A7) and (A9) in Appendix A).

Computations of the error S(*E*) were performed separately for specified durations *τ* of maximum precipitation. The value of the standard error of fit increased with increasing values of *τ* for both models. The smallest errors were obtained for the quantiles determined from the maximum precipitation values selected using the POT method and the PMAXTP model. An exception was the quantiles determined for the AMP values at the Chojnice station for duration *τ* equal to 720 and 4320 min. The errors of fit of the theoretical to the empirical distributions in the Bogdanowicz-Stachý model for precipitation values selected by the AMP method were on average 210% greater than those obtained with the PMAXTP model, and for the POT precipitation values, the errors were 300% greater. The most frequently selected most credible theoretical probability distribution for random samples of both AMP and POT maximum precipitation values, and for both the Bialystok and the Chojnice stations, was the WE distribution.

Figures 5–10 show a comparison of the functions for the probability of exceedance of maximum precipitations *P*AMP *<sup>τ</sup>* or *P*POT *<sup>τ</sup>* determined using the models, for Chojnice (Figures 5–7) and Bialystok (Figures 8–10). The plots contain density functions of probability distributions computed only for homogeneous observations of precipitation selected by the AMP and POT methods, in accordance with the results shown in Tables 2, 3 and 8. The diagrams show comparisons of: (i) the most credible probability functions for maximum

precipitation determined by the PMAXTP model for the AMP observation series (orange solid line) and for maximum precipitation selected by the POT method (blue solid line); (ii) upper limits of confidence intervals (orange and blue dotted lines); (iii) observations of AMP and POT maximum precipitation (orange and blue squares); and (iv) the probability function determined using the probabilistic Bogdanowicz-Stachý model (red solid line).

**Table 9.** Comparison of the PMAXTP and Bogdanowicz-Stachý methods for *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* , where *τ* = {5, 10, 15, 30, 45, 60, 90, 120, 180, 360, 720, 1080, 1440, 2160, 2880, 4320} (min), using the standard error of fit S(*E*). The comparison refers to the maximum precipitation values computed for the meteorological station in Chojnice (*P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* ) and in Bialystok (*P*POT *<sup>τ</sup>* ). Values in bold type are the smallest errors S(*E*) obtained separately for the Chojnice and Bialystok stations.


At the Chojnice station, for practically all of the analyzed durations of maximum precipitation, the quantile values from the Bogdanowicz-Stachý model are markedly higher than the observed precipitations and values of corresponding quantiles from the PMAXTP model, in relation to the maximum precipitations selected both by the AMP method (orange squares and solid line) and by the POT method (blue squares and solid line). The differences between the quantiles are particularly visible in the central region and in the region of the upper tails of the probability distributions. Similar maximum quantile values were obtained for precipitation with duration *τ* = {15, 30, 180, 1080} min. At Chojnice, the AMP values were described by the models GA and WE, while for description of the POT maximum precipitation values, the WE distribution was selected for short durations *τ*, and GA and LGA for medium and long durations.

At the Bialystok station, in the case of maximum precipitations with duration *τ* = {5, 45, 60, 90, 180} min (Figures 8 and 9), the quantile values determined using the Bogdanowicz-Stachý model (red solid line) are markedly higher than the corresponding quantiles obtained using the PMAXTP model for the maximum precipitations determined by the POT method (blue squares and solid line). Differences between quantiles are particularly visible in the central region and in the region of the upper tails of the probability distributions. The closest results for quantiles of POT maximum precipitations calculated using the PMAXTP method and from the Bogdanowicz-Stachý model were obtained for precipitation with duration *τ* = {720, 1080} min (Figure 9). For maximum precipitation with such durations, the most credible theoretical distribution was WE, while for short durations, *τ* = {5, 10} min, the respective distributions were LGA and GA. For maximum precipitation selected by the POT method with duration *τ* = {2160, ... , 4320} min, the Bogdanowicz-Stachý model returned markedly lower quantile values than the PMAXTP method.

**Figure 5.** Plots of functions of probability of exceedance for the random variables *P*AMP *<sup>τ</sup>* and/or *P*POT *<sup>τ</sup>* , where *τ* = {5, 15, 30, 120} min, for the most credible probability distributions, with indicated upper limits of quantile confidence intervals according to the PMAXTP method, compared with the model of Bogdanowicz and Stachý, for the Chojnice meteorological station.

**Figure 6.** Plots of functions of probability of exceedance for the random variables *P*AMP *<sup>τ</sup>* and/or *P*POT *<sup>τ</sup>* , where *τ* = {180, 360, 720, 1080} min, for the most credible probability distributions, with indication of upper limits of quantile confidence intervals according to the PMAXTP method, compared with the model of Bogdanowicz and Stachý, for the Chojnice meteorological station.

**Figure 7.** Plots of functions of probability of exceedance for the random variables *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* , where *τ* = {1440, 2160, 2880, 4320} min, for the most credible probability distributions, with indication of upper limits of quantile confidence intervals according to the PMAXTP method, compared with the model of Bogdanowicz and Stachý, for the Chojnice meteorological station.

**Figure 8.** Plots of functions of probability of exceedance for the random variables *P*POT *<sup>τ</sup>* where *τ* = {5, 10, 45, 60} min, for the most credible probability distributions, with indication of upper limits of quantile confidence intervals according to the PMAXTP method, compared with the model of Bogdanowicz and Stachý, for the Bialystok meteorological station.

**Figure 9.** Plots of functions of probability of exceedance for the random variables *P*POT *<sup>τ</sup>* where *τ* = {90, 180, 720, 1080} min, for the most credible probability distributions, with indication of upper limits of quantile confidence intervals according to the PMAXTP method, compared with the model of Bogdanowicz and Stachý, for the Bialystok meteorological station.

**Figure 10.** Plots of functions of probability of exceedance for the random variables *P*POT *<sup>τ</sup>* where *τ* = {2160, 2880, 4320} min, for the most credible probability distributions, with indication of upper limits of quantile confidence intervals according to the PMAXTP method, compared with the model of Bogdanowicz and Stachý, for the Bialystok meteorological station.

The final element of the verification of maximum precipitation values was a comparison of the estimated quantile error resulting from the randomness of the sample of maximum precipitations computed using the PMAXTP model for the random variables *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* at the meteorological station in Chojnice (Figures 11–13) and for *P*POT *<sup>τ</sup>* at the meteorological station in Bialystok (Figures 14–16).

The largest errors for values of maximum precipitation with high probabilities, such as 99.0 and 99.9, at the Chojnice station were observed for maximum precipitations selected using the AMP method (Figure 11 for *τ* = 5, Figure 12 for *τ* = {720, 1080}, Figure 13 for *τ* = {1440, . . . , 2160} min)—markedly higher errors for the AMP series than for the POT series at Chojnice. The largest errors for values of maximum precipitation with low probabilities, such as 0.01 and 0.001, were recorded for the Chojnice station (Figure 12 for *τ* = {720, 1080} and Figure 13 for *τ* = {1440, ... , 2160} min) for POT precipitations (markedly higher errors for the POT series than for the AMP series at Chojnice). The smallest differences in the quantile error in the entire range of theoretical occurrence of maximum precipitation were observed at Chojnice (Figure 13 for *τ* = {2880, 4320} min).

Calculations were made for 100 total rainfall measuring sites in Poland (Figure 17). Calculated characteristics of maximum rainfall totals, i.e., quantile values for *p*(%) ∈ {99.9, 99.5, 99, 98.5, 98, 95, 90, 80, 70, 60, 50, 40, 30, 20, 10, 5, 3, 2, 1, 0.5, 0.3, 0.2, 0.1, 0.05, 0.03, 0.02, 0.01} of a specified duration, *τ*(min) ∈ {5, 10, 15, 30, 45, 60, 90, 120, 180, 360, 720, 1080, 1440, 2160, 2880, 4320}, upper limits of the confidence interval and quantile errors were interpolated by the Thiessen Polygons (TP) method, which allowed for the assignment of certain areas for which measuring sites are representative as well as for the proportional division and distribution of sites within Poland. Higher resolution calculations can be achieved using Gaussian geostatistical simulation models [64] that accept any simple kriging model [65] or residual kriging model [66].

**Figure 11.** Comparison of estimated values of quantile error resulting from the randomness of the sample of maximum precipitations computed using the PMAXTP model for the random variables *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* with durations *τ* = {5, 15, 30, 120} min, for the Chojnice meteorological station.

**Figure 12.** Comparison of estimated values of quantile error resulting from the randomness of the sample of maximum precipitations computed using the PMAXTP method for the random variables *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* with durations *τ* = {180, 360, 720, 1080} min, for the Chojnice meteorological station.

**Figure 13.** Comparison of estimated values of quantile error resulting from the randomness of the sample of maximum precipitations computed using the PMAXTP method for the random variables *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* with durations *τ* = {1440, 2160, 2880, 4320} min, for the Chojnice meteorological station.

**Figure 14.** Comparison of estimated values of quantile error resulting from the randomness of the sample of maximum precipitations computed using the PMAXTP model for the random variable *P*POT *<sup>τ</sup>* with durations *τ* = {5, 10, 45, 60} min, for the Bialystok meteorological station.

**Figure 15.** Comparison of estimated values of quantile error resulting from the randomness of the sample of maximum precipitations computed using the PMAXTP model for the random variable *P*POT *<sup>τ</sup>* with durations *τ* = {90, 180, 720, 1080} min, for the Bialystok meteorological station.

**Figure 16.** Comparison of estimated values of quantile error resulting from the randomness of the sample of maximum precipitations computed using the PMAXTP model for the random variable *P*POT *<sup>τ</sup>* with durations *τ* = {2160, 2880, 4320} min, for the Bialystok meteorological station.

Interpolation also can be performed using the Inverse Distance Weighted (IDW) method, which uses a linearly weighted set of sampling points to determine mesh node values by using reverse weighted distance values. The weight is a function of the inverse distance, and the interpolated surface should be a variable surface depending on the position of the point [67]. An example of interpolating the maximum precipitation value *P*AMP *<sup>τ</sup>* with a duration of *τ* = 30 min with a probability *p* = 1% calculated using the IDW method is shown in Figure 18 (left part).

The IDW is a deterministic interpolation method because it is directly based on surrounding measured values. Another example is the set of geostatistical methods, such as the Kriging methods (right part of Figure 18), which include autocorrelation, which represents the statistical relationship between the measured points, thus providing a certain measure of reliability or accuracy of the forecast. The Kriging method is most suitable when one knows that there is spatial distance correlation or directional deviation in the data being analyzed.

**Figure 17.** Thiessen Polygons based on precipitation measurement sites in Poland.

**Figure 18.** Interpolation of maximum precipitations computed using the PMAXTP model for the random variable *P*AMP *τ* with durations *τ* = 30 min with probability of exceedance *p* = 1% using IDW method (**left** part) and kriging method (**right** part) for the Bialystok and Chojnice meteorological stations.

#### **5. Conclusions**

This paper described the PMAXTP model for a frequency analysis of maximum precipitation with a specified duration. It consists of two modules: statistical and computational. The first step selects values of maximum precipitation of a specified duration, which is conducted using two different methods: Annual Maximum Precipitation (AMP) and Peaks-Over-Threshold (POT). The advantage of the POT method is that it selects a larger number of observations of precipitation with the highest values in a given year, which leads to a better estimation of the characteristics of maximum precipitation with a specified duration and probability of exceedance. This is a significant issue in the design of drainage structures, particularly when they are at high risk of damage. The statistical module enables an analysis of the homogeneity of the series of measurements of maximum precipitation that serve as the input to the computational module, in which the mathematical models used for parameter estimation require a simple random sample, that is, one that satisfies the assumptions of independence and stationarity.

The computational module enables the selection of the best (the most credible) theoretical probability distribution by means of: (i) estimation of the parameters of four types of distributions belonging to the families gamma (GA), Weibull (WE), log-gamma (LGA), lognormal (LN), and Gumbel function (G); (ii) test of the hypothesis of goodness of fit of the theoretical probability distribution function with the empirical distribution using Pearson's χ<sup>2</sup> test; (iii) selection of the best-fitting function in each distribution type according to the criterion of minimum Kolmogorov distance; (iv) selection of the most credible distribution function from the set of best-fitting functions of various types; and (v) verification of the most credible distributions of precipitations *P*AMP *<sup>τ</sup>* and *P*POT *<sup>τ</sup>* using the single-dimensional tests DK−S, DA−D, DL−S, and DK.

The PMAXTP model was tested on data from two meteorological stations in Poland (Chojnice and Bialystok) representing two regions characterized by different spatial variability of maximum precipitation. The results were compared with those given by the Bogdanowicz-Stachý model—which to date has frequently been used in engineering practice in Poland—based on estimated values of the quantile error resulting from the randomness of the sample of maximum precipitation values computed for the tested stations.

In general, the errors of fit for the theoretical to the empirical distribution for the PMAXTP model were lower than the errors for the Bogdanowicz-Stachý model. The smallest errors were obtained for the quantiles determined on the basis of maximum precipitation POT using the PMAXTP model for both analyzed stations.

The following detailed conclusions may be drawn from the results:


**Author Contributions:** Conceptualization, B.O.-Z. and M.C.; methodology, B.O.-Z. and M.C.; software, M.C. and B.O.-Z.; validation, B.O.-Z., T.T. and J.A.; formal analysis, B.O.-Z. and M.C.; investigation, M.C. and B.O.-Z.; resources, Institute of Meteorology and Water Management–National Research Institute; data curation, Institute of Meteorology and Water Management–National Research Institute; writing—original draft preparation, M.C., B.O.-Z., T.T. and J.A.; writing—review and editing, M.C., B.O.-Z., T.T. and J.A.; visualization, M.C.; supervision, B.O.-Z. and T.T.; project administration, B.O.-Z.; funding acquisition, B.O.-Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** The authors acknowledge the financial and data support provided by the Polish Hydrological and Meteorological Service at the Institute of Meteorology and Water Management— National Research Institute.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A. The Density Function** *f***(***x***) and the Quantile Function** *xp*

The density function *f*(*x*) and the quantile function *xp* of the three-parameter GA distribution are written as [68]:

$$f(\mathbf{x}) = \frac{(\mathbf{x} - \boldsymbol{\varepsilon})^{\lambda - 1}}{\boldsymbol{\alpha}^{\lambda} \Gamma(\lambda)} \exp\left(-\frac{\mathbf{x} - \boldsymbol{\varepsilon}}{\boldsymbol{\alpha}}\right) \tag{A1}$$

$$
\omega\_p = \varepsilon + \mathfrak{a}t\_p(\lambda) \tag{A2}
$$

where Γ(*λ*) = <sup>∞</sup> <sup>0</sup> *t <sup>λ</sup>*−<sup>1</sup> exp(−*t*)*dt* is Euler's gamma function; *<sup>x</sup>* is an observation of the random variable *X*; *xp* is a quantile of the theoretical GA distribution; and *tp*(*λ*) is a quantile of the standardized gamma distribution, with probability of exceedance *p*.

The WE distribution is defined as [68]:

$$f(\mathbf{x}) \ = \frac{\lambda}{\alpha} \left( \frac{\mathbf{x} - \boldsymbol{\varepsilon}}{\alpha} \right)^{\lambda - 1} \exp \left[ - \left( \frac{\mathbf{x} - \boldsymbol{\varepsilon}}{\alpha} \right)^{\lambda} \right] \tag{A3}$$

$$\mathbf{x}\_{\mathcal{P}} = \mathfrak{a} [-\ln(1 - (1 - p))]^{\frac{1}{\Lambda}} + \mathfrak{e} \tag{A4}$$

The LGA distribution [69] is represented by the equations:

$$f(\mathbf{x}) = \frac{(\ln \mathbf{x} - \ln \varepsilon)^{\lambda - 1}}{\alpha^{\lambda} \Gamma(\lambda) \mathbf{x}} \exp \left( - \frac{\ln \mathbf{x} - \ln \varepsilon}{\alpha} \right) \tag{A5}$$

$$\mathfrak{x}\_{\mathcal{P}} = \varepsilon \exp\left[\mathfrak{a}t\_{\mathcal{P}}(\lambda)\right] \tag{A6}$$

The log-normal distribution (LN) [70] is represented as:

$$f(\mathbf{x}) = \frac{1}{(\mathbf{x} - \boldsymbol{\varepsilon})\alpha\sqrt{2\pi}} \exp\left[ -\frac{1}{2} \left( \frac{\ln(\mathbf{x} - \boldsymbol{\varepsilon}) - \mu}{\alpha} \right)^2 \right] \tag{A7}$$

$$\alpha\_p = \exp\left[\mu + \frac{a\sqrt{2}}{\text{erf}(2(1-p)-1)}\right] + \varepsilon \tag{A8}$$

where: erf( ... ) is the Gauss error function, and other symbols have the same meanings as above, except that *xp* denotes a quantile of the theoretical WE, LGA, and LN distributions, respectively.

The Gumbel distribution [71] is written as:

$$f(\mathbf{x}) = \frac{1}{\mathfrak{a}} \exp\left[ -\frac{\mathbf{x} - \boldsymbol{\mu}}{\mathfrak{a}} - \exp\left( -\frac{\mathbf{x} - \boldsymbol{\mu}}{\mathfrak{a}} \right) \right] \tag{A9}$$

$$\mathbf{x}\_{\mathcal{P}} = -\mathfrak{a}\,\,\ln[-\ln(1-p)] + \mu\tag{A10}$$

where *xp* is a quantile of the theoretical G distribution.

#### **Appendix B. The Goodness-of-Fit Tests**

The following are nonparametric goodness-of-fit tests used to test the goodness of fit of a mathematical model (theoretical distribution) with observations (empirical distribution). The Kolmogorov-Smirnov statistic DK−<sup>S</sup> [46]:

$$\text{D}\_{\text{K}-S} = \max\_{1 < i \le n} (\delta\_i)\_{\text{'}} \\ \text{gdzie } \therefore \delta\_i = \max \left[ \frac{i}{n} - F\_0(\mathbf{x}\_i; \boldsymbol{\theta}), \ F\_0(\mathbf{x}\_i; \boldsymbol{\theta}) - \frac{i-1}{n} \right] \tag{A11}$$

where *n* is the size of the random sample, and *F*<sup>0</sup> *xi*; ˆ *θ* is the distribution function of the theoretical probability distribution for the estimated parameter vector ˆ *θ*.

The Anderson-Darling statistic DA−<sup>D</sup> [48]:

$$\mathbf{D}\_{\mathbf{A}-\mathbf{D}} = -n - \frac{1}{n} \sum\_{i=1}^{n} \left\{ (2i - 1) \ln \mathbf{F}\_0 \left( \mathbf{x}\_i; \hat{\theta} \right) + (2n + 1 - 2i) \ln \left( 1 - \mathbf{F}\_0 \left( \mathbf{x}\_{n+1-i}; \hat{\theta} \right) \right) \right\} \tag{A12}$$

The Liao-Shimokawa statistic DL−<sup>S</sup> [49]:

$$\mathbf{D}\_{\rm L-S} = \frac{1}{\sqrt{n}} \sum\_{i=1}^{n} \frac{\max\left[\frac{i}{n} - F\_0(\mathbf{x}\_i; \boldsymbol{\theta}), F\_0(\mathbf{x}\_i; \boldsymbol{\theta}) - \frac{i-1}{n}\right]}{\sqrt{F\_0(\mathbf{x}\_i; \boldsymbol{\hat{\theta}}) \left[1 - F\_0(\mathbf{x}\_i; \boldsymbol{\hat{\theta}})\right]}} \tag{A13}$$

The Kuiper statistic DK [50]:

$$\mathbf{D}\_{\mathbf{K}} = \max\_{1 \le i \le n} \left( \hat{\delta}\_i^+ \right) + \max\_{1 \le i \le n} \left( \hat{\delta}\_i^- \right) \tag{A14}$$

$$\text{where } \delta\_i^+ = \max\left[\frac{i}{n} - F\_0(\mathbf{x}\_i; \theta)\right]; \delta\_i^- = \max\left[F\_0(\mathbf{x}\_i; \theta) - \frac{i-1}{n}\right].$$

#### **Appendix C. Formulas Used in the Probabilistic Model of Maximum Precipitation of Bogdanowicz and Stachý Model**

The Weibull probability distribution (extreme value type 3, EV3), *f*(*x*), and quantile of maximum precipitation *xp* are given as follows [1,2]:

$$f(\mathbf{x}) = \frac{\lambda}{\theta - \varepsilon} \left[ \frac{\mathbf{x} - \varepsilon}{\theta - \varepsilon} \right]^{\lambda - 1} \exp \left\{ - \left[ \frac{\mathbf{x} - \varepsilon}{\theta - \varepsilon} \right]^{\lambda} \right\} \tag{A15}$$

$$
\omega\_p = \varepsilon + \mathfrak{a} (-\ln p)^{\frac{1}{\lambda}} \tag{A16}
$$

where *ε* is the lowest bound; *ε*(*τ*) = 1.42*τ*0.33; *θ* is the quantile with probability of exceedance 1/e = 0.367 . . . ; *λ* is a shape parameter; and *α* = *θ* − *ε* is a scale parameter.

#### **References**

