Derivation of flood frequency curves through a bivariate rainfall distribution based on copula functions: application to an urban catchment in northern Italy's climate

Although the design event method is still a very popular approach for sizing hydraulic devices and for evaluating drainage system performances, its limitations have long been highlighted and debated (Adams & Howard 1986; Yen 1990). The main reason lies in assuming that the frequency of occurrence of the total depth of the design rainfall hyetograph equals that of the derived flood hydrograph. It is well known that the flood hydrograph depends on a number of factors, which are not accounted for by the rainfall statistics on which design event methods are based. Therefore, crucial characteristics, such as the catchment initial condition, the wet weather duration and the rainfall time pattern, are arbitrarily set. As a result, estimates derived by design event methods are always biased (Bacchi et al. 1993).

In order to overcome such deficiencies, continuous approaches should be preferred. In particular, the analytical-probabilistic methods, originally developed by Eagleson (1972), recently aroused researcher interest (Guo et al. 2012). In fact, these methods feature all the advantages of continuous approaches, even if they yield closed-form analytical expressions analogous to those of the design event methods.

According to the derived distribution theory, complete cumulative distribution functions (CDFs) of the runoff variables of interest can be defined from joint distribution functions (JDFs) of the rainfall variables (Adams & Papa 2000). Nevertheless, to have the possibility of analytically integrating the derived distributions, simplifying hypotheses must be adopted in the rainfall probabilistic model: usually exponential CDFs are utilized for the marginals and the mutual variable dependency is disregarded (Guo & Adams 1998; Balistrocchi et al. 2013).

When dealing with rainfall depth and wet weather duration JDFs, more complex CDFs are often needed to properly suit empirical marginal distributions, while the independence hypothesis is largely rejectable. Surprisingly, simplified JFDs produce results in better agreement with benchmarks, and more conservative, than those achieved by the more realistic ones (Adams & Papa 2000). Even if not totally explained, the reasons can lie in a compensation of errors in the simplified models and in drawbacks related to the JDF assessment by means of traditional techniques in the complete models.

Fortunately, copula functions (Nelsen 2006), recently introduced in hydrologic research (De Michele & Salvadori 2003; Favre et al. 2004), offer the opportunity to broaden the multivariate inference capability. Indeed, through copula approach, the assessment of the dependence structure relating to random variables can be carried out separately from those of marginal distributions (Genest & Favre 2007; Salvadori et al. 2007). As a consequence, the dependence structure analysis is no longer affected by marginal behaviours. Further, copula functions and marginal distributions belonging to different probability families, even complex, can be utilized to develop the JDF. The derivation procedure therefore becomes straightforward, eliminating additional sources of uncertainty, such as preliminary sample data transformations, and permitting an effective verification of model reliability by means of blanket test statistics (Genest et al. 2009).

Owing to these significant advances, the longstanding debate regarding the generalization of the return period concept, from the univariate case to the multivariate one, has been notably reopened. Traditionally, two methods based on simple logical expressions ‘AND’ ‘OR’, already utilized in conventional multivariate statistics but susceptible to copula interpretation, were utilized. More recently, Salvadori & De Michele (2010) proposed a new estimation method relying on Kendall function.

However, as can be seen in Gräler et al. (2013), where a review of approaches to develop multivariate design events is provided, all these methods lead to statistically different outcomes; moreover, a generally applicable solution cannot be suggested, so that the most suitable method should be selected in consideration of the specific hydraulic application at hand. Hydrologic event severity is actually related to hydraulic device performances: on the one hand, all events yielding identical performances can be associated with the same return period; on the other hand, the same event can be associated with different return periods depending on the type of device.

According to this device-targeted approach, the return period of a derived variable, exploited to express device performances, can be estimated by conventional univariate techniques and thus be associated with the input multivariate event. For instance, the return period of bivariate flood events, defined by using peak discharge and runoff volumes, was estimated with reference to routing reservoir performances by Requena et al. (2013) and Volpi & Fiori (2014), by using univariate statistics of maximum water level of the stored volume and maximum routed discharge, respectively.

Bearing in mind these encouraging results, herein the derivation of a flood frequency curve (FFC) from a JDF of storm volume and wet weather duration is developed for urban applications in medium-size catchments (area less than 100:200 ha) located in northern Italy. To do so, a hydrologic loss model and a runoff routing model must be established. As demonstrated by Candela et al. (2014) dealing with extended natural watersheds, a copula-based rainfall stochastic generator can be successfully coupled with a fully conceptual distributed model. Nevertheless, a lumped modelling based on hydrologic parameters operatively adopted in practical engineering was preferred in this study. Owing to rainfall–runoff transformation dynamics in urban catchments, this type of modelling was considered to be sufficiently reliable and more appealing for practitioners.

The objective of this paper is two-fold: first, quantifying the modelling improvement achievable by a complete JDF with respect to a simplified analytical-probabilistic model, in a case where the speculated compensation of errors does not occur; second, developing a device-targeted criterion to estimate the return period of bivariate rainfall events from univariate statistics of the runoff peak discharge.

Hereafter, copula analysis technique and JDF construction are briefly recalled in the first section, while special attention is given to the return period estimate in the second one, focusing on the existing method drawbacks and on the development of the proposed method. The continuous simulation model employed to derive the benchmarking FFC, available data and the test watershed are described in the following sections. JDF fitting results and the FFC derivation are then discussed. Conclusions are drawn in the last section.

In the following sub-sections, the three phases outlined above – (i) discretization procedure, (ii) underlying copula assessment and (iii) marginal distribution assessments – are discussed in more detail, justifying the selection of theoretical functions with reference to existing literature.

The separation of a continuous precipitation record into independent events can be carried out by applying two discretization thresholds: an interevent time definition (IETD) and a rainfall depth threshold (Balistrocchi et al. 2009). The first one represents the minimum dry weather period needed for two following hyetographs to be considered independent. Hence, if two rain bursts are detached by a dry weather period shorter than IETD, they are aggregated into a unique event, whose duration and depth are computed from the beginning of the first one to the end of the latter one.

The second parameter corresponds to the minimum rainfall depth that must be exceeded in order to have a rainfall relevant to analysis purposes. When this condition is not satisfied, the event is suppressed and the corresponding wet weather duration is assumed to be rainless and joined to the adjacent dry weather periods.

Discretization parameters strongly affect the derived sample statistics, namely the mean annual number of rainfall events and marginal parameters. Then, their values must be chosen very carefully. When dealing with FFC derivation, a suitable criterion is focusing on runoff discharge characteristics: the IETD can be estimated as the minimum time that avoids the hydrographs generated by two subsequent rainfalls to overlap, while the depth threshold can be identified with the initial abstraction (IA) of the catchment hydrological losses.

Exponents β and γ are dimensionless parameters that rule the distribution shape, while the denominators ζ (mm) and λ (h) are scale parameters. The lower limit IA in CDF (15) is set in accordance with the volume threshold employed in the rainfall discretization procedure. In fact, the shape parameter makes the distribution very versatile and leads to fittings as satisfactory as those achievable by using more complex models. The exponent less than one delineates a probability density function monotonically decreasing from the lower limit, where a vertical asymptote is present. The exponent greater than one shows that the probability density function has a finite mode and exhibits a right tail.

In order to quantify the impact on the whole model reliability of assuming exponential marginals, as in the most common analytical-probabilistic models, JDF with exponential marginals were fitted as well. Equations (15) and (16) consequently simplify, being the shape parameters β and γ unitary, and feature a finite mode in the lower limit, equal to the reciprocal of the scale parameter.

Despite the fundamental progress in understanding and overcoming the inherent limitations of the previous ones, estimate method demonstrates that it is affected by conceptual drawbacks and yields unacceptable estimates when applied to FFC derivation from the constructed JDF.

The transformation function relating to the input variables and the derived variable, called by some authors structure function (Volpi & Fiori 2014; Salvadori et al. 2015; Pappadà et al. 2016), must be a subjective function such that a sub-set of the input variable population is related to a single occurrence of the derived variable. This sub-set separates the sub-critical region from the super-critical one, by discriminating input events leading to lower or greater values of the derived variable, respectively. It is apparent that these regions are dichotomous and that infinite input events are associated with a unique return period.

Following this criterion, for instance, in Requena et al. (2013), copula simulation techniques were exploited to generate from a bivariate distribution of peak discharge and runoff volume a number of flood events, forcing a real-world routing reservoir. Hydraulic simulations were utilized as transformation function to obtain maximum water levels occurring in the storage volume during the routing process. Thus, flood events were associated with the return period of corresponding maximum water levels, estimated by means of their empirical frequencies.

Conversely, Volpi & Fiori (2014) analysed the same problem from a theoretical point of view, by assuming a constant inflow discharge and a linear behaviour for the reservoir. Under such simplifying hypotheses, an analytically closed-form transformation function relating inflow peak discharge and runoff volume to routed peak discharge can be derived. Such a function was exploited to delimitate sub-critical regions in the population of inflow flood variables for constant values of routed peak discharge. The bivariate return period estimate was therefore carried out with reference to non-exceedance probabilities, obtained by integrating copula density functions suggested in the literature on such regions.

Differently, peak discharge must be accounted for when dealing with open channel design or safety verification. This variable depends on the rainfall volume and the wet weather duration and can be derived from their JDF. Further, all rainfall events leading to an identical peak discharge can be associated with the same severity, or return period.

The integral in (30) must be computed numerically, unless the JDF is constructed by coupling exponential functions through the independence copula and Equation (27) is simplified by assuming the wet weather duration to be equal to the rainfall excess duration, as in Balistrocchi et al. (2013). Numerical computing is, however, facilitated by operating in the unitary square, which allows only proper integrals to be used.

Hydrological simulations were conducted to obtain the continuous discharge series generated by the test catchment, as a result of the observed rainfalls. A conceptual lumped model was developed, aiming at limiting discrepancies with respect to the one implemented in the derivation of the expression. Hence, the rainfall–discharge transformation process was again represented by the loss model (25), while the rainfall excess routing was executed by using a triangular instantaneous unitary hydrograph, with duration equal to t_c.

The JDF P_HD is herein constructed with reference to a 45-year long time series of observed rainfall depths, recorded every 30 minutes from 1949 to 1993 by the ITAS Pastori raingauge, located in Brescia, northern Italy, next to the transition between the Padan Plain and the southern foothills of the Alps. The rainfall regime is classified as sub-alpine, characterized by two maxima, in spring and autumn, and two minima, in winter and summer. During summer, short duration and high intensity storms separated by long dry weather periods occur, while longer but less intense precipitations are more common in the other seasons. The annual average rainfall depth is quite large with respect to other Italian climates and amounts to about 1.000 mm.

Bearing in mind small–medium size urban catchment applications, analyses referred to a synthetic test watershed characterized by this parameter set: catchment area A 100 ha, time of concentration t_c 20 minutes and hydrological losses expressed by an initial abstraction IA 5 mm and a runoff coefficient Φ 0.45, for the exceeding portion. The discretization of the continuous rainfall time series into independent rainfall events was hence conducted by using an IETD corresponding to the minimum suggested value of 3 h (Adams & Papa 2000) and a volume threshold equal to 5 mm.

The application of the discretization procedure to Brescia's rainfall time series yields a mean annual number ω of 48.2 independent rainfall events. This result can be considered reasonable in this rainfall regime, owing to the very small IETD adopted for the test catchment.

The maximum likelihood criterion was employed to fit the theoretical copula (8) to pseudo-observations, estimating a value of 1.41 for the dependence parameter θ. Although this value corresponds to a concordance coefficient τ_K equal to 0.29, a significant association is however shown. Indeed, when a goodness-of-fit test is performed according to statistics (12) by assuming that the dependence structure is given by the independence copula, that is, disregarding the association degree, the null hypothesis can be rejected for a p value (13) less than 0.1%.

Conversely, the null hypothesis that the underlying copula is expressed by copula (8) cannot be rejected for a p value greater than 80.0%. Moreover, the upper tail dependence coefficient λ_U (10) corresponding to θ estimate is 0.366, whereas the non-parametric estimate (14) yields a value of 0.346.

Marginals (15) and (16) were finally fitted to h and d samples by the maximum likelihood criterion, as well, leading to the following parameter set: 0.88, 1.18, 10.9 mm and 9.5 h for β, γ, ζ and λ, respectively. A β value less than one evidences a significant presence of rainfall events having very small depth, while a γ value greater than one shows that wet weather durations tend to be less asymmetric. The goodness-of-fit was verified by means of the confidence boundary test, revealing that these distributions cannot be rejected for conventional levels of significance of 5–10%. When marginals are simplified in exponential distributions, shape parameters are set equal to one, while scale parameter estimates are: 12.0 mm for ζ and 8.9 h for λ. However, such distributions can be rejected for a significance less than 1.0%.

As can be seen, a good agreement between CS FFC and GW FFC is evidenced. The largest residual, about 23%, is given by the peak estimate corresponding to the 45-year return period, that is the series length, and it can be considered acceptable. Instead, in the range of 5:20 years T_r, residuals are much smaller and the probabilistic model slightly overestimates CS FFC. This can be justified by considering that the rainfall series time step is comparable to t_c, so that simulated peak discharges are expected to be more attenuated than the real ones. However, GW FFC appears to be a little more conservative than benchmarking outcomes and its values match those suggested in the analysed rainfall regime for urban sewer design (a term of reference adopted by practitioners is about 100 l/(s ha) for T_r 10 years). This result definitively supports both the reliability of the JDF based on such functions and the device-targeted method to estimating the return period of bivariate rainfall events (30) herein developed.

On the contrary, in the other couplings, unacceptably large estimates are obtained. The worse result is provided by IE FFC, since almost double peak discharges than the benchmarking ones are derived. Hence, the most common assumptions on which analytical-probabilistic models commonly rely, in this situation, do not yield the compensation of errors occurring elsewhere. The most significant reliability improvement is due to a correct representation of the dependence structure rather than a better fitting of marginals. In fact, GE FFC is much closer to the CS FFC than IW FFC. This occurrence may be explained by considering that, even if rejected by test statistics, exponential distributions approximately suit the empirical marginal distributions, while the observed dependence structure is substantially different from the independence copula.

First, it is well known that, for a given T_r, the rainfall volume is expected to increase with respect to the wet weather duration. This behaviour is respected by isolines, which also provide a rainfall volume greater than zero for null durations. These trends are a consequence of the super-critical region delimitated by the method, which is mainly located in the lower-right corner of the unitary square. As can be seen in the example in Figure 1(d), in this approach, super-critical events are those featured by large depth but short duration, which exceed the hydrologic losses; the conceptual soundness of the proposed model is therefore further supported.

The extremely large peak discharge obtained from JDFs based on the independence copula can be explained, as well. In the independent copula, the copula density is uniformly distributed on ⁠, while, in the Gumbel–Hougaard copula, it concentrates on the main diagonal in the lower-left and in the upper-right corners. Thus, in the first case, the non-exceedance probability is underestimated, consequently leading to an excessively low return period.

When return periods are estimated by using and methods very similar outcomes are obtained. Isolines exhibit a decreasing trend of the rainfall volume with respect to the wet weather duration. In addition, since return period isolines correspond to JDF contour lines, a significant range of short durations next to the axis origin is excluded: the greater the return period, the larger this range is. On the contrary, the wet weather duration range must span from zero to infinite. Furthermore, in both methods, estimated return periods are evidently meaningless. More reasonable values are instead assessed by means of the method, although the rainfall volumes appear to be almost independent of the wet weather duration. This occurs because, in the region taken into consideration, survival copula contour lines are almost straight and parallel to wet weather duration axis.

Such huge discrepancies between and the others can be explained by the substantially different splitting of the bivariate population, by which they estimate the non-exceedance probability. The extreme sensitivity of the return period with respect to this estimate further accentuates such discrepancies. In complete agreement with Serinaldi (2015), ⁠, and methods can be adopted only if they match the failure mechanism of the analysed device.

For example, fewer discrepancies were found by Requena et al. (2013), when the estimate method based on routing reservoir performances is compared to and ⁠. In this application, peak discharge must actually decrease with the runoff volume when the maximum water level, or the routed peak discharge, is constant. Therefore, the failure mechanisms on which and rely are more similar to the device one. Conversely, as in the problem dealt with in this study, they do not provide suitable return period estimates.

In this paper, a FFC for a small–medium size urban catchment in a northern Italian rainfall regime has been derived, from a JDF of rainfall depths and wet weather durations. To do so, a Gumbel–Hougaard copula and Weibull marginals were combined and a bivariate return period estimate method was developed, by using a lumped conceptual scheme for the rainfall–runoff transformation and routing processes. The proposed model soundness is supported by a satisfactory agreement with benchmarking continuous simulations. Two main consequences can be drawn from this.

First, although a method of general applicability to estimating multivariate return periods cannot reasonably be established, it is however possible to delineate a common strategy, addressing the real-world application dynamics. In this circumstance, if device performances can be referred to a single constituent variable, the derived distribution theory allows the multivariate case to be traced back to the univariate one. In fact, with regard to the derived hydrologic variable, the population of the input random variable can effectively be split in the dichotomous regions, collecting sub-critical and super-critical events. Furthermore, it is important to point out that, in a multivariate framework, defining a single change in event for a given return period appears to be meaningless. When an estimate method is chosen, a single event corresponds to a unique return period but, on the contrary, infinite events share the same return period.

Second, the utilization of the copula approach in the rainfall JDF construction demonstrated to yield a substantial improvement in the overall model reliability. Differently from what is reported in the literature concerning the analytical-probabilistic modelling, JDFs properly accounting for the observed dependence structure of random variables, featuring a non-negligible association, performed better than those assuming their independence. Nonetheless, in the latter case, more conservative results are obtained. In any case, the model choice should be carried out carefully, balancing the definitive advantage of an analytical formulation and the real need for accuracy improvement following a more precise probabilistic model.

Future developments regarding the method herein proposed to estimating the multivariate return period are still desirable. On the one hand, the FFC derivation could be generalized towards natural watershed applications, whose scale and hydrological processes require more complex modelling. On the other hand, it could be attractive to more deeply understand, with regard to real-world test cases, whether the method can be suitably exploited to face other multivariate problems, or not. Indeed, a number of hydraulic design and verification problems, in which the hydrologic input would be more efficaciously outlined by a multivariate continuous approach, exists, namely landslide triggering risk assessment, routing reservoir design, spillway safety verification and flooding area delimitations.