Nonstationarities in the occurrence rates of flood events in Portuguese watersheds

An exploratory analysis on the variability of flood occurrence rates in 10 Portuguese watersheds is made, to ascertain if that variability is concurrent with the principle of stationarity. A peaks-over-threshold (POT) sampling technique is applied to 10 long series of mean daily streamflows and to 4 long series of daily rainfall in order to sample the times of occurrence (POT time data) of the peak values of those series. The kernel occurrence rate estimator, coupled with a bootstrap approach, was applied to the POT time data to obtain the time dependent estimated occurrence rate curves,λ̂(t), of floods and extreme rainfall events. The results of the analysis show that the occurrence of those events constitutes an inhomogeneous Poisson process, hence the occurrence rates are nonstationary. An attempt was made to assess whether the North Atlantic Oscillation (NAO) casted any influence on the occurrence rate of floods in the study area. Although further research is warranted, it was found that years with a less-than-average occurrence of floods tend to occur when the winter NAO is in the positive phase, and years with a higher occurrence of floods (more than twice the average) tend to occur when the winter NAO is in the negative phase. Although the number of analyzed watersheds and their uneven spatial distribution hinders the generalization of the findings to the country scale, the authors conclude that the mathematical formulation of the flood frequency models relying on stationarity commonly employed in Portugal should be revised in order to account for possible nonstationarities in the occurrence rates of such events.


Introduction
Nowadays, there seems to be a consensus among the scientific community that, due to climate change, there is an intensification of the hydrological cycle (Bates et al., 2008), assumably leading to more frequent and more intensive extreme hydrological phenomena, like floods and droughts.Milly et al. (2008) argue that such climate change undermines stationarity -a basic assumption that historically has assisted practice and research in the fields of hydrology and water resources management.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 a certain magnitude.
In Portugal, however, the generality of research on hydrological modeling relies on the principle of the stationarity of hydrological time series (Quintela and Portela, 2002), including recent academic works on flood hydrology (Portela and Dias, 2005;Delgado, 2007;Portela and Delgado, 2009), although general circulation models (GCMs) based on some scenarios of greenhouse gas emissions show an expected trend for an aggravation of extreme precipitation events in northern Portugal (Santos et al., 2002, p. 167).In such framework it is important to carry out research that could ascertain whether or not the hydrological time series that are used in the design of water resources systems exhibit definite signs of nonstationarity.If they do, the mathematical formalism applied to the planning and management of such structures must be revised, namely the statistical analyses that are Published by Copernicus Publications on behalf of the European Geosciences Union. A. T. Silva et al.: Nonstationarities in the occurrence rates of flood events in Portuguese watersheds usually employed in studies on extreme hydrological phenomena, such as extreme rainfall events and floods.
According to Ramos and Reis (2002), floods were the deadliest natural disasters during the 20th century in Portugal, followed by earthquakes.Those authors differentiate rain-induced floods in Portugal into two types: progressive floods, which are caused by prolonged heavy rainfall associated with a westerly zonal circulation that may persist for weeks; and, in contrast, flash floods, which are triggered by heavy and concentrated rainfall related to convective depressions (e.g.: active cold pools, depressions caused by interaction between polar and tropical circulations, and, very rarely, tropical depression).Progressive floods mainly affect larger watersheds such as the international Tejo (Tagus) and Douro river basins, which are heavily regulated with reservoirs and flood control systems, while flash floods mainly affect smaller watersheds.Due to the small size of most of the national watersheds (with times of concentration lower than one day), Portugal is particularly prone to flash flood events, which pose a great danger for human populations.
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).Studies carried out by Hurrell (1995) and Trigo et al. (2002a) have established links between the NAO phase and precipitation in western Europe.Although the NAO is evident throughout the year, its activity and impact on European surface climate is greater during the winter season (Osborn et al., 1999;Morán-Tejeda et al., 2011).There are also a number of studies on the influence of the NAO on precipitation and river flow in the western Iberian Peninsula in winter months (Rodríguez-Puebla et al., 2001;Corte-Real et al., 1998;Trigo et al., 2002bTrigo et al., , 2004;;Morán-Tejeda et al., 2011;Lorenzo-Lacruz et al., 2011), which have shown that when the NAO is in its negative phase precipitation and river flows tend to be above normal.Trigo et al. (2005) studied the influence of the NAO on monthly and seasonal precipitation on an area prone to landslides located north of the capital of Portugal, Lisbon, and concluded that the large inter-annual variability of winter precipitation in that region is largely modulated by the NAO mode.
The aim of the present study was to carry out an exploratory analysis on the variability of the streamflow regime in 10 mainland Portugal watersheds, in what concerns the occurrence rate of floods, and ascertain if that variability is concurrent with the hypothesis of nonstationarity.
This paper focuses primarily on mean daily streamflow data at 10 gauging stations geographically spread over mainland Portugal.Long series of daily rainfall records at rain gauges located in 4 of the studied watersheds were also utilized.The peaks-over-threshold (POT) sampling technique is applied to extract the dates (POT time data) and peak values (POT value data) of flood events and extreme rainfall events, respectively, from the streamflow and rainfall samples.
A nonparametric method for analyzing nonstationarities in the occurrence rate λ(t) of floods and extreme rainfall events was applied to the POT time data -the kernel occurrence rate estimation method, or kernel technique (Diggle, 1985;Mudelsee et al., 2003).
Given the aforementioned evidence of the role of the NAO in modulating rainfall and river flows in western Iberian watersheds, an exploratory analysis was conducted on the relationship between flood occurrence in the studied watersheds and the NAO.
Following this Introduction, in Sect.2, the utilized data sets as well as the methods used to sample POT data and to detect nonstationarities are presented.Section 3 presents and discusses the results of the study, and, finally, in Sect.4, the most relevant conclusions of the article are drawn and opportunities of future research are presented.

Streamflow and rainfall data and characteristics of the studied watersheds
As mentioned in Sect. 1, this paper focuses primarily on flood occurrence rates in 10 unregulated Portuguese catchments.For that purpose, two types of hydrological time series -identified in  Portela, 2007).It should be mentioned that the period of records and number of years shown in Table 1 refer to hydrological years, which in Portugal begin on 1 October.The daily rainfall data were utilized for the purpose of checking whether there are significant discrepancies between peaks in rainfall and streamflow, which would eventually suggest that the behaviour exhibited by the streamflows was under a significant anthropogenic influence.Such analysis allowed to consider that the flood regime in the studied watersheds is nearly pristine.
Figure 1b shows the location of the streamflow and rain gauging stations as well as the drainage divides of the        In panel (e), the vertical ticks show the points in time associated with the flood occurrences, the dotted diagonal line shows the uniform distribution and the solid diagonal lines show significance for a Kolmogorov-Smirnov statistic at 5 % and 1 %.catchments under study, whose areas are shown in Table 1.It should be commented, that the criteria for selection of data series resulted in a small number of studied watersheds which are, moreover, unevenly distributed in the Portuguese territory, as made evident in Fig. 1b: most of the studied watersheds are in the North of Portugal, only two are in the South and the central region is not represented.Consequently, it is clear that the results obtained from this data set should not be unquestionably generalized for the whole country.
As shown in Fig. 1a, in Portugal, the mean annual rainfall (MAR) varies from more than 2800 mm, in the northwestern region, to less than 400 mm, in the southern region, following a complex spatial pattern (N-S/W-E).Table 1 shows the mean annual flow depth, MAFD, and mean annual rainfall, MAR, of the streamflow and rainfall samples respectively.Those values show the diversity of hydrological regimes represented in the data set, the driest catchment being S6 -Monte da Ponte, in the South, and the wettest being S5 -Fragas da Torre, in the North.As for the represented hydrometeorological regimes, Fig. 1a and Table 1 show that, while the three rain gauging stations in the North are in very humid regions (R2, R3 and R4), the R1 station, in the South, is in a relatively dry region.

Peaks-over-threshold (POT) sampling and data
The peaks-over-threshold (POT) approach to hydrological frequency analysis, consists of applying a sampling technique to a time series that retains the peak values that exceed a given truncation level usually called base level or threshold (Lang et al., 1999).Todorovic and Zelenhasic (1970) and Todorovic (1978) introduced the marked point model for flood analysis that relies on the POT sampling technique.
Conventional analysis by the POT sampling technique consists of the following: from a continuous time series x t , with t running continuously from t 0 to t n , the i-th peak value x i corresponds to the largest value among those exceeding the threshold u, in the time interval [t 0 + t i , t 0 + t i + t], comprised in [t 0 , t n ].Formally, 16, 241-254, 2012 www.hydrol-earth-syst-sci.net/16/241/2012/If the time of occurrence of x i is denoted by T i ∈ [t 0 , t n ] and if there are m flood episodes of this kind, {T i , x i } m i=1 represents a marked point process for the variable X.
Usually, the peak values, x i , are supposed identically independently distributed (i.i.d.) variables.The time of occurrence of x i , that is T i , is associated with a Poisson process, with the number of occurrences in a given interval (e.g. a year) being a Poisson variate with parameter λ (Cunnane, 1979).In the context of POT sampling of hydrological variables, the time series discretized in daily intervals, such as those employed here, though not rigorously continuous, conform to the general requirements of a Poisson point process, as previously defined.
The selected peaks must meet the independence condition.Several criteria has been presented in the literature in order to verify this hypothesis (Lang et al., 1999).The criterion used in the present study, for the mean daily flow series, (NERC, 1975;Cunnane, 1979) determines that peaks should be separated in time by three times the time to peak and, furthermore, that the flow between two consecutive peaks should decrease below as much as two thirds of the first peak.Due to the relatively small areas of the watersheds being analyzed (Table 1), as most of the Portuguese watersheds, times of concentration lower than 1 day are expected, except for Quinta das Laranjeiras -S10.For that reason and given the daily time-step of the data, we considered that for the purpose of applying the forementioned independence criterion, the time to peak equals approximately one day at the analyzed gauging stations.
The independence criterion applied to the daily rainfall series is similar to the one applied to the mean daily flow series, the only difference being that the daily rainfall must decrease to zero between two peaks.
The selection of the threshold, u, is a procedure that involves a great level of subjectivity (Beguería, 2005).Lang et al. (1999) reviewed a number of systematic approaches to carry out this selection.Lang et al. (1999) remark, however, that there is no universal and unequivocal method for selecting the threshold, and that there is no unique threshold value that must be selected but rather a range of appropriate values.In the present study, in what concerns the mean daily flows, the threshold adopted equals 7 times the long term mean daily flow, or modulus, which according to Quintela (1984) provides a lower limit to identify the flood occurrences.Figure 3a shows the mean daily flows sample at Fragas da Torre (S5) and the corresponding threshold.The criterion adopted to determine the threshold, as applied to the daily rainfall records, considered that the mean number of over-threshold values should be similar to the one obtained for the mean daily flows.Such criterion resulted in a threshold equal to 12 times the mean daily rainfall.Table 1 shows the adopted threshold values for the mean daily streamflow and the daily rainfall series.
After the threshold u was defined, the correlation structure of the peaks or exceedances was analyzed, namely the lag-  one and lag-two autocorrelation coefficients, as exemplified in Fig. 3b.No significant serial correlation was found that could invalidate the application of the selected thresholds.Additionally, tests were carried out to evaluate the adequacy of the thresholds.Those tests focused on the mean number of over-threshold events in a year and on the mean exceedance above threshold as a function of variable threshold values, tests no. 1 and 2, respectively, as reviewed by Lang et al. (1999).The application of such tests did not invalidate the adopted threshold selection criterion, as it was found that the thresholds were such that a convenient overall number of events was sampled for the purpose of estimating λ(t), and they were generally in a domain where the mean exceedance above threshold is an approximately linear function of u.As an example, the results of the application of the tests to sample S5 -Fragas da Torre are presented in Fig. 3c  and d.
The inherent subjectivity regarding selection of the thresholds remains a source of ambiguity in statistical analyses based on POT data, despite the aforementioned efforts to validate that selection.Therefore it is important, in practice, to check if the results of such analyses are highly dependent on the threshold.For that purpose, the procedure suggested by Davison (2003, p.286), of imposing slight variations of threshold levels, around the adopted ones, and repeating the methods, did not cause any significant changes to the overall results presented in Sect.3, and to the final conclusions as well.

North Atlantic Oscillation (NAO) data
Traditionally, the NAO index has been defined as the difference in normalized surface pressure between Iceland (Stykkisholmur) and the Azores archipelago (Ponta Delgada).In recent decades, however, researchers have found that, during the winter season, stations located in the Iberian Peninsula could be used with some advantages over Ponta Delgada: Lisbon was used by Hurrell (1995) Trigo et al. (2005), although other definitions of the winter season could also be used (Osborn et al., 1999).Figure 3 shows the winter NAO indices used in this study along with a LOWESS (locally weighted regression) curve (Cleveland, 1979) fitted to the data using a smoothing parameter f = 0.1 in order to obtain an adequate inter-annual smoothing.

General remarks
The kernel technique is applied to the POT time data sampled from the original mean daily flow and daily rainfall time series.After the estimated flood and extreme rainfall occurrence rate curves, λ(t), are obtained, a pointwise bootstrap confidence band is constructed around them, allowing for a more rigorous interpretation of the results.
In the framework of this paper, λ(t) denotes the occurrence rate of flood events (Mudelsee, 2010, p. 249) when referring to mean daily flow, and the occurrence rate of extreme rainfall events when referring to daily rainfall.

Kernel occurrence rate estimation
The kernel technique is a nonparametric method developed by Diggle (1985) for smoothing point process data.For estimating the intensity of a point process such as the timedependent occurrence rate, λ(t), this technique may be formulated as: where K is the kernel function and h is the bandwidth.A Gaussian kernel was used as it can be efficiently calculated in Fourier space and yields a smooth estimated occurrence rate, λ(t) (Mudelsee et al., 2004(Mudelsee et al., , 2006)): The units of λ(t) are d −1 , i.e. the number of occurrences above threshold per day at a given point in time, t.However, to facilitate the interpretation of the results, λ(t) was multiplied by 365.25, such that, for a given instant, t, the λ(t) indicates the estimated number of occurrences above threshold per year (yr −1 ) The application of Eq. ( 2) may lead to a boundary bias near t 0 and t n consisting of an underestimation of λ(t) due to the nonexistence of data outside the interval [t 0 , t n ].This boundary effect can be reduced by generating pseudodata, i.e. pseudo extreme events, pT, outside of the observation interval, before estimating λ(t).The straightforward method of reflection was used to generate pseudodata (on the left side, for t<t 0 : , covering an amplitude of 3 times h before t 0 ; analogously on the right side, for t>t n ).Pseudodata generation is equivalent to the extrapolation of the empirical distribution of events near the boundaries, hence the estimation of λ(t) near the boundaries of the observation period should be analyzed with caution (Mudelsee et al., 2004;Mudelsee, 2010, p.251;Mudelsee, 2011).Considering T † as the original point data, T , augmented by the pseudodata, pT, and m † as the total number of points in T † , Eq. ( 2) can be rewritten as: Figure 4 shows the estimated occurrence rates λ(t) for the S8 -Ponte Juncais sample, with and without pseudodata generation, exemplifying the correction of the boundary bias via pseudodata generation.
The selection of the bandwidth, h, determines the bias and variance properties of the occurrence rate estimator λ(t): a too small h results in fewer data points that effectively contribute to the kernel estimation, which leads to a reduced bias and a high variance; on the other hand, a too large h leads to an oversmoothing of the estimator, resulting in a small variance and increased bias.The selection of the bandwidth can be seen as a compromise between those two cases.Furthermore, since there is a high seasonal variability of the hydrologic regime in Portugal and, as the objective is to describe inter-annual variability of λ(t), the bandwidth should be considerably higher than the year to avoid the effect of the seasonal variability.
In this study the selection of the bandwidth used a straightforward method: Silverman's rule of thumb (Silverman, 1986, p.48).This rule defines h as: h = 0.9 A −1/5 (5) with: where STD and IQR are, respectively, the standard deviation and the interquartile range of the POT time data.Silverman (1986, p. 48) comments that this method produces an adequate choice of bandwidth for many purposes.The bandwidths obtained range from 1141 days (sample S6) to 2157 days (sample S9), which is adequate for describing inter-annual variability, because it smooths out the withinthe-year seasonal variability.

Bootstrap confidence band
As previously described, a nonparametric method for estimating a varying point process intensity, the kernel method (Diggle, 1985), was used to characterize the occurrence rate of floods in the watersheds analyzed in this study.Point estimates may be difficult to interpret without some measurement of the uncertainties associated with those estimates.
For the purpose of quantifying the uncertainties associated with the results of Eq. ( 4), a pointwise confidence band was constructed around λ(t), by means of bootstrap simulations (Cowling et al., 1996;Mudelsee, 2011).This procedure consisted of drawing a set of n flood dates from the original POT data, augmented by pseudodata, with replacement, and calculating λ * (t) after Eq. ( 4), using the resampled data and the same bandwidth, h.The resampling-estimation procedure was repeated until 2000 estimated curves λ * (t) were obtained.Finally, an algorithm developed by Cowling et al. (1996) and used by Mudelsee et al. (2003Mudelsee et al. ( , 2004) ) and Mudelsee (2011), namely the percentile-t type confidence band was applied pointwise to the 2000 λ * (t) to construct a 90 % bootstrap confidence band around λ(t).

Preliminary data analysis
With the application of a POT sampling technique (Eq.1), two types of data were obtained: time data, {T i } m i=1 , which are the instants of occurrence of extreme events, and peak value data, {x i } m i=1 , which relates the magnitude of the extreme events themselves.
Although this study concerns primarily the forementioned time data, an exploratory analysis was made concerning the inter-annual variability of the peak values, {x (i)} m i=1 , for both the mean daily flow and the daily rainfall data.This analysis was performed by obtaining the series of running average exceedances over successive periods of 15 yr. Figure 5 shows those series, made dimensionless by their respective means, for both streamflow and rainfall data.In the horizontal axis of that figure, the year marks refer to the first trimester of the hydrologic year, which starts on 1 October of year k as previously mentioned.Although a visual analysis of Fig. 5 suggests that overall the samples do not exhibit trends or a significant inter-annual persistences above or below mean, linear relations were fitted to the running means and the significance of the slope parameter, β, was tested, and it was found that, out of the 10 series of running averages, only 3 of them have a slope that significantly differs from zero: S6, S8 and S9.It should be stressed that the beginning and the ending dates of the running averages vary according to the sampling period of each case study.
A preliminary analysis was also carried out concerning the parameters λ (occurrence rate) of the Poisson processes associated with the POT time data.Any counting process with stationary independent increments, i.e. the distribution of the number of events occurring in any interval depending only on the length of the interval, is a homogeneous Poisson process with parameter λ.With the introduction of a time dependence in the process, the Poisson parameter becomes a time function, λ(t).This time dependence is useful for studying nonstationarities in λ(t).When the process is homogeneous, λ is a constant and the variable w i = (T i − t 0 )/(t n − t o ) (normalized times of events above threshold) is distributed as order statistics of random sample from the uniform distribution in the interval [0, 1].A graphical analysis of that hypothesis can be performed by plotting the empirical distribution function of w i , F (w), and checking if there are any significant departures from the uniform distribution F (w) = w, 0≤w≤1 (Davison, 2003, p. 277-278).Figure 3e exemplifies the former graphical analysis applied to sample S5 -Fragas da Torre where the solid grey lines indicate the significance for a Kolmogorov-Smirnov statistic at levels 5 % and 1 %.Such figure suggests that the rate of the process may be non-uniform.It was verified that the generality of the samples exhibited significant departures from the uniform distribution.
An additional evaluation of whether there was any significant increase or decrease in the intensity of those processes focused on the number of extreme events per hydrologic year k, a discrete variable denoted by λ k .For each sample of both streamflow and rainfall data sets, λ k was smoothed by means of a LOWESS (locally weighted regression) curve (Cleveland, 1979), using a smoothing parameter f = 0.2, which was the smallest f that enabled an adequate smoothing of the data.The former analysis is represented in Fig. 6.
Figure 6 shows that there is some variability in the occurrence rate of extremes in the observed samples.Most notably, there seems to be an increase in the intensity of the Poisson process in the period ranging from the late 1950s to the late 1960s in all of the samples that cover this period, though in some samples it is more pronounced than in others.

Nonstationarity of flood occurrence rates
The techniques described in Sects.2.4.2 and 2.4.(grey area).The vertical ticks indicate the points in time when events occurred (POT time data).The code in the top left corner of each graph identifies the data series (Table 1).Fig. 7.Estimated flood and extreme rainfall event occurrence rate (black lines) with 95 % confidence intervals (grey area).The vertical ticks indicate the points in time when events occurred (POT time data).The code in the top left corner of each graph identifies the data series (Table 1).ig. 7.Estimated flood and extreme rainfall event occurrence rate (black lines) with 95 % confidence intervals grey area).The vertical ticks indicate the points in time when events occurred (POT time data).The code in he top left corner of each graph identifies the data series (Table 1).streamflow series and the daily rainfall series to obtain, respectively, estimated flood occurrence rates at the streamflow gauging stations, and, correspondingly, estimated extreme rainfall events occurrence rates at the rain gauging stations.
The results are presented in Fig. 7, which shows that the occurrence rates of extreme events in both mean daily flow and daily rainfall time series exhibit significant inter-annual variability.For example: in graph S2 of Fig. 7, the peak of λ(t) in the 1960s is significantly higher than the upper limit of the confidence band at 1990.
The results of Fig. 7 also show that there are some trends in the intensity of the inhomogeneous Poisson process that are exhibited consonantly among the analyzed mean daily streamflow samples, such as: (a) a peak in λ(t) in the early 1960s is visible in all the graphs with available data in those years, which indicates a higher frequency of extreme events; (b) in graphs S8, S9 there is a peak in the 1930s followed  The fact that such trends are visible in mean daily flow series from unregulated rivers that are geographically distributed (albeit unevenly) around the Portuguese territory (Fig. 1), suggests that those trends are not due to possible anthropogenic influences in the watersheds.The λ(t) estimates obtained for the rainfall data corroborate this hypothesis since, although the correspondence is not perfect, they exhibit some of the same trends that are visible in the flow data of the catchments under the influence of those particular rainfall gauging stations.The rainfall-flow influence relationships are visible in Fig. 1: the rainfall in station R1 has influence in streamflows at S6, R2 at S3, R3 at S4, and R4 at S2 and S5.As a detail of Fig. 7, Fig. 8 shows the λ(t) estimates for 2 of the previous rainfall-flow influence relationships without bootstrap confidence bands.
As mentioned in Sect.2.2, the methods were repeated for thresholds varying around the selected ones, and it was found that, while selecting higher thresholds lead to a lower λ, and vice-versa, the fluctuations of λ(t) around λ behave similarly.Figure 9 shows such fluctuations for threshold values of 6, 7 Winter NAO (NDJFM) index

S10
ber of floods in a hydrologic year, λ k plotted against the winter NAO (NDJFM) index of the same are identified with the sample number in the top left corner.
ples S1 to S10.Number of floods in a hydrologic year divided by the mean number, λ k /λ, of all ly flow samples plotted together against the winter NAO (NDJFM) index of year k.(the adopted ones), and 8 times the long term mean daily flow, or modular flow, q mod .

Relationship between flood occurrence and the North Atlantic Oscillation
The objective of this subsection is to ascertain whether or not the phase of the NAO has an influence on the occurrence rate of floods in the studied watersheds.By visually comparing Fig. 3 and Fig. 7 it is apparent that the peak in the estimated occurrences rate in the 1960s corresponds to a prolonged NAO negative phase in the same years.However, in the 1930s, where there are peaks in λ(t) in graphs S8 and S9 of Fig. 7, there is no similar negative phase in the winter NAO indices.
A quantitative analysis of the relationship between the winter NAO indices and the occurrence of floods was made on a discrete annual time basis.Figure 10 shows, for each streamflow sample, the number of floods per hydrologic year, λ k (discrete POT time data, Fig. 6), plotted against the winter NAO indices of the corresponding years.Such results suggest that years with positive NAO indices have a lower number of floods than years with negative NAO indices.Although that correlation does not seem to be particularly strong, Fig. 10 clearly shows that for every analyzed sample: (i) the majority of years without floods have positive NAO indices, and (ii) the years with the highest flood occurrence do not occur in positive NAO phases.
Figure 11 shows the same results as the previous analysis, but here λ k data from all the samples were made dimensionless by their respective mean values λ, and plotted together against the winter NAO indices.In that figure, box plots were drawn to represent the dispersion of the NAO indices for the values of λ k that are lower than the mean (λ k /λ<1); between one and two times the mean (1≤λ k /λ<2); and higher than two times the mean (λ k /λ≥2).The results clearly show that (i) years with a less-than-average number of floods tend to occur when the NAO is in a positive phase, and (ii) years with a higher number of floods (more that twice the average) tend to occur when the NAO is in a negative phase.
Although the relationship between the NAO and the occurrence of floods in Portuguese watersheds requires further investigation, the results of Figs. 10 and 11, and the comparison of Figs. 7 and 11, indicate that an increase in the rate of flood occurrence might be related to a decrease of NAO indices.That accordance is not strong enough to establish a cause-and-effect type of relationship, as one could expect, seeing that rainfall is strongly modulated by regional and local factors.However it does merit to be scoped in future research, namely on the possibility of establishing projections of the flood occurrence rates if reliable long-range forecasts of the NAO are made available (see Sutton and Allen, 1997;Rodwell et al., 1999)

Conclusions and future research
The current consensus on the effects of climate change on the hydrological cycle compromises the stationarity of hydrological time series.The objective of the research underlying this paper was to make an exploratory analysis on the variability of flood occurrence rates in 10 Portuguese watersheds, and to ascertain if that variability is concurrent with the principle of stationarity.
For that purpose, a peaks-over-threshold (POT) sampling technique was applied to 10 long series of mean daily streamflows and to 4 long series of daily rainfall in order to sample the times of occurrence (POT time data) of the peak values of those series.A preliminary analysis of the POT time data suggested that the occurrence rates of those events were nonuniform.The kernel occurrence rate estimator, coupled with a bootstrap approach, was applied to the POT time data to obtain the time dependent estimated occurrence rate curves, λ(t), of floods and extreme rainfall events.
The achieved results clearly show that, in the studied watersheds, the occurrence of floods constitutes an inhomogeneous Poisson process, hence the flood occurrence rates are nonstationary.There is a number of similarities in the behaviour of the estimated flood occurrence rates among the various samples, such as a peak in flood occurrence rates in the 1960s and the 1930s.Such peaks are also visible in the results of the application of the kernel technique to the daily rainfall data, from rain gauges located within 4 of the analyzed watersheds.Notwithstanding the unfeasibility of generalizing the results for the whole Portuguese territory, due to the small number and uneven spatial distribution of the samples, the similarities in the behaviour of λ(t) among different watersheds, and between rainfall and streamflow suggests that the observed trends are inherent to the natural inter-annual variation of the hydrological cycle, as opposed to potential anthropogenic influence on the catchments themselves.Furthermore, the observed trends in the estimates of λ(t) are robust against threshold selection in the POT sampling of extreme value data.
An attempt was made to assess whether or not the NAO index casted any influence on the flood occurrence rates in Portuguese watersheds.This was done by comparing the winter (NDJFM) NAO index with the number of floods within the hydrologic year, λ k .The results of that analysis show that, although the correlation is not particularly strong, the years with less floods tend to happen when the winter NAO is in its positive phase phase, and the years with a high flood occurrences tend to happen when the winter NAO is in a negative or neutral phase.This subject does, however, merit further research.
In conclusion, the results presented in this paper suggest that the mathematical formulation of the flood frequency models relying on stationarity, as those more commonly applied in Portugal, should be revised in view of accounting for possible nonstationarities in the occurrence rate of floods and extreme rainfall events.
Overall the obtained results lay the foundation for future research on nonstationary hydrological modelling of floods in Portuguese watersheds, such as a hybrid Poisson -extreme value distribution model which encompasses a nonparametric description of the time dependence via the inhomogeneous Poisson process and a parametric extreme value distribution.Having λ vary with time affects the flood quantile associated with a given annual exceedance probability, which is now allowed to change as a function of time, even if the peak exceedances over a threshold prove to be stationary.Such modelling of extreme hydrological phenomena, coupled with a potential dependence on climate patterns such as the NAO, could be of use to researchers and practitioners that deal with uncertainty in water resources systems planning and management, associated with hydrological extremes.
acknowledge the valuable comments and suggestions made by three anonymous referees and by Professor Félix Francés during the interactive discussion of this paper at HESSD.
Edited by: H. Cloke

Fig. 1 .
Fig. 1.Mainland Portugal.(a) Countour map of the mean annual rainfall.Source: Portuguese Environmental Agency (APA, http://sniamb.apambiente.pt/webatlas/),based on the climatological normals of reference, corresponding to the period 1931-1960, using a network of 334 stations(CNA, 1983).(b) Location of the rainfall and streamflow gauging stations, as identified by the codes presented in Table1.The shaded areas correspond to the catchments under study.

Fig. 2 .
Fig. 2. Peaks-over-threshold sampling technique applied to sample S5 -Fragas da Torre.(a) Mean daily flow sample; (b) autocorrelation function (ACF) of the exceedances with a 95 % confidence band; (c) mean number of exceedances in a year, (d) mean exceedance over the threshold with a 95 % confidence band; and(e) empirical distribution function of the times of the events.The selected threshold is identified by the dashed lines in panels (a), (c) and (d).In panel (e), the vertical ticks show the points in time associated with the flood occurrences, the dotted diagonal line shows the uniform distribution and the solid diagonal lines show significance for a Kolmogorov-Smirnov statistic at 5 % and 1 %. 18

Fig. 2 .
Fig. 2. Peaks-over-threshold sampling technique applied to sample S5 -Fragas da Torre.(a) Mean daily flow sample; (b) autocorrelation function (ACF) of the exceedances with a 95 % confidence band; (c) mean number of exceedances in a year, (d) mean exceedance over the threshold with a 95 % confidence band; and (e) empirical distribution function of the times of the events.The selected threshold is identified by the dashed lines in panels (a), (c) and (d).In panel (e), the vertical ticks show the points in time associated with the flood occurrences, the dotted diagonal line shows the uniform distribution and the solid diagonal lines show significance for a Kolmogorov-Smirnov statistic at 5 % and 1 %.

Fig. 3 .
Fig. 3. Winter NAO indices based on standardized sea-level pressures differences from November to March -NAO (NDJFM) -between Iceland and Gibraltar from 1900 to 2009 (the year corresponds to November of each NDJFM average).The solid black line is a smoothing LOWESS curve.

Fig. 4 .Fig. 4 .
Fig. 4. Estimated flood occurrence rate at S8 -Ponte Juncais, with (solid line) and without pseudodata generation; flood dates (represented by vertical lines) obtained from the POT time dat

Fig. 5 .Fig. 5 .Fig. 5 .
Fig. 5. Streamflow (top panel) and rainfall (bottom panel) data: running average exceedances above th over successive periods of (the previous) 15 years, and made dimensionless by the long term mean.

Fig. 10 .
Fig. 10.Number of floods in a hydrologic year, λ k plotted against the winter NAO (NDJFM) index of the same year.Graphs are identified with the sample number in the top left corner.

Fig. 11 .
Fig. 11.Samples S1 to S10.Number of floods in a hydrologic year divided by the mean number, λ k /λ, of all the mean daily flow samples plotted together against the winter NAO (NDJFM) index of year k.

Fig. 10 .
Fig. 10.Number of floods in a hydrologic year, λ k plotted against the winter NAO (NDJFM) index of the same year.Graphs are identified with the sample number in the top left corner.

Fig. 11 .
Fig. 11.Samples S1 to S10.Number of floods in a hydrologic year divided by the mean number, λ k /λ, of all the mean daily flow samples plotted together against the winter NAO (NDJFM) index of year k.

Table 1 .
Mean daily streamflow and daily rainfall data.Sample code, name and periods of records of the samples; mean annual flow depth, MAFD (mm), and area of the catchments; mean annual rainfall, MAR (mm); POT sampling threshold values, u (m 3 s −1 and mm, respectively).
Location of the rainfall and streamflow gauging stations, as identified by the codes presented in Table1.The shaded areas correspond to the catchments under study.
Jones et al. (1997)tion, andJones et al. (1997)used Gibraltar.In this work it was decided to use the Iceland-Gibraltar 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) and