Estimating Maximum Daily Precipitation in the Upper Vistula Basin , Poland

The aim of this study was to determine the best probability distributions for calculating the maximum annual daily precipitation with the specific probability of exceedance (Pmaxp%). The novelty of this study lies in using the peak-weighted root mean square error (PWRMSE), the root mean square error (RMSE), and the coefficient of determination (R2) for assessing the fit of empirical and theoretical distributions. The input data included maximum daily precipitation records collected in the years 1971–2014 at 51 rainfall stations from the Upper Vistula Basin, Southern Poland. The value of Pmaxp% was determined based on the following probability distributions of random variables: Pearson’s type III (PIII), Weibull’s (W), log-normal, generalized extreme value (GEV), and Gumbel’s (G). Our outcomes showed a lack of significant trends in the observation series of the investigated random variables for a majority of the rainfall stations in the Upper Vistula Basin. We found that the peak-weighted root mean square error (PWRMSE) method, a commonly used metric for quality assessment of rainfall-runoff models, is useful for identifying the statistical distributions of the best fit. In fact, our findings demonstrated the consistency of this approach with the RMSE goodness-of-fit metrics. We also identified the GEV distribution as recommended for calculating the maximum daily precipitation with the specific probability of exceedance in the catchments of the Upper Vistula Basin.


Introduction
The analysis of the maximum annual daily precipitation (P max ) is one of the crucial factors for the management of water resources in a catchment [1].By knowing the probability and frequency of such precipitation, decision-makers are able to mitigate the effects of floods.In fact, determining the maximum annual daily precipitation with a specific probability of exceedance (P maxp% ) is necessary to tackle the issues of excessive or deficient precipitation and, thus, to meet regional water needs.The estimation of P maxp% is also important for determining flood hazard zones, for designing specific hydraulic structures or systems [2] or for assessing the performance of sewage treatment plants under variable weather conditions [3].However, the growing atmospheric content of greenhouse gases and the resulting global warming we have been experiencing for a few decades significantly affected the precipitation structure (height and spatial variability), thus making the prediction of P maxp% increasingly difficult [4][5][6].
Based on a sufficiently long observation series of the maximum annual daily precipitation, it is possible to predict P maxp% by using probability distributions of random variables.In hydrological studies, P maxp% is commonly determined by using the following probability distribution functions: normal distribution, log-normal distribution, Pearson's type III, exponential function, Gumbel's distribution, generalized extreme value distribution (GEV), and Weibull's and Pareto's distributions [7][8][9][10][11][12][13]. Wdowikowski et al. [14] recommend the generalized exponential distribution for such a purpose.In [15,16], it is shown that hydrological records of typical length (some decades) may display a distorted picture of the actual distribution, suggesting that the Gumbel distribution is not an appropriate model for rainfall extremes.In addition, it is therein demonstrated that the extreme value distribution of type II (EV2) is a more consistent alternative.Selecting an appropriate form of probability distribution involves matching a proper theoretical function with an empirical distribution of a random variable.The selection of an inappropriate function for a region may result in underestimating or overestimating the maximum hydro-meteorological events with a specific probability of exceedance [17].
The application of model selection criteria within the field of frequency analysis of hydrological extremes is rare [18].The only methods are Akaike (AIC) or Schwartz (BIC) information criteria [19].Also, alternative statistical tests could be employed for assessing the performance of theoretical distributions, such as, for instance, the Anderson-Darling, which compares the whole range but gives more weight to the upper tail [20].The main disadvantage of these goodness-of-fit metrics is that they do not provide information on the accuracy limit and on the inclusion of atypical (outlying) values.Moreover, it is not possible to assess the value of the test statistics.On the other hand, methods commonly used for assessing the quality of rainfall-runoff hydrological models, such as the percentage error in peak flow, percentage error in volume, efficiency coefficient, peak-weighted root mean square error, the sum of absolute residuals, or the sum of squared residuals, seem promising alternatives.In fact, some of them, e.g., the efficiency coefficient, can be assessed against a scale indicating the quality of the results yielded by the models [21][22][23][24].Another option is a known but infrequently used goodness-of-fit metric based on the probability plot correlation coefficient [25,26].This approach allows for determining the strength of the relationship between the observed and calculated variables using statistical distributions, e.g., the Guilford's scale [27].It also shows the statistical significance of the resulting correlation coefficients, thus providing more information on the quality of fit of the probability distributions.
Previous studies focusing on the form of the probability distribution for P maxp% estimation only tend to achieve the best fit of the theoretical and empirical distributions and to describe the spatial variability of the forms of the probability distributions [28][29][30].The application of alternative methods (e.g., those used for rainfall-runoff modeling) for assessing the quality of fit of statistical and empirical distribution is not documented in the literature.Hence, the aim of this study was to determine the best probability distributions used for assessing P maxp% .The novelty of this work consists of identifying the best-fitting theoretical function by adopting the peak-weighted root mean square error (PWRMSE), commonly used for calibrating parameters in hydrological models.Additionally, the root mean square error (RMSE) and the coefficient of determination (R 2 ) were used as the goodness-of-fit metrics.Finally, the study aimed at determining the trend significance for the P max observation series and at estimating the density function for the stations showing a significant trend.

Description of the Study Area
The study area is the Upper Vistula Basin, which accounts for about 25% of the entire Vistula basin and about 15% of Poland's area, and occupies the southern part of the country.It covers part of the Carpathians, the Subcarpathian Valleys and Małopolska Uplands.The study area shows large variations in altitude, and this influences the height of precipitation and its extreme values [31].

Methods
The study involved the following elements: (i) preliminary analysis of the annual maximum precipitation time series, (ii) analysis of trends in the observed series, (iii) determination of the maximum annual daily precipitation with a specific probability of exceedance, (iv) selection of the best fit between theoretical and empirical distribution of random variables.

Preliminary Analysis of the Observation Series
The preliminary analysis of the annual maximum precipitation included the calculation of descriptive statistics.Specifically, we computed minimum (LPmax), mean (MPmax), and maximum (HPmax) values; measures of dispersion-standard deviation (s) and coefficient of variation (Cs); and measure of shape of the studied variate distribution-the coefficient of skewness (A).

Analysis of Trends in the Observed Series
Analysis of trends in the observed series was conducted by a Mann-Kendall test.This test was selected due to its advantages with respect to alternative tests.First, in order to apply the Mann-Kendall test, the data need not conform to any particular distribution.Second, the test exhibits low sensitivity to abrupt breaks due to an inhomogeneous time series [32].The significance level was set to α = 0.05.The H0 hypothesis of the Mann-Kendall test assumes a lack of a monotonic trend for the data, while the alternative H1 claims its existence.The S statistics of the Mann-Kendall was based on the following formula [33][34][35][36][37][38]:

Methods
The study involved the following elements: (i) preliminary analysis of the annual maximum precipitation time series, (ii) analysis of trends in the observed series, (iii) determination of the maximum annual daily precipitation with a specific probability of exceedance, (iv) selection of the best fit between theoretical and empirical distribution of random variables.

Preliminary Analysis of the Observation Series
The preliminary analysis of the annual maximum precipitation included the calculation of descriptive statistics.Specifically, we computed minimum (LP max ), mean (MP max ), and maximum (HP max ) values; measures of dispersion-standard deviation (s) and coefficient of variation (C s ); and measure of shape of the studied variate distribution-the coefficient of skewness (A).

Analysis of Trends in the Observed Series
Analysis of trends in the observed series was conducted by a Mann-Kendall test.This test was selected due to its advantages with respect to alternative tests.First, in order to apply the Mann-Kendall test, the data need not conform to any particular distribution.Second, the test exhibits low sensitivity to abrupt breaks due to an inhomogeneous time series [32].The significance level was set to α = 0.05.The H 0 hypothesis of the Mann-Kendall test assumes a lack of a monotonic trend for the data, while the alternative H 1 claims its existence.The S statistics of the Mann-Kendall was based on the following formula [33][34][35][36][37][38]: where: n-the number of elements of the time series, x j -observation at time j, x k -observation at time k.
On the basis of the standard statistics Z calculated according to the formula: Var(S)-variance S, determined on the basis of the formula: where: n-the number of elements of the time series.
If the value of Z is lower than the critical value of Z crit for the assumed significance level α, the hypothesis claiming the lack of a trend is acceptable.Otherwise, the H 0 hypothesis should be discarded in favor of the alternative hypothesis.The main assumption of the used Mann-Kendall test is the lack of autocorrelation in the data series.In the case of the analysis of the maximum annual daily precipitation, such dependencies may occur, which leads to an underestimation of the variance Var (S).Therefore, an adjustment for variance correction is included that is calculated only for data with significant autocorrelation in general [39]: where: n * s -the effective number of observations calculated as: where: ρ k -the lag k autocorrelation coefficients of the ranks of the observations.

Identification of Empirical Distributions with Kernel Estimators
For the rainfall stations showing a significant trend, a direct estimation of the unknown density function was made using the kernel estimators.The kernel estimation of the density function is a commonly used tool for analyzing the behavior of hydrological phenomena, including atmospheric precipitation, as evidenced by the works of numerous research teams [40][41][42].The outcomes allowed us to estimate the variability of the maximum annual daily precipitation.The estimators ( fh (x)) were determined based on the value of an n-element random sample according to the following formula [43,44]: where: n-sample size; h-smoothing parameter equated with so-called band width; K-kernel function; X i -i-th element of the sample.
The band width H was established as per Silverman method [45] and the kernel density function K was regarded as a Gaussian form of kernel [46].

Maximum Annual Daily Precipitation with a Specific Probability of Exceedance
The value of P maxp% was determined based on the following probability distributions: Pearson's type III (PIII), Weibull's (W), log-normal (LN), generalized extreme value (GEV), and Gumbel's (G).The rainfall stations showing a significant trend were excluded from the P maxp% calculations.Quantiles of maximum annual daily precipitation were calculated as per the following formulas [47,48]: Pearson's type III distribution: Weilbull's distribution: Log-normal distribution: Generalized extreme value (GEV) distribution: Gumbel's distribution: where: κ, λ, β-shape parameters; t(λ)-standardized variable; α-scale parameter; ε, ξ-location parameter; µ, σ-parameters of log-normal; p-probability of exceedance; u p -the standard normal variate of probability of exceedance p.
The values of the distribution parameters can be obtained using various methods, e.g., maximum likelihood, L-moments, and methods of moments.In this work, they were estimated by means of the maximum likelihood method.The maximum likelihood method is characterized by asymptotically unbiased and optimal (with the smallest variance) parameter estimators if the assumed estimation model is true [49].The L-moments method is a linear combiner of the sample.Typically, this estimation method gives comparable results to the maximum likelihood.The disadvantage of the L-moments method is that it cannot be applied to probability distributions that do not have explicit expressions per quantile.In addition, the formulas for the calculation of linear moments require attempts ordered in a non-trivial way, and this destroys the chronological order of events, thus hampering hydrological analysis of non-stationary processes.In the method of moments, the systematic error quickly increases with the degree of moments since the sample elements are raised to the power of the second and third when we want to estimate three parameters of the distribution.In addition, quantiles with a greater return period are underestimated when estimated by the method of moments [50].It should be mentioned that the distribution parameters are subject to uncertainty.As stated in [51], even when we have abundant, good-quality data to work with and a good model, our parameter estimates are still subject to a standard error.Although general guidance is available on how parameter uncertainty should be accounted for in probabilistic sensitivity analysis, there is no comprehensive guidance on the estimation of uncertainty in the parameters of the distributions used to represent stochastic uncertainty in statistical models [52].Therefore, to assess the consistency of the theoretical distribution with the empirical distribution, the Kolmogorov-Smirnov (K-S) test was adopted [53].If the K-S test indicated compliance for each analyzed theoretical distribution, then the best-fitted distribution goodness-of-fit measures were used.

Selection of the Theoretical Function Best Fitting The Empirical Distribution and Sensitivity to Outliers
The theoretical function best fitting the empirical distribution of the random variable was identified based on the metric used in the assessment, e.g., in hydrological modeling: Root mean square error (RMSE), coefficient of determination (R 2 ), and peak-weighted root mean square error (PWRMSE) The analyzed goodness-of-fit metrics are described by the following relationships [54][55][56][57]: where: n-size of the observation series; e i -difference between the observed and estimated value of the maximum daily precipitation for year i; O i -observed values for year i; P i -predicted values for year i; O-mean of observed values; P-mean of predicted values.
Finally, based on results derived from the analyzed metrics, we also recommended the most suitable form of the probability distribution for the determination of the maximum annual total daily precipitation with a specific probability of exceedance.The most suitable form of probability distribution was indicated using the rank method.The range of rank was from one to five (five probability distributions were analyzed).The rank no. 1 was for the best fitted distribution in a particular station, rank no. 2 for the second in order, up to rank no. 5 for the distribution with the poorest fit.Next, for all stations, the sum of ranks was computed.The best fitted distribution in the whole Upper Vistula Basin was the one with the lowest value of the rank sum.

Preliminary Analysis of the Annual Maximum Precipitation
Table 1 presents the values of the descriptive statistics for the annual maximum precipitation time series for all the stations.1, it is concluded that the coefficient of variation C s mostly remained below 40%, which indicates an average variability of P max in the time series.The lowest value of C s was for Szczawne (25%) and the highest for Kocierz Moszczanicki (61%).The values of the coefficient of skewness (A) were above zero in each time-series, which proved the right-sided asymmetry of the empirical distributions of the variates.This is due to the fact that for the analyzed time-series the majority of the series means observations are lower than the median values.

Analysis of Trends in the Observed Series
The observation series of the maximum annual daily precipitation were analyzed for significant monotonic trends with the Mann-Kendall tests.The findings are presented in Figure 2. Figure 3 presents a time series of the annual maxima daily precipitation for rainfall stations with a significant trend.

Analysis of Trends in the Observed Series
The observation series of the maximum annual daily precipitation were analyzed for significant monotonic trends with the Mann-Kendall tests.The findings are presented in Figure 2. Figure 3 presents a time series of the annual maxima daily precipitation for rainfall stations with a significant trend.A significant (only positive) trend indicating long-term change at the assumed significance level was found for seven out of 51 rainfall stations (that is, Rajcza, Radziszów, Książ Wielki, Wisłok Wielki, Radomyśl Wielki, Cisna, and Komańcza).Figure 3 presents the time series of the annual maxima daily precipitation for rainfall stations with significant trend along with the indication of generation mechanisms for the greatest values.The stations with a significant trend were located in the mountainous (Rajcza, Wisłok Wielki, Cisna, Komańcza) and highland (Radziszów, KsiążWielki, Radomyśl Wielki) parts of the Upper Vistula Basin.Also, it should be mentioned that Wisłok Wielki, Cisna, Komańcza, and Radomyśl Wielki are placed in the Bieszczady, which are located in the A significant (only positive) trend indicating long-term change at the assumed significance level was found for seven out of 51 rainfall stations (that is, Rajcza, Radziszów, Ksi ą ż Wielki, Wisłok Wielki, Radomyśl Wielki, Cisna, and Koma ńcza).Figure 3 presents the time series of the annual maxima daily precipitation for rainfall stations with significant trend along with the indication of generation mechanisms for the greatest values.The stations with a significant trend were located in the mountainous (Rajcza, Wisłok Wielki, Cisna, Koma ńcza) and highland (Radziszów, Ksi ą żWielki, Radomyśl Wielki) parts of the Upper Vistula Basin.Also, it should be mentioned that Wisłok Wielki, Cisna, Koma ńcza, and Radomyśl Wielki are placed in the Bieszczady, which are located in the mid-mountain part of the Carpathian.As can be seen in Figure 3d-g, most of the highest observations for these stations were for the years 1997 and 2010.Such years in Poland were characterized by extremely high precipitation, with the highest precipitation being observed in the south-east part of Poland.Such high precipitation was caused by incoming moist air masses from Ukraine and the Great Hungarian Plain.The other rainfall stations showed no significant trends in the observation series.This is consistent with studies on the behavior of the highest total daily precipitation in Southern Poland that reported irregular fluctuations rather than a trend [58].Niedźwiedź et al. [59], who also investigated the area of Southern Poland, obtained similar results usually indicating a lack of significant trends.Our findings are also consistent with trend results for maximum flows in the Upper Vistula Basin.As shown by studies conducted in this region over the past decades [60], these flows were characterized by the lack of any trend in their values.The relationship between the maximum annual daily precipitation and the maximum flows results from the physiographic characteristics of the Upper Vistula Basin that affect the rainfall-runoff transformation [61].
Ukraine and the Great Hungarian Plain.The other rainfall stations showed no significant trends in the observation series.This is consistent with studies on the behavior of the highest total daily precipitation in Southern Poland that reported irregular fluctuations rather than a trend [58].Niedźwiedź et al. [59], who also investigated the area of Southern Poland, obtained similar results usually indicating a lack of significant trends.Our findings are also consistent with trend results for maximum flows in the Upper Vistula Basin.As shown by studies conducted in this region over the past decades [60], these flows were characterized by the lack of any trend in their values.The relationship between the maximum annual daily precipitation and the maximum flows results from the physiographic characteristics of the Upper Vistula Basin that affect the rainfall-runoff transformation [61].

Identification of Empirical Distributions with Kernel Estimators
The statistical analysis was expanded by using kernel density estimators to evaluate the distribution of the maximum annual daily precipitation for the stations showing significant trends in the observation series.Results are presented in Figure 4.
The kernel estimation of the density function for the maximum annual daily precipitation revealed a multimodal nature of the density function for all cases exhibiting trends.This indicates a subpopulation within the investigated observation series, thus suggesting a change in the weather mechanisms that trigger intense rainfall, such as high content of water vapor, temperature

Identification of Empirical Distributions with Kernel Estimators
The statistical analysis was expanded by using kernel density estimators to evaluate the distribution of the maximum annual daily precipitation for the stations showing significant trends in the observation series.Results are presented in Figure 4.

Determination of the Maximum Annual Daily Precipitation with Specific Probability of Exceedance
The maximum annual daily precipitation with a specific probability of exceedance was calculated using four distributions.The calculations were carried out for the rainfall stations where the observation series did not show a significant trend.The Kolmogorov-Smirnov test for significance level α = 0.05 confirmed the consistency of the analyzed theoretical distributions with the empirical distributions of the random variables in all tested cases.Figure 5 presents the values of Pmax10% determined using the analyzed statistical distributions.Figure 6 presents the example Pmaxp% results obtained by using the theoretical distribution for Wisła-Malinka station.The kernel estimation of the density function for the maximum annual daily precipitation revealed a multimodal nature of the density function for all cases exhibiting trends.This indicates a subpopulation within the investigated observation series, thus suggesting a change in the weather mechanisms that trigger intense rainfall, such as high content of water vapor, temperature differences between the incoming and lingering air masses, time of low pressure persistence over the given area, thermodynamic equilibrium of the atmosphere and local conditions [62].This was also confirmed by significant trends of P max for the selected observation series.For six stations where multimodality was demonstrated on the basis of the density kernel analysis-Figure 4 (Rajcza, Radziszów, Wisłok Wielki, Radomysl Wielki, Cisna and Koma ńcza)-the circulation situation was analyzed, see Figure 3. Details on the circulation situation analysis were presented by Mły ński et al. [63].In 45.8% of the highest P max values, the situation corresponded to Nc (cyclonic low-pressure circulation with air advection from the north), NEc (cyclonic low-pressure circulation with air advection from the northeast) and Bc (trough, a low-pressure area, or an axis of a low-pressure trough with various advection vectors), but 37.5% of them were caused by Nc and NEc circulations.In the case of the rest of the P max values (smaller than the value of the highest values), they were caused by Cc (central cyclonic, center of low pressure -37% of all rainfall for the analyzed stations).The presented results show a different genesis for the formation of the maximum rainfall since the highest values are in relation to the other values.Additionally, it must be emphasized that extreme values are not necessarily clustered in recent years, and it cannot be said that the increasing trend could be attributed only to the occurrence of these extremes (the trends should also be attributed to an increase in the mid-range values).

Determination of the Maximum Annual Daily Precipitation with Specific Probability of Exceedance
The maximum annual daily precipitation with a specific probability of exceedance was calculated using four distributions.The calculations were carried out for the rainfall stations where the observation series did not show a significant trend.The Kolmogorov-Smirnov test for significance level α = 0.05 confirmed the consistency of the analyzed theoretical distributions with the empirical distributions of the random variables in all tested cases.Figure 5 presents the values of P max10% determined using the analyzed statistical distributions.Figure 6 presents the example P maxp% results obtained by using the theoretical distribution for Wisła-Malinka station.The values of Pmax10% and especially, quantiles of a smaller probability of exceedance are commonly used in hydrological engineering.This precipitation is considered reliable for designing rainwater draining systems in rural areas.Our calculations revealed that the highest values of Pmax10% were most commonly obtained from the log-normal distribution (58% of all rainfall stations).The lowest values of Pmax10% were usually derived from the Gumbel's distribution (84% of all rainfall stations).The highest Pmax10% quantile values obtained using the log-normal distribution are justified by the properties of such a model.In fact, the log-normal is one of the heavy-tailed distributions, which means that with the same order of the upper quantile, e.g., probability p ≤ 0.2, it generates much higher quantile values than other probability distributions [64].The GEV distribution shows similar properties; hence, it also yields high values for low quantiles, such as for the Zawoja station.Moreover, it is clear that quantile estimates of maximum precipitation for lower frequency have a higher difference between distributions than for higher frequency.Therefore, it is important to apply the goodness-of-fit metrics that are more sensitive to differences in quantiles with a low frequency (heavy-tailed distribution).The values of P max10% and especially, quantiles of a smaller probability of exceedance are commonly used in hydrological engineering.This precipitation is considered reliable for designing rainwater draining systems in rural areas.Our calculations revealed that the highest values of P max10% were most commonly obtained from the log-normal distribution (58% of all rainfall stations).The lowest values of P max10% were usually derived from the Gumbel's distribution (84% of all rainfall stations).The highest P max10% quantile values obtained using the log-normal distribution are justified by the properties of such a model.In fact, the log-normal is one of the heavy-tailed distributions, which means that with the same order of the upper quantile, e.g., probability p ≤ 0.2, it generates much higher quantile values than other probability distributions [64].The GEV distribution shows similar properties; hence, it also yields high values for low quantiles, such as for the Zawoja station.Moreover, it is clear that quantile estimates of maximum precipitation for lower frequency have a higher difference between distributions than for higher frequency.Therefore, it is important to apply the goodness-of-fit metrics that are more sensitive to differences in quantiles with a low frequency (heavy-tailed distribution).

Selection of the Best Fit between the Theoretical and Empirical Distribution of Random Variables
The analysis of fit between the theoretical and empirical distributions of the maximum annual daily precipitation was based on RMSE, R 2 , and the adapted PWRMSE goodness-of-fit metrics.The results of these calculations are shown in Table 2.
According to the PWRMSE goodness-of-fit metrics, the values presented in Table 2 identified the log-normal distribution as the most suitable probability distribution function for estimating the quantiles of P maxp% .This was confirmed for 41% of all rainfall stations.The other most suitable distributions according to this goodness-of-fit measures were GEV (29% of all rainfall stations), Pearson's type III (16%), and Weibull's (7%) and Gumbel's (7%).The RMSE goodness-of-fit metrics yielded the same pattern of results.Additionally, both the PWRMSE and the RMSE goodness-of-fit metrics indicated the same best-fitting theoretical distributions of the probability functions for all the analyzed stations.The value of PWRMSE was nearly always higher than that of RMSE (on average by 15%).This is due to the higher weight attributed to values above the average in the PWRMSE goodness-of-fit metrics, whereas, for the RMSE goodness-of-fit metrics, every value has the same weight.The values of the R 2 goodness-of-fit metrics revealed a very strong determination between the empirical data and the values derived from statistical distributions.The coefficient of determination R 2 identified GEV as the best-fitting distribution type (62% of all rainfall stations), followed by log-normal (18%), and Weibull's (11%), Gumbel's (7%), and Pearson's type III distributions (2%).Our study identified the log-normal and GEV functions (the same family of distributions) as the recommended statistical distributions for the determination of the P maxp% quantiles.This is evidenced by the results of inference carried out using PWRMSE for assessing the quality of the statistical distributions.In fact, this goodness-of-fit measure indicated the best fit of these distributions in the largest number of cases.While studies on the best-fit probability distributions for the estimation of P maxp% quantiles have not been conducted on a wider scale neither in Poland nor in Central Europe and there is no detailed research in this study area in southern Poland, a report by Cebulska [65] also indicated the log-normal and Weibull's distributions as the best for the Orava-Nowy Targ Basin (the Carpathian part of the Upper Vistula Basin).Another study [66] suggested Weibull's and log-gamma distribution as the best-fitting ones for central Poland.The research carried out by Wdowikowski et al. [14] confirmed the generalized exponential distribution as the best measure for estimating P maxp% in the Upper Odra Basin.Villarini [67] claimed the GEV function to be the best distribution to estimate the maximum annual daily precipitation with a specific probability of exceedance in the countries of Central Europe due to its upper tail.Despite such studies, it should be noted that the form of the probability distributions (and particularly their parameters) are closely related to the areas (physical, geographical, and meteorological characteristics affecting the precipitation) for which the maximum daily precipitation with a specific probability of exceedance is determined [28].

Conclusions
Our calculations showed a lack of significant trends in the observation series of the investigated random variables for a majority of rainfall stations in the Upper Vistula Basin (44 out of 51).The multimodality of the empirical density function for the rainfall stations with significant trends confirmed the change in meteorological mechanisms in the surveyed multi-year period that affected the maximum daily precipitation in the Upper Vistula Basin.We found that the peak-weighted root mean square error (PWRMSE) goodness-of-fit metrics, commonly used in the quality assessment of rainfall-runoff models, may be used for identifying the best-fit statistical distributions.This was also confirmed by the root mean square error (RMSE) goodness-of-fit metrics.We also identified the log-normal and generalized extreme value GEV distributions as suitable for calculating the maximum daily precipitation with a specific probability of exceedance in the catchments of the Upper Vistula Basin (according to the PWRMSE goodness-of-fit metrics, the log-normal was the best for 41% and the GEV for 29% of all the stations; according to the RMSE goodness-of-fit metrics, consistent results were found; according to coefficient of determination R 2 , the GEV was optimal for 62% and the log-normal for 18% of all the stations).Since the obtained results from PWRMSE goodness-of-fit metrics, with regard to the best-fitted distribution, were similar to RMSE and R 2 , it was concluded that the analyzed methods can be used interchangeably.In addition, it was concluded that PWRMSE can be used as the goodness-of-fit metric instead of other goodness-of-fit measures.
The input data were the maximum total daily precipitation recorded in the years 1971-2014 at 51 rainfall stations in the Upper Vistula Basin.The data were obtained from the public database of the Institute of Meteorology and Water Management, National Research Institute in Warsaw.The Upper Vistula Basin and rainfall stations are depicted in Figure1.Atmosphere 2018, 9, x FOR PEER REVIEW 3 of 18 stations in the Upper Vistula Basin.The data were obtained from the public database of the Institute of Meteorology and Water Management, National Research Institute in Warsaw.The Upper Vistula Basin and rainfall stations are depicted in Figure 1.

Figure 1 .
Figure 1.The digital elevation model of Upper Vistula Basin and rainfall stations location (black markers) on the map of Poland.

Figure 1 .
Figure 1.The digital elevation model of Upper Vistula Basin and rainfall stations location (black markers) on the map of Poland.

Figure 2 .
Figure 2. Results of the Z statistics for the Mann-Kendall test.

Figure 2 .
Figure 2. Results of the Z statistics for the Mann-Kendall test.

18 Figure 5 .
Figure 5. Values of Pmax10% determined for the rainfall stations using the analyzed statistical distributions.

Figure 5 .
Figure 5.Values of P max10% determined for the rainfall stations using the analyzed statistical distributions.

Figure 5 .
Figure 5. Values of Pmax10% determined for the rainfall stations using the analyzed statistical distributions.

Table 1 .
The characteristics of the annual maximum precipitation in the analyzed period.
LP max -minimum values of precipitation, MP max -mean values of precipitation, HP max -maximum values of precipitation, C s -coefficient of variation, s-standard deviation, A-the coefficient of skewness Based on Table

Table 2 .
Values of the analyzed goodness-of-fit measures matching the theoretical distributions.-peak-weighted root mean square error, RMSE-the root mean square error, R 2 -coefficient of determination, PIII-Pearson's type III distribution, W-Weibull distribution, GEV-generalized extreme value distribution, G-Gumbel distribution. PWRMSE