Dealing with uncertainty in the probability of overtopping of a flood mitigation dam

In recent years, copula multivariate functions were used to model, probabilistically, the most important variables of flood events: discharge peak, flood volume and duration. However, in most of the cases, the sampling uncertainty, from which small-sized samples suffer, is neglected. In this paper, considering a real reservoir controlled by a dam as a case study, we apply a structure-based approach to estimate the probability of reaching specific reservoir levels, taking into account the key components of an event (flood peak, volume, hydrograph shape) and of the reservoir (rating curve, volume–water depth relation). Additionally, we improve information about the peaks from historical data and reports through a Bayesian framework, allowing the incorporation of supplementary knowledge from different sources and its associated error. As it is seen here, the extra information can result in a very different inferred parameter set and consequently this is reflected as a strong variability of the reservoir level, associated with a given return period. Most importantly, the sampling uncertainty is accounted for in both cases (single-site and multi-site with historical information scenarios), and Monte Carlo confidence intervals for the maximum water level are calculated. It is shown that water levels of specific return periods in a lot of cases overlap, thus making risk assessment, without providing confidence intervals, deceiving.

An important application of this multivariate analysis is the determination of the risk of failure of a hydraulic structure.De Michele et al. (2005) were the first to check the adequacy of a dam's spillway under a bivariate hydrological load, followed by Requena et al. (2013), while Volpi and Fiori (2014) formalised the idea that the return period of a failure of a structure depends on the structure of interest, and therefore the interaction between the hydrological loads and the structure should be taken into consideration by fixing a "structure-based" return period.
In the same conceptual framework, Serinaldi (2016) suggested that the choice between a univariate and multivariate risk assessment should not be based on whether one or the other overestimate/underestimate the risk but rather on the operational criteria of the problem, or simply on what is the mechanism of failure.
Copulas are functions that combine marginal distributions with the joint cumulative distribution; therefore, the latter is only indirectly affected by the choice of the marginals.So, the practical problem of identification and estimation of the joint distribution is handled from two non-interwinding aspects: the dependence structure of the set of variables and the marginal distributions.
In the majority of the studies, the communication of the sampling uncertainty -an integral component in a univariate framework -is overlooked in a multivariate case.Serinaldi (2013) studied the effect of sample size on the confidence bands of the probability of exceedance curves of a joint peak-volume event and showed that in small and medium sample sizes these curves largely overlap.Similarly, Zhang et al. (2015) implemented a Bayesian inference approach to account for the uncertainty of parameter estimation and for the occurrence of a drought, coming to the conclusion that the 95 % confidence interval of a 20-year event can include the expected values of 10 to 50-year events.
In order to account for the sampling uncertainty of multivariate cases, where a variable of interest can be expressed as a function of one or more variables, Serinaldi (2016) proposed a Monte Carlo procedure.He also underlined the importance of including confidence intervals when providing point estimates of the variable of interest, which is even more necessary in the multivariate frequency analysis, where the unknown dependence structure contributes to the uncertainty.
The sampling uncertainty in a joint peak-volume event was quantified by Dung et al. (2015) who used two bootstrapping methods -one developed by the authors and the second by Serinaldi (2013) -and concluded, as the previous, that the model selection and parameter estimation methods are of minor significance in uncertainty estimation with respect to sampling uncertainty, even in relatively large sample sizes.They suggested that efforts should be focused on the expansion of the data set in order to achieve a reduction of uncertainty.
The data expansion can be temporal, spatial and causal (Merz and Blöschl, 2008), thus enriching the available evidence with information from neighbouring basins, previous periods and by comprehension of the flood-generating mechanisms.In the past, many researchers (Parent and Bernier, 2003;Reis Jr. and Stedinger, 2005;Gaume et al., 2010;Halbert et al., 2016;Parkes and Demeritt, 2016;Viglione et al., 2013) have dealt with the extension of the available data using information on paleo-floods, historical flood reports, marks of the river stage during important flood events, expert judgement, etc. with the aim of reducing the range of uncertainty bands or simply to reach a more realistic design value.Bayesian inference allows the integration of information from different sources and their associated uncertainty and errors and provides a means of conveying hydrological reasoning in a mathematical context.
In this research, we validate a methodology of flood risk assessment in a real case study, where risk is expressed in terms of probability of exceeding a given reservoir level in an online flood mitigation dam.We consider this level as a function of flood peak, volume and hydrograph shape and, consequently, multivariate modelling is implemented with the use of copulas.The characteristics of the reservoir -also a function of the level -are synthesised in the rating curve and the volume-level curve.The main scope is to integrate the associated sampling uncertainty and to build confidence intervals for each water level through Monte Carlo simulations.Furthermore, we incorporate additional information on the peaks in a Bayesian framework and we examine its effect on the distributions, their confidence intervals, as well as the ones of the reservoir-level frequency curve.

Case study and data
We have focused our interest on the Panaro catchment -an important influent of the Po river in northern Italy.In particular, the watershed under investigation is composed of the Panaro river itself, the Scoltenna and the Leo tributaries with an outlet upstream of the Panaro dam (Fig. 1) and occupies an area of 876 km 2 .The Panaro has its source at Monte Cimone (2165 m a.s.l.) and flows into the Po at Bondeno; it takes its name from the valley of Montespecchio after converging with the Leo and Scoltenna streams that constitute the upper part of the river network.The hydrographic network of the watershed shows a low degree of hierarchy, indicating an evolving state which is also evident by the existence of torrential dynamic phenomena (Autorità di bacino del fiume Po, 2006).
The influence of snowfall is negligible due to the modest land elevation and the majority of rainfall events occurring seasonally (September-April).The average precipitation height ranges between 700 and 2000 mm yr −1 (Autorità di bacino del fiume Po, 2006).
The basin's permeability is low, and therefore the runoff is influenced little by water infiltration.In fact, the study basin consists mostly of sandstones and silicatic alternating sequences (44 % of total area) and marls and clay (34 %).
The Panaro dam is a concrete gravity dam (150 m in length), located near the city of Modena and constructed for flood mitigation purposes.The hydraulic system consists of two reservoirs, a principal on the river course and a secondary at the right river bank, and a series of levees that enclose them.The crest of the principal levees is at 44.85 m a.s.l.The reservoirs can hold in total 23.66 hm 3 up to the spillway's crest at 41.1 m a.s.l.There are also nine discharge outlets at Hydrol.Earth Syst.Sci., 21, 2497Sci., 21, -2507Sci., 21, , 2017 www.hydrol-earth-syst-sci.net/21/2497/2017/ the bottom of the dam that ensure constant flow to the downstream valley.
The available flood data included a 52-year discharge series (1936-1943, 1945 and 1946 were missing) with an hourly time interval from the Bomporto station located downstream near the current location of the dam.The hydrological characteristics of the study basin are briefly presented in Table 1.
Additional data included the annual peaks of the missing years from the same station (Servizio Idrografico Italiano, 1939, 1953) and recent annual peaks from upstream stations, after consulting the annual hydrological reports of ARPA -Emilia Romagna published in its website (http://www.arpae.it), in particular, for 2003 from Pievepelago, for 2004, 2005 and 2015 from Ponte Samone, and for 2006 to 2014 from Spilamberto.
In a report about natural disaster risks in the city of Modena (Nora and Ghinoi, 2009), it is also stated that the most disastrous flood events of the 20th century happened during the last 40 years (1966, 1969, 1972 and 1973).In November 1966, the flooded area from the Panaro covered 9400 ha; in September 1972, it covered 2540 ha; and in September 1973, it covered 6000 ha (Nora and Ghinoi, 2009).

Data regionalisation
In order to rescale the flood information from subcatchments and from the downstream station, depending on the area, the following scale function was used: (1) where Q(A 1 ) is the rescaled discharge, Q(A 2 ) a known discharge and m a regional-scale exponent.De Michele and Rosso (2002) clustered basins with similar flood generation mechanisms and checked the homogeneity of the grouped regions.The study area was located in northwestern Italy and included the Panaro watershed.The proposed scale exponent m for this region is 0.772 with a standard deviation of 0.072.In our case, the rescaling regarded different locations of the same basin, although in theory neighbouring basins could have been used (e.g. the Secchia basin), but they did not add additional information here.

Incorporating additional data
Thomas Bayes' theorem expresses how an individual's degree of belief can change after the presence of new evidence.Bayes' theorem can be formulated as where π(θ ) is the prior density distribution of the parameters θ , p(θ|D) is the posterior distribution after the introduc- 1449.9 t c (h) 10.0 tion of the observed information D and l(D|θ ) is the likelihood of the data.The denominator serves only as a normalisation constant to ensure unity of the area under p(θ |D), so the equality sign can be substituted with the proportionality sign.This integral cannot be solved analytically, so for its computation Monte Carlo Markov chain algorithms are employed.In each Markov chain, the aim is the maximisation of the logarithm of the unnormalised joint posterior distribution starting from an initial value and proceeding iteratively, in order to arrive at each target distribution (Statisticat, LLC, 2016).
In a Bayesian framework, the model's parameters are handled as stochastic variables in order to incorporate the uncertainty of their values (Ouarda and El-Adlouni, 2011).In the present case, the model's parameters were the parameter of the peak marginal distribution and the scale exponent.We used noninformative prior distribution for the marginal and a normal prior distribution for the exponent m N (0.772, 0.072).We integrated a perception threshold X P -a value which only in k number of years in a historic period of h years was exceeded; here, it is set as at about 1000 m 3 s −1 , a value thought to be exceeded only once in a historic period of 117 years since flood reports indicate that during the early years of the 1900s, when systematic records were nonexistent, the flood events were of less significance in comparison with the events occurring in the 1970s.Additionally, we have introduced the uncertainty of the scale exponent m.Therefore, the likelihood function of the data is set as where is the binomial coefficient, s is the number of different sites of the recorded flood peaks, n i is the number of recorded peaks for each site and y ij are the annual peaks from the different sites.
The Bayesian inference was conducted in R with the package LaplacesDemon (Statisticat, LLC, 2016) and the MCMC algorithm utilised was the Componentwise Hit-And-Run Metropolis.The logarithm of the posterior distribution to be maximised is the sum of the logarithm of the likelihood and the logarithm of the priors: where µ and σ are the mean and shape parameters of the peak distribution.

Copula and marginal inference
Copulas are functions that describe the dependence structure between variables independently of the choice of marginal distributions.The joint distribution functions and the marginals are linked by Sklar's theorem (Sklar, 1959): for all x ∈ R d , where the F i , i = 1, . . ., d are the marginals of F and C is the copula function.
Copulas provide a powerful tool for the statistical modelling of multivariate data: for a theoretical introduction, see Nelsen (2006), Joe (2014) and Durante and Sempi (2015); for a practical engineering approach, see Genest and Favre (2007), Salvadori et al. (2007) and Salvadori and De Michele (2007).
The application of copula functions has facilitated overcoming some inadequacies of traditional multivariate distributions such as that the marginals must derive from the same distribution family and their parameters may define the dependence structure between the variables (Salvadori et al., 2007).
The degree of relation between pairs of variables was examined by measures of association.These include Kendall's τ , Spearman's ρ S (which expresses the existence or absence of concordance) and Pearson's ρ P (which expresses linear dependence).For the observed discharge/volume pairs, these were equal to 0.58, 0.77 and 0.81, respectively, and indicate strong dependence.
In the absence of a long sample, the copulas that fit the data can be numerous and goodness-of-fit tests cannot distinguish between them (Serinaldi, 2013).Since inferring the "correct" copula model is not the aim of this research and since this endeavour at this point can be futile, given the available data set, the final choice was based partly on the preference of previous published research of the Gumbel distribution, including conference proceedings by Balistrocchi et al. (2014) who fitted the Gumbel on peaks obtained from a peak-over-threshold method on the same discharge time series.In the present case, both the Gaussian and the Gumbel-Hougaard one-parameter copulas passed the goodness-offit tests (Cramér-von Mises and Kolmogorov-Smirnov) and demonstrated the smallest Akaike weights -or else the probability that the chosen model is the most apt among the tested ones (Burnham and Anderson, 2004).However, we thought that if tail dependence exists, Gumbel would be more appropriate (belonging to the extreme value copula family), as the Gaussian has no tail dependence.We recall the Gumbel-Hougaard copula as where u, v are the pseudo-observations and θ the copula parameter.
The existence of tail dependence between peak and volume was also implied by some historical evidence.Many significant events in Italy occur when a frontal perturbation, generated by the cold high masses coming from the North Atlantic Ocean or the Arctic Ocean, meets Mediterranean southward warm fronts.Depending on the persistence of the south and north currents, the generated front begins to develop, covering a large area (e.g. 10 4 km 2 ).Inside this warm front, the energy content is very high.This causes local convective phenomena enhanced by orographic effects.So, thunderstorms can appear locally producing rainfall whose values can surpass one-third of the mean annual in 24-30 h.In the vicinity of the local thunderstorms, the rainfall is moderately high, producing large soil saturation and increasing, significantly, the contribution to the groundwater.This kind of rainfall event produces not only maximum observed peaks of flood in many rivers of small (< 100 km 2 ) and medium (< 2000 km 2 ) sizes but also the largest observed volumes associated with the persistence of the global event.This is the case, for example, for the flood in Florence and Triveneto on 4 November 1966, in Valtellina on 18-25 July 1987, along the Tanaro on 5-6 November 1994, along the Po in Piedmont on 17-21 October 2000, etc.
Unfortunately, tail dependence estimators such as the ones of Frahm et al. (2005) and Schmidt and Stadtmüller (2006) can be biased and susceptible to high uncertainty even in large sample sizes (Serinaldi et al., 2015); thus, their use in this case is discouraged.
Regarding the choice of the marginal distributions, we preferred distributions that were more parsimonious, thus reducing the additional statistical uncertainty introduced by an extra parameter, following the logic of Occam's razor, and this provided a nice visual fit.The differences between the corrected Akaike information criterion (AIC), Bayesian information criterion (BIC) and Akaike-weighted values were not sufficient to make a safe distinction between the models.The peaks were modelled with the inverse Gaussian distribution (two parameters instead of three of the generalised extreme value (GEV) distribution) and their corresponding volumes with the one-parameter Rayleigh.It is, however, imperative to note that there was no clear indication of overall performance superiority of the chosen distributions.
The parameters of the inferred distributions (copula and marginals) are presented in Table 2.

Hydrograph selection
The shape of the "design hydrograph" is often considered an important factor in the design procedure and is related to the spatial and temporal rain distribution as well as the Hydrol.Earth Syst.Sci., 21, 2497-2507, 2017 www.hydrol-earth-syst-sci.net/21/2497/2017/ basin's shape and behaviour (Singh, 1997).Therefore, typical hydrographs were determined from the annual maxima flood events extracted from the available time series.In particular, the events' hydrographs were clustered according to their characteristics and utilising the methodology proposed by Dyck and Peschke (1995), which suggests the normalisation of the hydrograph (after the removal of the baseflow) by where Q max is the hydrograph's peak, q base the base flow and t max the time to peak, starting from the rising limb of the flood event.
Consequently, the normalised peak equals 1 at time 1.All normalised hydrographs were extended to a common duration (for comparison purposes) and cluster analysis with the Ward method and Euclidean distances was implemented (Aronica et al., 2012).

Structure-based risk analysis
Volpi and Fiori (2014) associated the return period with the structure of interest by relating the structure design parameter to the hydrological load through the function Z = g(X).Consequently the structure-oriented return period of a value z takes the form where µ T is the mean interarrival time between two consecutive occurrences of z (in our case, µ T = 1 year), F Z is the probability distribution function of the derived variable Z, which in this case is the reservoir level and X ≡ (Q, V , shape, etc.).Here, the structure function is very complex, since the reservoir level is a function of the spillway's rating curve and the flood's natural variables, and therefore the whole analysis must be based on Monte Carlo simulations.This adds to the computational burden, specifically when dealing with the quantification of the uncertainty.

Uncertainty estimation
In order to account for sampling uncertainty and to estimate the confidence intervals, the following Monte Carlo procedure was implemented, originally proposed by Serinaldi (2016).
1. Estimate the parameter θ of the copula for the observed sample as well as the parameter of the flood volume distribution.
2. Simulate B bivariate samples of size n (equal to the number of years of the observed sample) using the estimated copula parameter and then transform into volume using the estimated parameter of the marginal.
3. Calculate the copula parameter θ and the volume marginal parameter for each sample with the same estimation method used for the observed sample.
4. Simulate B bivariate samples of size M with the copula parameter θ estimated in the previous step.
5. Transform the samples from the unit interval to discharge and volume using the estimated marginal parameters.Generate B × M hydrographs with an assigned peak, volume and shape and route them to calculate the reservoir level and the frequency curves of all B samples.
6. Build the confidence intervals of the reservoir level frequency curves.
In the present research, B was set equal to 10 000 and M equal to 1000.The confidence intervals of the peaks' marginal distribution parameters have been estimated in a Bayesian framework, as stated previously, in order to incorporate the additional knowledge and to account for the scaling uncertainty.
The parameter uncertainty of additional distributions that fit the data could be introduced in the procedure, leading to larger confidence intervals.However, in this case, only the parameter uncertainty of the inferred models was of interest.Maximum water level (m) q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q qq q q q q qq q q q q q q q q q q q q q q Generated (clusters=4) Generated (clusters=1) Observed

Results and discussion
Initially, we have clustered the hydrograph shapes into four characteristic groups.After simulating 10 000 peak-volume pairs from the inferred distributions, we assigned to each one a specific hydrograph shape (respecting their frequency of occurrence).Then, we denormalised and routed the hydrographs; we repeated the same procedure but after clustering into only one group, thus considering a global mean hydrograph.The characteristic shapes are depicted in Fig. 2a, along with the mean shape (Fig. 2b).The level frequency curve showed small differences between the two cases (Fig. 3) and,  as we shall see later, this difference is negligible in comparison with the estimated uncertainty.So, we have proceeded with the uncertainty assessment considering one mean hydrograph.
We have implemented the Bayesian framework on the peaks extracted from the systematic discharge series recorded at the Bomporto station, adding also the uncertainty of the scaling exponent of the regionalisation relation.In the second scenario, we also included recently recorded annual peaks from other hydrometric stations of the same basin, mentioned previously, as well as information from flood re- ports, while integrating the uncertainty of the scaling exponent.As it can be seen in Fig. 4, when ignoring the additional information, the estimate of discharge peak is lower for a given return period.We focus our attention on relatively small return periods since any extrapolation beyond the available time period is subject to great uncertainty.Indicatively, a peak with a univariate return period of 50 years can increase by 18 %, exceeding the confidence intervals of the fitted distribution.This occurs because during the last 10 years crucial flood events appeared in the area.In Table 2, we note that in the second case the distribution's mean is bigger and the shape parameter is smaller with a significantly lower standard deviation.
The 95 % confidence interval of both of the peak distributions can be wide (Fig. 4); e.g. for a univariate return period of 50 years, it can span from 665 to 941 m 3 s −1 (29 % difference) and from 837 to 1092 m 3 s −1 (23 % difference) for the first and second cases.
Similarly, in the case of the flood volume, the 95 % confidence interval for a univariate return period of 50 years can span from 95.5 to 120 hm 3 (20 % difference) (Fig. 5).
The confidence intervals of the parameters of the inferred distributions are presented in Table 2.
The results of the increased peaks are reflected also on the frequency curve of the maximum water level (MWL).The return period here corresponds to a water level, so it is considered structure based, since the level is a function of the structural and operational characteristics of the dam, among others.As it can be seen (Fig. 6), the MWL is significantly lower in the case of no extra information, especially for greater return periods.For a return period of 50 years, the average MWL can differ by 1.2 m -a considerable magnitude in terms of volume and when considering that the safety margin above the spillway's crest is in some cases 1 m.In Fig. 7, the highest density regions depict a sort of confidence interval of the MWL for specific return periods for each case (single-site information and multi-site with historical information cases).These regions can be defined as the smallest areas in the sample space with a certain probability coverage and they have the advantage of displaying multimodal distributions; thus, they may consist of disjoint subsets (Hyndman, 1996).They are particularly suitable in multivariate cases or for asymmetrical distributions.In the highdensity region (HDR) box plots, the mode (horizontal line) substitutes the median and the darker region corresponds to a probability coverage of 50 %, the lighter to a coverage of 95 % and the points outside to the data beyond the 95 % probability.
The span of the highest density regions slightly decreases as more information is introduced.However, this decrease in uncertainty seems small and we cannot come to a conclusion, whereas the extra information has contributed to a systematic uncertainty reduction.An increase in the simulation size could lead to a slightly different picture; however, the added computational burden is prohibitive and anyhow, as Salvadori et al. (2015) stated, a clear rule of thumb regarding the simulation size does not exist.In any case, the size is considered large enough for some safe conclusions, especially for the smaller return periods.Within a certain return period, the parameter uncertainty can lead to substantial MWL variations, e.g. for 20 years (in the case of extra information), the span of the MWL with a density of 50 and 95 % is of 0.8 and 2.3 m, respectively, which correspond to huge volume differences.These can have devastating effects not only in the case of overtopping but also because the remaining water can cause bank failure due to piping.These spans could increase for larger return periods, where the uncertainty is bound to get vaster.
As the results of previous studies suggested (Serinaldi, 2013;Dung et al., 2015;Zhang et al., 2015;Serinaldi, 2016), the regions of the return periods can overlap; in this case, the 95 % confidence interval of an event of 20 years can include, marginally, the expected values of events of 10 to 50 years (e.g.single-site scenario).For the multi-site with the extra information scenario, the overlapping region is smaller.
In Fig. 8, the highest density regions (50, 75 and 95 %) are depicted in a two-dimensional plane (discharge-volume) that correspond to a return period of 50 years.For events that result in a MWL with a specific return period, the variation of discharge and volume can be huge even when looking in the smaller density regions (e.g.50 %).For example, in the 95 % region, the discharge can assume values with a univariate return period from 1 to 50 years (for the multi-site scenario), which is a strong indicator of the nonlinearity of the problem.The same applies for the flood volume.

Conclusions
This analysis focuses on the uncertainty introduced when calculating the probability of exceeding specific water levels in a flood control reservoir, which is a result of the parameter Hydrol.Earth Syst.Sci., 21, 2497Sci., 21, -2507Sci., 21, , 2017 www.hydrol-earth-syst-sci.net/21/2497/2017/ uncertainty of the marginals of the hydrological variables, as well as the copula multivariate function, due to the small size that characterises, in most cases, a hydrological sample.Therefore, we attempted to quantify this uncertainty, without aiming our attention to copula/marginals inference.Instead, we studied the effect of additional flood information not only on the distribution parameters but also on the uncertainty range in a Bayesian framework that, among others, permits the consideration of errors from different sources.
The extra flood data that included additional peaks from different hydrometric stations led to a peak distribution with bigger mean and smaller shape parameters and thus to elevated peaks, since the data include flood events of recent years that exceed the events of the historical data series in magnitude.Consequently, including the additional information translates into a general bigger estimate of the peaks, which is also reflected on the MWLs, as the peak is a driving factor of the routing process.
The uncertainty range of discharge and volume is considerable and affects, along with the copula parameter, the MWL.The variations in the MWL for the same structurebased return period correspond to significant variation in the stored water volume.Most importantly, the return period of a specific water level cannot be determined with certainty because the return periods of the events overlap.Naturally, the range of discharge and volume values for a given structurebased return period is very ample due to the wide range of the parameters of the inferred distributions.
A clear observation of whether uncertainty is systematically reduced with the introduction of additional information cannot be made here.Nonetheless, a Bayesian framework allows a certain degree of transparency (Parkes and Demeritt, 2016).Incorporating knowledge about water levels during historic events, e.g. at the locality of Navicello for the 1783 and 1842 floods, could result in a more significant change in the uncertainty range, as past research has shown.However, one must consider also the great amount of error involved in these data in the Bayesian framework.
As a general remark, one can deduce that the process of risk estimation is inherently crippled by uncertainty that can be quantified or at least approximated.Any attempt to obscure this uncertainty could create a false notion about its existence in a multivariate problem with eventual implications in dam safety.

Figure 2 .
Figure 2. Characteristic normalised hydrograph shapes (a) with a certain probability of occurrence and (b) "mean" normalised hydrograph shape.

Figure 3 .
Figure 3. Frequency curves of maximum water level of synthetic hydrographs for four clusters and one cluster and corresponding levels of observed hydrographs.

Figure 4 .
Figure 4. Flood frequency curves with 95 % confidence limits for the single-site data and the multi-site data with the extra information case.Observed peaks are also plotted with the Gringorten plotting position.

Figure 5 .
Figure 5. Flood volume curve with 95 % confidence limits.Observed flood volumes are also plotted with the Gringorten plotting position.

Figure 6 .
Figure 6.MWL frequency curves with 95 % confidence limits for the single-site data and the multi-site data with the extra information case.

Figure 7 .
Figure 7. High-density region box plots of MWL for return periods of 10, 20 and 50 years for (a) single-site data and (b) multi-site data with historical information.

Figure 8 .
Figure 8.The 50th, 75th and 95th percentiles of the kernel density areas of MWLs with a return period of 50 years on the discharge-volume plane for single-site data (a) and multi-site data with historical information (b).

Table 1 .
Main hydrological characteristics of the Panaro watershed (area A, main stream length L, minimum elevation H min , elevation drop H , time of concentration t c ).

Table 2 .
Estimated parameters of the inferred distributions and their confidence interval (95 %).