Stoch Environ Res Risk Assess DOI 10.1007/s00477-013-0813-z

ORIGINAL PAPER

On peaks-over-threshold modeling of floods with zero-inflated Poisson arrivals under stationarity and nonstationarity Artur Tiago Silva • Maria Manuela Portela Mauro Naghettini

•

Ó Springer-Verlag Berlin Heidelberg 2013

Abstract The peaks-over-threshold (POT) model with Poisson arrivals and generalized Pareto (GP) distributed exceedances remains a popular and useful tool for modelling hydrologic extremes. The use of the Poisson–GP model for flood frequency analysis requires the validation of the hypothesis that the distribution of the annual number of flood events may be described by a Poisson distribution. Such hypothesis is not always valid in practical applications. The present study concerns the use of an alternative distribution for modelling the annual number of floods— the zero-inflated Poisson (ZIP) distribution with two parameters. A ZIP–GP model for flood frequency analysis is proposed. This model is less restrictive than the Poisson– GP model since it allows for a more accurate description of the occurrence process in a POT framework if the fraction of years with no exceedances is significantly higher than the theoretical mass at zero of the Poisson distribution. Furthermore, a nonstationary model (NSZIP–GP) is presented, in which the parameters of the ZIP are allowed to change in time as a function of a covariate, which, even for stationary peak magnitudes, affects the annual maximum flood quantiles with a given non-exceedance probability. Applications of the ZIP–GP model to flood data from Northern Portugal and the evaluation of its performance relative to the Poisson–GP model, including assessments of quantile uncertainty, are presented. An illustrative application of the NSZIP–GP model, using the North Atlantic Oscillation as a covariate is also presented. The

A. T. Silva (&)  M. M. Portela CEHIDRO, Instituto Superior Te´cnico, Lisbon, Portugal e-mail: [email protected] M. Naghettini Federal University of Minas Gerais, Belo Horizonte, Brazil

applications of both models include assessment of quantile uncertainty. Keywords Flood frequency analysis  Peaks-overthreshold  Zero-inflated Poisson  Nonstationarity

1 Introduction The peaks-over-threshold (POT) method is an established sampling technique in flood frequency analysis. It differs from the more commonly used annual maximum series (AMS) method, since it samples all the peak values above a given threshold, thus allowing for a more comprehensive and representative selection of extreme values. The peaks that constitute a POT sample may be modeled by a dualdomain-model that is able to describe both the magnitudes and the times of occurrence of the selected events. The aim of the method is to derive the distribution of annual floods from a combination of the distribution of the flood peaks with the distribution of the annual number of floods (Lang et al. 1999). The generalized Pareto (GP) distribution is commonly adopted to model the exceedance magnitudes in a POT framework (Davison and Smith 1990; Rosbjerg et al. 1992; Madsen et al. 1997; Lang et al. 1999; Silva et al. 2012a), although alternative distributions have been proposed in the technical literature, such as the gamma (Zelenhasic 1970), the Weibull (Ekanayake and Cruise 1993), the lognormal (Rosbjerg and Rasmussen 1991), and the generalized logistic (Bhunya et al. 2012) distributions. Justification for the use of the GP is given in extreme value theory (EVT): the second theorem of EVT, i.e., the Pickands–Balkema–de Haan theorem (Balkema and de Haan 1974; Pickands 1975) defines that the asymptotic tail

123

Stoch Environ Res Risk Assess

distribution of a random variable X is of the generalized Pareto family of distributions. The GP has the cumulative distribution function (CDF) (  xlj1 j 6¼ 0; x [ l; r ; GðxÞ ¼ 1  1  j xl ð1Þ 1  exp  r ; j ¼ 0; x [ l; where j is the shape parameter, r is the scale parameter, and l is the location parameter equal to the threshold. The GP reduces to the Exponential distribution as a special case (j = 0). Shane and Lynn (1964) proposed that the arrival of peaks may be described by a homogeneous Poisson process with intensity k (given by the mean annual number of threshold exceedances), i.e., that the annual number of peaks follows a Poisson distribution. Under this assumption, the CDF of annual maximum floods, F(x) (Davison and Smith 1990) is given by FðxÞ ¼ expfk½1  GðxÞg:

ð2Þ

By combining (1) and (2), one obtains n  (  j1 o exp k 1  j xl ; j 6¼ 0; x [ l; FðxÞ ¼   xlr  j ¼ 0; x [ l; exp k exp  r ; ð3Þ which is the generalized extreme value (GEV) distribution with a Poisson–GP parametrization. Equation (3) reduces to the extreme value type-1 (EV1) or Gumbel distribution for j = 0. This correspondence is in agreement with the first theorem of extreme value theory, i.e. the Fisher– Tippett–Gnedenko theorem (Fisher and Tippett 1928; Gnedenko 1943), according to which the limiting distribution of the annual maxima of a random variable is one of the three upper-tail types of a GEV. Given the above, the Poisson–GP setup for POT modeling remains a popular and useful alternative to common statistical analysis of AMS. However, that setup requires the validation of the Poisson assumption for the overthreshold arrivals. Cunnane (1979) and NERC (1975, pp. 197–198) propose a test for the verification of such assumption based on the dispersion index (ratio between variance and mean) of the observed annual number of over-threshold exceedances. Those studies include an application of the dispersion index test to streamflow data from 26 gauging stations in Great Britain and show that, generally, the annual number of peaks in those samples did not conform to the Poisson assumption, and that the departures from it could be explained by variances significantly higher than the means. Cunnane (1979) and NERC (1975, pp. 199–201) proposed an alternate discrete distribution— the negative binomial (NB)—which allows for higher variances, but conclude their study with the perception that

123

the use of the NB did not seem to offer any definite improvement. Ben-Zvi (1991) analysed the goodness of fit of the Poisson and NB distributions to data from eight rivers in Israel and found that the acceptance rate of the NB was larger than that of the Poisson. ¨ no¨z and Bayazit (2001) compared flood estimates O based on the Poisson, binomial and NB distributions for count data, combined with the Exponential distribution for the peak magnitudes and concluded that such estimates, as well as their sampling variance, were nearly identical. In this paper an alternative to the Poisson, binomial, and NB distributions is proposed for modeling the annual overthreshold exceedance arrival counts—the zero-inflated Poisson (ZIP) distribution. According to Lambert (1992), this distribution was introduced by Singh (1963). The ZIP distribution can be combined with the GP to result in a new ZIP–GP model for flood frequency analysis. The proposed ZIP–GP model is based on the stationarity assumption—a basic assumption that historically has assisted practice and research in the fields of hydrology— i.e., its parameters are invariant in time and so the annual maximum flood quantile with a given non-exceedance distribution is also invariant in time. Recently, however, there seems to be a consensus among the scientific community that, due to climate change, there is an intensification of the hydrological cycle, assumably leading to more frequent and more intensive extreme hydrological events, like floods and droughts (Katz 2013). Milly et al. (2008) argue that such climate change undermines stationarity. Clarke (2007) questions the widespread assumption of stationarity in hydrological practices and argues that the next few decades should see an increase in understanding the processes causing climate changes and variability, not only for the purpose of forecasting the development of such changes, but also for predicting the frequency of occurrence of events of certain magnitudes. On an alternate viewpoint, Koutsoyiannis (2011) and Lins and Cohn (2011) contend that hydrologic change can coexist with stationarity through long-term persistence, and encourage the adoption of a Hurst-Kolmogorov analysis framework to explicitly consider such changes while maintaining the stationarity assumption. There are a number of studies linking changes in hydrological extremes to changes in teleconnection indices such as the North Atlantic Oscillation (NAO) and the El Nin˜o-Southern Oscillation (ENSO) (Trigo et al. 2005; Lee and Ouarda 2010; Grimm and Tedeschi 2009; Delgado et al. 2012; Silva et al. 2012b; Lo´pez and France´s 2013). In this context, if hydrological data that support the planning and management of water resources systems show significant signs of change, the mathematical formalism of the statistical models commonly used for hydrological

Stoch Environ Res Risk Assess

frequency analysis should be reviewed in order to account for the dynamic behaviour of the data. In recent years, studies on the nonstationarity of hydrometeorological extremes have become a topic of high scientific interest (Renard et al. 2006; Hundecha et al. 2008; Villarini et al. 2009, 2011; Katz 2010, 2013; Delgado et al. 2010; Siliverstovs et al. 2010; Galiatsatou and Prinos 2011; Kwon et al. 2011; Salas and Obeysekera 2013; Tramblay et al. 2013). Given the aforementioned emergence of nonstationary assessment of the frequency of hydrological extremes, a nonstationary zero-inflated Poisson–generalized Pareto (NSZIP–GP) model for flood frequency analysis is proposed, in which the two parameters of the ZIP are allowed to vary in time as a function of a covariate, by means of a ZIP regression (Lambert 1992). Under a POT approach, changes in the underlying distribution of peak occurrences alter the annual maximum flood quantile, even for stationary peak exceedances.

in a year. If the peak magnitudes XOT are assumed to be independent of their times of occurrence, (6) becomes FðxÞ ¼ PðY ¼ 0Þ þ

1 X

fPðY ¼ yÞ½GðxÞy g;

ð7Þ

y¼1

where G(x) is the PDF of the magnitude of the overthreshold peaks. Under the stationarity assumption, the parameters / and k are considered constant in time. Moreover, the probability distribution of the exceedances over a threshold u, G(x), is also taken as invariable in time. After algebraic expansion outlined in Appendix 2, (7) becomes FðxÞ ¼ / þ ð1  /Þ expfk½1  GðxÞg

ð8Þ

and if we take G(x) from (1), (8) becomes FðxÞ n  8   1 o < / þ ð1  /Þ exp k 1  j xl j ; r ¼   xl : / þ ð1  /Þ exp k exp  r ;

j 6¼ 0; x[ l; j ¼ 0; x[ l; ð9Þ

2 Model formulation 2.1 The stationary model, ZIP–GP The simplest distribution for modeling count data, such as the annual number of over-threshold exceedances, is the Poisson distribution with parameter k and probability mass function (PMF) given by y

PðY ¼ yÞ ¼ expðkÞ

k ; y!

ð4Þ

where k is estimated by the sample mean. The ZIP distribution (Lambert 1992) is a two-component mixture model combining a point mass at zero with the Poisson distribution (4), resulting in the following PMF ( PðY ¼ 0Þ ¼ / þ ð1  /Þ expðkÞ; ð5Þ y PðY ¼ yÞ ¼ ð1  /Þ expðkÞ ky! ; y 2 N1 ; where the parameter / is the probability inflation at y = 0. The two parameters of the ZIP distribution, k an /, can be estimated via maximum likelihood (see Appendix 1). Under a POT sampling scheme of streamflow data, if the peaks above a threshold are i.i.d., the CDF of the annual maxima, XAM, can be formally expressed as FðxÞ ¼ Pð XAM  xÞ ¼ PðY ¼ 0Þ þ

1 X y¼1

" P

y \

# \ ð XOTk  xÞ ðY ¼ yÞ ;

ð6Þ

k¼1

where the random variable Y is the annual arrival counts and XOTk is the magnitude of the kth over-threshold peak

and the corresponding quantile functions, xF, 8 n h  ij o < l þ r 1   1 log F/ ; j¼ 6 0; x [ l; j k h  1/ i xF ¼ F/ : l  r log  1 log ; j ¼ 0; x [ l: 1/ k ð10Þ Equations (9) and (10) define the ZIP–GP model for flood frequency analysis. 2.2 The nonstationary model, NSZIP–GP Lambert (1992) introduced the zero-inflated Poisson regression model, which enables the estimation (via maximum likelihood) of the parameters of the ZIP as a function of one or more covariates. Having the parameters of the ZIP vary in time as a function of a single covariate, z, with an annual time step, enables the introduction of nonstationarity in the flood frequency model outlined in Sect. 2.1, leading to the following frequency model for annual maximum floods Fðx; zÞ ¼ /ðzÞ þ ½1  /ðzÞ expfkðzÞ½1  GðxÞg;

ð11Þ

where k(z) is the log link function log kðzÞ ¼ k0 þ k1 z , kðzÞ ¼ expðk0 þ k1 zÞ; and /(z) is the logit link function logit /ðzÞ ¼ /0 þ /1 z , /ðzÞ ¼

expð/0 þ /1 zÞ : 1 þ expð/0 þ /1 zÞ

123

Stoch Environ Res Risk Assess

The forementioned link functions are suggested by Lambert (1992), as the log k(z) is the canonical link that linearizes Poisson means [ensures that k(z) [ 0], and the logit /(z) is the canonical link that linearizes Bernoulli probabilities of success (ensures that /ðzÞ 2 ð0; 1Þ). The Nonstationary zero-inflated Poisson–generalized Pareto, NSZIP–GP, model for floods is given by

( Fðx; zÞ ¼

n   j1 o /ðzÞ þ ½1  /ðzÞ exp kðzÞ 1  j xl ;   xlr  /ðzÞ þ ½1  /ðzÞ exp kðzÞ exp  r ;

and the corresponding quantile functions, xF,z, xF;z 8 n h  ij o 1 > log F/ðzÞ ; < l þ rj 1   kðzÞ 1/ðzÞ ¼ h   i > : l  r log  1 log F/ðzÞ ; kðzÞ 1/ðzÞ

j 6¼ 0; x [ l; j ¼ 0; x [ l: ð13Þ

It should be mentioned that nonstationarity could be introduced not only in the count data but also in the parameters of G(x) (Coles 2001; Katz 2013). 2.3 Quantile uncertainty To fully evaluate the performance of the ZIP–GP and NSZIP–GP models, confidence intervals were constructed for evaluating the accuracy of the quantile estimation. Asymptotic theory of sampling distributions shows that the quantiles XF are approximately normal with mean x^F and variance VF (Coles 2001, pp. 31–33). This result enables the construction of approximate 100(1 - a)% quantile confidence intervals. The asymptotic quantile variance VF of a k-parameter (h ¼ ðh1 ; . . .; hk Þ) distribution model can be approximated by VF ¼ h> Rh;

ð14Þ

where R is the limiting variance–covariance matrix of the parameter estimators and h is the gradient of the parameter estimators, obtained by numerical differentiation,  oxF  h¼ ; h ¼ ðk; /; r; jÞ; ohi F;h^ for the ZIP–GP model, and  oxF;z  h¼ ; h ¼ ðk0 ; k1 ; /0 ; /1 ; r; jÞ; ohi F;z;h^

123

for the NSZIP–GP model as, in this case, the quantiles are a function of the non-exceedance probability, F, and of the covariate z. In Eq. (14) h> denotes the transpose of h: For estimating R; it is convenient to separate the parameters of the distribution of the peak occurrence process from the parameters of the distribution of peak magnitudes and assume that the sample properties of the peak magnitudes are

j 6¼ 0; x [ l; j ¼ 0; x [ l;

ð12Þ

independent of the occurrence process. Under the former assumption, for the ZIP–GP model, we have h ¼ ðhZIP ; hGP Þ; with CovðhZIPi ; hGPj Þ ¼ 0; hence the matrix R becomes

RZIP 02;2 R¼ ; 02;2 RGP where 0m,n in a m-by-n null matrix, and RZIP and RGP are the asymptotic variance–covariance matrices of the parameters of the ZIP and GP, respectively, which can both be obtained numerically by inverting the Hessian of the respective loglikelihood functions at the point of maximum likelihood. Analogously, for the NSZIP–GP model, we have h ¼ ðhNSZIP ; hGP Þ; with CovðhNSZIPi ; hGPj Þ ¼ 0 and R becomes

04;2 RNSZIP R¼ ; 02;4 RGP where RNSZIP is the asymptotic variance–covariance matrix of the NZSZIP, and can also be obtained numerically following the methods given by Lambert (1992).

3 Data 3.1 Streamflow data The suitability of the ZIP distribution to describe the annual number of over-threshold peaks was studied using mean daily streamflow data from 8 gauging stations in Northern Portugal. The locations of the gauging stations are shown in Fig. 1. Table 1 shows the main characteristics of the streamflow samples. The samples in Table 1 constitute part of the data set used by Silva et al. (2012b). The data were acquired by the Portuguese Water Institute, INAG, and made available via the SNIRH database (Sistema Nacional de Informac¸o˜es de Recursos Hı´dricos, http://snirh.pt), which has high data quality standards and is the main source of Portuguese hydrological and

Stoch Environ Res Risk Assess Fig. 1 Map of Northern Portugal. Location of the gauging stations presented in Table 1 and the catchments under study

Table 1 Mean daily streamflow data Code

Name

Period of records (no. of years)

MAFD (mm)

Catchment area (km2)

S1

Castelo Bom

1957/1958–2003/2004 (47)

348

897

S2

Castro Daire

1945/1946–2003/2004 (59)

738

291

S3

Cunhas

1949/1950–2005/2006 (57)

824

338

S4

Ermida Corgo

1956/1957–2005/2006 (50)

888

291

S5

Fragas da Torre

1946/1947–2005/2006 (60)

992

660

S6

Ponte Juncais

1918/1919–1974/2075 (57)

492

604

S7

Ponte Sta. Clara-Da˜o

1921/1922–1972/2073 (52)

455

177

S8

Quinta das Laranjeiras

1942/1943–2005/2006 (64)

246

3,464

Sample code, name and periods of records of the samples; mean annual flow depth, MAFD (mm), and area of the catchments

hydrometeorological data used by researchers and practitioners of water resources engineering and science. The mean daily streamflow samples were selected based on two main criteria: (1) the rivers must have no significant regulation that could influence the watersheds’ response to floods, and (2) the data series’ length should be at least 30 years. A few missing values were filled in using corroborated methods for Portugal (Pulido-Calvo and Portela 2007). It should be mentioned that the period of records and the number of years shown in Table 1 refer to hydrological years, which in Portugal begin on 1 October. 3.2 North Atlantic Oscillation (NAO) data The North Atlantic Oscillation (NAO) is a prominent and recurrent pattern in climate variability of the Northern Hemisphere, which refers to a redistribution of atmospheric masses between the Arctic and the subtropical Atlantic

(Hurrell et al. 2003). Silva et al. (2012b) found a relationship between flood occurrence rates and the winter NAO index in the catchments presented in Table 1: the years with less floods tend to happen when the winter NAO is in a positive phase, and the years with a high number of flood occurrences tend to happen when the winter NAO is in a negative or neutral phase. Due to the previous observations, the winter NAO was used as a covariate in the application of the NSZIP–GP model. The utilized NAO data consists of the IcelandGibraltar index developed by the Climate Research Unit (Jones et al. 1997 ; http://www.cru.uea.ac.uk/cru/data/ nao/), based on instrumental pressure measurements in Gibraltar and SW Iceland back to 1821. Also, the winter season has been defined as November to March (NDJFM) in accordance with Jones et al. (1997), Trigo et al. (2005) and Silva et al. (2012b), although other definitions of the winter season could also be used (Osborn et al. 1999).

123

Stoch Environ Res Risk Assess

"

4 Applications CI100ð1aÞ% ðdÞ ¼ 4.1 Relative performance of the ZIP In this section, the suitability of the Poisson and ZIP distributions is investigated, regarding their capabilities of describing the random variable Y, annual number of overthreshold exceedances, in the studied watersheds. The selection of the threshold remains the most subjective aspect of POT modeling and there is no universally accepted rule for that selection. In fact, different thresholds lead to obtaining different samples of both the peak magnitudes and their times of occurrence. For that reason, rather than evaluating the suitability of the Poisson and ZIP distributions based on a single threshold, the methods presented in this section were applied to threshold values varying quasi-continuously from the lower annual maximum in each sample (so as to insure that there is at least one year i such that yi = 0) to the fifth largest independent peak. This interval contains ranges of thresholds which may not be valid to use in a POT framework: too low thresholds might lead to the sampling of dependent peaks and/or peaks that do not correspond to extreme events; and, on the other hand, too high thresholds lead to reduced sample sizes of peak magnitudes which are insufficient for making a proper estimation of parameters over the time scales of the analyzed samples. It is assumed that the range of valid thresholds for POT modeling is somewhere in the middle. In order to guarantee the independence of the peaks, it was imposed that the selected ones must be separated in time by a minimum of three days, and that the flow between two peaks should decrease below as much as two thirds of the earlier peak of the two, following a suggestion by NERC (1975) and Cunnane (1979). As an alternative, the methods proposed by Claps and Laio (2003) could be utilized. The suitability of the Poisson distribution was assessed via the dispersion index test proposed by Cunnane (1979). The dispersion index statistic, d, is the ratio between the sample variance and sample mean d¼

s2y y

ð15Þ

:

If Y* Poisson(k), then E(Y) = Var(Y) = k, hence the dispersion index should take values close to 1 if the data adjusts to the Poisson distribution. Furthermore, if Y* Poisson(k), then (Cunnane 1979; Beguerı´a 2005; Naghettini and Pinto 2007) N X ðyi  yÞ2 i¼1

y

 v2N1 ;

ð16Þ

where N is the number of years in the sample. From the result (16), 100(1 - a)% confidence intervals, CI100(1-a)%, for d can be constructed as follows

123

Fv1 2

N1

a 2

N1

;

 # Fv1 1  a2 2 N1

N1

;

ð17Þ

denotes the inverse of the chi-squared distriwhere Fv1 2 N1 bution with N - 1 degrees of freedom. The Poisson hypothesis should not be rejected if the estimated dispersion index is bounded by the limits of the confidence interval. Figure 2a, b shows the dispersion index test applied to the studied streamflow samples using increasing threshold values, u. In order to plot the dispersion index curves of all the samples on the same graphs, in Fig. 2a the abscissae are the standardized thresholds obtained by dividing u by the long term means, or modulus, qmod, while in Fig. 2b the abscissae are the mean annual number of occurrences. The figure suggests that the Poisson assumption should be rejected over a wide range of thresholds in most samples at a 5 % confidence level. Also, as the d statistic exceeds the upper limit of the confidence interval, the rejection is due to overdispersion (a larger variance relative to the mean). Given these results, it does not seem likely that the annual number of overthreshold floods in the North of Portugal conforms to the Poisson distribution. Furthermore, it is apparent that the results obtained here are similar to the ones obtained by NERC (1975) and Cunnane (1979), in Great Britain. Figure 2c shows the difference between the fraction of years with no exceedances (the emprirical mass at zero), p0, and the theoretical mass at zero of the Poisson distribution, PPoi(Y = 0), for increasing standardized thresholds. That figure shows that in all the catchments under analysis the Poisson distribution underestimated the number of years with zero exceedances for most threshold values. The performance of the ZIP distribution, relative to the Poisson distribution, was evaluated by the Akaike Information Criterion (Akaike 1974), given by ^ AIC ¼ 2k  2‘ðhÞ;

ð18Þ

^ is the logwhere k is the number of parameters and ‘ðhÞ likelihood function at the point of maximum likelihood. The model that minimizes the AIC is the one with the best performance (Galiatsatou and Prinos 2011). Figure 2d shows the difference in the AIC of the Poisson, AICPoi, and of the ZIP, AICZIP, distribution fits over increasing standardized thresholds. The figure suggests that, generally, the ZIP provides a better fit for the y samples for most thresholds (as AICPoi [ AICZIP) in the catchments under analysis. 4.2 Stationary flood frequency analysis using the ZIP– GP The following is an example of flood frequency analysis using the ZIP–GP model proposed in Sect. 2.1. The

3.5

(b)

10

15

20

1.5

2.0

2.5

3.0 5

S5 S6

0.5

Dispersion index 0

S3 S4

S7 S8

S1 S2

0.0

0.0

S1 S2

1.0

2.5 2.0 1.5 1.0 0.5

Dispersion index

3.0

(a)

3.5

Stoch Environ Res Risk Assess

25

0

1

S3 S4 2

S5 S6

3

4

S7 S8 5

6

y

u q mod

(c)

(d)

40 30 20

AICPoi − AICZIP

50

S1 S2 S3 S4 S5 S6 S7 S8

0

0.00

10

0.10

0.20

0.30

S1 S2 S3 S4 S5 S6 S7 S8

0

5

10

15

20

25

u q mod

0

5

10

15

20

25

u q mod

Fig. 2 a, b Dispersion index (continuous lines) and 95 % confidence intervals (dotted lines) plotted against a the standardized thresholds; b the mean annual number of occurrences. The confidence intervals correspond to the color attributed to each sample. c Difference

between the empirical and Poisson masses at zero, p0 and PPoi(Y = 0), respectively. d Difference in AIC of the Poisson and ZIP distribution fits, plotted against the standardized thresholds

application used the data from the gauging station S2— Castro Daire (Table 1). The threshold, u0, was selected using the mean residual life (MRL) or mean excess plot technique described by Davison and Smith (1990) and Coles (2001). The mean of the threshold excesses over each threshold u, for increasing values of u, provides an empirical estimate of the MRL. The threshold uo is then selected, such that the MRL plot is approximately linear for u [ u0. Confidence intervals, based on the normality assumption of the sample means can be added to the plot. Furthermore, once a candidate uo is defined, theoretical MRL functions can be established using the parameters of the GP (including the Exponential special case) distribution fitted to the POT data resulting from uo (Davison and Smith 1990) and Coles (2001),

The parameters of the GP were estimated by maximum likelihood estimation using the package ismev (Coles and Stephenson 2006) in R (R Core Team 2012). The threshold u0 = 60 m3/s (u0 =qmod ¼ 8:832; y ¼ 2:288) was selected based on the MRL technique. Figure 3 shows the MRL plot, complete with 95 % confidence bands and GP and Exponential MRL functions for u0 = 60 m3/s. The figure shows that the confidence bands can accommodate a linear, but not constant, MRL for u [ 60 m3/s. Table 2 shows the parameter estimates and standard errors of the Exponential and GP distributions. In Table 2, NSGP stands for nonstationary generalized Pareto as it will be clarified in Sect. 4.3. Since we are dealing with nested models, i.e., one model (GP) reduces to the other (Exponential) when certain parameters are fixed(j = 0), an asymptotic likelihood ratio test (LRT) may be applied to compare the fit of the two models. For comparing the fits of the alternate model, M1 with k parameters (hM1 ¼ ðh1 ; . . .; hk Þ), and the null model,   M0 with k - q parameters (hM0 ¼ h1 ; . . .; hkq Þ; the LRT

MRL ðuju0 ; r; jÞ ( rjðuu0 Þ ; 1þj ¼ r;

k [  1; u [ u0 ; r  uj [ 0; k ¼ 0; u [ u0 ; r [ 0:

ð19Þ

123

50

Stoch Environ Res Risk Assess

Parameter estimate (st. error)

Max. loglikelihood

p Value

AIC

30

Model Poisson k

2.2881 (0.1969)

-135.1

272.29

k

3.1522 (0.2927)

-120.6

\10-7

245.26

/

0.2740 (0.0633) -111.8

\10-3

231.51

ZIP

20 0

10

3 Mean residual life (m s)

40

Table 3 Parameter estimates and standard (st.) errors for the Poisson, ZIP and NSZIP (winter NAO covariate) models fitted to the sample of annual number of over-threshold occurrences, Y, for u0 = 60 m3/s; maximized log-likelihood and p value for likelihood ratio test comparing a model (alternate) with the one (null) immediately above; AIC

NSZIP k0

Empirical MRL 95% confidence bands GP MRL (u0 = 60) Exponential MRL (u0 = 60) 50

100

150

1.1094 (0.0980)

k1

-0.2782 (0.0938)

/0

-1.4932 (0.4421)

/1

0.8977 (0.3854)

u (m s) 3

Fig. 3 Empirical mean residual life plot with 95 % confidence bands; theoretical Exponential and GP mean residual life

Table 2 Parameter estimates and standard (st.) errors for the Exponential, GP and NSGP (winter NAO covariate) distributions fitted to the peak magnitudes above the threshold uo = 60 m3/s; maximized log-likelihood and p value for likelihood ratio test comparing a model (alternate) with the one (null) immediately above; AIC Model

Parameter estimate (st. error)

Max. loglikelihood

p Value

AIC

Exponential l = u0 r GP

-641.6

1285.1

42.6193 (3.6680)

l = u0 r

53.5390 (6.0275)

j

0.2510 (0.0758)

-638.5

0.013

1280.9

-638.0

0.343

1282.0

NSGP l = u0 r0

4.0085 (0.1203)

r1

-0.0519 (0.0507)

j

0.2931 (0.0901)

uses the asymptotic result (Casella and Berger 1990, pp. 488–492; Coles 2001, p. 35; Davison 2003, pp. 126–127) n o 2 ‘ðh^M1 Þ  ‘ðh^M0 Þ  v2q ; ð20Þ where ‘ðh^M1 Þ and ‘ðh^M0 Þ are the maximized log-likelihoods of the alternate and null models, respectively. Table 2 shows that the p value (0.013) of an LRT, where

123

the null model is the Exponential and the alternative model is the GP, suggests that the Exponential distribution should be rejected in favor of the GP at the 5 % significance level. That table also shows that AICExponential [ AICGP. For the selected threshold, the Poisson assumption is rejected, as d = 2.427 and CI95% ðdÞ ¼ ½0:670; 1:395 (Fig. 2a), and AICPoi [ AICZIP (Fig. 2d). Table 3 shows the parameter estimates and corresponding standard errors of the Poisson and ZIP models. Figure 4 shows the empirical PMF and CDF of the annual number of occurrences, Y, along with the fitted Poisson and ZIP models. Here it seems that the ZIP provides the best fit for the data. Table 3 also shows that the p value (\10-7) of an LRT, where the null model is the Poisson and the alternative model is the ZIP, suggests that the Poisson distribution should be rejected in favor of the ZIP. We also note that AICPoi [ AICZIP (Fig. 2d; Table 3). Given the difference in AIC, the outcomes of the dispersion index test, the LRT, and the plots of Fig. 4, there seems to be a strong evidence of the superiority of the ZIP distribution to describe the annual number of threshold exceedances in the present case study. Figure 5a shows the Exponential and GP distributions fitted to the XOT data. As suggested by the downward slope of the MRL plot in Fig. 3, and by the results of Table 2, the peak magnitudes do not seem to be exponentially distributed but rather GP distributed with a positive shape parameter j. In Fig. 5b, the Poisson–GP and ZIP–GP models are compared against the observed XAM sample. It is apparent that the two models coincide at the upper tail, both in quantile estimates and variance (denoted by the perfect agreement of the confidence envelopes and of the quantile curves), despite the ZIP– GP having one more parameter. Overall, however, the ZIP– GP model provides a better representation of the observed data, especially for lower quantiles, namely for F \ 0.75.

(b)

0.6 0.4 0.2

Distribution

0.8

Empirical Poisson ZIP

Empirical Poisson ZIP

0.0

Mass

(a)

0.00 0.05 0.10 0.15 0.20 0.25 0.30

Fig. 4 a PMF and b CDF of the empirical distribution of the annual number of occurrences and of the fitted Poisson and ZIP distributions

1.0

Stoch Environ Res Risk Assess

0

2

4

6

8

10

0

2

4

250 200 150

Annual maxima (m 3 s)

300 250 200 150

0

100

AMS sample Poisson−GP ZIP−GP Threshold

1

2

3

4

5

6

0.1

0.75

10

2

0.999

8 6 2

4

1 0 −1

Winter NAO index

0.99

Y (observed) λ(z) φ(z) λ (stationary) φ (stationary)

200 150 100

0.95

F Gumbel probability paper

(b)

Winter NAO

0

−2

50

Mean daily flow (m 3 s)

250

− Log(1 − G(x))

MDF Threshold

10

100

(b) POT Sample Exponential Generalized Pareto

50

Peaks over threshold (m 3 s)

(a)

0

(a)

8

50

Fig. 5 a POT sample; Exponential and GP fits and 95 % confidence envelopes (dashed curves). b AMS sample; Poisson–GP and ZIP– GP fits and 95 % confidence envelopes (dashed curves)

6

y

300

y

1950 1960 1970 1980 1990 2000

Year

−2

−1

0

1

2

Winter NAO index, z

Fig. 6 a Mean daily flows, MDF, at Castro Daire (S2), Winter NAO (NDJFM) index plotted on the 1st of January of the corresponding winter. b Annual number of over-threshold peak occurrences plotted

against the winter NAO index, z, of the corresponding year; NSZIP parameters as functions of z, k(z) and /(z); ZIP parameters, k and / (constant)

4.3 Nonstationary flood frequency analysis using the NSZIP–GP

application also used streamflow data from gauging station S2—Castro Daire (Table 1) and the NAO data presented in Sect. 3.2. Figure 6a shows the mean daily flows at Castro Daire and the winter NAO index in the period.

The following is an example of flood frequency analysis using the NSZIP–GP model presented in Sect. 2.2. The

123

Stoch Environ Res Risk Assess

(b) F=0.999 Threshold

0

50 −2

−1

0

1

Winter NAO index, z

123

z=−1.5 z=0 z=1.5 Threshold

50

Annual maxima (m 3 s)

F=0.900 F=0.990

100 150 200 250 300 350

time (Coles 2001, p. 106). Table 2 shows the p value (0.343) of an LRT test where the alternative model is the NSGP and the null model is the stationary GP, which determines that the stationarity of the peak magnitudes should not be rejected at the 5 % significance level. Furthermore, that table shows that AICNSGP [ AICGP. The previous results suggest that the distribution of the peak magnitudes does not vary as a function of the winter NAO, and that the assumption of their stationarity is not invalidated. In a POT model for flood frequency analysis, such as the ones given by Eqs. (3), (9), and (12), if the parameters of the annual number of occurrences model are allowed to change in time, then the probability of the annual maxima also varies even if the peak exceedances may be considered stationary. Figure 7 shows the results of the application of NSZIP–GP model for flood frequency analysis. Figure 7a shows observed annual maximum floods and flood quantiles with a non-exceedance probability, F, of 0.900, 0.990 and 0.999, as a function of the winter NAO index, z. Figure 7b shows the NSZIP–GP CDF of the annual maximum flood magnitudes for three values of z (-1.5, 0, and 1.5). The confidence intervals in Fig. 7 were estimated using the methods described in Sect. 2.3. The results of Fig. 7 show that, for a fixed nonexceedance probability, F, the annual maximum flood quantile decreases as the winter NAO index increases. This relationship seems to be in accordance with the observations, as can be seen in Fig. 7a, especially in the higher observations. Figure 7b shows that the probability of there being no exceedances above threshold, i.e., P (XAM \ u0) = P (Y = 0), increases as the winter NAO increases, which also seems to be in accordance with the observation of Fig. 7a, where the annual maximum floods lower than the threshold generally occur when the winter NAO is in its positive phase.

0

(a) Annual maxima (m 3 s)

Fig. 7 a Annual maximum flood magnitudes (circles) plotted against the corresponding winter NAO index; estimated quantiles of the NSZIP–GP model for 3 nonexceedance probabilities (continuous lines) and 95 % confidence envelopes (dashed lines). b NSZIP–GP model for 3 values of the winter NAO index (continuous lines) and 95 % confidence envelopes (dashed lines)

100 150 200 250 300 350

The application of the NSZIP–GP model used the threshold of u0 = 60 m3/s, selected in Sect. 2.1. The parameters of the NSZIP were estimated using the zeroinfl() function (Zeileis et al. 2008) in R. The estimated parameters and corresponding standard errors are shown in Table 3. That table also shows the p value (\10-3) of an LRT where the null model is the ZIP and the alternative model is the NSZIP, which suggests that the NSZIP model fits the data better than both the Poisson and (stationary) ZIP distributions. Table 3 also shows that the NSZIP is the model with the lowest AIC. Figure 6b shows the observed annual number of occurrences plotted against the winter NAO index of the corresponding year along with the parameter link functions k(z) and /(z) of the NSZIP–GP model (Eq. 12) as functions of the winter NAO. For the sake of comparison, Fig. 6b also shows the constant parameters k and / of the stationary ZIP model. The results of Fig. 6b and of Table 3 pertaining to the NSZIP model show that there is a statistically significant relationship between the number of over threshold peak floods in a year and the winter NAO index of the same year. That relationship is characterized by the number of occurrences of over-threshold floods in the catchment being inversely related to the winter NAO index. This result is in accordance with the findings by Silva et al. (2012b). In order to assess the stationarity assumption of the peak magnitudes a nonstationary generalized Pareto model, NSGP, was fitted to the peak magnitudes using the winter NAO as a covariate via a log link function applied to the scale parameter r, log r(z) = r0 ? r1 z , using the methods described in Coles (2001) and Katz (2013). The shape parameter, j, was kept constant as extreme value shape parameters are difficult to estimate with precision and it is usually unrealistic to allow them to vary with

2

0.1

0.75

0.95

0.99

F Gumbel probability paper

0.999

Stoch Environ Res Risk Assess

5 Conclusions This paper introduces the ZIP distribution as a valid alternative to the Poisson distribution for modeling the annual arrival counts of over-threshold peaks, Y, in a POT approach for flood frequency analysis. That distribution consists of combining a point mass at zero with the Poisson distribution and is able to adequately model Y when the observed fraction of zero counts is significantly higher than the theoretical mass at zero of the Poisson distribution. Using streamflow data from 8 gauging stations located in the North of Portugal, it was found that the assumption that Y is described by a Poisson distribution should be rejected over a wide range of possible thresholds, and that the departure from the Poisson is characterized by the variance of Y being significantly larger than its mean mean. It was also found that the Poisson consistently underestimates the probability of Y = 0. The ZIP distribution was fitted to the same data and it was found to provide a better fit than the Poisson. A ZIP–GP model was proposed where the ZIP model for Y is combined with the GP for threshold exceedance magnitudes, under a POT framework. The performance of the ZIP–GP was tested and compared to the commonly used Poisson–GP model based on an application to daily streamflow data at one of the gauges. In that application it was found that the ZIP–GP provided a better fit to the observed annual maxima than the Poisson–GP model. Moreover, despite the ZIP having one more parameter, there is no increment in uncertainty for higher quantiles of the ZIP–GP model. It was also found that the for high nonexceedance probabilities the two models coincide both in quantile estimates and asymptotic variance. These results suggests that an accurate description of the random nature of the occurrence process is more important for less extreme (more frequent) floods, while it bears little important when dealing with very rare floods. Hence, if the modeler is concerned only with upper quantiles, the Poisson distribution may still be used, even if the overall adherence to data is not confirmed by a statistical test. The application described herein also demonstrated how nonstationarity may be introduced in the proposed model for annual maxima via the parameters of the ZIP distribution, thereby obtaining the NSZIP–GP, a nonstationary model which is useful for modeling extremes when the occurrence rate of peaks changes over time. The threshold exceedance peak magnitudes, modeled by the GP, may be assumed stationary (constant GP parameters) or not (varying GP parameters). It is therefore possible to obtain a nonstationary model for annual maxima where nonstationarities are modeled separately for the arrival counts and the peak magnitudes.

The methods presented in this paper are a contribution to researchers and practitioners of flood risk assessment as it enables: (i) the application of the POT approach in cases where the Poisson assumption is rejected; and (ii) to consider a dynamic planning and operation paradigm when dealing with uncertainty associated with hydrological extremes under nonstationarity, with the advantage of being able to model the occurrence process of such extremes and their peak magnitudes separately. Acknowledgments The research was supported by the Portuguese Science and Technology Foundation (FCT) through a scolarship for A. T. Silva (Grant SFRH/BD/86522/2012). The authors also which to thank two anonymous reviewers for their valuable and insightful comments.

Appendix 1: Log-likelihood function of the ZIP and parameter estimation The PMF of the ZIP (5) can be written as PðY ¼ yÞ ¼ ½/ þ ð1  /Þ expðkÞIf0g ðyÞ  ky 1If0g ðyÞ  ð1  /Þ expðkÞ ; y!

ð21Þ

where the indicator function I{0}(y) is equal to one when y = 0 and zero otherwise. The likelihood function becomes ( n Y Lðk; /jyÞ ¼ ½/ þ ð1  /Þ expðkÞIf0g ðyi Þ i¼1



k yi  ð1  /Þ expðkÞ yi !

1If0g ðyi Þ )

ð22Þ ;

and the log-likelihood function ( n X ‘ðk; /jyÞ ¼ If0g ðyi Þ log½/ þ ð1  /Þ expðkÞ i¼1

 ) ky i þ 1  If0g ðyi Þ log ð1  /Þ expðkÞ : yi ! 



ð23Þ ^ kg; ^ estimates are obtained by maximizing (23). The f/;

Appendix 2: Expansion of (7) Here the expansion of (7), resulting in (8) is shown. From (5), PðY ¼ 0Þ ¼ / þ ð1  /Þ expðkÞ and y PðY ¼ yÞ ¼ ð1  /Þ expðkÞ ky! for y [ 0. Hence

123

Stoch Environ Res Risk Assess

FðxÞ ¼ / þ ð1  /Þ expðkÞ  1 X ky ð1  /Þ expðkÞ ½GðxÞy þ y! y¼1 ¼ / þ ð1  /Þ expðkÞ  1 X ½kGðxÞy þ ð1  /Þ expðkÞ y! y¼1 ! 1 X ½kGðxÞy ¼ / þ ð1  /Þ expðkÞ 1 þ y! y¼1

 y 1 X ½kGðxÞ ¼ / þ ð1  /Þ expðkÞ y! y¼0 ¼ / þ ð1  /Þ expðkÞ exp½kGðxÞ ¼ / þ ð1  /Þ expfk½1  GðxÞg ð24Þ

References Akaike H (1974) A new look at the statistical model identification. IEEE Trans Autom Control 19(6):716–723 Balkema AA, de Haan L (1974) Residual life time at great age. Ann Probab 2(5):792–804 Beguerı´a S (2005) Uncertainties in partial duration series modelling of extremes related to the choice of the threshold value. J Hydrol 303(1):215–230 Ben-Zvi A (1991) Observed advantage for negative binomial over Poisson distribution in partial duration series. Stoch Hydrol Hydraul 5(2):135–146 Bhunya PK, Singh RD, Berndtsson R, Panda SN (2012) Flood analysis using generalized logistic models in partial duration series. J Hydrol 420:59–71 Casella G, Berger RL (1990) Statistical inference, vol 70. Duxbury Press, Belmont Claps P, Laio F (2003) Can continuous streamflow data support flood frequency analysis? An alternative to the partial duration series approach. Water Resour Res 39(8):1216 Clarke RT (2007) Hydrological prediction in a non-stationary world. Hydrol Earth Syst Sci 11(1):408–414. doi:10.5194/hess-11-4082007 Coles S (2001) An introduction to statistical modeling of extreme values. Springer, London Coles S, Stephenson A (2006) Ismev: an introduction to statistical modeling of extreme values. Original S functions by Stuart Coles and R port and R documentation files by Alec Stephenson Cunnane C (1979) A note on the Poisson assumption in partial duration series models. Water Resour Res 15(2):489–494 Davison AC (2003) Statistical models. Cambridge University Press, Cambridge Davison AC, Smith RL (1990) Models for exceedances over high thresholds. J Roy Stat Soc B Met 52(3):393–442 Delgado JM, Apel H, Merz B (2010) Flood trends and variability in the Mekong river. Hydrol Earth Syst Sci 14(3):407–418 Delgado JM, Merz B, Apel H (2012) A climate-flood link for the lower Mekong River. Hydrol Earth Syst Sci 16(5):1533–1541 Ekanayake ST, Cruise JF (1993) Comparisons of Weibull-and exponential-based partial duration stochastic flood models. Stoch Hydrol Hydraul 7(4):283–297

123

Fisher RA, Tippett LHC (1928) Limiting forms of the frequency distribution of the largest or smallest member of a sample. In: Math Prof Cambridge, Cambridge Univ Press, vol 24, pp 180–190 Galiatsatou P, Prinos P (2011) Modeling non-stationary extreme waves using a point process approach and wavelets. Stoch Environ Res Risk Assess 25(2):165–183 Gnedenko B (1943) Sur la distribution limite du terme maximum d’une se´rie ale´atoire. Ann Math 44(3):423–453 Grimm AM, Tedeschi RG (2009) ENSO and extreme rainfall events in South America. J Clim 22(7):1589–1609 Hundecha Y, St-Hilaire A, Ouarda TBMJ, El Adlouni S, Gachon P (2008) A nonstationary extreme value analysis for the assessment of changes in extreme annual wind speed over the Gulf of St. Lawrence, Canada. J Appl Meteorol Clim 47(11):2745–2759 Hurrell JW, Kushnir Y, Ottersen G, Visbeck M (2003) An overview of the North Atlantic oscillation. In: Hurrell J, Kushnir Y, Ottersen G, Visbeck M (eds) The North Atlantic Oscillation: climatic significance and environmental impact. Geophysical Monograph, vol 134. AGU American Geophysical Union, pp 1–36 Jones PD, Jonsson T, Wheeler D (1997) Extension to the North Atlantic Oscillation using early instrumental pressure observations from Gibraltar and south-west Iceland. Int J Climatol 17(13):1433–1450 Katz RW (2010) Statistics of extremes in climate change. Clim Change 100(1):71–76 Katz RW (2013) Statistical methods for nonstationary extremes. In: Extremes in a changing vlimate. Springer, New York, pp 15–37 Koutsoyiannis D (2011) Hurst–Kolmogorov dynamics and uncertainty. J Am Water Resour Assoc 47(3):481–495 Kwon HH, Sivakumar B, Moon YI, Kim BS (2011) Assessment of change in design flood frequency under climate change using a multivariate downscaling model and a precipitation-runoff model. Stoch Environ Res Risk Assess 25(4):567–581 Lambert D (1992) Zero-inflated Poisson regression, with an application to defects in manufacturing. Technometrics 34(1):1–14 Lang M, Ouarda TBMJ, Bobe´e B (1999) Towards operational guidelines for over-threshold modeling. J Hydrol 225(3–4):103–117 Lee T, Ouarda TBMJ (2010) Long-term prediction of precipitation and hydrologic extremes with nonstationary oscillation processes. J Geophys Res Atmos 115:D13 Lins HF, Cohn TA (2011) Stationarity: wanted dead or alive? J Am Water Resour Assoc 47(3):475–480 Lo´pez J, France´s F (2013) Non-stationary flood frequency analysis in continental Spanish rivers, using climate and reservoir indices as external covariates. Hydrol Earth Syst Sci Discuss 10:3103–3142 Madsen H, Rasmussen PF, Rosbjerg D (1997) Comparison of annual maximum series and partial duration series methods for modeling extreme hydrologic events: 1. At-site modeling. Water Resour Res 33(4):747–757 Milly PCD, Betancourt J, Falkenmark M, Hirsch RM, Kundzewicz Z, Lettenmaier D, Stouffer R (2008) Stationarity is dead: whither water management? Science 319(5863):573 Naghettini M, Pinto EJA (2007) Hidrologia estatı´stica. CPRM, Belo Horizonte NERC (1975) Flood Studies Report, vol I. Natural Environment Research Council, London ¨ no¨z B, Bayazit M (2001) Effect of the occurrence process of the O peaks over threshold on the flood estimates. J Hydrol 244(1):86–96 Osborn TJ, Briffa KR, Tett SFB, Jones PD, Trigo RM (1999) Evaluation of the North Atlantic Oscillation as simulated by a coupled climate model. Clim Dyn 15(9):685–702 Pickands J (1975) Statistical inference using extreme order statistics. Ann Stat 3(1):119–131 Pulido-Calvo I, Portela MM (2007) Application of neural approaches to one-step daily flow forecasting in Portuguese watersheds. J Hydrol 332(1–2):1–15

Stoch Environ Res Risk Assess R Core Team (2012) R: a language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. ISBN: 3-900051-07-0 Renard B, Lang M, Bois P (2006) Statistical analysis of extreme events in a non-stationary context via a Bayesian framework: case study with peak-over-threshold data. Stoch Environ Res Risk Assess 21(2):97–112 Rosbjerg D, Rasmussen F (1991) Modelling of exceedances in partial duration series. In: Proceeding of the international hydrology and water resources symposium, vol 3. National Conference Publication—Institute of Engineers, pp 755–760 Rosbjerg D, Madsen H, Rasmussen PF (1992) Prediction in partial duration series with generalized Pareto-distributed exceedances. Water Resour Res 28(11):3001–3010 Salas JD, Obeysekera J (2013) Revisiting the concepts of return period and risk for nonstationary hydrologic extreme events. J Hydrol Eng Shane RM, Lynn WR (1964) Mathematical model for flood risk evaluation. J Hydraul Eng Div 90(HY 6):4119–4122 ¨ tsch R, Kemfert C, Jaeger CC, Haas A, Kremers H Siliverstovs B, O (2010) Climate change and modelling of extreme temperatures in Switzerland. Stoch Environ Res Risk Assess 24(2):311–326 Silva AT, Portela MM, Baez J, Naghettini M (2012a) Construction of confidence intervals for extreme rainfall quantiles. Risk Anal VIII:167–293

Silva AT, Portela MM, Naghettini M (2012b) Nonstationarities in the occurrence rates of flood events in Portuguese watersheds. Hydrol Earth Syst Sci 16(5):241–254. doi:10.5194/hessd-88609-2011 Singh SN (1963) A note on inflated Poisson distribution. J Indian Stat Assoc 1(3):140–144 Tramblay Y, Neppel L, Carreau J, Najib K (2013) Non-stationary frequency analysis of heavy rainfall events in southern France. Hydrol Sci J 58(2):280–294 Trigo RM, Zezere JL, Rodrigues ML, Trigo IF (2005) The influence of the North Atlantic Oscillation on rainfall triggering of landslides near Lisbon. Nat Hazards 36(3):331–354 Villarini G, Smith JA, Serinaldi F, Bales J, Bates PD, Krajewski WF (2009) Flood frequency analysis for nonstationary annual peak records in an urban drainage basin. Adv Wat Resour 32(8):1255–1266 Villarini G, Smith JA, Serinaldi F, Ntelekos AA (2011) Analyses of seasonal and annual maximum daily discharge records for central Europe. J Hydrol 399(3):299–312 Zeileis A, Kleiber C, Jackman S (2008) Regression models for count data in R. J Stat Softw 27(8):1–25 Zelenhasic E (1970) Theoretical probability distributions for flood peaks. Hydrology paper 42. Colorado State University, Fort Collins

123