Depth–Duration–Frequency Relationship Model of Extreme Precipitation in Flood Risk Assessment in the Upper Vistula Basin

: The Upper Vistula Basin is a ﬂood-prone region in the summer season (May–October) due to intensive rainfall. From the point of view of water management, it is particularly important to assess the variability in this main factor of ﬂood risk, as well as to establish the depth–duration– frequency (DDF) relationship for maximum precipitation, this having not yet been derived for the region. The analysis of a 68-year (1951–2018) data series of summer maximum precipitation collected by 11 meteorological stations showed the series’ stationarity, which supports the conclusion that there is no increase in the risk of rainfall ﬂoods due to the intensiﬁcation of extreme precipitation. A new approach is proposed for the determination of the DDF relationship, where the best-ﬁtted distribution for each station is selected from among the set of candidate distributions, instead of adopting one ﬁxed distribution for all stations. This approach increases the accuracy of the DDF relationships for individual stations as compared to the commonly used approach. In particular, the traditionally used Gumbel distribution turns out to be not well ﬁtted to the investigated data series, and the advantage of the recently popular GEV distribution is not signiﬁcant.


Introduction
Establishing the depth-duration-frequency (DDF) relationship, or an equivalent intensity-duration-frequency (IDF) relationship for the maximum annual or seasonal rainfall is one of the basic tasks when it comes to water management. For almost a hundred years, the DDF (IDF) curves have been a useful tool for determining how the return levels of maximum rainfall depths (intensities) vary with duration over a range of return periods [1]. The initial formulas were semi-empirical [2,3], while, with the development of statistical modeling of hydrological time series, the general formula for the IDF relationship was introduced by Koutsoyiannis et al. [4], in accordance with the probabilistic theory of hydrological maxima analysis. The basic DDF relation and the underlying structure of the DDF curves have been widely applied and discussed in the literature with respect to maximum precipitation [5][6][7], as well as maximum flows [8][9][10]. One of the modifications to the traditional DDF method is multi-scaling based on the multifractal properties of temporal rainfall [6,11,12]. The assumption about stationary conditions during the DDF (IDF) derivation procedure has also been considered [13][14][15]. This issue is investigated here in relation to the analyzed data, while the general discussion is beyond the scope of this article.
To the best of the author's knowledge, the DDF or IDF relationship for the Upper Vistula Basin has not been determined so far. It is a region in which hydrological processes are highly dynamic and has a high risk of flooding in the summer season, from May to October, when the maximum floods are caused by intense rainfall. The advantages of a seasonal approach to modeling peak flows in Poland and other countries with similar hydrological The predominance of summer rainfall floods can be easily seen. The greatest floods in the observation period on the Upper Vistula occurred in the years 1958, 1960, 1962, 1970, 1972, 1997, 2010, and 2014; see, e.g., [29]. The most catastrophic in terms of extent, flood volume, force, and flood losses were the floods in the years 1997 and 2010. All the above-mentioned flood events were caused by heavy rains lasting several days, the greatest amounts of which occur within 3-5 days [30][31][32].
It is commonly suggested that the large floods in recent years are caused or at least reinforced by human activities, such as regulation and hydro-technical structures in riverbeds, changes in land cover and land use, and growing urbanisation increasing the area of impermeable surfaces. The proofs of water management impacts, such as regularization of riverbeds resulting in shortening (by cutting meanders) and narrowing (by levies and cross dikes construction), on flood formation processes are evident; see, e.g., [33,34]. They can accelerate and magnify the flood peaks. However, this research is out of the scope of this study. A statistical evaluation of the role of anthropogenic factors will result here from a comparison of trends in peak flows and trends in rainfall totals for 1, 3, and 5 days.
The main goal of this study is to answer the research question, whether the main driving climatic force of dominating floods in the Upper Vistula Basin, i.e., the maximum rainfall totals in 1, 3, and 5 days, shows significant temporal changes which can be treated as an increasing climatic flood risk in the region. The detailed objective is to determine the DDF relationship for the maximum rainfall of selected stations from the Upper Vistula Basin. The approach proposed here of choosing the best-fitted probability distribution for particular data series, instead of the common practice of adopting one fixed distribution for all stations in the region, seems appropriate given the conditions of a mountainous area with varying altitudes and sloping terrain, and with a small number of stations mainly measuring rainfall in the valleys. The inference about the spatial variability of quantiles in such terrain is still a challenge for researchers and practitioners.
The paper is organized as follows. After the introduction to the topic, the characteristics of the research area and the measurement data collected by 11 stations for maximum precipitation in the summer season for the period 1951-2018 are presented. The The predominance of summer rainfall floods can be easily seen. The greatest floods in the observation period on the Upper Vistula occurred in the years 1958, 1960, 1962, 1970, 1972, 1997, 2010, and 2014; see, e.g., [29]. The most catastrophic in terms of extent, flood volume, force, and flood losses were the floods in the years 1997 and 2010. All the above-mentioned flood events were caused by heavy rains lasting several days, the greatest amounts of which occur within 3-5 days [30][31][32].
It is commonly suggested that the large floods in recent years are caused or at least reinforced by human activities, such as regulation and hydro-technical structures in riverbeds, changes in land cover and land use, and growing urbanisation increasing the area of impermeable surfaces. The proofs of water management impacts, such as regularization of riverbeds resulting in shortening (by cutting meanders) and narrowing (by levies and cross dikes construction), on flood formation processes are evident; see, e.g., [33,34]. They can accelerate and magnify the flood peaks. However, this research is out of the scope of this study. A statistical evaluation of the role of anthropogenic factors will result here from a comparison of trends in peak flows and trends in rainfall totals for 1, 3, and 5 days.
The main goal of this study is to answer the research question, whether the main driving climatic force of dominating floods in the Upper Vistula Basin, i.e., the maximum rainfall totals in 1, 3, and 5 days, shows significant temporal changes which can be treated as an increasing climatic flood risk in the region. The detailed objective is to determine the DDF relationship for the maximum rainfall of selected stations from the Upper Vistula Basin. The approach proposed here of choosing the best-fitted probability distribution for particular data series, instead of the common practice of adopting one fixed distribution for all stations in the region, seems appropriate given the conditions of a mountainous area with varying altitudes and sloping terrain, and with a small number of stations mainly measuring rainfall in the valleys. The inference about the spatial variability of quantiles in such terrain is still a challenge for researchers and practitioners.
The paper is organized as follows. After the introduction to the topic, the characteristics of the research area and the measurement data collected by 11 stations for maximum precipitation in the summer season for the period 1951-2018 are presented. The existing DDF (IDF) method is summarized, and its modification is proposed in the next section. The following section presents the results of studies on the stationarity of the analyzed data series, followed by the results of the construction of the depth-duration-frequency relation obtained with the new method. Discussion of the results and a comparison of the classical and new approaches are presented in the next section. Finally, the studies are summarized in the conclusion.
All calculations in this paper were made with the use of Fortran and Python software supported by Excel software.

Study Area and Data Sets
The research area covers the western part of the Upper Vistula River Basin. The total area of the Upper Vistula Basin is 33,458 km 2 (including 1,952.4 km 2 outside the Polish territory) and refers to the 399 km long section from the source to the mouth of the river San. The studied region covers the area of the Polish Western Carpathians and part of the Podkarpacie region. According to the commonly used Koppen climate classification system [35], the upland parts of the basin have a predominantly humid continental climate with mild summers and precipitation all year round. As the height above sea level increases, the sub-arctic climate begins to dominate with harsh winters, no dry season, and cool summers, and in the highest mountain parts a polar tundra climate prevails [36,37].
The research focuses on the Upper Vistula Basin, as this area has a high flood potential. The main threat here is summer floods caused by intense, long-lasting rainfall, while the floods in the Middle Vistula and the Lower Vistula mainly indicate the regime of snowmelt floods or mixed snowmelt floods; see, e.g., [19,38]. The share of the annual runoff from the Upper Vistula Basin in the runoff balance of the Vistula is about 1.7 times greater than the share of the Upper Vistula Basin area in the Vistula Basin area. The voivodeships located in the basin (especiallyŚląskie and Małopolskie) have the highest population density index in Poland; therefore, the flood risk potential, i.e., the number of people and their property, as well as economic potential, is higher than the national average. Flood losses and damages recorded in the Upper Vistula Basin are significant. In 1953-2006, the sum of flood losses in the Upper Vistula water region accounted for nearly 50% of all flood losses in Poland, while the area of this region covers only 15% of the country's area [39].
The article analyzes the series of daily precipitation during the summer half-year (i.e., 1 May-31 October) from the period 1951-2018 collected by 11 meteorological stations. The selected stations have data series of daily precipitation for a continuous period of 68 years. The rainfall data from the measurement network of the Institute of Meteorology and Water Management-National Research Institute (IMGW-PIB) are available freely at [28].
In Table 1 the stations are listed according to their location within the main river and tributary basins. Left (L) and right (R) tributaries of the Vistula are indicated in brackets. Next, the geographical coordinates and altitude are presented, as well as the seasonal mean, along with the maximum value of daily precipitation in the summer half-year for 1951-2018. The locations of the stations is shown in Figure 2.   Poland, a Central European country, is located on the border of the collision of the sea air masses from the Atlantic Ocean (flowing from the west) and the continental air masses from Asia (from the east). Moist air masses from the Baltic Sea reach Poland from the north, but the inflow of hot and dry air from the southern part of Europe is inhibited by mountain ranges. The general spatial variability of annual precipitation is latitudinal. The mean annual amount of precipitation increases from north to south with the general Poland, a Central European country, is located on the border of the collision of the sea air masses from the Atlantic Ocean (flowing from the west) and the continental air masses from Asia (from the east). Moist air masses from the Baltic Sea reach Poland from the north, but the inflow of hot and dry air from the southern part of Europe is inhibited by mountain ranges. The general spatial variability of annual precipitation is latitudinal. The mean annual amount of precipitation increases from north to south with the general slope of the terrain. The lowest totals (c.a. 500 mm) are observed in the central part of the country, the highest (more than 1000 mm) in the high mountains (Tatra) in the Upper Vistula Basin.
The Upper Vistula Basin is situated in the area of three geographical regions. These are: the mountain belt (Karpaty), the Podkarpackie valleys belt, and the highland belt. Especially in the belt of the mountains and valleys, in the right-bank part of the basin, a large variety of terrain is conducive to the occurrence of orographic and convection rainfall occurring alone or strengthening the frontal rainfall.
For each year in the period 1951-2018, the maximum depth of precipitation D in the summer half-year was extracted for three duration periods d for d = 1, 3, and 5 days, giving a 68-year-long observation series D(d).

Formulation of the DDF Curve
The basic requirement that the statistical model of maximum precipitation should meet is the stochastic ordering of distribution functions F(x; d) for different durations d. We say that the random variables X d 1 , . . . , X d n are stochastically ordered if, for each probability F, the quantiles of the order F (x(d i , F)) of these variables satisfy the inequality: In other words, one can say that the distribution of X d j is stochastically lower than the distribution of X d k (for all j, k = 1, . . . , n and j < k), i.e., F x; d j ≥ F(x; d k ) for all x ∈ R + and all j, k; j < k.
Stochastic ordering of maximum precipitation distributions of different duration is related to the need to ensure the logical consistency of the hydrological design, i.e., the expected increase in precipitation, along with its duration, at a fixed F value. A similar problem occurs in other issues of determining reliable hydrological characteristics when their duration counts.
The estimation procedure was based on the method proposed by Javelle et al. [8,9] which uses the mean flows exceeded in consecutive d days to describe the peak parts of a flood hydrograph. This procedure meets the requirement of statistical ordering of probability distributions for growing duration d.
The relationship implemented when approaching the problem of DDF (and IDF) estimation is the formula given by Equation (2), also known as the Montana curve [40]: is the return level of the summer (May-Oct) maximum precipitation depth for duration d and return period T, which is related to the probability of non-exceedance F of the event according to Equation (3); see, e.g., [41]: Since a(T) in Equation (2) doesn't depend on d, the family of DDF curves defined by Equation (2) is parallel for different T values, and the shape of the curves is determined by the denominator. The η parameter affects the slope of the straight part of the DDF curves, while the θ parameter affects the point of curvature change [42]. Both parameters characterize the dynamics of the extreme precipitation process as a function of the duration d and are related to the climate [5].
Let D(d, i) denote the i-th element in the non-descending ordered N-element series of maximum seasonal precipitation for the duration of d, i.e., D(d, 1 : N) ≤ . . . ≤ D(d, N : N). It is easy to notice that D(d, i) is an empirical equivalent of the model given by Equation (2). Let x(d, i) denote the scaled values of D(d, i) according to Equation (4): Note that index i approximately corresponds here to the empirical return period. According to the procedure described in [8,9], the optimum values of climatic parameters of θ and η are determined by minimizing the dispersion ξ of the scaled ordered observations: where x(i) is the mean of the scaled observation values: For the optimal parameters of θ and η, the series x(i) corresponds to a(T) from the model (Equation (2)). By fitting the probability distribution to the series of mean scaled observations x(i) for i = 1, . . . , N, the theoretical quantiles of any order can be determined, and then, by virtue of Equation (2), the theoretical quantile D(d, T) corresponds to any duration d.

Distributions of Maximum Precipitation
For the modeling of annual or seasonal maximum precipitation, probability distributions with no upper bound on the domain and non-negative skewness are usually used; see, e.g., [43]. The set of eight distributions commonly used for the modeling of extreme precipitation are investigated here for the scaled maximum precipitation depth X. They are three-parameter distributions: generalized extreme value (GEV), Pearson 3 (P3), Weibull (We3), log-normal (LN3), and their two-parameter counterparts, i.e., Gumbel (Gu) (a special case of GEV), gamma (Ga), Weibull (We2), and log-normal (LN2). The cumulative distribution functions (CDF) and the formulas for the quantile of order F are presented in Table 2, where: ε ≥ 0 is a location parameter, α > 0 is a scale parameter, k > 0 is a shape parameter, and random variable x ≥ ε ≥ 0. Table 2. Cumulative distribution functions used in the article, along with the quantile formulas (the formulas are from [41], except for GEV [44]).

Distribution
Cumulative Distribution Function (CDF)

Estimation Procedure
The methods for identifying probability distributions of maximum precipitation used in the article require the homogeneity and stationarity of analyzed time series [41]. Here, the Mann-Kendall test [46,47] was used for trend detection in the 68-year series of summer maximum precipitation D(d) for d = 1, 3, and 5 days, respectively. The Mann-Kendall test is a non-parametric (distribution-free) test commonly used in the analysis of environmental, climatic, or hydrological data [48][49][50]. It is based on a rank correlation test for the observed values and their order in time and indicates the significance of a monotonic trend. However, the series of maximum precipitation may also show a more abrupt change at a given time. Thus, the Pettitt test was performed to investigate the homogeneity of the data series and to detect a change point. This non-parametric rank test is widely presented in the literature [51][52][53].
In the next step, for each duration d = 1, 3, and 5 days, the observation series of summer maximum precipitation is ordered in a non-decreasing order as D(d, i) for i = 1, . . . , 68, and then scaled according to Equation (4). Afterwards, the optimum values of parameters θ and η are calculated from Equations (5) and (6).
To determine the quantile of any order for the series of the mean scaled observation x(i) (i = 1, . . . ., 68), the probability distribution to which the series is subject was assumed. Instead of the commonly used approach of adopting one fixed distribution for all considered stations, a new approach was proposed in this step, namely, the selection of the best-fitting distribution from the set of alternative distributions. The eight two-and three-parameter distributions were investigated here ( Table 2). The maximum likelihood method (MLM) [54] was used for the estimation of distribution parameters, and three criteria were used for distribution selection: 1.
The AIC criterion (Akaike information criterion) is based on the likelihood distribution function maximized with respect to the parameter values [55]; 2.
The KS (Kolmogorov-Smirnov) goodness of fit test [56], where the differences between the theoretical and empirical probabilities of an ordered sample are examined. The KS statistics were also used here for the test of goodness of fit of the individual distribution to the series of the mean of scaled observations x(i).

3.
The PCC (Pearson's correlation coefficient) [57] was applied here in the investigation of the correlation between the theoretical and empirical quantiles of the same order.
Regarding the decision rules of the distribution selection criteria 1-3, the probability distribution with the smallest value of the AIC and the KS criteria denotes the distribution (model) best-fitted to the empirical data among the candidate distributions, while in the case of the PCC, the distribution with the highest value of the selection statistic represents the best fit.
For the assumed distribution type of the series x(i), the theoretical quantiles of any order were determined, as described in Section 3.1. The flowchart of the estimation procedure is illustrated in Figure 3.

Stationarity of the Series of Seasonal Maximum Precipitation
As the procedure for constructing DDF curves for extreme precipitation assumes the homogeneity and stationarity of the data series, for each of the 11 meteorological stations, the observation series ( ) were investigated using the Mann-Kendall test for trend detection and the Pettitt test for change point detection considering the duration of 1, 3,

Stationarity of the Series of Seasonal Maximum Precipitation
As the procedure for constructing DDF curves for extreme precipitation assumes the homogeneity and stationarity of the data series, for each of the 11 meteorological stations, the observation series D(d) were investigated using the Mann-Kendall test for trend detection and the Pettitt test for change point detection considering the duration of 1, 3, and 5 days, successively. The p-values of both tests are presented in Table 3. The p-values of the Mann-Kendall test vary from 0.111 to 0.916, and for the Pettitt test from 0.167 to 0.969-these values clearly being higher than the assumed significance level α = 0.05. Therefore, it can be assumed that there is no trend in the studied series of maximum precipitation depths and that the data series are stationary.

DDF Curve Parameter Estimation
In order to determine the depth-duration-frequency relationship for the maximum summer (May-Oct) precipitation for each of the 11 analyzed meteorological stations, the procedure described in Section 3 was performed. The θ and η parameters affect the shape of the DDF curves. Their estimated values are summarized in Table 4, while an example plot of DDF curves for various fixed T values and durations of 1, 3, and 5 days is illustrated in Figure 4 for the Bielsko Biała and Węglówka stations under the assumption that each data series follows the GEV distribution.
Water 2021, 13, x FOR PEER REVIEW 9 of 20 illustrated in Figure 4 for the Bielsko Biała and Węglówka stations under the assumption that each data series follows the GEV distribution.    The values of the θ parameter vary from 0.305 for Katowice to 1.795 for the Węglówka station and indicate a high variability of the point of curvature change. Meanwhile, the values of the η parameter show a much smaller range and vary from 0.704 for the Rycerka Górna station to 1.054 for the Węglówka station. This indicates relatively little variation in the slope of the straight part of the DDF curves. In the case of both parameters, the highest values are for the Węglówka station, which may be due to an orographic effect.

Quantile Estimation
Instead of typically adopting one common distribution for all stations to model the series of the mean of scaled observation x(i), here, eight hypothetical distributions were successively assumed (see Table 2). The MLM method was applied to the estimation of the distribution parameters. For each distribution, the Kolmogorov-Smirnov goodnessof-fit test was performed, which indicated the acceptance of the null hypothesis in all cases. Using the formulas listed in the last column of Table 2, theoretical quantiles of the x(i) series were calculated for each of the distributions, followed by the determination of theoretical quantiles D(d, T) corresponding to any duration d (Equation (2)). Finally, the above theoretical quantiles were compared with the empirical quantiles of the observation series for summer (May-Oct) maximum precipitation of adequate durations and plotted in Figure 5 for all 11 stations.
From the above plots, one can see a regularity related to the location of individual stations, as stations are listed according to their location within the main river and tributary basins. The first two stations, Skoczów and Bielsko Biała, show a fairly good agreement of the empirical and theoretical quantiles for all distributions (except the We2) in their main mass of probability, and a worse match in the range of high probabilities, i.e., in the tail of the distributions. Meanwhile, passing through successive stations, namely, Katowice, Rycerka Górna, Węglówka, the matching in the tail of the distributions improves in all duration periods. For the next stations: Kraków, Kasprowy Wierch, Szaflary, as well as for Harkabuz station located on the lee side of the mountain range, the compatibility of quantiles is high in the entire range of probability variability. A slight deterioration of the tail fit occurs for the last two stations, i.e., Białka Tatrzańska (for two-parameter distributions) and Tarnów (for two-and three-parameter distributions). For all stations and for all duration periods, the worst fit shows the We2 distribution, for which poor agreement is observed in the range of both low and high probabilities.
In visual assessment, the two-parameter distributions in general yield poorer fits of theoretical and empirical quantiles than their three-parameter counterparts. In particular, the Gumbel distribution, traditionally used to derive the DDF (or IDF) relationship, turns out to be inferior for the modeling of summer maximum precipitation in the Upper Vistula Basin relative to the recently proposed GEV distribution.  From the above plots, one can see a regularity related to the location of individual stations, as stations are listed according to their location within the main river and tributary basins. The first two stations, Skoczów and Bielsko Biała, show a fairly good agreement of the empirical and theoretical quantiles for all distributions (except the We2) in their main mass of probability, and a worse match in the range of high probabilities, i.e., in the tail of the distributions. Meanwhile, passing through successive stations, namely, Katowice, Rycerka Górna, Węglówka, the matching in the tail of the distributions improves in all duration periods. For the next stations: Kraków, Kasprowy Wierch, Szaflary, as well as for Harkabuz station located on the lee side of the mountain range, the compatibility of quantiles is high in the entire range of probability variability. A slight deterioration of the tail fit occurs for the last two stations, i.e., Białka Tatrzańska (for two-parameter distributions) and Tarnów (for two-and three-parameter distributions). For all stations and for all duration periods, the worst fit shows the We2 distribution, for which poor agreement is observed in the range of both low and high probabilities.
In visual assessment, the two-parameter distributions in general yield poorer fits of theoretical and empirical quantiles than their three-parameter counterparts. In particular, the Gumbel distribution, traditionally used to derive the DDF (or IDF) relationship, turns

Selection of the Best Fit Distribution
To determine which of the eight alternative distributions is best to determine the depth-duration-frequency relationship for summer (May-Oct) maximum precipitation of selected stations in the Upper Vistula Basin, three criteria of distribution selection were used (Section 3.3) and the results are presented in Table 5. For each criterion of distribution selection, the individual distributions were assigned fitting ranks (1-the best distribution; 8-the worst), and their sum is shown in the fifth column. The best distribution among two-and three-parameter distributions, respectively, was selected according to the sum of ranks. The example quantiles of the order F = 0.99 for duration d = 3 days are given in parentheses for the best-fitted distributions.
The results do not show clear rules related to the location of individual stations within the river basins or the altitude above sea level. The indications for the best distribution vary depending on the selection criterion. The KS and PCC criteria select the three-parameter distributions as the best-fitting for 10 stations. The exception is the Katowice station, where the gamma and Gumbel distributions are superior, according to the KS and PCC criteria, respectively. The PCC criterion most often indicates the GEV distribution (eight times), while the KS indicates various distributions. Meanwhile, the AIC criterion gives different results. The three-parameter distribution is selected only for five stations, and most often it is Weibull3 (three times), while, for six stations, the two-parameter LN2 distribution is the best fit. This is related to the fact that the AIC is defined according to the principle of parsimony and in the case of a similar quality of fitting for all candidate distributions, the AIC favors two-parameter over three-parameter distributions. Summing up the ranks, various three-parameter distributions are selected for eight stations (most often GEV, i.e., three times) and two-parameter distributions for three stations (most often LN2, i.e., two times). In the category of three-parameter distributions, there is no clearly dominant distribution; the most frequently indicated is GEV, but it fits best for only four out of the 11 stations. Meanwhile, in the category of two-parameter distributions, the LN2 distribution clearly dominates, as it is the best for eight out of 11 stations. Table 5. The results of the distribution selection criteria for DDF curves of summer (May-Oct) maximum precipitation for the 11 stations analyzed in the article, along with the quantiles of the order F = 0.99 for duration d = 3 days (given in brackets in (mm)).

Station
Best The values of the quantiles corresponding to the best-fitted distributions according to the various criteria do not show significant variability for most stations, as can also be seen in Figure 5. However, with the exception of two stations, the high-order quantiles (as here the 0.99 quantile) determined under the assumption of the GEV distribution are higher than under the assumption of other considered distributions. The differences in the presented quantiles are due to the worse or better adaptation of the selected distribution types to the observation data. The high values of design quantile estimates may, in fact, be a consequence of a better fit of the GEV distribution to the data series.

Discussion
In this article, the statistical model of a depth-duration-frequency relation for extreme precipitation was proposed. The obtained DDF relation can serve as a factor of flood risk in the Upper Vistula Basin. According to reports, in this area floods occur mainly in the summer season and are caused by intense rainfall, including the millennium flood in 1997, as well as floods in 2001 and 2010 [58][59][60]. The model proposed here is based on the method developed by Javelle [8,9] for peak parts of flood hydrographs; however, it has been modified by introducing the selection of the best-fitting distribution to the observation data series instead of the common practice of adopting one fixed distribution for all stations in the region. Both the first and second approaches may be debatable; the choice depends on the purpose of the research. In practical applications, it is more efficient to use one fixed distribution, due to the simplicity of calculations and for a comparison of the results, as well as the possibility of using the DDF relation for the ungauged (uncontrolled) parts of the region. However, the use of one fixed distribution may lead to an underestimation or overestimation of the maximum precipitation quantiles for some stations. For example, the GEV distribution, which is a heavy-tailed distribution, can provide overestimated values. It has been proven that the GEV distribution fits only partially with the maximum daily rainfall in the Upper Vistula region, where the LN distribution is a significant complement [61], and the GEV poorly matches the Polish hydrological data [19,62]. The DDF relationships of extreme precipitation intensities are a set of the most commonly used tools in water resources engineering for the planning, design, and operation of water resources projects. In particular, the estimates of high quantiles of maximum precipitation and maximum flows (most often the quantile of 0.99) serve as design values for the construction of hydrotechnical structures, which affects investment costs. Overestimated values of design quantiles increase investment costs, while underestimated values of quantiles may result in incomplete protection against the effects of extreme hydro-meteorological phenomena. The approach of one fixed distribution chosen among a set of candidates was successively developed and used in design practice in Poland [63,64]. The statistical model covers all Polish territory except the areas of the Sudetes and the Carpathian Mountains due to their complex rainfall field with an insufficient number of gauging stations (mainly located in the valleys) to perform the proper statistical inference. Therefore, the approach proposed in this article of selecting the best distribution among alternative distributions is advantageous, as it provides a more accurate estimate of the DDF (or IDF) relationship for maximum precipitation in the studied area.
Testing the homogeneity and stationarity of a maximum precipitation series is often overlooked in studies on the design of DDF or IDF curves, and it is assumed that this condition is met [5,6,42,65]. The results of the homogeneity of the series of seasonal maximum precipitation of Polish stations investigated here are consistent with the literature, since, unlike in other European countries, no statistically significant trends in both seasonal and annual precipitation totals were found in Poland and Slovenia [66][67][68]. Only Młyński et al. [27] found a growing significant trend in the annual maximum daily precipitation for the period 1971-2014 (i.e., a period shorter than the one treated of here) at the level of 0.05 for 12% of 51 meteorological stations surveyed in the Upper Vistula region. Referring to the maximum flows, the trends in daily flow maxima on the Vistula and some of its tributaries in the seasonal timeframe were studied in [68]. The research did not reveal any significant changes in the high waters regime. The Mann test (not presented here) confirmed also that there are no statistically significant trends in summer peaks flow in the period 1951-2018 at the hydrological stations on the Upper Vistula reach. Therefore, based on the hydrometeorological data available at present, it can be concluded that flood risk does not change with time, and the impact of human pressures cannot be assessed (yet) as significant.
Comparing the empirical quantiles from the observations with the theoretical quantiles for the proposed eight distributions shows consistency in the main mass of probabilities for all stations and all durations ( Figure 5). However, for some stations, a poor match in the high probability range is observed, i.e., for the stations Skoczów, Bielsko Biała, Katowice, and Tarnów. This is related to the type of distribution itself, as well as to the estimation method used, since the MLM concentrates on the main mass of probability (i.e., far from the upper quantiles). A possible future approach would be to use the two-component extreme value (TCEV) probability distribution, which is based on the hypothesis that the extreme values do not come from the same population as the ordinary ones [69][70][71], or to use a probability mixture model, which is based on the aggregation of the probabilities of non-exceedance of a constant rainfall value from the candidate distributions [72].
The results presented in the Q-Q plots in Figure 5 show that, in general, the threeparameter distributions seem to be adequate; however, there is no universal distribution for all stations. Both the GEV and LN3 distributions yield high quantiles agreement, although P3 and We3 are only slightly worse. Among the two-parameter distributions, LN2 stands out as the best, followed by Gu and Ga, while We2 is definitely the worst. Thus, the obtained results do not coincide with the commonly used approach of adopting one fixed probability distribution for all meteorological stations in a given region-typically this is the Gumbel distribution [4,5,73,74]. The assumption of GEV distribution found in some studies [6,7,42,75] seems to be more reasonable.
The choice of the best-fitting distribution for individual stations depends on the selection criterion. The results in Table 5 indicate that, regardless of the adopted criterion, the regional distribution of DDF curves by the type of best-fitted distribution is difficult to determine. The only criterion which shows the dominance of one distribution (GEV) is the PCC; however, here, also, no regional consistency is observed. For example, Szaflary and Białka Tatrzańska stations are located closest to each other among all considered stations, but individual criteria (including PCC) showed a different probability distribution for each station. The situation is similar with the closely located stations of Skoczów and Bielsko Biała, where the AIC and KS criteria indicated different best-fitted distributions. Another group includes the Katowice, Kraków, and Tarnów stations located in the northern part of the studied area and for them the KS and PCC criteria yield different distributions. Katowice is the only station where, according to each of the criteria, two-parameter distributions are best suited, among them the Ga and Gu distributions, not selected at any other station. The specificity of Katowice is related to the location of the meteorological station at the airport, where strong winds may affect the accuracy of precipitation measurements. Moreover, the air quality in the city of Katowice is poor (high pollution). In order to draw conclusions about the regional distribution of the best-fitted probability distributions, it would be necessary to increase the number of analyzed stations. Even then, however, it may not be possible due to the high diversity of the terrain.

Conclusions
Studies on the determination of the depth-duration-frequency (DDF or equivalently IDF) relationship are particularly important in the region of the Upper Vistula, as this area is exposed to floods of rainfall origin in the summer season. The analysis of maximum annual or seasonal precipitation as a determining factor of floods is usually carried out by regional methods. This approach is very useful in hydrological practice and research. However, despite the advantages of the method, its application in an orographically diverse terrain, as characterises a large part of the analysed area, the regional approach is a big challenge, especially when we have at our disposal only a limited number of stations with observations taken over sufficiently long periods. As precipitation in mountain regions is strongly influenced by the local conditions, general findings on the suitability of regionalization methods for such an area can only be obtained by analysing the data from a large number of stations at the same time. In addition, the stations should represent all climatic floors and expositions of the slopes. That is why the results presented here cannot be extended to greater areas, especially in the mountainous part of the Basin, and can serve as preliminary information about distributions and quantiles.
The research reported here on the statistical model of the depth-duration-frequency relation for summer (May-Oct) extreme precipitation measurements from the 11 stations in the Upper Vistula Basin provides the following conclusions:

1.
For the 11 weather stations analyzed in the article, the series of the maximum sums of 1-, 3-, and 5-day precipitation are stationary. This allows for the conclusion that there is no basis for stating that the risk of flooding from rainfall in the summer season (May-Oct) increases. However, it does not exclude the occurrence of significant floods in the coming years.

2.
The statistical DDF model proposed in the article provides a more accurate estimation of the DDF curves for individual stations than the commonly used approach of one distribution for all surveyed stations. However, the choice of approach depends on the purpose of the research. When the precision of the DDF relationship is important, the approach proposed in the article is recommended. 3.
The Gumbel distribution that is traditionally used for the construction of DDF curves turned out to be inferior in relation to the distributions proposed in the article; in particular, it turned out to be worse than the two-parameter log-normal distribution. 4.
The three-parameter distributions show a better fit to the seasonal (May-Oct) maximum precipitation in the Upper Vistula Basin than their two-parameter counterparts.
The GEV distribution turned out to be the best, but its advantage over other distributions is not significant.

5.
The estimates of the upper quantiles of maximum precipitation totals in 1, 3, and 5 days (i.e., quantiles) are smaller in the lowlands and greater in the highlands. 6.
Due to the insufficient number of stations and the highly diverse terrain, the regional distribution of DDF curves cannot be determined.
Further research will concern the extension of the observation data set with data from additional meteorological stations (though the period of observations will be shorter) and the preparation of maps showing altitude, slopes exposure, land cover, and land use of the study area, as well as data on circulation types. We will then attempt to determine the regions and the regional DDF curves with the option of selecting the best-fitting distribution. We are also planning research based on Monte Carlo simulations. This would provide a comprehensive theoretical and practical tool for determining the DDF (IDF) relationship for maximum precipitations in the Upper Vistula Basin region.