Bivariate Return Period for Design Hyetograph and Relationship with T-Year Design Flood Peak

This study focuses on the return period evaluation for design hyetographs, which is usually estimated by adopting a univariate statistical approach. Joint Return Period (JRP) and copula-based multivariate analysis are used in this work to better define T-year synthetic rainfall patterns which can be used as input for design flood peak estimation by means of hydrological simulation involving rainfall-runoff (RR) models. Specifically, a T-year Design Hyetograph (DH) is assumed to be characterized by its peak H, at the chosen time resolution ∆t, and by the total rainfall height W, cumulated on its critical duration dCrit, which has been a priori fixed. As stated in technical literature, the choice of the expression for JRP depends on which event is deemed as critical for the investigated system; the most important cases are: (i) all the variables must exceed a certain magnitude to achieve critical conditions; or (ii) at least one variable must be greater than a threshold; or (iii) critical conditions are induced by all the events with a joint Cumulative Density Function (CDF) overcoming an assigned probability threshold. Once the expression for JRP was chosen, the relationship among multivariate T-year design hyetographs and T-year design flood peak was investigated for a basin located in Calabria region (southern Italy). Specifically, for the selected case study, a summary diagram was obtained as final result, which allows the main characteristics of T-year DHs to be estimated, considering both the univariate and the copula based multivariate analysis, and the associated T-year design flood peaks obtained through the simulation with a RR model.


Introduction
Many hydrological issues require the use of design rainfall models in order to evaluate the effects of intense precipitation events over a basin or urban area.With this aim, the most adopted practice is to generate Design Hyetographs (DHs), which are consistent with the rainfall time series of the study area.
In general, a T-year DH is characterized by several variables: (i) the rainfall peak H, cumulated on a time resolution ∆t; (ii) the critical duration d Crit ; (iii) the total rainfall height W, cumulated on d Crit ; (iv) the hyetograph shape and in particular the peak position in the temporal pattern.
A common approach, often adopted by technicians, is represented by the univariate analysis: H and W are evaluated from a T-year Amount-Duration-Frequency (ADF) curve, in which the d Crit to be adopted can be estimated by analyzing historical flood events, or it can be assumed equal to the time of concentration of the watershed [1,2] or, as a further alternative, it is set to the duration that leads to maximum peak discharge after the application of a rainfall-runoff model.In general, most hydrological analyses that are aimed at determining peak discharge for hazard or risk assessments, especially in ungauged catchments, often require the estimation of catchment response time parameters as primary input, and large errors can be ascribed to errors in the estimation of these crucial characteristics [3].One of the time parameters most frequently used to express catchment response time is the already mentioned time of concentration.In acknowledging this, it is worthy to note that a universally accepted working definition of this parameter is currently lacking and several definitions can be found in the technical literature along with related estimation procedures [4][5][6][7].In this paper, in accordance with the following theoretical definitions often used in the technical literature [8], we intend the time of concentration as the time from the end of rainfall excess to the time of the end of direct runoff.
However, all of the above listed quantities related to a DH are characterized by an intrinsic variability, and as such a multivariate approach should be more suitable.For this reason, in this paper the authors carried out a multivariate analysis by using copula functions [9], with the aim of assessing the impact of different DHs in the estimation of design flood peak.Applications of copulas in geosciences, and in particular for hydrological studies, are receiving increasing attention in scientific and technical literature (see for example [10] and references herein).The main and well-known advantage of using copula functions is constituted by the opportunity to split the analyses of marginal distributions (which can be different from one variable to another) and of the dependence structure among variables.
In this work a bivariate analysis is discussed: the critical duration d Crit is fixed equal to the time of concentration of the chosen watershed, and the random variables in the bivariate analysis are the rainfall peak H(cumulated on ∆t = 20 min) and the net volume W n = W − H.The latter variable is chosen to eliminate the apparent correlation between H and W [10].The basin of Corace at Grascio (located in Calabria region, southern Italy) was selected as case study.In this analysis, all the DHs are assumed to have the rainfall peak H in correspondence of d Crit /2, while W n is uniformly distributed on all the remaining time intervals with a resolution ∆t = 20 min.
Another point that is crucial for the estimation of DHs is the definition of the Joint Return Period (JRP) involving the considered random variables.In fact, different expressions can be obtained depending on which event is deemed as critical for the investigated system [11,12]; the most important cases are: (i) all the variables must exceed a certain magnitude to achieve critical conditions; or (ii) at least one variable must be greater than a threshold; or (iii) critical conditions are induced by all the events with a joint Cumulative Density Function (CDF) overcoming an assigned probability threshold.Once the specific expression for JRP was chosen for the selected case study, the authors provided as a final result a summary diagram, which allows for: (1) the identification of the main characteristics of T-year DHs; and (2) the analysis of the relationship among T-year DHs and T-year design flood peaks, obtained through a simple lumped Rainfall-Runoff (RR) model.
The paper is organized as follows: Section 2 provides a brief overview of copula functions suitable for a multivariate analysis; Section 3 regards the definition of a bivariate return period, while the selected study area is described in Section 4. Section 5 presents the main results and conclusions are drawn in Section 6.

Overview of Copula Functions
Copula functions constitute an efficient tool for multivariate analysis: their use allows for the analysis of global dependence among the variables, with marginal distributions studied separately (i.e., each distribution can assume a different mathematical expression).In this paper, a brief overview of 2D copula theory is provided.For further details, in particular about the extension to multivariate cases, see Nelsen [9].
Let F H (h) and F W n (w n ) be the marginal Cumulative Density Functions (CDFs) for H and W n ; on the basis on Sklar's theorem [13] the joint CDF can be written as: Water 2017, 9, 673 3 of 13 where U H = F H (h)and U W n = F W n (w n ) are standard uniform random variables, and C is a 2D copula function, characterized by the following properties: H , u W n , u H and u (1) For hydrological applications, Archimedean copula functions are mainly used, which are defined as: where ϕ(.) is the generator, a continuous function strictly decreasing and convex from I = [0, 1] to [0, ϕ(0)].In this work, two particular kinds of Archimedean copulas were considered: Gumbel-Hougaard with ϕ(t) = (− ln t) ϑ and ϑ ≥ 1, and Frank Copula ϑ−1 and ϑ > 0.
In both expressions, the parameter ϑ is an indicator of the dependence between the variables.In particular, in Equation ( 3) U H , U W n are independent for ϑ = 1 and positively correlated for ϑ > 1, while in Equation (4) U H , U W n are independent for ϑ = 1, positively correlated for 0 < ϑ < 1, and negatively correlated for ϑ > 1.
Evaluation of ϑ can be carried out by using the Canonical Maximum Likelihood (CML) method [14][15][16][17].After parameter estimation, the choice of the copula function that best fits the sample can be made by considering the following C-measure, named as Kendall function [16,18,19]: which represents the CDF of Probability plots (comparing empirical and theoretical CDFs of a copula function) and Q-Q plots (in which K C (t) is plotted against standard uniform quantiles U(0, 1)) can then be used to assess the performance of the adopted copula.

Bivariate Return Period Definition
Under the hypothesis of stationary process herein considered, in technical literature the return period T (expressed in years) is commonly calculated as: where P represents the probability of observing realizations exceeding reference threshold values, and µ (expressed in year) is the average inter-arrival time between two consecutive realizations.When annual maxima are considered, µ is equal to 1 year.If a univariate analysis is carried out, by considering, for example, the variable H, then Equation ( 6) becomes: Water 2017, 9, 673 4 of 13 and it is clear that the critical event [H > h] uniquely defines the denominator 1 − F H (h).
Conversely, when a bivariate analysis is considered, the attention is focused on both variables, in this case H and W n , and the critical event can be defined in several ways [11,[20][21][22][23]: indicated as an OR event, for which the corresponding return period T OR (h, w n ) is: • [H > h ∩ W n > w n ], indicated as an AND event, for which the corresponding return period T AND (h, w n ) is: , indicated as a COND event, for which the corresponding return period T COND (h| w n ) is: , indicated as a KEN event because the probability of this event can be computed by using the Kendall function.The associated Kendall return period T KEN (h, w n ) is: from which: It can be demonstrated [2,18] for an assigned pair (h, w n ) that: As stated in Shiau [11] and Serinaldi [12], the choice for the expression of T depends on which event is critical for the investigated system.If both variables must exceed specific values to achieve critical conditions, then T AND (.) should be adopted.On the contrary, if it is enough that only one variable is greater than a threshold, then T OR (.) is more suitable.If the exceedance of a variable is conditioned on the exceedance of the other one, then T COND (.) should be preferred.T KEN (.) should be used when critical conditions are induced by all the events (h, w n ) with a joint CDF greater than an assigned threshold (see Figure 1 in Serinaldi [12] and Figure 2 in Gräler et al. [22] for further graphical details).
Moreover, for a fixed expression written for the bivariate return period (Equations ( 8)-( 11)), it is clear that an infinite set of points (h, w n ) belonging to a contour line on the OHW n coordinate system is associated to a particular value of T.Although all the events on a contour line have the same return period, some combinations (h, w n ) can greatly differ in terms of their magnitudes; concerning T OR (.) (adopted for the selected study area), at the edge of a contour line one variable tends to its marginal value, while the other indefinitely increases, and the pairs (h, w n ) differ also in terms of joint probability density function (PDF) evaluated along the contour line [24].In this context, users can consider the most-likely design realization that maximizes the joint PDF [25], or an ensemble of design Water 2017, 9, 673 5 of 13 events along a portion of the contour line representing an assigned confidence band, in order to reflect the variability within the set [24].The catchment response has been modelled using a lumped and event-based RR model [26,32] composed by the Soil Conservation Service Curve Number (SCS-CN) method [33] for rainfall excess estimation, and by Nash's model [34] for runoff estimation at basin outlet.Although the SCS-CN method is a popular (due its relative simplicity and its reliance on a limited number of parameters, resulting in a robust tool for those catchments that are partially or poorly gauged) and ubiquitous means for estimating storm runoff from a given rainfall event, i.e., as an alternative to empirical or physically-based infiltration models like Green-Ampt or Horton's equations, it should be noted that its usage at the sub-daily time scale is controversial [35].The major conflict is that the classic SCS-CN procedure allows for the estimation of cumulative rainfall excess (and thus cumulative losses) with rainstorms, but does not take time dimension into account in its equations.Several attempts have been made in order to overcome the conceptual limitations and inconsistencies of the original approach, for example, making explicit the concept of soil moisture accounting in a such a manner that is more suited for continuous simulations of rainfall-runoff processes [36,37] or proposing a mixed approach like the CN4GA, which aims to include the event scale correct information from the SCS-CN into the Green-Ampt infiltration equation [38].Nevertheless, the SCS-CN method is currently indicated by the River Basin Authority of Calabria region as a part of the procedure for design flood estimation within the study region and thus was used in the following sections of the paper.

Results and Discussion
Figure 2 shows the good comparison (on EV1 probabilistic plot) between the frequency distributions of observed samples and a generated series of spatially averaged values of 20-min annual maxima H. Figure 3 illustrates empirical ADF values and theoretical ADF curves, both obtained from the generated 500 years of synthetic data, for return periods T equal to 50, 100, and 200 years.In Figure 3, R , is the rainfall height relative to a duration d and a return period T. As regards the theoretical expression of ADF curves, a power formula was used [26], which is commonly adopted in Italy.

Data and Materials
The case study is the basin of the Corace River at Grascio, located in Calabria region (southern Italy, Figure 1).The watershed area is 177.34 km 2 , its mean elevation is 822 m above sea level, the length of the longest drainage path is 43.84 km, and the mean annual precipitation is 1173.46mm.The gauging network (Table 1) comprises five rain gauges (Tiriolo, Albi, Taverna-Ciricilla, Parenti, and Catanzaro) and one stream gauge (Grascio).The authors considered the time series of daily and 20-min rainfall height, annual maximum peak discharges, which are available for the periods 1916-2016, 1990-2016, and 1925-1970, respectively.The time series of rainfall data were also used to reconstruct the sample sequences of the volume W associated to each 20-min rainfall annual maximum.All the sample sizes are indicated in Table 1.The time series of the average daily and 20-min rainfall data over the basin, which are of interest for application, were obtained by using the technique of Voronoi polygons.Then, the two-stage rainfall generator described in Biondi and De Luca [26] was used in order to obtain a 500-year time series of 20-min rainfall height.In detail, as the size of 20-min data is usually shorter than the size of daily data, firstly a daily rainfall generator, whose parameters are estimated on the longer daily sequence, is used and then a downscaling scheme, developed for southern Italy [27] and calibrated on shorter fine-scale records, is applied in order to disaggregate the generated daily data into 20-min rainfall values.From each generated year of data, the corresponding annual maxima for H and W were computed: W was estimated as the maximum aggregated rainfall (comprising H) along the moving window of duration d Crit , which was set equal to the time of concentration of the basin, 7.33 hours.This value was taken from the VAPI (VAlutazione delle PIene, Italian acronym of Flood Estimation) report for the Calabria region [28], in which the times of concentration for some basins, including the Corace river, were derived from available rainfall-runoff events data set.Then W n was derived as W n = W − H.Many papers have highlighted the importance of rainfall scenarios with a high peak (and consequently have used H) in order to obtain the maximum discharge.As an example, Alfieri et al. [29] focused on the application of five different design hyetographs: with an assigned Instantaneous Unit Hydrograph (IUH), the maximum discharge was obtained with the well-known Chicago hyetograph [2].Moreover, by considering the so-called rational method [30], which is a worldwide procedure for flood design, a Rectangular Design Hyetograph (RDH) is usually adopted: the scenario with a duration d 1 and a rainfall intensity I 1 derived from an Intensity-Duration-Frequency (IDF) curve provides a discharge greater than that which is estimated by using a RDH with a duration d 2 (with d 2 > d 1 ) and I 2 derived from an IDF curve (with I 2 < I 1 ), and it is well-known that the total rainfall volume W 1 is less than W 2 [29,31].Consequently, the choice of focusing attention on annual maxima H and then selecting W comprising H appears more suitable in order to obtain the maximum discharge, with respect to considering annual maxima of W and then evaluating the associated peaks (which could not be equal to H).
The catchment response has been modelled using a lumped and event-based RR model [26,32] composed by the Soil Conservation Service Curve Number (SCS-CN) method [33] for rainfall excess estimation, and by Nash's model [34] for runoff estimation at basin outlet.Although the SCS-CN method is a popular (due its relative simplicity and its reliance on a limited number of parameters, resulting in a robust tool for those catchments that are partially or poorly gauged) and ubiquitous means for estimating storm runoff from a given rainfall event, i.e., as an alternative to empirical or physically-based infiltration models like Green-Ampt or Horton's equations, it should be noted that its usage at the sub-daily time scale is controversial [35].The major conflict is that the classic SCS-CN procedure allows for the estimation of cumulative rainfall excess (and thus cumulative losses) with rainstorms, but does not take time dimension into account in its equations.Several attempts have been made in order to overcome the conceptual limitations and inconsistencies of the original approach, for example, making explicit the concept of soil moisture accounting in a such a manner that is more suited for continuous simulations of rainfall-runoff processes [36,37] or proposing a mixed approach like the CN4GA, which aims to include the event scale correct information from the SCS-CN into the Green-Ampt infiltration equation [38].Nevertheless, the SCS-CN method is currently indicated by the River Basin Authority of Calabria region as a part of the procedure for design flood estimation within the study region and thus was used in the following sections of the paper.

Results and Discussion
Figure 2 shows the good comparison (on EV1 probabilistic plot) between the frequency distributions of observed samples and a generated series of spatially averaged values of 20-min annual maxima H. Figure 3 illustrates empirical ADF values and theoretical ADF curves, both obtained from the generated 500 years of synthetic data, for return periods T equal to 50, 100, and 200 years.In Figure 3, R d,T is the rainfall height relative to a duration d and a return period T. As regards the theoretical expression of ADF curves, a power formula was used [26], which is commonly adopted in Italy.

Results and Discussion
Figure 2 shows the good comparison (on EV1 probabilistic plot) between the frequency distributions of observed samples and a generated series of spatially averaged values of 20-min annual maxima H. Figure 3 illustrates empirical ADF values and theoretical ADF curves, both obtained from the generated 500 years of synthetic data, for return periods T equal to 50, 100, and 200 years.In Figure 3,   T  d R , is the rainfall height relative to a duration d and a return period T. As regards the theoretical expression of ADF curves, a power formula was used [26], which is commonly adopted in Italy.Focusing on the marginal CDFs F H (h) and F W n (w n ), the TCEV (Two Component Extreme Value, [39]) and the Lognormal distributions were tested.The fitting of the empirical CDFs, resulting from 500-year synthetic data, is reported in Figures 4 and 5    Concerning the choice of copula function, Gumbel expression (Equation ( 3)) provided the best results.From application of the CML method,  = 1.643 was obtained; Figure 6 highlights the good performance of Gumbel copula, in terms of Kendall function KC(-), on probabilistic and Q-Q plots.Moreover, Figures 7 and 8 show the scatter plots of empirical pairs (derived from 500-year synthetic data) with 10,000 pairs generated with the Monte Carlo method by the assumed bivariate distribution, in terms of marginal CDFs and random variables, respectively.Concerning the choice of copula function, Gumbel expression (Equation ( 3)) provided the best results.From application of the CML method,  = 1.643 was obtained; Figure 6 highlights the good performance of Gumbel copula, in terms of Kendall function KC(-), on probabilistic and Q-Q plots.Moreover, Figures 7 and 8 show the scatter plots of empirical pairs (derived from 500-year synthetic data) with 10,000 pairs generated with the Monte Carlo method by the assumed bivariate distribution, in terms of marginal CDFs and random variables, respectively.Concerning the choice of copula function, Gumbel expression (Equation ( 3)) provided the best results.From application of the CML method, ϑ = 1.643 was obtained; Figure 6 highlights the good performance of Gumbel copula, in terms of Kendall function K C (-), on probabilistic and Q-Q plots.Moreover, Figures 7 and 8 show the scatter plots of empirical pairs (derived from 500-year synthetic data) with 10,000 pairs generated with the Monte Carlo method by the assumed bivariate distribution, in terms of marginal CDFs and random variables, respectively.
) derived from 500-year synthetic data, and 10,000 pairs of CDF generated with the Monte Carlo method by the assumed bivariate distribution.
For the investigated basin, a statistical analysis of the sample annual maximum peak discharges provided the following quantiles for the return periods T = 50, 100, and 200 years: Q50 = 523 m 3 /s, Q100 = 629 m 3 /s, and Q200 = 735 m 3 /s.
As a final step of this study, the RR model, introduced in Section 4, was calibrated with the available time series of rainfall and discharge data according to a procedure placed in the context of Bayesian inference described in Biondi and De Luca [32].Specifically, the assessment of posterior distributions of model parameters were derived on the basis of the closeness of simulated hydrological signatures obtained with different parameter sets to those derived from actual at-site observations.The hydrological signatures considered to restrict hydrological model parameters are the first three L-moments of annual streamflow maxima; the best fitting set of parameters is considered in the following.
) derived from 500-year synthetic data, and 10,000 pairs of CDF generated with the Monte Carlo method by the assumed bivariate distribution.
For the investigated basin, a statistical analysis of the sample annual maximum peak discharges provided the following quantiles for the return periods T = 50, 100, and 200 years: Q50 = 523 m 3 /s, Q100 = 629 m 3 /s, and Q200 = 735 m 3 /s.
As a final step of this study, the RR model, introduced in Section 4, was calibrated with the available time series of rainfall and discharge data according to a procedure placed in the context of Bayesian inference described in Biondi and De Luca [32].Specifically, the assessment of posterior distributions of model parameters were derived on the basis of the closeness of simulated hydrological signatures obtained with different parameter sets to those derived from actual at-site observations.The hydrological signatures considered to restrict hydrological model parameters are the first three L-moments of annual streamflow maxima; the best fitting set of parameters is considered in the following.For the investigated basin, a statistical analysis of the sample annual maximum peak discharges provided the following quantiles for the return periods T = 50, 100, and 200 years: Q 50 = 523 m 3 /s, Q 100 = 629 m 3 /s, and Q 200 = 735 m 3 /s.
As a final step of this study, the RR model, introduced in Section 4, was calibrated with the available time series of rainfall and discharge data according to a procedure placed in the context of Bayesian inference described in Biondi and De Luca [32].Specifically, the assessment of posterior distributions of model parameters were derived on the basis of the closeness of simulated hydrological signatures obtained with different parameter sets to those derived from actual at-site observations.The hydrological signatures considered to restrict hydrological model parameters are the first three L-moments of annual streamflow maxima; the best fitting set of parameters is considered in the following.Through the RR model, it was possible to determine all the pairs   h w n , which provide the flood quantiles Q50, Q100, and Q200.For each QT, an envelope curve was obtained and compared with the contour lines associated to the bivariate return periods OR T = 50, 100, and 200 years, and with the pairs   h w n , related to the classical univariate approach (Figure 9).From the analysis of the obtained results, it is evident that with the hypothesis of Crit d equal to the time of concentration of the basin and with the adopted shape for the hyetograph: Comparison among pairs of (w n , h) derived from 500-year synthetic data, and 10,000 pairs of (w n , h) generated with the Monte Carlo method by the assumed bivariate distribution.
Through the RR model, it was possible to determine all the pairs (w n , h) which provide the flood quantiles Q 50 , Q 100 , and Q 200 .For each Q T , an envelope curve was obtained and compared with the contour lines associated to the bivariate return periods T OR = 50, 100, and 200 years, and with the pairs (w n , h) related to the classical univariate approach (Figure 9).Through the RR model, it was possible to determine all the pairs   h w n , which provide the flood quantiles Q50, Q100, and Q200.For each QT, an envelope curve was obtained and compared with the contour lines associated to the bivariate return periods OR T = 50, 100, and 200 years, and with the pairs   h w n , related to the classical univariate approach (Figure 9).From the analysis of the obtained results, it is evident that with the hypothesis of Crit d equal to the time of concentration of the basin and with the adopted shape for the hyetograph: From the analysis of the obtained results, it is evident that with the hypothesis of d Crit equal to the time of concentration of the basin and with the adopted shape for the hyetograph: • the design hyetograph from the univariate analysis (based on the ADF curves represented in Figure 3) should correspond to a return period greater than 200 years in order to get, by using the chosen RR model, a peak discharge equal to Q 50 ; • the bivariate analysis allows for obtaining Q 50 , Q 100 , and Q 200 with design hyetographs whose pairs (w n , h) can be associated to return periods T OR equal to the T of the peak discharges and, in any case, considerably less than those obtained from the univariate analysis; • on the basis of Equation ( 13), the adoption of other forms for JRP (T AND , T KEN , T COND ) would imply T values for the DH greater than those obtained with T OR .
The diagram in Figure 9 clearly constitutes a powerful tool that is aimed at determining the relationship between T-year DHs and the associated T-year design flood peaks for an assigned basin.Other shapes for the hyetograph (preserving the peak at d Crit /2) were also tested, but for the selected case study the obtained results were not dissimilar to the one represented in Figure 9.

Conclusions
The paper highlighted the importance of deriving DHs using multivariate analysis, which is in general more informative for design purposes than a univariate one when hydrological events involve correlated random variables.
In this topic, copula functions constitute an efficient tool, because they allow for marginal distributions to be studied separately (i.e., each distribution can assume a different mathematical expression) and then for their global dependence to be analyzed.
However, the choice of the most appropriate JRP depends on which event is critical for the investigated system (all the variables must exceed a certain magnitude, at least one variable must be greater than a threshold, and so on).Once the expression for JRP is fixed (T OR in this study), an ensemble of infinite critical combinations is related to a particular T value.For the selected basin, a very useful diagram was obtained as a result, which allows the relationship between the T-year DHs and the associated T-year design flood peaks to be evaluated for the adopted RR model.
The main outcomes of the study can be summarized as follows: as expected, the common univariate approach, based on the ADF curves, implies very high T-values for the DHs, compared with the return period of the peak discharges, thus confirming the poor consistency of the classical iso-frequency assumption between rainfall and peak flow; -the bivariate analysis conducted on the pairs of random variables (rainfall peak, net volume), while the critical duration is kept constant, allows for the desired T-year flood quantile to be obtained with at least one pair of the design hyetograph characteristics associated to the same return period, and, in any case, significantly less than those obtained from the univariate analysis.

Figure 1 .
Figure 1.The basin of Corace at Grascio: the green dots represent the rain gauges, while the red dot is the stream gauge.

Figure 1 .
Figure 1.The basin of Corace at Grascio: the green dots represent the rain gauges, while the red dot is the stream gauge.

Figure 2 .
Figure 2. EV1 probability plot of 20-min annual maxima H at basin scale: comparison between observed (estimated by using the technique of Voronoi polygons) and generated series for Corace at Grascio.

Figure 2 .
Figure 2. EV1 probability plot of 20-min annual maxima H at basin scale: comparison between observed (estimated by using the technique of Voronoi polygons) and generated series for Corace at Grascio.

Figure 2 .
Figure 2. EV1 probability plot of 20-min annual maxima H at basin scale: comparison between observed (estimated by using the technique of Voronoi polygons) and generated series for Corace at Grascio.

Figure 4 .
Figure 4. EV1 probabilistic plot: comparison between empirical Cumulative Density Function (CDF) and Two Component Extreme Value (TCEV) distribution for both H (left) and n W (right).

Figure 5 .
Figure 5. Lognormal probabilistic plot: comparison between empirical CDF and Lognormal distribution for both H (left) and n W (right).

Figure 4 .
Figure 4. EV1 probabilistic plot: comparison between empirical Cumulative Density Function (CDF) and Two Component Extreme Value (TCEV) distribution for both H (left) and W n (right).

Figure 4 .
Figure 4. EV1 probabilistic plot: comparison between empirical Cumulative Density Function (CDF) and Two Component Extreme Value (TCEV) distribution for both H (left) and n W (right).

Figure 5 .
Figure 5. Lognormal probabilistic plot: comparison between empirical CDF and Lognormal distribution for both H (left) and n W (right).

Figure 5 .
Figure 5. Lognormal probabilistic plot: comparison between empirical CDF and Lognormal distribution for both H (left) and W n (right).

Figure 6 .
Figure 6.Performance of Gumbel copula, in terms of Kendall function KC(-) (which represents the CDF of copula function), on probabilistic plot (left) and Q-Q plot (right).

Figure 6 . 13 Figure 6 .
Figure 6.Performance of Gumbel copula, in terms of Kendall function K C (-) (which represents the CDF of copula function), on probabilistic plot (left) and Q-Q plot (right).

Figure 7 .
Figure 7.Comparison among pairs of CDF (F W n (w n ), F H (h)) derived from 500-year synthetic data, and 10,000 pairs of CDF generated with the Monte Carlo method by the assumed bivariate distribution.

Figure 9 .
Figure 9.Comparison between univariate (square boxes) and bivariate analysis (solid lines) for the definition of the design hyetograph.The dotted lines correspond to all the pairs   h w n , providing quantiles Q50, Q100, and Q200 of discharge when fed into the adopted rainfall-runoff (RR) model.

Figure 8 .
Figure 8.Comparison among pairs of (w n , h) derived from 500-year synthetic data, and 10,000 pairs of (w n , h) generated with the Monte Carlo method by the assumed bivariate distribution.

Figure 9 .
Figure 9.Comparison between univariate (square boxes) and bivariate analysis (solid lines) for the definition of the design hyetograph.The dotted lines correspond to all the pairs   h w n , providing quantiles Q50, Q100, and Q200 of discharge when fed into the adopted rainfall-runoff (RR) model.

Figure 9 .
Figure9.Comparison between univariate (square boxes) and bivariate analysis (solid lines) for the definition of the design hyetograph.The dotted lines correspond to all the pairs (w n , h) providing quantiles Q 50 , Q 100 , and Q 200 of discharge when fed into the adopted rainfall-runoff (RR) model.

Table 1 .
Station network referred to Corace at Grascio.

Table 1 .
Station network referred to Corace at Grascio.
on EV1 and Lognormal probabilistic plot, respectively.It is clear the TCEV reproduces better than Lognormal distribution the empirical CDF for both H and W n , and thus it was chosen for bothF H (h) and F W n (w n ). ) and the Lognormal distributions were tested.The fitting of the empirical CDFs, resulting from 500-year synthetic data, is reported in Figures4 and 5on EV1 and Lognormal probabilistic plot, respectively.It is clear the TCEV reproduces better than Lognormal distribution the empirical CDF for both H and n W , and thus it was chosen for both   ) and the Lognormal distributions were tested.The fitting of the empirical CDFs, resulting from 500-year synthetic data, is reported in Figures4 and 5on EV1 and Lognormal probabilistic plot, respectively.It is clear the TCEV reproduces better than Lognormal distribution the empirical CDF for both H and n W , and thus it was chosen for both  