Effect of climate change and variability on extreme rainfall intensity – frequency – duration relationships : a case study of Melbourne

The increased frequency and magnitude of extreme rainfall events due to anthropogenic climate change, and decadal and multi-decadal climate variability question the stationary climate assumption. The possible violation of stationarity in climate can cause erroneous estimation of design rainfalls derived from extreme rainfall frequency analysis. This may result in significant consequences for infrastructure and flood protection projects since design rainfalls are essential input for design of these projects. Therefore, there is a need to conduct frequency analysis of extreme rainfall events in the context of non-stationarity, when nonstationarity is present in extreme rainfall events. A methodology consisting of threshold selection, extreme rainfall data (peaks over threshold data) construction, trend and nonstationarity analysis, and stationary and non-stationary generalised Pareto distribution (GPD) models was developed in this paper to investigate trends and non-stationarity in extreme rainfall events, and potential impacts of climate change and variability on intensity–frequency–duration (IFD) relationships. The methodology developed was successfully implemented using rainfall data from an observation station in Melbourne (Australia) for storm durations ranging from 6 min to 72 h. Although statistically significant trends were detected in extreme rainfall data for storm durations of 30 min, 3 h and 48 h, statistical non-stationarity tests and non-stationary GPD models did not indicate non-stationarity for these storm durations and other storm durations. It was also found that the stationary GPD models were capable of fitting extreme rainfall data for all storm durations. Furthermore, the IFD analysis showed that urban flash flood producing hourly rainfall intensities have increased over time.


Introduction
Over the last 100 years, global surface temperature has increased approximately by 0.75 • C, and this warming cannot be explained by natural variability alone (IPCC, 2007).IPCC (2007) stated that excessive greenhouse gas emissions due to human activities are the main reason for current global warming.Increasing frequency and magnitude of extreme weather events is one of the main concerns caused by global warming.Increases in extreme rainfall frequency and magnitude have already been recorded in many regions of the world (Mueller and Pfister, 2011;Dourte et al., 2013;Bürger et al., 2014;Jena et al., 2014), even in some regions where the mean rainfall has shown decreasing trends (Tryhon and DeGaetano, 2011).Moreover, the magnitude and frequency of extreme rainfall events are very likely to increase in the future due to global warming (IPCC, 2007).
Increased frequency and magnitude of extreme rainfall events questions the stationary climate assumption (i.e. the statistical properties of the rainfall do not change over time), which is an underlying assumption of frequency analysis of extreme rainfalls.Khaliq et al. (2006) explained that the classical notions of probability of exceedance and return period are no longer valid under non-stationarity.The possible violation of stationarity in climate increases concerns amongst hydrologists and water resources engineers about the accuracy of design rainfalls, which are derived from frequency analysis of extreme rainfall events under the stationary climate assumption.Erroneous selection of design rainfalls can cause significant problems for water infrastructure projects and flood mitigation works, since the design rainfalls are an important input for design of these projects.Therefore, there is a need to conduct frequency analysis of extreme rainfall events under the context of the non-stationarity.Sugahara et al. (2009) carried out a frequency analysis of extreme daily rainfalls in the city of São Paulo using data over the period of 1933-2005.They considered nonstationarity in frequency analysis through introducing time dependency to the parameters of generalised Pareto distribution (GPD), which is one of the widely used distributions in frequency analysis of extreme values.Park et al. (2011) developed non-stationary generalised extreme value (GEV) distribution (another commonly used extreme value distribution) models for frequency analysis of extreme rainfalls in Korea considering non-stationarity similar to Sugahara et al. (2009).Tramblay et al. (2013) performed non-stationary heavy rainfall (it should be noted that "heavy" rainfall used here as same as "extreme" rainfall in Sugahara et al., 2009) analysis using daily rainfall data of the period 1958-2008 in France.They incorporated the climatic covariates into the generalised Pareto distribution parameters to consider nonstationarity.
There are very few studies, which investigated extreme rainfall frequency analysis in the context of non-stationarity in Australia.Jakob et al. (2011a, b) investigated the potential effects of climate change and variability on rainfall intensityfrequency-duration (IFD) relationships in Australia, considering possible non-stationarity of extreme rainfall data in design rainfall estimates.Yilmaz and Perera (2014) developed stationary and non-stationary GEV models using a single station in Melbourne considering data for storm durations ranging from 6 min to 72 h to construct IFD curves through frequency analysis.They investigated the advantages of nonstationary models over stationary ones using graphical tests.
In this paper, we aim to investigate extreme rainfall nonstationarity through trend analysis, non-stationarity tests and non-stationary GPD models (NSGPD).The extreme rainfall trend analysis was performed using data from a rainfall station in Melbourne considering storm durations of 6 and 30 min and 1, 2, 3, 6, 12, 24, 48 and 72 h.Trend analysis was used to determine whether the extreme rainfall series have a general increase or decrease over time.However, trends do not necessarily mean non-stationarity.The mean and variance of extreme rainfall data series may not change over time (i.e.stationarity), despite the presence of trends in extreme rainfall data series (Wang et al., 2006).Therefore, further analysis should be conducted to check whether the detected trends may correspond to extreme rainfall non-stationarity.Non-stationarity analysis of the extreme rainfall data was further carried out using statistical non-stationarity tests and NSGPD models in this study.
Potential effects of climate change on the IFD relationship were investigated through GPD models in this study following the stationarity analysis.Expected rainfall intensities for return periods of 2, 5, 10, 20, 50 and 100 years were derived and compared for two time slices : 1925-1966 (i.e. cooler period) and 1967-2010 (warmer period) after selecting 1967 as the change point based on the findings of Yilmaz and Perera (2014).Yilmaz and Perera (2014) conducted the change point analysis for extreme rainfall data for storm durations ranging from 6 min to 72 h in Melbourne, and stated that the year 1966 is the change point.Moreover, Jones (2012) listed the period 1910-1967 as stationary and 1968-2010 as nonstationary according to the observed minimum and maximum temperature and rainfall data in south-eastern Australia (which includes the Melbourne region).Therefore, the entire data set was divided into two periods (i.e. 1925-1966 and 1967-2010) and the IFD information was generated for the two periods to understand whether there are any changes in rainfall intensities between these cooler and warmer periods.
Changes in rainfall intensities (i.e.IFD information) over time can occur due to both climate change and natural climate modes (i.e.natural climate variability).The El Niño-Southern Oscillation (ENSO) with El Niño and La Niña phases (Verdon et al., 2004), the Indian Ocean Dipole (IOD) (Ashok et al., 2003), the Southern Annual Mode (SAM) (Meneghini et al., 2007), and the Interdecadal Pacific Oscillation (IPO) (Verdon-Kidd and Kiem, 2009) were expressed as significant climate modes, which have an influence on the precipitation variability in Victoria (Australia), which includes the Melbourne region.IPO affects the precipitation variability in Victoria itself; it also modulates the association between ENSO and Australian climate (Power et al., 1999;Kiem et al., 2003;Micevski et al., 2006).ENSO and Australian climate relationship was strong in particular during the IPO negative phases (i.e associated with wetter conditions).Moreover, Kiem et al. (2003) stated that La Niña events, which were increased during the negative IPO phases, are the primary driver for flood risk in Australia.It can be seen from the above studies that there is a need to investigate the IPO and extreme rainfall relationship due to its direct effects on Australian rainfall as well as effects of IPO on ENSO, which has a strong link to Australian rainfall.The effects of IPO on extreme rainfalls were investigated in this study through extreme rainfall IFD analysis during IPO negative and positive phases.Salinger (2005) and Dai (2013) defined time periods of IPO negative and positive phases as 1947-1976 and 1977-1998, respectively.Therefore, extreme rainfall IFD analysis was performed for these two periods to explain the relationship between IPO and extreme rainfalls.It should be noted that potential effects of climate change on design rainfall intensities (IFD information) were investigated through GPD models developed for 1925-1966 and 1967-2010 time periods, whereas the IPO and extreme rainfall relationship was investigated with GPD models for the periods of IPO negative (1947IPO negative ( -1976) ) and positive  phases.
As mentioned earlier in this section, there are very limited studies in the literature investigating IFD relationships in Australia considering non-stationarity of extreme rainfall data (e.g.Jakob et al., 2011a, b;Yilmaz and Perera, 2014).However, Jacob et al. (2011a, b) did not develop nonstationary extreme rainfall models to investigate their performances over stationary models, as is done in this study.Although Yilmaz and Perera (2014) developed non-stationary models for the same study area as in this study, they simply used annual maxima as extreme rainfall input to the stationary and non-stationary models.Several studies recommended the use of peaks-over-threshold (POT) data (derived by selecting values over a certain threshold) instead of annual maxima as extreme rainfall data input to frequency analysis (e.g.Re and Barros, 2009;Tramblay et al., 2013), since the POT approach results in larger data sets, leading to more accurate parameter estimations of extreme value distribution.Therefore, this study used POT data to develop stationary and non-stationary GPD models.

Study area and data
The city of Melbourne, Australia, was selected as the case study area.Data of the Melbourne Regional Office rainfall station (site no: 086071; lat 37.81 • S, long 144.97 • E) were provided by the Bureau of Meteorology in Australia.This station was selected for the study since it has long rainfall records, which are essential for trend and extreme rainfall IFD analysis.The approximate location of the station is shown in Fig. 1.
Six-minute pluviometer data are available from April 1873 to December 2010 at the Melbourne Regional Office station.These data were used to generate rainfall data for storm durations including 30 min and 1, 2, 3, 6 and 12 h.Furthermore, daily rainfall data starting from April 1855 are available at the Melbourne Regional Office station.Daily rainfall data were used to produce 48 and 72 h rainfall data.Although daily rainfall record is complete, there are missing periods in the 6 min data record.Missing periods in the 6 min data record were from January 1874 to July 1877 and from July 1914 to December 1924.Therefore, rainfall data over the period 1925-2010 from both sources (i.e. 6 min and daily) were used for all storm durations in this study.

Methodology
The methodology of this study consists of the following four steps.
1. Extreme rainfall data were constructed based on the POT approach after selection of suitable thresholds for storms of different durations.2. Trend analysis of POT data of all storm durations was carried out using non-parametric tests.Then, stationarity analysis was performed for the same data sets using statistical non-stationarity tests and non-stationary GPD models.
3. Stationary GPD models were developed, and design rainfall estimates were derived for standard return periods considering two time slices (1925-1966 and 1967-2010) in order to investigate potential effects of climate change on design rainfall intensities (extreme rainfall IFD information).
4. Stationary GPD models were constructed to obtain design rainfall intensities for IPO negative  and positive  phases to investigate the IPO and extreme rainfall relationship.

Threshold selection and extreme rainfall data set construction
The first step of the extreme rainfall frequency analysis is to construct the extreme rainfall data set.There are two widely used approaches to construct such data sets: block maxima and POT, also called the partial duration series approach (Thompson et al., 2009;Lang et al., 1999).In the block maxima approach, a sequence of maximum values is taken from blocks or periods of equal length, such as daily peak rainfall amount over an entire year or season.On the other hand, rainfall values that exceed a certain threshold are selected in the POT approach.Although the block maxima approach is the commonly used method due to its simplicity, it has a very important shortcoming in that it uses only one value from each block (Sugahara et al., 2009).This may cause the loss of some important information, and also smaller sample sizes, which affect the accuracy of the parameter estimates.Moreover, the POT method has the advantage of investigation of changes in number of events per year as well as magnitude (Jakob et al., 2011a).Due to the above-mentioned reasons, the POT approach is recommended for frequency analysis of extreme events (Re and Barros, 2009;Tramblay et al., 2013).It should be noted that "extreme rainfall data" and "POT data" terminology has been used interchangeably in the rest of the paper.Despite the above-mentioned advantages of the POT method over the block maxima approach, the POT approach is prone to producing dependent data.Data independency is an underlying requirement for use of extreme value distributions in frequency analysis.Therefore, the data dependency was removed in this study from the POT data of all storm durations through the method recommended by Jakob et al. (2011a).If there is a cluster of POT events, they recommended that the POT values 24 h prior to and after the peak rainfall event be removed from the data set.For example, if a peak rainfall value in a cluster of POT data is selected for 9 November 2013, rainfall values over the threshold on 8 and 10 November 2013 are not considered in the POT data set.None of the POT data sets (after the application of the method by Jakob et al., 2011a) showed dependency, even at the 0.1 significance level, according to the autocorrelation test as explained in Chiew and Siriwardena (2005).
The critical step in the construction of POT data is the selection of the appropriate threshold value.Researchers have proposed several procedures for selecting the thresholds, but a general and objective method is yet to emerge (Lang et al., 1999;Coles, 2001;Katz et al., 2005).The threshold selection task is a compromise between bias and variance.If the threshold is too low, the asymptotic arguments underlying the derivation of the GPD model are violated.On the other hand, too high a threshold will result in fewer excesses (i.e.rainfall values above threshold) to estimate the shape and scale parameter, leading to high variance.Therefore, in the threshold selection, it should be considered whether the limiting model provides a sufficiently good approximation compared with the variance of the parameter estimate (Coles, 2001;Katz et al., 2005).Beguería et al. (2011), Coles (2001), and Lang et al. (1999) recommended the mean residual plots to select the threshold.The mean residual plot indicates the relationship between mean excesses (i.e.mean of values above the threshold) and various thresholds.Mean excess is a linear function of threshold in GPD (Coles, 2001).Therefore, threshold value should be selected from the domain, where the mean residual plot shows linearity (i.e.linearity between mean excess and threshold) (Hu, 2013).The exact threshold value can be determined from the linear domain in such a way that on average 1.65-3.0extreme events per year are selected (e.g.Jakob et al., 2011a;Cunnane, 1973).This study adopted the mean residual plot method for selection of appropriate thresholds for all storm durations.

Trend and non-stationarity tests
Trend tests can be broadly grouped into two categories: parametric and non-parametric methods.Non-parametric tests are more appropriate for non-normally distributed and censored hydrometeorological time series data (Bouza-Deano et al., 2008).However, data independency is still a requirement of these tests.Mann-Kendall (MK) and Spearman's rho (SR) are non-parametric rank-based trend tests, which are commonly used for trend detection of hydrometeorological data (Yue et al., 2002).Formulation and details of the MK and SR tests can be found in Kundzewicz and Robson (2000).MK and SR tests were applied to POT data sets of all storm durations (6 and 30 min, and 1, 2, 3, 6, 12, 24, 48 and 72 h) over the period of 1925-2010 after applying the autocorrelation test as explained earlier.
Trend tests are used to determine whether the time series data has a general increase or decrease in trend.However, increasing or decreasing trends do not guarantee non-stationarity even if they are statistically significant.Therefore, it is useful to conduct further analysis in order to investigate non-stationarity of the data sets.In this study, three statistical tests, namely augmented Dickey-Fuller (ADF), Kwiatkowski-Phillips-Schmidt-Shin (KPSS) and Phillips-Perron (PP), were employed to investigate the non-stationarity in extreme rainfall data.These tests were selected due to their proven capability in hydrological studies (Wang et al., 2005;Wang et al., 2006;Yoo, 2007).Sen and Niedzielski (2010) and van Gelder et al. (2007) explain the details of these tests.Non-stationarity of data is the null hypothesis of ADF and PP tests, whereas the null hypothesis of the KPSS test is stationarity of the data series.Tests were performed at 0.05 significance level in this study.Whenever the significance level is higher than the p value (probability) of the test statistic, the null hypothesis is rejected.

Stationary GPD models
Several studies recommended the use of GPD for frequency analysis of POT data (e.g.Beguería et al., 2011).Therefore, GPD is used in this study to derive the extreme rainfall IFD relationships.GPD is a flexible, long-tailed distribution defined by shape (γ ) and scale (σ ) parameters.Equation (1) shows the cumulative distribution function of GPD.It should be noted that the stationary GPD model corresponds to conventional GPD models with constant shape and scale parameters.
The scale parameter (σ in Eq. 1) characterises the spread of distribution, whereas the shape parameter (γ in Eq. 1) characterises the tail features (Sugahara et al., 2009).Rainfall intensities in millimetres per hour for different return periods (2, 5, 10, 20, 50 and 100 years in this study) are calculated using the inverse cumulative distribution function.Details of GPD can be found in Sugahara et al. (2009), Coles (2001) and Rao and Hamed (2000).
There are different approaches such as maximum likelihood and L moments to estimate the parameters of GPD.In this study, the L moments method was used to estimate GPD parameters since it is less affected by data variability and outliers (Borujeni and Sulaiman, 2009).Hosking (1990) described the details of the L moments method.
Goodness of fit of the stationary GPD models was determined using the graphical diagnostics and statistical tests.The probability (P-P) and the quantile (Q-Q) plots are common diagnostic graphs.In P-P plot, the x axis shows empirical cumulative distribution function (CDF) values, whereas the y axis shows theoretical CDF values.In Q-Q plot, the x axis includes input (observed) data values, whereas the y axis shows the theoretical (fitted) distribution quantiles calculated by where F −1 (x) is inverse CDF, F n (x) is empirical CDF, and n is sample size.
Close distribution of the points of probability and quantile plots around the unit diagonal indicates a successful fit.Probability and quantile plots explain similar information; however, different pairs of data are used in probability and quantile plots.It is beneficial to use both plots to assess the goodness of fit, since one plot can show a very good fit while the other can show a poor fit.Coles ( 2001) explains the details of the diagnostic graphs.When probability and quantile plots show different results, statistical tests are useful to determine adequacy of the fit.
In addition to diagnostic graphs, Kolmogorov-Smirnov (KS), Anderson-Darling (AD), and chi-square (CS) statistical tests were used in this study to check the goodness of fit.These tests were used in the past hydrological applications of extreme value analysis (Laio, 2004;Salarpour et al., 2012).They are used to determine whether a sample comes from a hypothesised continuous distribution (GPD in this study).Null hypothesis (H 0 ) of the tests is "data follow the specified distribution".If the test statistic is larger than the critical value at the specified significance level, then the alternative hypothesis (H A ), which is "data do not follow GPD", is accepted (Yilmaz and Perera, 2014).Details of these tests can be seen in Di Baldassarre et al. ( 2009) and Salarpour et al. (2012).
As explained in Sect. 1, the extreme rainfall data of all storm durations were fitted to the stationary GPD models for 1925-1966 and 1967-2010 periods to investigate the climate change effects on IFD information, and for 1947-1976 (IPO negative phase) and 1977-1998 (IPO positive phase) to investigate the IPO and extreme rainfall relationship.

Non-stationary GPD (NSGPD) models
NSGPD models were used along with statistical nonstationarity tests in this study to identify whether the de-tected trends based on MK and SR tests correspond to nonstationarity.If it is proven that extreme rainfall data show non-stationarity over time, it is preferable to use NSGPD models instead of stationary GPD models.Non-stationary GPD models can be developed through the incorporation of non-stationarity features (i.e.time dependency or climate covariates) into the scale parameter of the stationary GPD model in Eq. ( 1) (Coles, 2001;Khaliq et al., 2006).Thus, the scale parameter is not constant and varies with time in non-stationary models.It is also possible to incorporate the non-stationarity into the shape parameter.However, it is very difficult to estimate the shape parameter of the extreme values distribution with precision when it is time dependent, and therefore it is not realistic to attempt to estimate the shape parameter as a smooth function of time (Coles, 2001).
In this study, two types of non-stationary GPD models were developed with parameters as explained below: In the above models, β 0 , β 1 and β 2 modify the scale parameters of NSGPD models.It should be noted that the exponential function has been adopted to introduce time dependency into the scale parameter to ensure its positivity.There are other functions which result in positive scale parameters; however an exponential function was used in this study, since it was recommended by some studies (e.g.Furrer et al., 2010) in the literature.NSGPD1 and NSGPD2 were applied to POT data of all storm durations in this study over the two periods (1925-1966 and 1967-2010).
The maximum likelihood method was used for parameter estimation of NSGPD models because of its suitability for incorporating non-stationary features into the distribution parameters as covariates (Sugahara et al., 2009).Shang et al. (2011) explain the details of the maximum likelihood method.
Superiority of the NSGPD models over the stationary GPD models were investigated through the deviance statistic test.Let M 0 and M 1 be the stationary and the non-stationary models, respectively, such that M 0 ⊂ M 1 .The deviance test is used to compare the superiority of M 1 over M 0 using the log-likelihood difference (D) using the following equation (Coles, 2001;El Adlouni et al., 2007): where l 1 (M 1 ) and l 0 (M 0 ) denote the maximised loglikelihood under models M 1 and M 0 , respectively.The test of the validity of one model against the other (in this case M 1 against M 0 ) is based on the probability distribution of D, which is approximated by chi-square distribution.For instance, consider comparing NSGPD1 model with three parameters (i.e.β 0 , β 1 , γ ) denoted by M 1 in Eq. ( 3) with stationary GPD model with two constant parameters (i.e.  3).Under the null hypothesis, the statistic D is approximately chi-square-distributed with 1 degree of freedom (the degree of freedom is decided based on difference between the number of parameters of M 0 and M 1 models).Stationary GPD model (M 0 ) should be rejected in favour of NSGPD1 (M 1 ) if D > c α , where c α is the (1 − α) quantile of a chi-square distribution at the significance level of α.Large values of D suggest that model M 1 explains substantially more of the variation in the data than M 0 .More detailed information about the deviance statistic test can be found in Coles ( 2001), Tramblay et al. (2013) and Beguería et al. (2011).Superiority of M 1 is evidence of non-stationarity of extreme rainfall data.In this case, NSGPD models should be used to generate rainfall intensity estimates.

Threshold selection
The thresholds for all storm durations were selected using the mean residual plots based on the linearity of data in these plots, as explained in Sect.3.1.A range of different threshold values in the linear domain of the mean residual plots were tested to select the final threshold so that the number of extreme rainfall events per year is in the range of 1.65-3.0events (Cunnane, 1973;Jakob et al., 2011a).For example, thresholds of 3.6 and 9.8 mm were selected for 6 min and 1 h storm durations, respectively, using the mean residual plots as shown in Fig. 2. Selected threshold values for all other storm durations are listed in Table 1.

Trend and non-stationarity tests results
Table 2 summarises the results of the trend tests.The trend tests (i.e.MK and SR) showed that extreme rainfall data of 30 min, 3 and 48 h exhibited statistically significant increasing trends at different significance levels.The 30 min data set showed the most significant data trend according to both MK and SR tests.It should be noted that only SR test indicated statistically significant trend for 48 h data set.Data sets of all other storm durations except 6 h also showed increasing trends; however these trends are not statistically significant even at the 0.1 significance level.
Trends in number of POT events per year were also investigated in this study.It was found that there is an increasing trend in the number of POT events for storm durations less than or equal to 2 h, whereas the number of POT events per year for storm durations greater than 2 h showed decreasing trends.However, none of these trends were statistically significant even at the 0.1 significance level.Furthermore, the ADF, KPSS and PP non-stationarity tests did not indicate non-stationarity in any of the extreme rainfall data sets.

NSGPD models
Non-stationary models (NSGPD1 and NSGPD2) for all storm durations were developed for 1925-1966 and 1967-2010 time slices.The deviance statistic test showed that there was no evidence that any of the non-stationary models outperformed their counterpart stationary models.For example, for the extreme rainfall data set of 3 h storm duration over period 1967-2010, the maximised log-likelihood of stationary GPD model (M 0 in Eq. 3) is 189.4,whereas the maximised log-likelihood values of NSGPD1 and NSGPD2 (M 1 in Eq. 3) are 189.3 and 189.1, respectively.D, calculated by Eq. ( 3), is smaller than c α for both non-stationary cases (NSGPD1 and NSGPD2).Therefore, it can be stated that non-stationary models do not outperform stationary models for these data sets.This is the case for all other storm durations (including the durations, in which extreme rainfall data showed statistically significant increasing trends) in both time periods (i.e. 1925-1966 and 1967-2010).As explained in Sect.4.2, the statistical non-stationarity tests (i.e.ADF, KPSS and PP) also showed that there was no evidence for non-stationarity of extreme rainfall data sets used in this study.Therefore, the stationary GPD models were used for the frequency analysis of extreme rainfall data sets to compare rainfall intensity estimates.

Stationary GPD models
POT data were used in stationary GPD models for two different pairs of periods, 1925-1966 and 1967-2010 (to investigate the effects of climate change), -IPO negative  and positive  phases (to investigate the IPO and extreme rainfall relationship), to compute IFD information under stationary conditions.This section explains the results of stationary GPD models over the periods of 1925-1966and 1967-2010, whereas Sect. 4.5 , whereas Sect. 4.5 shows the results of GPD models developed for the IPO analysis.
The graphical diagnostic and statistical tests showed that all extreme data sets (for all storm durations) were successfully fitted with the stationary GPD models.Figure 3 shows examples of the diagnostic graphs (i.e.probability and quantile plots) of stationary GPD models for the extreme rainfall data of 6 min, 3 h and 24 h storm durations over the 1925-1966 period.Table 3 indicates the results of the stationary GPD analysis (i.e.rainfall intensity estimates), whereas Fig. 4 illustrates the same information graphically for all storm durations.
Primary findings of the stationary GPD analysis are listed below: -Rainfall intensity estimates of the stationary GPD models over the period 1925-1966 were larger than those estimates of the period 1967-2010 for all storm durations equal or greater than 24 h (i.e. 24, 48 and 72 h) except 24 h storm duration of 2-year return period.
-For return periods less than or equal to 10 years, rainfall intensity estimates of sub-daily storm durations for the period of 1967-2010 were larger than those estimates of the 1925-1966 period.
-For the return periods above 10 years, the majority of hourly rainfall intensity estimates over the period 1967-2010 were larger than those estimates for the period of 1925-1966.
It is possible to conclude then that urban flash-floodproducing (i.e.flooding occurring in less than 6 h of rain; Hapuarachchi et al., 2011) hourly rainfall intensities have increased over time (i.e. from 1925-1966 to 1967-2010), with minor exceptions (i.e. 1 h storm durations of return periods above 10 years, and 2 h storm duration of 50-and 100-year return periods).It should be noted that 90 % confidence limits of rainfall intensity estimates were also calculated, but they are not shown in Fig. 4 in order to remove the clutter in the plots.

IPO analysis
The relationship between extreme rainfall data and IPO was investigated through IFD analysis for the periods of IPO negative (1947IPO negative ( -1976) ) and positive  phases.Results of the IPO analysis are shown in Table 4 and Fig. 5.
Results of the IPO analysis (Table 4 and Fig. 5) can be summarised as follows: -The rainfall intensities of storm durations equal to or greater than 24 h (24, 48 and 72 h) for all return periods during the IPO negative phase were larger than the corresponding rainfall intensities during the IPO positive phase.
-The rainfall intensities of all storm durations for the return periods greater than or equal to 20 years (i.e. 20, 50 and 100 years) during the IPO negative phase exhibited larger values relative to those rainfall intensities for the IPO positive phase, as can be seen Table 4 and Fig Table 3. Rainfall intensity (mm h −1 ) estimates derived from stationary GPD models over the periods 1925-1966 and 1967-2010.Duration/Return 2 years 5 years 10 years 20 years 50 years 100 years period 1925-1966 1967-2010 1925-1966 1967-2010 1925-1966 1967-2010 1925-1966 1967-2010 1925-1966 1967-2010 1925-1966 1967 -Rainfall intensities of storm durations below 3 h for the return periods less than or equal to 10 years (i.e. 2, 5 and 10 years) during the IPO negative phase were lower than those design rainfall intensities for the positive phase.This was also the case for the rainfall intensity estimates for storm durations between 3 and 12 h for return periods of 2 and 5 years.
In summary, increases in rainfall intensities were observed during the IPO negative phase for storms with long durations and high return periods, which is consistent with the literature (Kiem et al., 2003).In other words, the IPO negative phase can be the driver of higher rainfall intensities for long durations and high return periods.However, the trends in extreme rainfall data and differences in rainfall intensities for short storm durations and return periods cannot be explained with the IPO influence.In this study, only the relationship of IPO and extreme rainfall was investigated since the literature indicated IPO to be a very influential climate mode on extreme rainfall events in Victoria.However, there is a need to examine relationships between extreme rainfalls and other climate modes to correctly identify the primary driver for the extreme rainfall trends and differences in rainfall intensity estimates.Also, it is necessary to conduct similar analysis using data of other stations to assess the findings of this study.

Climate change and extreme rainfalls
Anthropogenic climate change may be the reason for the findings of this study (differences in rainfall intensity estimates over time and detected trends).Anthropogenic climate change can impact not only the extreme rainfalls directly, but also the dynamics of key climate modes.Climate change causes increases in intensity and frequency of extreme rainfalls, since the atmosphere can hold more water vapour in a warmer climate (Chu et al., 2014).The increase in rainfall extremes is larger than changes in mean rainfall in a warmer climate, because extreme precipitation relates to increases in moisture content of atmosphere (Kharin and Zwier, 2005).
Some studies (e.g.Murphy and Timbal, 2008;CSIRO, 2010) on rainfall changes in south-eastern Australia stated that, although there is no clear evidence to attribute rainfall change directly to the anthropogenic climate change, it still cannot be ignored.Rainfall changes are linked at least in part to the climate change in south-eastern Australia.Nevertheless, it is very difficult to attribute extreme rainfall trends and rainfall intensity differences to anthropogenic climate change due to the limited historical data records and strong effects of natural climate variability (Westra et al., 2010).Further analysis to investigate the reasons of the extreme rainfall trends and design rainfall intensity differences is beyond the scope of this paper.

Conclusions
A methodology consisting of threshold selection, extreme rainfall data (peaks over threshold data) construction, trend and non-stationarity tests, and stationary and non-stationary generalised Pareto distribution (GPD) models was developed in this paper to investigate the potential effects of climate change and variability on extreme rainfalls and intensityfrequency-duration (IFD) relationships.The methodology developed was successfully implemented using extreme rainfall data of a single observation station in Melbourne (Australia).The same methodology can be adopted for other stations in order to develop studies of larger spatial scale by analysing data of multiple stations.Major findings and conclusions of this study are as follows: -Statistically significant extreme rainfall (in millimetres) trends were detected for storm durations of 30 min, 3 h and 48 h, considering the data from 1925 to 2010.
-Statistically insignificant increasing trends in the number of POT events were found for storm durations less than or equal to 2 h, whereas statistically insignificant decreasing trends were detected in the number of POT events per year for storm durations greater than 2 h.
-Despite to the presence of trends in extreme rainfall data for above storm durations (i.e. 30 min, 3 h and 48 h), there was no evidence of non-stationarity according to statistical non-stationarity tests and non-stationary GPD models.The developed non-stationary GPD models did not show any advantage over the stationary models.
-The stationary GPD models were capable of fitting extreme rainfall data for all storm durations according to the graphical and statistical tests.
-Urban flash-flood-producing hourly rainfall intensities have increased within the time periods 1925-1966 and 1967-2010.-Analysis on relationship between the Interdecadal Pacific Oscillation (IPO) and extreme rainfalls showed that the IPO could be responsible for higher rainfall intensities for long durations and high return periods.On the other hand, the IPO cannot be shown as a driver for the trends in extreme rainfall data and differences in rainfall intensities for short storm durations and return periods.
It should be noted that this study used data from a single station to demonstrate the methodology for future studies.It is not realistic to extrapolate the findings of this study for larger spatial scales such as even the entire Melbourne metropolitan area without further analysis using rainfall data from multiple observation stations within the area.It is recommended that the methodology developed in this study be applied using data from multiple stations for larger spatial scales.It is also recommended that similar analysis of this study be conducted for future time periods using future rainfall data derived from climate models, since several studies highlighted very likely increases in intensity and frequency of extreme rainfalls in future.

Figure 1 .
Figure 1.Approximate location of the Melbourne Regional Office rainfall station.
et al.: Climate change and variability on rainfall IFD relationships γ ) denoted by M 0 in Eq. (

Figure 3 .
Figure 3. Goodness of fit for extreme rainfall data of (a) 6 min, (b) 3 h and (c) 24 h over the 1925-1966 period.

Figure 5 .
Figure 5. Graphical illustration of results of IPO negative-and positive-phase analysis.

Table 1 .
Threshold values obtained from the mean residual plot.