Climate change and its influence on design rainfall at-site in New South Wales State, Australia

Climate change has increasingly led to more frequent floods, droughts and heatwaves, and has affected many sectors including water resources. Understanding the future changes, particularly on regional scales, is important for the development of adaptation and mitigation measures (Homsi et al. 2020). Obviously, the increasing frequency and intensity of heavy rainfall is an expected result of climate change at many locations around the globe (Fotovatikhah et al. 2018; Shamshirband et al. 2020). Therefore, there has been considerable attention and research on the design rainfall based on the stationary and non-stationary frequency analysis approaches around the world (e.g. Akpan & Okoro 2013; Ball et al. 2016; Yilmaz et al. 2017; Hajani & Rahman 2018; Lawin et al. 2019). Design rainfall as the form of intensity–frequency–duration (IFD) rainfall data is defined as the rainfall depth associated with a given average return interval and duration (Al Mamoon et al. 2014). The establishment of IFD curve relationship was done as early as 1931 and 1932 (Sherman 1931; Bernard 1932). The rainfall IFD relationship is commonly required for the planning and design of various water resource projects (El-Sayed 2011).

The recommended IFD data for Australia are provided in Australian Rainfall and Runoff (ARR). The ARR1987 selecting Log Pearson Type 3 (LP3) distribution to derive IFD curves (I. E. Aust. 1987). In 2013, IFD curves have been estimated selecting generalized extreme value (GEV) distribution as part of the ARR by the Australian Bureau of Meteorology (BOM). The main assumption in estimating the IFD curves is that the rainfall data at a specific site are based on the stationary method (Ball et al. 2016; Hajani & Rahman 2018). This means that climatic changes have a negligible effect on the IFD curves. However, it has been proved by several recent studies that climatic change does exist, and it has considerable effect on rainfall data (Commonwealth Scientific and Industrial Research Organisation (CSIRO) 2014; BOM 2015). It is thus necessary to conduct further research to assess the impact of climate change on IFD relationship in Australia.

Many researchers have specifically addressed trends in annual, seasonal, monthly, daily and extreme rainfall data (Cai & Cowan 2009; Jakob et al. 2011; Li et al. 2012; Hajani et al. 2017). Most of these studies have found trends in rainfall data, the common findings being that the northwest had experienced an increase in rainfall over the last 50 years, while much of eastern Australia and the far southwest have experienced a decrease. The effect of climate change on rainfall data in Australia has been documented; for example, Yilmaz & Perera (2013) investigated heavy rainfall in the IFD relationship adopting the rainfall data from the Melbourne Regional Office station located in Melbourne, Australia. Their results showed that the IFD relationships vary over the time period and increase in sub-hourly design rainfall intensities. In another study, Bates et al. (2015) showed that in eastern Australia, rainfall frequency estimates derived assuming stationarity lie within the uncertainty bounds derived using the non-stationary method.

The objective of this study is to evaluate the possible impact of climate change on design rainfall using stationary and non-stationary GEV distributions. In addition, the IFD curves obtained in the current study (i.e. a particular station depends on stationary and non-stationary methods) have been compared with the Australian national IFD values (i.e. regional IFD curves data depend on the stationary method). It is expected that the outcomes of this study will provide valuable information regarding the design rainfall at-site location under change of a climate regime to conduct the frequency analysis of extreme rainfall events when non-stationarity is present in extreme rainfall events.

This study selects New South Wales (NSW) State in Australia as the study area. The NSW state has a very variable climate. The north-east of the state is dominated by summer rainfall and relatively dry winters, while the south of the state is dominated by regular rainfall during the winter (NSW Government 2016). The mean annual rainfall in NSW varies between 150 and 500 mm. Among all the Australian states, NSW has the highest average annual flood damage cost of about $150 million (BOM 2016). According to the Köppen climate classification, the major part of NSW, west of the Great Dividing Range, has an arid to semi-arid climate (EPA 1997).

Initially, 10 pluviograph stations with 6-min temporal resolution rainfall data were collected over eastern NSW as shown in Figure 1. The extract annual maximum (AM) rainfalls of 6-min duration were then aggregated into six sub-hourly durations (6, 12, 18, 24, 30 and 48 min), six sub-daily durations (1, 2, 3, 6, 8 and 12 h) and three daily durations (1, 2 and 3 days) for each of the selected pluviograph stations. There were numerous missing periods in the rainfall data, totalling 5% gaps in the AM rainfall series. It should be noted that the threshold record length of 55 years (i.e. study period of 1960–2015) was arbitrarily selected, and this should be long enough to characterize changes in rainfall. The gaps in the rainfall data of a specific station were filled by regression analysis with a nearby station where there was no data gap.

Finally, data from Sydney-066062 (Observatory Hill) pluviograph station as shown in Figure 1 were retained. Sydney-066062 station was chosen for estimating the stationary and non-stationary at-site IFD curves because the AM rainfall data exhibit statistically significant trends at the 10 and 5% significance levels for most of the selected rainfall durations. The rainfall data have been obtained from the database of the Australian BOM (BOM 2016). The primary rainfall data were recorded at 6-min intervals. A computer code in FORTRAN language was developed to extract AM rainfall events of six sub-hourly durations (6, 12, 18, 24, 30 and 48 min), six sub-daily durations (1, 2, 3, 6, 8 and 12 h) and three daily durations (1, 2 and 3 days). The standard regional IFD curve data in ARR for Sydney-066062 pluviograph station (i.e. ARR87 and ARR13 IFD curve data) were obtained from the BOM IFD calculator (BOM 2016).

The current IFD data in Australia (as part of its national IFD data, BOM 2013) are based on a stationarity assumption, which assumes that the climate is stationary. However, as noted in the section ‘Introduction’, numerous studies have found that rainfall has been changing. Many rainfall stations in Australia show either a positive or negative trend in their AM rainfall event data. In this study, the Mann–Kendall (MK) test (Mann 1945; Kendall 1975) has been adopted to analyse the trends in AM rainfall data for the 10 pluviograph stations shown in Figure 1 at 1, 5 and 10% significance levels. Also, three goodness-of-fit tests (i.e. Kolmogorov–Smirnov, Anderson–Darling (AD) and Chi-Square tests) have been adopted to assess the goodness of fit of the GEV distribution.

To tackle the non-stationarity, which is identified in the AM rainfall data, the non-stationary GEV distribution has been used. The reason for employing GEV distribution is that it has been widely adopted in rainfall frequency analysis in previous studies (e.g. Elsebaie 2011; Haddad & Rahman 2014; Ball et al. 2016). Also, this distribution is often fitted to environmental time-series of extreme values such as AM rainfall data. The distributional parameters in GEV distribution, such as the location and scale, are varied as a function of time or with relevant covariates.

The maximum-likelihood estimation (MLE) method has been used in this study to estimate the parameters of the GEV distribution for both the stationary and non-stationary cases. This approach has been gaining popularity for IFD estimation; hence, it has been applied here to the station exhibiting significant trends found by the MK test (Mann 1945; Kendall 1975). A brief description of each of these statistical techniques used in this paper is presented below.

Many studies require long-term rainfall data; therefore, a test must be conducted to check the consistency of the rainfall record. This is necessary because over a time period, it may happen that there will be some obstructions (trees and buildings) that may have emerged after the installation of a gauge or station location might have changed, or the observational procedure might have changed (Singh 1994). The inconsistency of the rainfall record can be checked by one of the most common and widely accepted methods (i.e. the double mass curve method).

The theory of the double mass curve method was proposed by Merriam (1937). It assumes that the graph of the mean cumulation of one quantity of the individual station against the mean cumulation of another quantity for a group of stations during the same period will show a straight line as long as the data are proportional. The slope of the line will represent the constant of proportionality between the quantities. A break in the slope of the double mass curve means that a change in the constant of proportionality between the two variables has occurred or perhaps that the proportionality is not a constant at all rates of cumulation. Two methods of adjustments are presented below:

• 1.

  A plot of accumulated AM rainfall at a station (being tested) vs. the accumulated values of the AM rainfall of the group of stations is prepared.

• 2.

  A break in the slope of this plot indicates a change in the AM rainfall of the station (being tested).

In detecting trends, the non-parametric methods, such as MK and Spearman's Rho (SR) tests, are widely used (e.g. Laz et al. 2014; Sadeghi & Hazbavi 2015; Hajani et al. 2017). A number of researchers have also adopted parametric tests (e.g. linear regression and Sen's slope) in detecting trends in rainfall data (e.g. Reza et al. 2011; Loveridge & Rahman 2013; She et al. 2016). The advantages of the non-parametric tests are that they are robust with respect to non-normality, nonlinearity, missing values, serial dependency and outliers (Villarini et al. 2011; Ishak et al. 2013).

The null hypothesis is rejected at a significance level (α) if ⁠, where is the value of the standard normal distribution with an exceedance probability of α/2.

In this study, three goodness-of-fit tests are used to assess how well a given distribution fits the rainfall data series of a given duration. The detail of each test is provided below.

The GEV distribution has been adopted in many applications in hydrology. In the United Kingdom, the GEV was recommended for at-site flood frequency analysis (Natural Environment Research Council (NERC) 1975) and in the United States, for rainfall frequency dataset (Willeke et al. 1995). In practice, it is widely used in modelling a wide variety of natural extremes, including maxima of rainfall, floods, winds, temperatures, and others. The GEV distribution is very attractive because its inverse has a closed form, and parameters are easily estimated (Martins & Stedinger 2000). The GEV distribution has three parameters (location: μ, scale: σ and shape: ξ), and it is used as an approximation to model the maxima of long sequences of random variables. The GEV distribution is equivalent to a Gumbel (ξ = 0), Fréchet (ξ > 0) or Weibull (ξ < 0) distribution depending on the shape parameter. The non-stationary frequency analysis based on the GEV distribution was undertaken on the case of dependency on time (t). The non-stationary behaviour is considered in the μ and σ parameters of the GEV distribution, where the ξ parameter is treated as a constant (i.e. does not change with time) as recommended previously (Cannon 2010).

In this study, four models have been employed with the emphasis on comparing the three non-stationary models dependent on time (t) to the stationary case in order to select the best model.

The models are:

• 1.

  Stationary model: the μ, σ and ξ parameters are constants, i.e. parameters being independent of the time (t).

• 2.

  Non-stationary model: the μ parameter linearly dependent on the time (t)

• 3.

  Non-stationary model: the σ parameter linearly dependent on the time (t)

• 4.

  Non-stationary model: the μ and σ parameters linearly dependent on the time (t)

Thus, the model with lower values of AIC and BIC was considered to be the ‘optimum model’ in representing the AM rainfall data series.

The results are presented in the following six sub-sections. The first five sub-sections focus on checking the consistency and the trends in AM rainfall intensity data at the local level, parameters of stationary and non-stationary IFD curves and estimating stationary and non-stationary IFD curves at the selected station. The last sub-section focuses on estimating IFD curves for different future periods and explores differences between stationary and non-stationary IFD curves and with ARR stationary IFD curves.

To test for the consistency in rainfall data for the study period of 1966–2015, the double mass curve was adopted (described in the section ‘Double mass curve method’). Figure 2 presents the cumulative plot of the double mass curve of AM rainfall for the station being tested (i.e. Sydney-066062) vs. the cumulative plot of the double mass curve of AM rainfall provided by the nine rainfall stations (shown in Figure 1) accompanied by the best linear fit. It was found that the break in the slope of the curve occurs around 1988. The change in slope in Figure 2 shows the effect of climate change on the rainfall data. The outcome of analysis is that correction factor is 1.084. Also, it is found that the average AM rainfall before and after applying the double mass curve for Sydney-066062 station is 9.21 and 9.73 mm, respectively.

From the results of the MK statistics test shown in Table 1, it can be observed that climate change is influencing nearly all extreme rainfall events for the 10 rainfall stations (shown in Figure 1). For some specific rainfall durations, the rainfall data exhibit statistically significant trends, which appear to be linked to rainfall extremes increasing/decreasing over these durations. Some regional changes in Australian rainfall have been linked to human-induced climate change and those from climate models forced with increasing greenhouse gases.

As can be observed in Table 1, the MK test for Sydney-066062 station showed that extreme rainfall data exhibit statistically significant increasing/decreasing trends at the 10 and 5% significance levels for most of the selected rainfall durations (here, 30 min to 1 day). It can also be observed in Table 1 that trends are not significant even at the 0.1 significance level for any of the 15 rainfall durations used in this study. Therefore, the presence of a statistically significant increasing/decreasing trend in AM rainfall data violates the basic assumption of stationary IDF curves. The rainfall data for Sydney-066062 station have been adopted for analysing the stationary and non-stationary IFD curves. In addition, Figure 3 shows that among the 10 selected stations, Sydney-066062 station has a maximum value of the average rainfall of 9.42 mm. The AM rainfall for Sydney-066062 station was found to range between 4.50 mm for 1983 and 18.60 mm for 1966. Meanwhile the average AM rainfall values for all the 10 stations were found equal to 8.68 mm for the study period. Among all the 10 stations, it was found that the highest maximum rainfall of 32.11 mm occurred on 1 June 1978, at Colo-061211 station, while the lowest rainfall value of 0.80 mm occurred on 13 July 1978 at Wallacia-067029 station.

Depending on the AM rainfall data series of the selected station (i.e. Sydney-066062) for the 15 rainfall durations, a total of 15 datasets (i.e. 1 × 15) were available for the goodness-of-fit testing. The results of goodness-of-fit tests, in Table 2, show that out of the 15 rainfall durations, the three goodness-of-fit tests (i.e. KS and AD and χ² tests) have not rejected GEV distribution in any of the 15 cases at the three significance levels (10, 5 and 1%).

The four models (i.e. Model 1 for the stationary IFD case and Models 2–4 for the non-stationary IFD case) were applied to AM rainfall data series of the selected pluviograph station for each of the 15 durations. This station was selected as it showed evidence of non-stationarity, i.e. significant trends in AM rainfall data series as discussed in the section ‘Consistency check for rainfall data’. The non-stationary behaviour was considered in the μ and σ parameters of the GEV distribution, whereas the ξ parameter was assumed to be stationary/constant. The use of the GEV models with covariates allowed examination of the changes in rainfall quantiles in future.

The parameter (ξ) which estimates for the selected station based on the MLE approach over the 55-year study period for 15 rainfall durations is presented in Table 2. The results in Table 2 show that the estimated values of the shape parameters were positive for the four GEV models (i.e. Model 1 for the stationary IFD case and Models 2–4 for the non-stationary IFD case). Furthermore, the evaluation of the standard error values of the shape parameters did not show any swinging between positive and negative values, which indicates that the AM rainfall data series would likely follow the Fréchet family of the GEV distribution.

The difference between the GEV models was demonstrated by comparing the estimated AIC and BIC values for each of the nested models with those of the stationary model. To achieve this, the minimized nllh functions for the candidate models were computed, and the results are presented in Table 3. For each candidate model, Table 4 lists the estimated AIC and BIC values along with the number of parameters in individual models. As shown in Table 5, the non-stationary models (i.e. Models 2–4) produce nearly similar values of AIC, which implies that the models with the linear trend in the location parameter and in the log-transformed scale parameters are nearly similar to the stationary model (i.e. Model 1). Depending on the results of AIC and BIC tests, the non-stationary models of the GEV distribution show a minor difference in fitting the AM rainfall data series.

The second-degree polynomial based on Equation (14) has been used to derive the relationship of the stationary IFD curves for the seven rainfall quantiles (1, 2, 5, 10, 20, 50 and 100 years) for the selected station. For stationary IFD curves, Figure 4 shows R² values for the fitted second-degree polynomial based on the maximum-likelihood parameter of the IFD curves for all the seven return periods. The R² values range between 0.978 and 0.999 (average = 0.992). These R² values are close to 1.00 and hence represent an acceptable fit.

For the selected station, the 95% lower and 95% upper confidence limits (CLs) of the rainfall quantiles for seven adopted return periods (i.e. 1, 2, 5, 10, 20, 50 and 100 years) are shown in Figure 5. The CLs presented in Figure 5 show that the uncertainty in rainfall estimates decreases with increasing return periods (except 2- and 100-year return periods). It should be noted here that the CL is mainly governed by the standard error values, which is inversely proportional to the sample size, i.e. to reduce the CLs, data length should be higher. Figure 5 shows that uncertainty in IFD data is not significant for the selected station. In general, it can be observed that the uncertainties for hourly and daily durations are higher than sub-hourly durations. The results of uncertainty in rainfall data in the current study are different from those of Yilmaz & Perera (2013), who investigated heavy rainfall trends and IFD relationship for a single rainfall station in Melbourne in Australia (covering the period of 1925–2010) and found that uncertainties for short durations (i.e. sub-hourly durations) were higher than long durations (i.e. hourly and daily durations). However, they found that there was no evidence to state that increase in design rainfall intensity was statistically significant. On the other hand, the results are similar to those of Jakob et al. (2011) who found a high degree of uncertainty at sub-daily durations in the south-eastern part of Australia.

Using the three non-stationary models (i.e. Models 2, 3 and 4), which have been described in the section ‘Methods’, the effect of climate change on IFD curves is investigated. The trend considered in the models is a simple one, i.e. as a function of time only. However, in the current study, all four models are adopted to estimate the stationary and non-stationary IFD quintiles.

Once the parameters of the IFD curves are estimated, the non-stationarity assumption is tested using three models (i.e. Models 2, 3 and 4) as discussed in the section ‘Methods’. For estimating the future IFD curves, by taking time as 2030, Figure 6 shows that when the μ parameter is a function of time (Model 2), a more consistent sets of IFD curves are achieved compared with Model 3 (standard deviation parameter is a functions of time) and Model 4 (both mean parameter and standard deviation parameter are function of time). For non-stationary IFD curves, Figure 6 shows that there is a large difference between the stationary and non-stationary IFD curves for long duration compared with short duration, and this confirms the outcome of the MK test. Figure 6 shows R² values of the fitted second-degree polynomial IFD curves for all the seven return periods. The R² values range between 0.979 and 0.999 (average = 0.993) for the selected station.

Depending on Models 2, 3 and 4 in the section ‘Non-stationary GEV distribution’, three future periods 2030, 2060 and 2100 have been used to simulate future changes in IFD curves and compared with those of the stationary methods depending on Model 1 in the section ‘Non-stationary GEV distribution’. Figures 7 and 8 show the stationary and non-stationary IFD curves for both 5- and 100-year return periods. For the selected station, generally Figures 7 and 8 show that a decrease in design rainfall depths is found for all three models. It should be noted here that the trend test for the selected station showed significant negative trends for most of the durations. Also, Figures 7 and 8 show that when the μ parameter is a function of time (Model 2), a more consistent sets of IFD curves are achieved compared with Model 3 (σ parameter is a function of time) and Model 4 (both μ and σ parameters are functions of time).

The new IFD curves at-site estimated in the current study based on both the stationary and non-stationary methods were compared with the standard regional IFD curves (i.e. ARR87 and ARR13 IFD curves), recommended by ARR, (I. E. Aust. 1987; Ball et al. 2016) to examine the expected degree of variation between the at-site and regional IFD data. The comparison is illustrated for two return periods (5 and 100 years) and three years (2030, 2060 and 2100) as shown in Figures 9 and 10. Since Model 2 (when μ parameter is a function of time) seems to be the best non-stationary model (as shown in Figures 7 and 8), comparison is done between the results of Model 2 and stationary IFD curves (i.e. Model 1) with ARR IFD curves.

In the stationary method, Figures 9 and 10 show that there is a greater difference between the new IFD curves and both ARR87 and ARR13 IFD curves for the 5-year return period compared with the 100-year return period. In addition, the results show that the ARR87 and ARR13 IFD curves are generally higher than the at-site IFD curves derived in the current study. The median difference between at-site IFD curves and regional ARR-recommended IFD curves is in the range of 11–18%. It should be noted that ARR87 adopted LP3 distribution to derive IFD curves and ARR13 adopted GEV distribution; however, GEV distribution has been adopted in the current study.

In the non-stationary method, it can be observed in Table 6 that the absolute median RE values in percentage (defined by Equation (15)) for the 5- and 100-year return periods range between 9 and 142% with ARR87 IFD data, and between 10 and 83% with ARR13 IFD data. In general, the highest absolute median RE values are found in Models 3 and 4, while the lowest absolute median RE values are found in Model 2.

The results of the non-stationary IFD curves of the current study confirm the results of the MK test. In addition, Table 7 shows the absolute median RE values (%) between the stationary and non-stationary IFD curves at three future times (i.e. 2030, 2060 and 2100). It can be observed in Table 7 that the absolute median RE values range between 1 and 44%, with a mean RE value of 8%, and with a standard deviation of RE values of 7% (based on year 2030). Also, the results in Table 7 show that as moving further in the future (in year), the absolute median RE values generally increase. For example, the median RE values range between 1 and 15% (based on future period 2030), while the absolute median RE values range between 4 and 25% (based on future period 2060) and the absolute median RE values range between 5 and 44% (based on future period 2100). Moreover, over the three future periods, the average highest absolute median RE value is found between the stationary and non-stationary approaches for Model 4, while the average lowest absolute median RE value is found for Model 2. Furthermore, as the return period (in years) increases, the median RE values increase between the non-stationary and stationary GEV models.

In Australia, the current IFD curves as recommended by the ARR are based on the assumption that AM rainfall data series are stationary, i.e. the IFD curves were derived assuming that any trends in rainfall data would have a negligible effect on the IFD curves. Moreover, the ARR IFD curves (i.e. ARR87 and ARR13 IFD curves) have been developed based on the regional frequency analysis technique to a rainfall dataset consisting of many stations in the region. The main challenge in the evaluation of the impact of climate change on design rainfall is that large uncertainty exists in the results. Therefore, in the current study, the stationary and non-stationary IFD curves have been developed using a particular case in NSW state in Australia to examine the expected degree of variation between the at-site and the standard regional AAR IFD curves.

From 10 candidate rainfall stations in NSW, Sydney-066062 pluviograph station was selected as it showed significant trends. The possible impact of climate change on design rainfall was investigated using stationary and non-stationary GEV distribution. The MLE method was used to estimate the parameters of the GEV model. In the adopted approach, the parameters of the GEV model such as the location and scale parameters were varied as a function of time or with relevant covariates representing climate change. The IFD curves resulting from both stationary and non-stationary approaches were compared with both ARR87 and ARR13 IFD curves.

Based on the three goodness-of-fit tests (i.e. KS test, AD test and Chi-Square test), it was found that GEV distribution fits the annual maximum rainfall data for the selected station. Also, it was found that the consideration of the location parameter of the GEV distribution is dependent on climatic change, while both scale and shape parameters are treated as a constant, giving more reliable results for estimating non-stationary IFD curves.

Also, it was found that when the future period increases (i.e. 2030, 2060 and 2100), the relative difference between the stationary and non-stationary IFD curves generally increased. A decrease in design IFD curves has been found for all the non-stationary models adopted in the current study. Under the stationary condition, model parameters used for simulations are basically tuned for historic data, thus being insensitive to future climate change. Due to the strong dependency of non-stationary model parameters on climate conditions, parameters are misrepresented in the stationary condition resulting in unrealistically low future rainfall data.

The comparison of the new IFD curves based on the stationary and non-stationary approaches with ARR IFD curves illustrates that there is a better match between the ARR IFD curves and the new stationary IFD curves compared with the non-stationary IFD curves. The methodology of the current study should be extended to more rainfall stations across Australia to have a better understanding of the possible differences in IFD data under the effect of climate change at-site and regional analyses. In addition, with the same distribution (i.e. GEV distribution) selected in this study or by adopting different distribution such as LP3, a similar study can be conducted in the future to evaluate the effects of climate variability such as El Niño Southern Oscillation, Interdecadal Pacific Ocean and Southern Annual Mode in Australia's rainfall data.

The author acknowledge the Australian Bureau of Meteorology for providing the pluviograph data.

Data cannot be made publicly available; readers should contact the corresponding author for details.