Regional Flood Frequency Analysis of the Sava River in South-Eastern Europe

: Regional ﬂood frequency analysis (RFFA) is a powerful method for interrogating hydrological series since it combines observational time series from several sites within a region to estimate risk-relevant statistical parameters with higher accuracy than from single-site series. Since RFFA extreme value estimates depend on the shape of the selected distribution of the data-generating stochastic process, there is need for a suitable goodness-of-distributional-ﬁt measure in order to optimally utilize given data. Here we present a novel, least-squares-based measure to select the optimal ﬁt from a set of ﬁve distributions, namely Generalized Extreme Value (GEV), Generalized Logistic, Gumbel, Log-Normal Type III and Log-Pearson Type III. The ﬁt metric is applied to annual maximum discharge series from six hydrological stations along the Sava River in South-eastern Europe, spanning the years 1961 to 2020. Results reveal that (1) the Sava River basin can be assessed as hydrologically homogeneous and (2) the GEV distribution provides typically the best ﬁt. We offer hydrological-meteorological insights into the differences among the six stations. For the period studied, almost all stations exhibit statistically insigniﬁcant trends, which renders the conclusions about ﬂood risk as relevant for hydrological sciences and the design of regional ﬂood protection infrastructure.


Introduction
Floods are among the most impactful natural disasters and can cause significant human casualties and economic damage [1,2]. In Europe, numerous large rivers have been affected by extreme flood events since the 1990s, causing billions of Euros in damages [3][4][5]. Such impacts are expected to increase in the future partly due to climate change, and partly due to exposure of larger populations and economic growth [6,7].
The State of the Global Climate 2020 report [21] (WMO, 2021) marked that year as exceptional for worldwide flooding and other hydro-meteorological extremes. Such major events are expected to become more common in the future due to the increase of water vapour in the atmosphere (i.e., the Clausius-Clapeyron equation) leading to more the question (3) after the study's novelty and limitations and the practical consequences that can be drawn for decision makers.

Research Area
The Sava River basin, with a total area of around 97,000 km 2 , is the largest river catchment in South-eastern Europe and the second-largest sub-basin of the Danube River basin (Figure 1). It is shared by Slovenia, Croatia, Bosnia and Herzegovina, and Serbia; and a small part of 0.2% of the basin lies within Albania [50]. The Sava River originates in Slovenia, where two mountainous streams, Sava Dolinka and Sava Bohinjka, merge together near the town of Radovljica. The main tributaries of the Sava River are in downstream order the rivers Kolpa, Una, Vrbas, Bosna, Drina and Kolubara. Between the Una and the Drina, the majority of water comes from the south. Una, Vrbas, Bosna and Drina jointly contribute a discharge of 1149 m 3 /s to the Sava River, which is as much as 68% of the flow of the Sava at its confluence with the Danube [51]. Drina is the largest tributary with 371 m 3 /s or 22% of the flow of the Sava at its confluence with the Danube [51]. The length of the Sava River is 945 km to its confluence with the Danube River in Belgrade. Around 8 million people live in the Sava River basin [36]. The headwaters of the Sava River have a mountainous climate that transitions downstream to a temperate continental climate [52]. The temperate continental climate characterizes the right tributaries originating in the Dinaric mountains, whereas a moderate continental climate features on the left side of the basin, including the The headwaters of the Sava River have a mountainous climate that transitions downstream to a temperate continental climate [52]. The temperate continental climate characterizes the right tributaries originating in the Dinaric mountains, whereas a moderate continental climate features on the left side of the basin, including the Pannonian plain [50]. The average annual precipitation in the Slovenian part of the Sava River ranges from 2210 mm representing the basin area of the Radovljica gauging station [32] to 1570 mm (Čatež gauging station) [53]. Over the Croatian part of the basin, the average annual precipitation is 855 mm [51]. The highest monthly amount of rainfall is usually recorded during June (91 mm), and the lowest typically in February (42 mm). For the Serbian part of the basin (S. Mitrovica station), the average is 577 mm, with the highest mean in June (84 mm) and the lowest in February (38 mm) [36].

Data
We utilized a 60-year (1961 to 2020) data set of daily discharges for six stations located on the Sava River ( Figure 1). This data set was collected from three different national water authorities within the investigated river basin (namely the Slovenian Environment Agency; Meteorological and Hydrological Service, Croatia; and Republic Hydrometeorological Service of Serbia) ( Table 1). The assumption of independence of extreme flows [56,57] needed for statistical analysis was met using the annual maximum series (AMS) of river discharges. AMS is a data type that is most frequently used for FFA as it is the most intuitive example of the block maxima method of the extreme value theory [31,58]. In FFA, a block usually represents a year, and so the highest discharge value within each year is selected [39]. One of the advantages of the AMS method is that it is straightforward to derive the distribution of the AMS values. Descriptive summary statistics of the AMS data are shown in Table 2.

Methods
We applied the RFFA method combined with L-moments estimation [43] to the AMS data from six gauging stations along the Sava River. The advantage of RFFA is that it can be applied for flood risk estimation in situations where for a specific site the recording period is too short for a reliable estimation (or even where data are not available). The advantage of the L-moments technique (as compared with maximum likelihood estimation) is that it yields robust parameter estimates for various candidate distributions. In a statistical sense, the robustness of a method means that it yields reliable results even if assumptions underlying an analysis (such as the distributional shape) are violated [59]. The homogeneity of a region can then be assessed on the basis of a comparison test of the L-moments parameter estimates for the different gauging stations, supported by a goodness-of-fit test and L-moments diagrams [60,61].
The principle assumptions employed for the RFFA are that the observed data have been generated by a stationary random process [62,63]. One of the advantages of the AMS data type is that it allows statistical science to derive the distributional shape of the AMS values. This distribution is called Generalized Extreme Value (GEV) distribution, which emerges as the limit form of the maximum of many independent and identically distributed random elements in a time block (such as one year), and this statistical work has been performed in the late 1920s [64,65]. Furthermore, the AMS method provides flood event data that can be considered as independent and the frequency distributions generally conform to theoretical distributions [66]. For recent methodological developments in RFFA and its applications, see, for example, contributions about Canada [67] or Australia [68].
There are several guidelines issued by national hydrological services related to the choice of the distribution used in FFA. In the United States, the Log-Pearson Type III (LPIII) distribution should be used for FFA according to the United States Water Resources Council [69]; in the United Kingdom, usage of the GEV or the Generalized Logistic (GLO) distributions is advised according to the Flood Estimation Handbook [70]; in Poland, usage of the Pearson Type III (PIII) distribution is recommended [71]; and in Australia, it is recommended to compare the LPIII, GEV and Generalized Pareto estimation results before selecting a distribution [72].
For the rivers within the Pannonian basin, previous studies indicate that the best-fitting regional distribution shape could be lognormal [31]. However, further research conducted at individual stations showed that other distributions should be considered. For example, for the Bogojevo station on the Danube River, the PIII distribution was presented as the best fit [42]. A similar study conducted for the basin of the Tisza River (a tributary of the Danube River) also found that the best-fitting distribution is PIII [41]. Regarding the Sava River, Morlot et al. [42] demonstrated that PIII distribution is the best-fitting distribution for the Orljava-Pleternica gauging station (Croatia) as well as for the S. Mitrovica gauging station (Serbia). These results for the Sava River were further confirmed [36]. So far, in the Sava River basin and in the countries that share this basin, there are no official guidelines for FFA. Here, we applied the five most commonly used distributions in this region, namely: GEV, GLO, Gumbel, Log-Normal Type III (LNIII) and LPIII. The Appendix A provides the cumulative distribution functions (CDFs) for these shapes.

L-Moments Estimation for RFFA
The L-moments method was used to estimate distribution parameters. L-moments are a linear combination of the ranked observations [73]. What distinguishes L-moments from "conventional" statistical moments is their ability to describe a wider range of distributions and that they are more robust to the presence of outliers and, hence, less susceptible to estimation bias [43]. The most commonly employed method to calculate the L-moments is via probability weighted moments (PWM) [35,63]: where b r represents the rth-order PWM, while F X (x) is the CDF of X. Further, b i for i = 0, 1, 2, 3, represents the unbiased sample estimators of the first four PWM as described by Hosking and Wallis [43]: In these equations, X (j) indicates the ranked AMS values, that is, X (n) is the largest and X (1) is the smallest value. Furthermore, the first four L-moments [43] are given as follows: That means, the first four L-moments follow straightforwardly from the PWM sample estimates. Another advantage of the usage of the L-moments approach is that these moments provide directly interpretable values for the mean or location (λ 1 ), scale (λ 2 ), shape (λ 3 ) and kurtosis (λ 4 ) of the distribution functions, a fact that is particularly useful in RFFA. Lastly, the L-moments ratios are calculated as follows [43]: For RFFA, the L-moments approach should be applied in four steps [43], comprising: (1) screening of data, (2) identifying homogeneous regions, (3) choosing a regional frequency distribution and (4) estimating the parameters of the frequency distribution.
When determining the regional distribution and homogeneity of a region, the L-moments ratio diagram [66] should be constructed. The close association between LCs and LCk has been shown for the most frequently used probability distributions [43]. If a region is homogeneous, then the locations tend to group together in the L-moments ratio diagram. Closeness of the sample versions of LCs and LCk to a theoretical line of distributions (in case of threeparameter distributions), or of the sample mean to the theoretical points (in case of twoparameter distributions), determines the best-fitting parameter distribution describing the data from multiple sites [35]. If in the L-moments ratio diagram a point is located near the curve that represents a candidate distribution, that distribution can be considered as an acceptable selection for the regional distribution. The distance between a point and curve can be considered as measure of the goodness-of-fit. It is important to highlight that this method is dependent on the homogeneity of the regional data [74].
For further confirmation and for determination of the regional distribution, a statistical test of spatial homogeneity of the selected region should be conducted [43]. This test is based on a comparison of (1) the between-site variation in sample-LCv and (2) the expected variation for a homogeneous region. This test is calculated as follows: In these equations, N is the number of sites; n i is the record length at site i; τ (i) , τ 3 (i) , τ 4 (i) are the sample values of LCv, LCs and LCk at site i, respectively; and τ R , τ R 3 , τ R 4 are the regional sample averages of LCv, LCs and LCk, respectively. For a region to be considered as homogeneous, the test statistic values, V i , for i = 1, 2, 3, should all be less than unity. If 1 ≤ V i < 2, then a region is possibly homogeneous and if V i ≥ 2, then the region can be assumed as heterogeneous [31,73].
In order to test whether a particular distributional shape provides a good fit to data, and furthermore in order to compare the fit quality of various candidate distributions, we employed two test measures.
The first test measure is the conventionally employed Z-statistic [43]. It is based on the comparison of sample-and population-LCk values for the different distributions: In this equation, the acronym DIST refers to a particular distribution, τ DIST 4 is the population-LCk value for a particular candidate distribution, τ R 4 is the regional sample average of the LCk values and σ 4 is the regional sample standard deviation of the LCk values [31]. The test is performed by means of a comparison with a quantile of the standard normal distribution. We adopted a significance level of 5%, this means, a candidate distribution is considered suitable if |Z| < 1.96.
A second test measure is introduced here since this paper aims not only to find the optimal distributional shape for the data from the Sava River basin, but also to contribute data-analytical tools for hydrology. As with the Kolmogorov-Smirnov test [75,76], our metric employs the empirical CDF of the AMS values. More specifically, let n be the sample size; X (j) , j = 1, . . . , n, be the size-sorted AMS values (as in Section 4.1); and the empirical CDF, F n (x), be defined as the number of X (j) ≤ x divided by n; then F n (x) is a step-like curve in dependence on X (j) , which goes from 0 to 1 and makes steps of size 1/n at each X (j) value. The smaller the distance is between a certain theoretical candidate distribution, F theo (x), the better the fit. The distance can be measured by means of the differences, While the Kolmogorov-Smirnov test utilizes max{w(j)}, we employed a least-squares criterion, Our metric allows all data points to contribute to the distance measurement since we require here that a suitable theoretical distribution should supply a good fit over the full range of data values. The factor (n − p) −1 is used to correct for the number of employed parameters, p, of a theoretical distribution. (That means, it is "easier" to fit a distribution with many parameters to data than with few, and Ockham's Razor [56] indicates to prefer descriptions with fewer parameters.)

Trend Estimation
To test the stationarity assumption for the RFFA analysis, we employed linear ordinary least-squares regression [56,63] to determine trends in the AMS data. Let {T(k), X(k)} n k=1 denote the time series, where T(k) is time and X(k) is AMS. The linear model is given by where a is the intercept parameter, b the slope parameter and X noise (k) the noise component. The minimization of the least-squares sum, can be straightforwardly achieved by setting the first derivatives of SSQ (with respect to a and b) equal to zero, which yields two linear equations (called normal equations) that can be easily solved for the parameter estimates of a and b [63]. To assess the significance of a trend, an uncertainty measure (standard error) for the slope parameter, b, is indispensable. This study employed moving block bootstrap resampling [63] to obtain slope standard errors that are reliable. Moreover, the tests on X noise (k) can indicate the presence of (1) non-normal distributions or (2) autocorrelation. The validity of the bootstrap for obtaining reliable uncertainty measures in regression analysis has been shown by means of statistical simulation experiments with artificial time series [63].

Regional Flood Frequency Analysis
The principle assumptions employed for the RFFA and the fitting of candidate probability distributions are that the observed data have been generated by a stationary random process [62,63]. Therefore, the first step was the calculation of the basic statistics of the AMS data (mean, median, maximum, standard deviation, skewness and kurtosis) for the six gauging stations ( Table 2). These basic statistics were then used to define the selected frequency distributions for deriving the probability of exceeding a certain discharge value [77].
Hosking and Wallis [43] highlighted that there are little gains in accuracy of quantile or return level estimates when using regions that contain a large number of stations (above, say, 20) and that it is preferable to have a region with fewer stations but longer at-site records. In this study we used the six stations on the Sava River, which mainly have 60 years of data ( Table 2). The positive skewness values (Table 2) indicate that the AMS data for the stations show deviations from a normal shape (in the form of right-skewness), as is to be expected for a block maxima series [63]. This result also justifies usage of the bootstrap procedure for trend estimation (Section 4.2).

Trend Analysis
Linear trend analysis of the AMS at the six gauging stations along the Sava River showed mixed results (Figure 2, Table 3). AMS at the Radovljica andČatež gauging stations in Slovenia demonstrate increasing but statistically insignificant trends (i.e., the absolute value of the slope estimate is less than its standard error). A statistically insignificant increase of AMS was also observed at the Zagreb gauging station during the study period. Also shown are linear trends (red lines); slope estimates are given in Table 3.
In general, negative AMS trends in the eastern part of the basin may be related to decreased snowmelt contributions from the Dinaric mountains, which represent the main water input for the right-hand tributaries. This influences the discharge regime of the Sava River downstream of Jasenovac [55]. Negative trends observed in AMS have also been attributed to (1) regional warming and (2) ocean-atmosphere modes of climate variability [2], the combination of which may further have caused a decline in precipitation and an increase in evaporation in the basin [22]. Comparable results of insignificantly decreasing river discharge trends were previously reported for the Danube and Sava Rivers [36,78].  Table 3. In general, negative AMS trends in the eastern part of the basin may be related to decreased snowmelt contributions from the Dinaric mountains, which represent the main water input for the right-hand tributaries. This influences the discharge regime of the Sava River downstream of Jasenovac [55]. Negative trends observed in AMS have also been attributed to (1) regional warming and (2) ocean-atmosphere modes of climate variability [2], the combination of which may further have caused a decline in precipitation and an increase in evaporation in the basin [22]. Comparable results of insignificantly decreasing river discharge trends were previously reported for the Danube and Sava Rivers [36,78].
At the two most downstream gauging stations (Županja and S. Mitrovica), maximum annual discharges are decreasing, likely as a consequence of a considerably larger catchment area for these stations and because large tributaries join the Sava River downstream. Consequently, the river discharge regime is more stable and maximum discharges are less influenced by severe and spatially limited precipitation. Presented results are consistent with other studies [22,79,80], which project decreasing discharge for major rivers in central and southern Europe. Furthermore, river modelling indicates a decreasing discharge of up to 40% for the Sava River at its lower course (S. Mitrovica station) [81].

Homogeneity Assessment and Regional Flood Frequency Analysis
The numerical results of the three homogeneity measures V 1 , V 2 and V 3 ( Table 4) indicate that after Hosking and Wallis [43] we can consider the Sava River basin as acceptably homogeneous. Since all measures are less than unity, no further inspection of the data or of the cross-correlation between sites was necessary. That means, that the AMS data from the Sava River basin can be meaningfully applied for RFFA and also for a subsequent quantile estimation for gauged and ungauged rivers within the basin. Table 4. Homogeneity test of Sava River basin.
0.14 0.13 0.14 0.04 0.07 0.077 Since the basin homogeneity could be confirmed, the goodness-of-fit test and the L-moments diagram were used for the determination of the best-fitting distribution for RFFA ( Figure 3).
Owing to the robustness of the L-moments, station sample points can be expected to be randomly distributed above and below the theoretical line of an appropriate distribution. The results (Figure 3) show that there is no clear variation of LCk and LCs within the Sava River basin. Most data points display a tendency to group around the GEV and GLO distribution functions. In such cases, the mean value can be used as a reliable indicator of the best-fitting distribution [82]. Our results appear to indicate that the GEV distribution is the best-fitting distributional shape in the case of the Sava River basin.
For the purpose of identifying a common regional distribution that can be applied for ungauged rivers within the basin, the Z-statistic was calculated. The results reveal that all five distributions satisfy the test criterion |Z| < 1.96, namely GEV (Z = −0.06), GLO (Z = −0.20), Gumbel (Z = 0.28), LNIII (Z = 0.38) and LPIII (Z = 0.27). Since the absolute value of Z for the GEV distribution is the lowest, this is another indication that this distribution is valid for the Sava River basin as a whole. These findings are consistent with an FFA for Slovenia, where the GEV was identified as the best-fitting regional distribution [83]. Furthermore, results show that data from the gauging stations in the Sava River basin can be used for RFFA, and this application increases the reliability of the quantile estimates for both ungauged and gauged tributaries.  Owing to the robustness of the L-moments, station sample points can be expected to be randomly distributed above and below the theoretical line of an appropriate distribution. The results (Figure 3) show that there is no clear variation of LCk and LCs within the Sava River basin. Most data points display a tendency to group around the GEV and GLO distribution functions. In such cases, the mean value can be used as a reliable indicator of the best-fitting distribution [82]. Our results appear to indicate that the GEV distribution is the best-fitting distributional shape in the case of the Sava River basin.
For the purpose of identifying a common regional distribution that can be applied for ungauged rivers within the basin, the Z-statistic was calculated. The results reveal that all five distributions satisfy the test criterion | | < 1.96, namely GEV (Z = −0.06), GLO (Z = −0.20), Gumbel (Z = 0.28), LNIII (Z = 0.38) and LPIII (Z = 0.27). Since the absolute value of Z for the GEV distribution is the lowest, this is another indication that this distribution is valid for the Sava River basin as a whole. These findings are consistent with an FFA for Slovenia, where the GEV was identified as the best-fitting regional distribution [83]. Furthermore, results show that data from the gauging stations in the In order to examine the goodness-of-fit of the selected distribution to the empirical data, the CDFs were calculated for the analysed hydrological stations (Figure 4).
From Figure 4 it is notable that at the highest values (i.e., floods), several candidate distributions are in good agreement with the empirical data. However, the GEV distribution is in good agreement for the entire range of empirical data and therefore delivers the best goodness-of-fit measures (Table 5).
Owing to the superiority of the distributional fit (Figure 4), the GEV is used to estimate the risk curves [63], that is, the return level in dependence on the return period ( Figure 5). It is important to point out the inherent uncertainties (expressed as 95% confidence bounds) in the estimated return levels, which are due to the sampling variability, that is, a finite data set of noise-affected observations [31]. As expected, the range of uncertainty increases for longer return periods [29,31,56,63]. Despite the associated uncertainties of the estimated risk curves, the GEV distributions with their 95% confidence bounds can provide useful quantitative information for the construction of water infrastructure. For example, a precautionary guideline could be to apply the upper 95% confidence bound for a return period of 100 years. Mudelsee [63] recommends that the return period should be of the same order as the length of a series because the uncertainties increase with the return period ( Figure 5). For the Sava River with record lengths of 60 years (Table 1), an upper bound estimate for a return period of 100 years may be suitable, whereas a bound of 500 years is certainly not. Sava River basin can be used for RFFA, and this application increases the reliability of the quantile estimates for both ungauged and gauged tributaries. In order to examine the goodness-of-fit of the selected distribution to the empirical data, the CDFs were calculated for the analysed hydrological stations (Figure 4). From Figure 4 it is notable that at the highest values (i.e., floods), several candidate distributions are in good agreement with the empirical data. However, the GEV distribution is in good agreement for the entire range of empirical data and therefore delivers the best goodness-of-fit measures (Table 5).   One of the main purposes of FFA is to determine the quantiles (return levels) in the extreme upper end (i.e., flood levels) of the best-fitting distribution. This was achieved by estimating the parameters of the overall best-fitting (i.e., the GEV) distribution for the selected stations. The regional shape parameter is the average of the at-site shape One of the main purposes of FFA is to determine the quantiles (return levels) in the extreme upper end (i.e., flood levels) of the best-fitting distribution. This was achieved by estimating the parameters of the overall best-fitting (i.e., the GEV) distribution for the selected stations. The regional shape parameter is the average of the at-site shape parameters [43]. The return discharge levels for the various stations follow then from the regional shape parameter plus the other at-site GEV parameter estimates ( Figure 6). parameters [43]. The return discharge levels for the various stations follow then from the regional shape parameter plus the other at-site GEV parameter estimates ( Figure 6). One of the advantages for using RFFA is that the estimates at a single site can be increased by applying data from other sites within homogeneous region that are confirmed to have similar frequency distribution [43,84]. Application of the L-moments technique is appropriate for increasing the sample size for the RFFA. The L-moments approach simultaneously uses data from several homogeneous basins in a hydrological analysis [85]. Our study shows that this technique is an effective approach for discharge estimation of flood peaks related to the selected return periods in basins with missing data or relatively short time span of data. One of the advantages for using RFFA is that the estimates at a single site can be increased by applying data from other sites within homogeneous region that are confirmed to have similar frequency distribution [43,84]. Application of the L-moments technique is appropriate for increasing the sample size for the RFFA. The L-moments approach simultaneously uses data from several homogeneous basins in a hydrological analysis [85]. Our study shows that this technique is an effective approach for discharge estimation of flood peaks related to the selected return periods in basins with missing data or relatively short time span of data.

Conclusions
We employed the FFA and RFFA methods using various theoretical distributions fitted to AMS discharge series from six hydrological stations in the Sava River basin. Our first aim was to determine the best-fitting distributional shape for the various datasets.
Using long-term series (1961-2020) with few missing data enabled a robust statistical inference. The homogeneity test confirmed that the investigated area can be considered as homogeneous, so no further subdivision into different basin parts was necessary. Our second aim was to compute regional frequencies of the AMS (i.e., flood risk parameters). For this purpose, quantiles and return levels were estimated at a confidence level of 1 − α = 95%. The results show that the considered distributions (GEV, GLO, Gumbel, LNIII and LPIII) are for most of the stations in good agreement with the test requirement (|Z|<1.96). Overall, the GEV distribution delivered the smallest absolute Z-values, almost uniformly the smallest u-values, and also the regional mean confirms that the GEV is the best-fitting option. Hence, we conclude that the GEV distribution provides the best fit for the Sava River basin.
A trend analysis showed that although the majority of stations (four out of six) exhibit increasing trends, few of the trend estimates are statistically significant. The observed variations between positive and negative trends are thought to be mainly a consequence of modified conditions in the upstream river sections. In such cases, the size of the basin and the degree of human interventions make it difficult to quantify the relationships of climate elements (i.e., the temperature rise) and discharge properties [86]. Model projections suggest a decreasing trend in thickness and duration of the snow cover in the Alpine and Dinaric area in the twenty-first century [55]. These changes must also be taken into consideration for adaptation planning [87]. Particularly important is the observed change at the relatively low elevations (600 to 1300 m) [55].
The results of this study provide helpful information for designers of the various hydraulic structures (levees, bridges, spillways, etc.) and also for policymakers in need of a decision-support tool for flood estimation along the Sava River. The results also enable more robust estimates of design discharges corresponding to different return periods for the design of the hydraulic structures, which could lead to the reduction of the risk of the failure of these structures and to reduce the negative impacts and environmental damages from future flooding.
We provide a short outlook from a practical flood risk protection viewpoint. On the one hand, the underlying datasets are (1) limited in length, (2) limited in time resolution and (3) plagued by measurement errors (e.g., from the water stage-runoff calibrations). This situation leads inevitably to statistical uncertainties of made inferences [63], such as the return levels ( Figures 5 and 6). On the other hand, however, the data constitute the best what is, and what will be, available and we believe that the trend analyses and the RFFA methodology (Section 4) are quite capable of taking the data quality into account and provide reliable uncertainty measures. An interesting option for future short-term flood risk protection is usage of modern instrumental, wide-range measurement data, which have the potential to contribute flood information at high spatial resolution. This has been attempted recently for Bangladesh with the help of satellite data [88] and Iran with the help of remote-sensing data [89]. Finally, in a rational world, where science informs politicians-what decisions should follow from the risk curves such as for the return period ( Figures 5 and 6)? Answer: it depends, namely on the financial resources and the political instruments at the disposal of the decision makers. If money is available and also the will to protect from damages, these curves may help to find out at which geographical place a good return (in terms of prevented losses) can be expected from a made investment. Scientists should be ready to respond if asked for advice by decision makers.

Acknowledgments:
The discharge data used in the present study were obtained from the Slovenian Environment Agency (ARSO), Croatian Meteorological and Hydrological Service (DHMZ) and Republic Hydrometeorological Service of Serbia (RHMZ). We thank Professor Rob Wilby (Loughborough University, United Kingdom) for a critical read of the manuscript. We are grateful for all three detailed and constructive reviews received on the initial submission.

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

Appendix A. Mathematical Formulas for the Distribution Functions
We briefly write down the CDFs for the five theoretical distributions fitted to the data. If there exist several variants of a CDF, then we only consider those with nonzero shape parameter and right-skewness. Note that there exist various ways of expressing the formulas and various parameter notations. More details, such as general descriptions and numerical approximations, can be found in standard textbooks [43,63,[90][91][92].
The GLO distribution has the following CDF: where the parameters satisfy the following conditions: -∞ < µ < ∞ and σ > 0. The Gumbel distribution has the following CDF: where the parameters satisfy the following conditions: −∞ < µ < ∞ and σ > 0.
The LNIII distribution has no explicit analytical form for the CDF. However, the probability density function, f X (x), where is given by where ln(·) is the natural logarithm of an argument and the parameters satisfy following conditions: −∞ < µ < ∞, x > µ and σ > 0.