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

Abstract

Regional flood 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-fit measure in order to optimally utilize given data. Here we present a novel, least-squares-based measure to select the optimal fit from a set of five distributions, namely Generalized Extreme Value (GEV), Generalized Logistic, Gumbel, Log-Normal Type III and Log-Pearson Type III. The fit 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 fit. We offer hydrological-meteorological insights into the differences among the six stations. For the period studied, almost all stations exhibit statistically insignificant trends, which renders the conclusions about flood risk as relevant for hydrological sciences and the design of regional flood protection infrastructure.

1. 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 following list of works on observed river floods during the past decades and centuries may help the reader to put the current analysis into a wider geographical context: [4,8,9,10,11,12,13,14,15,16,17,18,19,20]. These papers include analyses of meteorological and demographic causes of extreme river flooding effects and also provide additional costs estimates.

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 intense precipitation events and, as a consequence, to more floods [22,23,24,25]. This has been confirmed by high-resolution regional climate modelling for Europe [26,27]. Moreover, expansion of cities, construction of roads and bridges across rivers, can all affect floodplain hydraulics and thereby change flood risks for exposed populations and assets [4,28]. Thus, for flood risk management purposes it is important that periodic assessments of rivers are conducted, especially regarding long-term discharge patterns. This can be achieved via flood frequency analysis (FFA), which can offer important information regarding river hydrological behaviour, also in changing environment [29,30,31,32].

For the determination of flood magnitude with a given return period, the most commonly applied statistical data-based approach is FFA. The standard method of FFA fits mathematical functions to available data and extrapolates the tails of the distribution to estimate the probability and magnitude of extreme events [33]. When sufficient data are available at a single gauging station, an extreme value distribution is typically fitted to the observed data, indicating an at-site FFA or local model [34]. However, for stations with incomplete or short data series, as well as for ungauged sites, data from comparably situated or nearby gauging stations can be used. This approach is known as regional flood frequency analysis (RFFA) [34]. Reliable estimates of magnitude and frequency of floods are essential for the design of different hydraulic structures such as spillways, levees, bridges, dams, irrigation ditches, bridges and urban drainage systems [35].

In recent years, many at-site FFA of different hydrological variables have been conducted for Europe and other parts of the world [22,29,35,36,37].

However, there is growing interest in the use of RFFA across Europe as evidenced by research in Germany [38], Norway [39], Poland [40], the Pannonian basin [31], the Tisza River basin [41] and the Danube River basin [42]. This is because accuracy can be improved by drawing on the information content of multiple sites [43]. Moreover, RFFA techniques offer opportunities for further development and application in decision-making [35,44].

Floods in Serbia are mostly caused by either (1) flow of humid air masses originating from the Atlantic Ocean, or (2) rapid snow melt from the Dinaric mountains associated with föhn winds in early spring. In addition, we note that the formation of floods in further downstream parts of the study area is influenced by erosion processes and the presence of a dense river network (note that the average river length per area in Serbia is 747 m/km²), see for example [28]. In the spring of 2014, Croatia, Bosnia and Herzegovina, and Serbia, were impacted by a flood event unparalleled in the hydrological records of the region [45,46]. The major cause of this flood was a Vb-type cyclonic weather system [5] that remained stationary over South-eastern Europe for several days causing rainfall accumulations of up to 300 mm in three days [45,47]. In Serbia alone, these floods affected the territories of 24 municipalities and caused EUR 1.525 billion in damages [48]. Before the 2014 event, there were other widespread floods such as in 1933, 1937, 1940, 1947, 1964, 1974 and 2007 [49]. The most remarkable of these occurred in 1974, when historical record water levels and discharges were exceeded at all stations downstream of the Drina River confluence. These extremes were exceeded by the 2014 flood event.

Given that extreme floods in the Sava River are expected to occur more frequently with climate change, reliable FFA is essential for engineering practice. Accordingly, this paper concentrates on the estimation of extreme river discharges as follows. First, the research area is described with emphasis on the Sava River basin (Section 2). Next, we present the data (Section 3) and methods (Section 4) applied. In the discussion (Section 5), we compare our results with those obtained in similar studies. In particular, we address the following questions. (1) What is the appropriate regional frequency distribution for selected gauging stations along the Sava River? (2) Which frequency distribution should be used for RFFA for ungauged rivers within the Sava River basin? This paper is, to the best of our knowledge, the first RFFA for the Sava River basin. Hence, we also address the question (3) after the study’s novelty and limitations and the practical consequences that can be drawn for decision makers.

2. Research Area

The Sava River basin, with a total area of around 97,000 km², 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³/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³/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 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].

The two rivers Sava Dolinka and Sava Bohinjka (Radovljica station) in Slovenia have an Alpine nival-pluvial regime [54]. Towards the Čatež gauging station, the river regime changes to an Alpine pluvial-nival type [54]. Further downstream, the river regime changes into the Peripannonian pluvial-nival (Zagreb and Jasenovac) and Pannonian pluvial-nival types (Županja and S. Mitrovica) [55].

3. 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. Geographical characteristics for observed stations.

Station	Country	Elevation (m a.s.l.)	Location (River-km)	Basin Area (km2)
Radovljica	Slovenia	408	901	908
Čatež	Slovenia	137	737	10186
Zagreb	Croatia	112	664	12450
Jasenovac	Croatia	87	501	38953
Županja	Croatia	76	262	62891
S. Mitrovica	Serbia	72	139	87966

).

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. Descriptive statistics for the AMS data from the six gauging stations on the Sava River; STD, standard deviation.

Station	n	Mean (m3/s)	Median (m3/s)	Maximum * (m3/s)	STD (m3/s)	Skewness †	Kurtosis †
Radovljica	60	445	433	809	149	0.60	−0.39
Čatež	60	1744	1727	3811	588	0.95	1.49
Zagreb	60	1664	1643	3005	441	0.59	0.07
Jasenovac ‡	56	2078	2089	2814	302	0.36	−0.69
Županja	60	2915	2843	5317	573	1.26	3.66
S. Mitrovica	60	4027	3823	6420	910	0.56	0.49

Notes: * Maximum of the daily values from the interval 1961–2020; ^† a normal distribution has skewness zero (dimensionless) and kurtosis zero (dimensionless); ^‡ missing years (1992–1995) are due to the civil war in former Yugoslavia.

.

4. 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.

4.1. 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]:

$$b_r=E\left\{X{\left[F_X\left(x\right)\right]}^r\right\},$$

(1)

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]:

$$b_0=n^{-1} {\displaystyle \sum }_{j=1}^n{X}_{\left(j\right)},$$

(2)

$$b_1=n^{-1}{\displaystyle \sum }_{j=2}^n X_{\left(j\right)}\frac{\left(j-1\right)}{\left(n-1\right)},$$

(3)

$$b_2=n^{-1} {\displaystyle \sum }_{j=3}^n X_{\left(j\right)} \left(j-1\right)\left(j-2\right)/\left[\left(n-1\right)\left(n-2\right)\right],$$

(4)

$$b_3=n^{-1} {\displaystyle \sum }_{j=4}^n X_{\left(j\right)}\frac{\left(j-1\right)\left(j-2\right)\left(j-3\right)}{\left[\left(n-1\right)\left(n-2\right)\left(n-3\right)\right]}.$$

(5)

In these equations, X_(j) indicates the ranked AMS values, that is, X_(n) is the largest and X₍₁₎ is the smallest value. Furthermore, the first four L-moments [43] are given as follows:

$${\lambda }_1=b_0,$$

(6)

$${\lambda }_2=2b_1-b_0,$$

(7)

$${\lambda }_3=6b_2-6b_1+b_0,$$

(8)

$${\lambda }_4=20b_3-30b_2+12b_1-b_0.$$

(9)

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 (λ₁), scale (λ₂), shape (λ₃) and kurtosis (λ₄) of the distribution functions, a fact that is particularly useful in RFFA. Lastly, the L-moments ratios are calculated as follows [43]:

$$LCv={\tau }_2={\lambda }_2/{\lambda }_1,$$

(10)

$$LCs={\tau }_3={\lambda }_3/{\lambda }_2,$$

(11)

$$LCk={\tau }_4={\lambda }_4/{\lambda }_2.$$

(12)

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 three-parameter distributions), or of the sample mean to the theoretical points (in case of two-parameter 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:

$$V_1={\left\{{{\displaystyle \sum }}_{i=1}^N{n}_i·{({\tau }^{\left(i\right)}-{\tau }^R)}^2/{{\displaystyle \sum }}_{i=1}^N{n}_i\right\}}^{1/2},$$

(13)

$$V_2={{\displaystyle \sum }}_{i=1}^N{n}_i·{\left\{{({\tau }^{\left(i\right)}-{\tau }^R)}^2+{({\tau }_3^{\left(i\right)}-{\tau }_3^R)}^2\right\}}^{1/2}/{{\displaystyle \sum }}_{i=1}^N{n}_i,$$

(14)

$$V_3={{\displaystyle \sum }}_{i=1}^N{n}_i·{\left\{{({\tau }_3^{\left(i\right)}·{\tau }_3^R)}^2+{({\tau }_4^{\left(i\right)}-{\tau }_4^R)}^2\right\}}^{1/2}/{{\displaystyle \sum }}_{i=1}^N{n}_i.$$

(15)

In these equations, N is the number of sites; n_i is the record length at site i; τ⁽ⁱ⁾, τ₃⁽ⁱ⁾, τ₄⁽ⁱ⁾ are the sample values of LCv, LCs and LCk at site i, respectively; and ${\tau }^R$, ${\tau }_3^R$, ${\tau }_4^R$ 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:

$$Z=({\tau }_4^{\mathrm{DIST}}-{\tau }_4^R)/{\sigma }_4.$$

(16)

In this equation, the acronym DIST refers to a particular distribution, ${\tau }_4^{\mathrm{DIST}}$ is the population-LCk value for a particular candidate distribution, ${\tau }_4^R$ is the regional sample average of the LCk values and 𝜎₄ 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 $\left|Z\right|<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,

$$w\left(j\right) = F_n\left(x = X_{\left(j\right)}\right)-F_{\mathrm{theo}}\left(x=X_{\left(j\right)}\right).$$

(17)

While the Kolmogorov-Smirnov test utilizes max{w(j)}, we employed a least-squares criterion,

$$u={(n-p)}^{-1} {\displaystyle \sum }_{j=1}^n {[w\left(j\right)]}^2.$$

(18)

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)⁻¹ 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.)

4.2. 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\left(k\right), X\left(k\right)\}}_{k=1}^n$ denote the time series, where T(k) is time and X(k) is AMS. The linear model is given by

$$X\left(k\right)=a+b ·T\left(k\right)+X_{\mathrm{noise} }\left(k\right),$$

(19)

where a is the intercept parameter, b the slope parameter and $X_{\mathrm{noise}}$(k) the noise component. The minimization of the least-squares sum,

$$SSQ\left(a,b\right)={{\displaystyle \sum }}_{k=1}^n {[X\left(k\right)-a-b ·T(k)]}^2,$$

(20)

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_{\mathrm{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].

5. Results and Discussion

5.1. 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).

5.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. Slope estimates for the AMS data.

Station	$\mathit{b}$ *
Radovljica	0.5 ± 1.1
Čatež	3.3 ± 4.6
Zagreb	1.9 ± 4.5
Jasenovac	3.1 ± 2.4
Županja	−0.7 ± 4.2
S. Mitrovica	−9.6 ± 6.4

Note: * Estimate ± standard error (obtained from 2000 moving block bootstrap resamplings). Unit is m³/s per year.

). 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.

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].

5.3. Homogeneity Assessment and Regional Flood Frequency Analysis

The numerical results of the three homogeneity measures V₁, V₂ and V₃ (

Table 4. Homogeneity test of Sava River basin.

${\mathit{\tau }}^{\mathit{R}}$	${\mathit{\tau }}_3^{\mathit{R}}$	${\mathit{\tau }}_4^{\mathit{R}}$	${\mathit{V}}_{1 }$	${\mathit{V}}_2$	${\mathit{V}}_3$
0.14	0.13	0.14	0.04	0.07	0.077

) 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.

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 $\left|Z\right|<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.

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. Results of goodness-of-distributional-fit measures u (Equation (18); dimensionless) of candidate distributions to observed station data.

Distribution	Radovljica	Čatež	Zagreb	Jasenovac	Županja	S. Mitrovica
GEV	0.3580	0.3516	0.3576	0.3280	0.3516	0.3575
GLO	0.3664	0.3607	0.3754	0.3674	0.3607	0.3785
Gumbel	0.3597	0.3611	0.3589	0.3399	0.3621	0.3612
LNIII	0.3575	0.3518	0.3578	0.3575	0.3515	0.3580
LPIII	0.3827	0.3598	0.3697	0.3798	0.3498	0.3666

Note: Bold type denotes best fit result (i.e., smallest u).

).

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.

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).

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.

6. 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 (Figure 5 and Figure 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 (Figure 5 and Figure 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.