Use of observed scaled daily storm profiles in a copula based rainfall disaggregation model

Advances in Water Resources 45 (2012) 26–36

Contents lists available at SciVerse ScienceDirect

Advances in Water Resources journal homepage: www.elsevier.com/locate/advwatres

Use of observed scaled daily storm profiles in a copula based rainfall disaggregation model Yeboah Gyasi-Agyei ⇑ Centre for Railway Engineering, Faculty of Sciences, Engineering & Health, CQUniversity, Rockhampton QLD 4702, Australia

a r t i c l e

i n f o

Article history: Available online 4 December 2011 Keywords: Stochastic Rainfall Disaggregation Copula Storm profile Daily

a b s t r a c t It has been demonstrated in this paper that observed scaled daily storm profiles (OSDSPs) can be used at a different climatic site to disaggregate daily rainfall data into fine timescales by following some selection rules. For a given daily rainfall depth, the normal copula model predicts the total daily wet periods’ duration, capturing the seasonal variation of the copula and marginal lognormal model parameters by a first harmonic Fourier series. An OSDSP type, defined according to whether the first and/or last periods of the day are wet, is selected based on whether the day is an isolated wet day or in a cluster of wet days. For two consecutive wet Julian days, a transition probability is used to ascertain whether rainfall will cross midnight and spread to the next Julian day. A pool of the selected OSDSP type having the same, or close, numerical value of the copula predicted total daily wet periods’ duration is created. From the pool is selected one OSDSP that has the closest value of lag-1 autocorrelation of the natural logarithm of the wet period depths to a simulated value conditioned on the total daily wet periods’ duration. Multiplication of the OSDSP’s period depths by the daily rainfall depth gives the complete daily storm profile. The applicability of the OSDSP based disaggregation model is demonstrated by using two case study sites and OSDSPs from Australian State capital cities situated in different climates. Using the capital cities’ OSDSPs independently, and also together, predicted quite well the monthly gross rainfall statistics down to 6 min. Ó 2011 Published by Elsevier Ltd.

1. Introduction 1.1. Background In many countries, including the USA and Australia, long records of daily rainfall for over 100 years of a dense network are available. The Australian SILO Data Drill (www.longpaddock.qld.gov.au/silo/) and the Australian Water Availability Project (www.bom.gov.au/ jsp/awap) have established facilities for extracting daily rainfall data, at 5 km spatial resolution over the entire continent dating back to 1889, from an archive of interpolated rainfall and climate surfaces. A similar archive of rainfall data is managed by the USA’s National Climatic Data Centre. A major advantage of daily rainfall data is the preservation of seasonal and annual trends and easy identification of climate change phenomena. However, environmental processes such as erosion and pollutant transport, and the planning, design, and management of hydraulic structures in urban areas require very fine timescale rainfall data. Fine timescale rainfall data are, however, limited worldwide as a result of the high cost and low reliability associated with monitoring such data. This paper contributes to disaggregation modelling of the long-term daily rain⇑ Tel.: +61 7 49309977; fax: +61 7 49306984. E-mail address: [email protected] 0309-1708/$ - see front matter Ó 2011 Published by Elsevier Ltd. doi:10.1016/j.advwatres.2011.11.003

fall to a very fine timescale, such as 6 min. Continuous rainfall–runoff modelling will be facilitated by the availability of such fine timescale rainfall data. Rainfall data disaggregation has a long history starting with simplified approaches (e.g., [4,15]) that were based on distribution functions of the storm event properties, namely starting times, volume, intensity and duration. Recent approaches have been based on variants of the rectangular pulse models (e.g., [2,7,9,11,18,13]) and fractal and multifractal models (e.g., [8,19]). Serinaldi [26] has recently provided an extensive review of the later group of disaggregation models. A daily rainfall storm profile is multidimensional with properties of start time, duration (or total wet periods’ duration), depth, peak intensity and time to peak. These properties are stochastic in nature and have different values for each wet day. Copulas have the ability to join multivariate distribution functions to their one-dimensional marginal distribution function based on the Sklar theorem (e.g., [16]). The pair-wise dependence of these variables has traditionally been modelled using bivariate distributions such as the extreme-value distribution, with the drawback being the assumption that the variables follow the same parametric univariate distributions [6]. This limitation is avoided by copula models where the dependence between the variables is independent of the marginal distributions which could be different for the variables. While copula models

Y. Gyasi-Agyei / Advances in Water Resources 45 (2012) 26–36

have been available for years, particularly in the field of actuarial and financial applications, they are now making fast progress into hydrology. Notable pioneering applications of copulas in hydrology include De Michele and Salvadori [3], Salvadori and De Michele [22] and Favre et al. [5] with the initial focus on frequency analysis. Since then copulas have been applied in diverse areas of hydrology. The International Commission on Statistical Hydrology of the International Association of Hydrological Sciences (ICSH-IAHS), ‘‘Reference collection initiative on Copula function’’ group keeps record of almost all papers on copula applications in hydrology, indicating over 30 papers per year in international journals since 2007 (www.stahy.org). Serinaldi [25] and Bardossy and Pegram [1] demonstrated the use of copula models for multisite daily rainfall modelling. Zhang and Singh [29,30] have derived rainfall intensity–frequency–duration curves using copula models. A stochastic design storm generator based on copulas and mass curves was presented by Vandenberghe et al. [27]. Gyasi-Agyei [14] developed a copula based continuous daily rainfall disaggregation model which uses a vine pair copula construction to model the dependence structure between total daily rainfall depth, total wet periods’ duration and the maximum proportional depth of a wet period. This paper is a further development of his rainfall disaggregation model as explained in the section. 1.2. The basic concepts of the disaggregation model This paper attempts to greatly simplify the modelling process presented by Gyasi-Agyei [14] by firstly modelling only the dependence structure between daily rainfall depth and the total wet periods’ duration. Secondly, some rules are followed to select an observed scaled daily storm profile (OSDSP) from different climatic regions, thus eliminating the need for the binary wet–dry sequence and the AR(1) depth process sub-models. In this way the natural rainfall process is mimicked as much as possible and the number of model parameters is considerably reduced. Fig. 1 shows a definition sketch of an observed daily storm profile consisting of 3 bursts of durations di (hours) separated by dry periods. For the sampling timescale ts (hours) there are nL (=10 in the sketch) wet periods of depths ri (mm) summing up to the total daily rainfall depth R (mm). The scaled depths, or wet period proPL portional depths wi are defined as wi = ri/R, noting that ni¼1 wi ¼ 1, and the distribution of the wis within the 24-h period gives the OSDSP. The total daily wet periods’ duration L (hours) is the sum of the burst durations (d1 + d2 + d3) which is equal to nL ⁄ ts. Variables R and L are treated as random and their dependence structure is captured by a bivariate copula. Below are the steps for disaggregating a daily rainfall amount R (mm):

rainfall depth (mm)

(a) estimate L for a given value of R from a bivariate copula as explained in Section 2; (b) estimate the number of wet periods nL as the nearest integer of (L/ts + 0.5);

rmax = r3 timescale

r5

r10

st start time

d1

(c) select an OSDSP type and determine the start time st within the wet day as explained in Section 3; (d) create a pool of OSDSPs of the selected type and of the same, or very close, numerical value of nL or L; (e) select from the pool one OSDSP that has the closest value of lag-1 autocorrelation of the natural logarithm of the wis to a simulated value conditioned on L as explained in Section 3; and (f) multiply the wis of the selected OSDSP by R to obtain the complete daily storm profile, which is shifted to start at st if it is of type SSP1 as described in Section 3. 1.3. Data Six minutes timescale rainfall data at 0.01 mm depth resolution for the 8 Australian State capital cities obtained from the Bureau of Meteorology (BOM), Australia, are used to develop the OSDSPs. The capital cities are situated in different climatic regions as shown in Fig. 2. Some capital cities (e.g., Perth and Adelaide) experience maximum rainfall in winter, and for others the maximum occurs in summer (e.g., Darwin and Brisbane). For some capital cities (e.g., Canberra and Melbourne) there is not much seasonal variation in monthly rainfall averages (http://www.bom.gov.au/climate/ averages/tables/). Verification of the model is carried out using 1min timescale continuous rainfall data, aggregated to six minutes timescale to conform to the BOM data set, collected at the Gregory Erosion Control Experimental field trials site [10,12] between January 1998 and December 2009. This 12-year continuous data set was collected by a tipping bucket of tipping depth resolution of 0.18 mm. Fig. 2 shows the location of the Gregory field trials site and the nearest two BOM stations (Emerald and Dingo) within 70 km radius with 6 min timescale, 0.01 mm depth resolution, rainfall data. The Emerald (035147) and Dingo (035025) rainfall data were combined into one data set (E–D) and used to develop the regional model parameters. As indicated in Table 1, a high proportion of the BOM rainfall data are embedded in months with missing wet days. However, the model formulation presented allows the use of every good quality wet day data irrespective of missing data for the other days within the same month. This is a major advantage over other models in the literature that discard monthly data with missing values. For each State capital city, the data were grouped into 13 classes based on L quantiles of E–D [14]) as shown in Table 1. The rainfall data were then analysed at nine aggregation levels (0.1, 0.2, 0.5, 1, 2, 4, 8, 12, and 24 hours) which were carefully chosen so that 24 h (a day) is an integral multiple of each one of them. False events defined as daily rainfall amounts less than 0.1 mm were not included in the analysis. Note that, due to the short record length, the Gregory data was grouped into only 4 L-quantile classes (Table 1). An overview of copula modelling philosophy and parameter estimation procedure is first presented, followed by detailed analysis of the OSDSPs from the Australian State capital cities. The disaggregation model is then verified through two case studies in the penultimate section before providing the summary and conclusions.

2. Copula modelling of the dependence structure between R and L ts = sampling

r1

27

d2

d3

time (hours)

Fig. 1. Definition sketch of an observed daily storm profile.

The copula modelling described here is limited to the 2-dimensional (2D) case, and higher dimensional cases can be found in Joe [16], Nelson [20] and Salvadori et al. [23]. In this paper the 2 random variables to be modelled by bivariate copulas are the total daily rainfall depth R (mm) and the total wet periods’ duration L (hours). Following the Sklar theorem which is the basis of copula modelling, the joint cumulative distribution function (cdf), F, of the random vector X = (R, L) is given as:

28

Y. Gyasi-Agyei / Advances in Water Resources 45 (2012) 26–36

Fig. 2. Location of Australian State capital cities and the case study sites.

FðxÞ ¼ PrðX 6 xÞ ¼ CfF R ðrÞ; F L ðlÞg

ð1Þ

where FR(r) and FL(l) are the continuous marginal cdfs of the random variables R and L and C is a unique copula function that maps the 2dimensional uniform distributions to a single dimension uniform distribution expressed as [0, 1]2 ? [0, 1]. Assigning the probabilities u = FR(r) and v = FL(l) for the quantiles r and l, the marginal distributions transform quantiles (r and l) into the unit [0, 1] interval having the standard uniform distribution U(0, 1). This implies that the quantiles can be estimated by the inverse of the marginal distribu1 tions as r ¼ F 1 R ðuÞ and l ¼ F L ðv Þ leading to the Sklar inversion theorem given as:

n o 1 Fðr; lÞ ¼ F F 1 R ðuÞ; F L ðv Þ ¼ Cðu; v Þ

ð2Þ

For a given pair of data sets of length n, Genest and Favre [6] recommended estimating the empirical pair values (ui, vi) by their ranks, Si and Ti, defined as:

ui ¼

Si ; nþ1

vi ¼

Ti nþ1

ð3Þ

Since the random variables R and L are considered continuous, discontinuities and ties in the rainfall data as a result of the sampling process needed to be removed first. Following Gyasi-Agyei [14] and Vandenberghe et al. [28], ties were removed and the rainfall data made continuous by applying a random noise sampled from a uniform distribution. For the random variable R which has a resolution of 0.01 mm, the random noise rai added to each ri was sampled from rai  U (0.005, 0.005) before applying a proportional adjustment to ensure that the total amount added to R matches Ra drawn from the same distribution as rai. In the case of the random variable L which has a resolution of 0.1 h (6 min), the added random noise La was drawn from La  U (0.05, 0.05). A pair-wise (R–L) scatter plot of the empirical rank distribution revealed that, for all months, the dependence structure is nearly elliptical with a strong positive correlation. For the E–D data set the Kendall’s Tau correlation coefficient ranged between 0.31 and 0.62 and the Spearman’s rank ranged between 0.46 and 0.82. Seven single parameter bivariate copulas were considered: Ali–Milhail– Haq (AMH), Gumbel–Hougaard (GH), Clayton (CL), Frank (FR), Plackett (PL), normal (NO) and t (with different degrees of free-

dom). The goodness-of-fit of the copulas considered were assessed using the R copula package function gofCopula which is based on the Cramer–von Mises Sn statistic [24,21,17]. The maximum likelihood method was used to estimate the parameters with the number of multiplier iterations set to 10,000. From Table 2 it is observed that the elliptical copulas (normal and t) generally performed better than the other copulas considered, and at the 5% significance level both normal and t10 copulas cannot be rejected for all months except for November which is acceptable at the 1% significance level. The normal copula (NO) having the highest monthly average value was selected to model the dependence between R and L for all months. On a monthly basis the histograms of the marginals of the random variables R and L are positively skewed, so the following 2parameter family distributions were selected for testing: Pareto, Gumbel, Weibull, gamma and lognormal. For the marginal distribution fitting, the Anderson–Darling goodness-of-fit test was used to select the best distribution. The lognormal, with scale (l = meanlog) and shape (r = sdlog) parameters estimated by the maximum likelihood method, emerged as the best distribution for both R (monthly p-values ranging between 0.08 and 0.82) and L (monthly p-values ranging between 0.09 and 0.96), and therefore it cannot be rejected at the 5% significance level. The lognormal shape parameter r appeared to have minimal monthly variation, so it was considered to be the same for all months. A first harmonic Fourier series was used to capture the seasonal variation of the normal copula (hm) and the lognormal scale (lm) parameters given as:

hm ; or



lm ¼ a0 þ a1 cos m

   2p 2p þ a2 sin m 12 12

ð4Þ

where a0, a1 and a2 are constants to be calibrated and m is the monthly index (1 for January and 12 for December). For the marginals, a0, a1, a2 and r were calibrated together by maximising the sum of the log-likelihood functions of all the months. Table 3 provides the calibrated values and Fig. 3 shows with solid lines the fitted first harmonic Fourier series. Seasonal variability of the lognormal marginals and the normal copula model parameters is apparent. By fitting the first harmonic Fourier series, the number of model parameters is considerably reduced from 60 to just 11, 8 for the marginals and 3 for the normal copula.

29

Y. Gyasi-Agyei / Advances in Water Resources 45 (2012) 26–36 Table 1 Distribution of the number of wet days according to total daily wet periods’ duration (L). Station

Year

ID No.

066062

086071

040223

009021

023000

094029

Name

SYD

MEL

BNE

PER

ADE

HOB

CAN

From to

1913 2006

1910 2006

1949 2000

1961 2006

1910 1979

1911 2005

1937 2006

E–D L quantile Class

Probability

070014

014015

Gregory site

DAR

035147 035025 E–D

1953 2006

1963 2006

1998 2009

Number of wet days in a quantile class

GRE L quantile

Max (h)

1 0.10 0.48 2 0.20 0.81 3 0.30 1.19 4 0.40 1.68 5 0.50 2.37 6 0.60 3.23 7 0.70 4.38 8 0.80 6.05 9 0.90 9.00 10 0.95 12.13 11 0.97 14.82 12 0.99 19.57 13 1.00 24.00 Total number of wet days % in a month with missing wet days data

1210 1058 902 945 1072 1020 946 976 985 541 317 246 110 10,328 44

1448 1531 1435 1557 1734 1603 1465 1316 1068 417 176 96 28 13,874 51

549 536 461 478 523 478 440 382 368 199 82 81 36 4613 60

354 246 232 270 345 368 384 397 394 218 98 115 35 3456 70

The conditional simulation procedure for sampling the quantile l of random variable L for a given quantile r of the random variable R is as follows (e.g., Salvadori et al. [23]): (a) the probability u = FR(r) is calculated from the lognormal marginal distribution FR () for R with the monthly parameters obtained from Eq. (4) or Fig. 3); (b) independent probability s is sampled from the standard uniform distribution U(0, 1); (c) quantiles x = U1(u) and y = U1(s) are calculated, where U1() is the inverse of the standard normal distribution U  N(0, 1); qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi (d) y is transformed as y ¼ hm  x þ 1  h2m  y where hm is the normal copula dependence parameter between R and L for month m obtained from Eq. (4) or Fig. 3; (e) the probability v is estimated as v = U(y); and 1 (f) the quantile l is given by l ¼ F 1 L ðv Þ, where F L ðÞ is the inverse of the lognormal marginal distribution for L with the monthly parameters obtained from Eq. (4) or Fig. 3. 3. Analysis of observed daily storm profiles The question being addressed here is: can an OSDSP from one site be used at another site of different climatic regions? If the answer is ‘yes’, then there is a large pool of OSDSPs in Australia, and indeed around the world, to sample from without any numerical modelling of this facet of the daily rainfall disaggregation model. All that would be required then is a copula model for the dependence structure between R and L as presented in the preceding section. To address this question the daily data for each station were analysed based on the 13 L-quantile classes of E–D (Table 1). The properties being investigated here are the dry probability DP (a signature for the wet and dry sequences), the lag-1 autocorrelation qw of the natural logarithm of wet period depths ris for each wet day. Since ri is assumed to follow a lognormal distribution [14], the statistic qw was obtained by fitting the lag-1 autoregressive model to the natural logarithm of the ris for each wet day. 3.1. Diurnal effect and types of daily storm profiles Fig. 4 shows the superposition of Emerald–Dingo (E–D) L-quantile class DP curves with those of Sydney and Adelaide, and also compares

GRE

919 937 872 904 1049 998 903 787 599 181 42 22 2 8215 49

1239 1146 1052 1186 1342 1173 1081 1063 1012 479 197 156 54 11,180 70

448 394 378 391 430 437 492 427 421 230 95 78 28 4249 48

545 510 497 410 452 409 383 346 320 164 69 66 29 4200 71

296 296 295 296 296 295 296 295 296 148 59 59 30 2957 45

Probability

Max (h)

No.

0.3 0.6 0.9 1.0

0.26 0.90 3.68 24.00

259 258 258 87

863 0

the three south eastern seaboard stations (Sydney, Brisbane and Hobart). The average absolute percentage errors of the DP curves of the State capital cities from those of Emerald–Dingo are: SYD (12.4), MEL (11.7), BNE (14.0), PER (11.3), ADE (20.1), HOB (12.4), CAN (7.3) and DAR (15.4). It is observed that the DP curves for the south eastern capital cities are very similar but they are different from those of the Emerald–Dingo region. Canberra DP curves are the closest to the Emerald–Dingo ones and the Adelaide DP curves have the maximum average absolute percentage error. Diurnal effects largely explain the differences in the DP curves of the various rainfall stations. For example, a 4-h storm starting between midnight and 8 am, or between noon and 8 pm will give a 12-h DP value of 0.5, but the value will be zero if the starting time is between just after 8 am and just before 12 noon as a result of the aggregation process. Hence sites experiencing a lot of early morning or late afternoon rainfall will tend to have higher aggregated DP values compared to sites associated with uniform storm start times during the day, or occurring during late mornings. In general, sites experiencing prolonged rainfall with high intermittency will have lower DP values. Therefore, there should be rules in selecting OSDSPs of a site for use at a different climatic site. An OSDSP will fall into one of the following types: (a) SSP1 – the start and end times are within the second and last but one periods, i.e., no rainfall during the first or last periods of the wet day; (b) SSP2 – there is no rainfall during the first period but there is rainfall during the last period, indicating rainfall crossing midnight and spreading to the next Julian day; (c) SSP3 – there is rainfall during the first period but there is no rainfall during the last period, indicating rainfall on the previous Julian day extending past midnight; and (d) SSP4 – there is rainfall during both the first and the last periods, indicating rainfall on the previous and current Julian days extending past midnight. It needs to be underlined that there are the same numbers of SSP2s as SSP3s for each station as they occur concurrently. An SSP4 can occur after an SSP2 provided there are at least 3 consecutive Julian wet days. The last wet day within a cluster of wet days should experience an SSP3. All isolated wet days will experience an SSP1, and so will some wet days within a cluster. Fig. 5 shows the proportions of the various OSDSP types in the data sets. There is a

Y. Gyasi-Agyei / Advances in Water Resources 45 (2012) 26–36

the TP is estimated using Eq. (5) and the L value obtained from the bivariate copula. If a value sampled from a standard uniform distribution is less than or equal to the TP, then an SSP2 is selected, otherwise it is an SSP1. All wet days following a previous wet day of SSP2 or SSP4, and at the end of a cluster of wet days, should be an SSP3, where a wet day follows an SSP2 or SSP4 and it is not at the end of a cluster of wet days, then an SSP3 or SSP4 are both possible. The next Julian wet day following an SSP3 should be an SSP1 if it is at the end of a cluster of wet days, or an SSP2 if it is in the middle of a wet cluster depending on the relative values of TP and the value sampled from the standard uniform distribution. It is observed that the start time of an SSP1 can be approximated by a uniform distribution between 0.1 (second period) and 24-L, the available start time, represented as st  U(0.1, 24-L). For the other storm profile types, the start and end times are the same as those of the selected OSDSP in accordance with their definition.

Table 2 p-Values of the Cramer–von Mises Sn statistic used to assess the goodness-of-fit of the selected copulas. Month

GH

FR

PL

NO

t5

t10

1 2 3 4 5 6 7 8 9 10 11 12 Average

0.000 0.003 0.502 0.091 0.003 0.165 0.079 0.016 0.078 0.002 0.000 0.002 0.078

0.010 0.489 0.035 0.008 0.472 0.008 0.379 0.224 0.227 0.195 0.014 0.208 0.189

0.036 0.189 0.111 0.064 0.009 0.031 0.162 0.057 0.533 0.294 0.015 0.336 0.153

0.282 0.174 0.389 0.329 0.220 0.129 0.278 0.151 0.598 0.498 0.024 0.436 0.292

0.167 0.014 0.375 0.610 0.009 0.145 0.169 0.032 0.618 0.173 0.004 0.285 0.217

0.315 0.081 0.445 0.460 0.075 0.148 0.253 0.100 0.698 0.401 0.016 0.433 0.285

Copulas: Gumbel–Hougaard (GH), Frank (FR), Plackett (PL), normal (NO) and t; Ali– Milhail–Haq (AMH) – all p-values equal 0; Clayton (CL) – all p-values less than 0.004 except November which gave 0.026.

3.2. Pool creation and final selection All OSDSPs of the selected type (SSP1, SSP2, SSP3, or SSP4) having L values equal, or close, to that obtained from the bivariate copula are put in a pool. For an SSP1, the storm start time of each in the pool is shifted (Fig. 6) to the start time st sampled from U(0.1, 24-L) and, if any of the wis is lost (falls outside the 24-h period) as a result of the shifting, it is discarded from the pool. Where the number of pooled OSDSPs is less than a set minimum (say 10), additional ones with L ± k  ts are included in the pool, where k is an integer starting from 1 and increased by 1 at a time until the target number is reached. This is only likely to happen for large L values, say greater than 12 h. Fig. 7(a) shows the distribution of the qw statistic estimated for each wet day for Sydney. Also shown on this plot are the 13 Lquantile class mean values which exhibit an increasing trend with the random variable L. An example distribution of the qw statistic within L-quantile class 8 of Sydney is shown in Fig. 7(b). While the class distribution type appears to be invariant, its parameters vary significantly from site to site as demonstrated by the qqplot of Sydney against the other stations for class 8. Clearly, the distribution parameters of the Emerald–Dingo data set are different from the other sites and therefore the qw statistic should be taken into account in the selection of an OSDSP from the pool. This is confirmed by the p-values (Table 4) of the two-sample KS-test used to check whether the empirical cumulative distributions of the same class of Sydney and the other stations are from the same continuous distribution. Most of the classes for Hobart, Melbourne and Brisbane cannot be rejected at the 1% significance level as coming from the same distribution as Sydney. None of the Emerald–Dingo and Darwin classes passed at the 1% significance level. The simulation model for the qw statistic proposed by GyasiAgyei [14] is used to sample it for a given value of L for the

Table 3 Emerald–Dingo (E–D) first harmonic Fourier series parameters. Parameter

a0

a1

a2

R – meanlog (l) L – meanlog (l) R–L normal copula, h

1.531 0.824 0.659

0.148 0.206 0.111

0.097 0.050 0.020

a0, a1 and a2 are the coefficients defined in Eq. (4).

considerable variation in the proportion of the different storm profile types from site to site, and, in general, this variation increases with L. The Emerald–Dingo region experiencing the highest proportions of SSP2/SSP4 indicates a higher degree of diurnal effects with rainfall overlapping two consecutive Julian wet days. Since the main trigger for an SSP3 or SSP4 is for an SSP2 to occur within a cluster of wet days, it suffices to model the probability of a transition from SSP2 to SSP3/SSP4 given two consecutive Julian wet days. This is essentially the probability of an SSP2 occurring given two consecutive Julian wet days. Such a transition probability (TP) is displayed in Fig. 5(c) which shows a positive linear trend with L. Again, the transition probability is highest for the Emerald–Dingo data set indicating a higher degree of diurnal effects. Adelaide rainfall is the least influenced by diurnal effects. The straight line fitted to the Emerald–Dingo values is given as:

TP ¼ 0:23 þ 0:029L

ð5Þ

(a) 1.3

R

2

4

6

month

8

10

12

σ = 1.118

(b) 2

L 4

6

month

8

10

12

0.7 0.6

parameter θ

0.8

0.8 0.9 1.0 1.1

parameter μ

1.6 1.5 1.4

parameter μ

1.7

σ = 1.404

0.5 0.6 0.7

1.8

If the daily rainfall to be disaggregated is isolated, i.e., separated by dry days, then automatically it is an SSP1. For a daily rainfall beginning a cluster of at least two consecutive Julian wet days,

0.5

30

(c) 2

R-L 4

6

8

10

12

month

Fig. 3. Calibrated lognormal (l, r) and normal copula parameters (h): observed monthly values in open circles and fitted first harmonic Fourier series in solid lines, Emerald– Dingo.

31

0.1

1

10

1.0 0.8 0.4

0.6

(c)

0.0

dry probability (-)

SY D BNE HOB

0.2

1.0 0.8 0.6 0.4

dry probability (-)

(b)

0.0

0.2

0.4

0.6

(a)

E-D ADE

0.2

0.8

E-D SY D

0.0

dry probability (-)

1.0

Y. Gyasi-Agyei / Advances in Water Resources 45 (2012) 26–36

0.1

timescale (hours)

1

10

0.1

timescale (hours)

1

10

timescale (hours)

Fig. 4. Comparison of Emerald–Dingo (E–D) class dry probability curves with (a) Sydney (SYD), (b) Adelaide (ADE) and (c) of Sydney (SYD), Brisbane (BNE) and Hobart (HOB): the topmost curve represents the driest (class 1) and the bottommost curve the wettest (class 13) with classes 4, 6, 8, 9, 10, 11, 12 in between.

0

time (hours) Fig. 6. An example of shifted daily storm profile: selected profile (dashed line) shifted to a new start time (solid line).

needed for sites close to the eastern seaboard where its variation may not be significant as there are many OSDSPs to choose from, thus reducing the number of parameters by a third.

ð6Þ

For the first L-quantile class, qtw is sampled from the normal distribution N(l, r2⁄) with the parameters derived as:

4. Verification of the disaggregation model The Gregory site case study is presented in this section to show how available long records of daily data can be disaggregated. Normally the site of interest may not have fine timescale rainfall records to establish the required model parameters. It is suggested that nearby sites with the required timescale rainfall records be used. In the example presented here, Emerald and Dingo BOM stations are the nearest, being within 70 km on either side of the case study site (Fig. 2). Because of the short record length, the two data sets were combined into one (E–D) and used to estimate the 19 parameters (3 for the normal copula, 8 for the R and L marginal lognormal distributions, 2 for the transition probability Eq. (5)), and 6

5

10

15

20

mean wet duration (hours), L

25

0

5

10

15

20

mean wet duration (hours), L

25

0.8 0.6

0.4 0.3 0.2

(b)

(c)

0.4

0.5

0.6 0

HOB CAN DAR E-D

0.1

0.2 0.1

(a)

SYD MEL BNE PER ADE

0.0

CAN DAR E-D

proportion of SSP4

PER ADE HOB

0.3

0.4

SYD MEL BNE

0.0

proportion of SSP2

0.5

where l and r are obtained from the fitted polynomial curves of Eq. (6). The qw statistic is then obtained as (1  qtw ). Hence, for each wet day with a given value of L, the OSDSP with the least absolute difference value of its qw statistic from the one sampled from the lognormal distribution is picked from the pool to represent the wet day. The wis of the selected OSDSP are then multiplied by R to obtain the complete daily storm profile. It needs to be underlined that the incorporation of the qw statistic in the final selection may not be

HOB CAN DAR E-D fitted

SYD MEL BNE PER ADE

0.2

ð7Þ

0.0

l ¼ exp½l þ r2 =2 r2 ¼ exp½2l þ r2 ½expðr2 Þ  1

24

old st

new st

probability: SSP2 to SSP3/SSP 4

lðLÞ ¼ 0614  1:398log10 ðLÞ  0:302log10 ðLÞ2 rðLÞ ¼ 0:704 þ 0:518log10 ðLÞ  0:432log10 ðLÞ2

rainfall depth (mm)

Emerald–Dingo data set. In order to make all values of qw positive for a right skewed distribution fitting, the qw values were transformed as qtw ¼ ð1  qw Þ which is basically a reflection along the axis qw = 1. The Anderson–Darling goodness-of-fit test was adopted to select a common marginal distribution for the 13 Lquantile classes from the Normal, Weibull, gamma and lognormal distributions using the maximum likelihood method to estimate the parameters. Here too the lognormal distribution emerged as the best with 11 L-quantile class p-values greater than 0.05, with only the first L-quantile class failing at the 1% significance level but passing at the 5% significance level for the normal distribution. The variation of the lognormal parameters (l = meanlog and r = sdlog) with the random variable L class mean were fitted with the following polynomials:

0

5

10

15

20

25

mean wet duration (hours), L

Fig. 5. Variation of proportions of (a) SSP2 and (b) SSP4 with mean wet duration and (c) given two consecutive Julian wet days, the transition probability of rainfall to cross midnight of first day and spread to the next day for the different stations: straight line is fitted to E–D values.

0.5

station quantile

60

class 8

-0.5

0

10

-0.4

wet duration (hours), L

0.0

0.2

0.4

0.6

0.8

1.0

HOB CAN DAR E-D

MEL BNE PER ADE

20

40

frequency 1

SYD: class 8

(c)

0.0

100

(b)

80

0.5 0.0 -0.5

SY D indiv idual class mean

-1.0

lag-1 auto-cor. coef.

(a)

1.0

Y. Gyasi-Agyei / Advances in Water Resources 45 (2012) 26–36 1.0

32

-0.5

0.0

lag-1 auto-cor. coef.

0.5

1.0

SYD quantile

Fig. 7. Lag-1 autocorrelation (qw statistic): (a) distribution for Sydney, (b) L-quantile class 8 histogram for Sydney and (c) qqplot of Sydney versus the other stations for Lquantile class 8.

Table 4 KS-test p-values comparing Sydney empirical cumulative distribution function of the qw statistic with the other stations for the different L-quantile classes. 2

3

4

5

6

7

8

9

10

11

12

13

0.004 0.000 0.258 0.169 0.682 0.003 0.000 0.002

0.698 0.000 0.322 0.598 0.311 0.119 0.000 0.000

0.250 0.000 0.273 0.213 0.791 0.000 0.000 0.000

0.042 0.000 0.046 0.000 0.273 0.014 0.000 0.000

0.009 0.006 0.000 0.000 0.370 0.028 0.000 0.000

0.016 0.009 0.000 0.000 0.265 0.043 0.000 0.000

0.042 0.011 0.000 0.000 0.365 0.005 0.000 0.000

0.033 0.839 0.000 0.000 0.029 0.035 0.000 0.000

0.000 0.769 0.000 0.000 0.041 0.000 0.000 0.000

0.383 0.152 0.684 0.003 0.366 0.000 0.001 0.000

0.302 0.444 0.665 0.009 0.403 0.000 0.000 0.000

0.188 0.133 0.126 0.058 0.699 0.000 0.000 0.000

0.028 0.346 0.006 0.936 0.525 0.052 0.007 0.000

0.0

0.2

0.4

0.6

0.8

rank distribution, R

0.05 1.0

0.1

0.5

5.0

50.0

daily rainfall depth (mm), R

5.0 0.5 1.0 2.0

wet duration (hours), L

(b)

GRE: Dec

0.1 0.2

1.00

5.00

E-D: Dec

0.20

wet duration (hours), L

0.8 0.6 0.4 0.2

(a)

0.0

rank distribution, L

E-D: Dec

20.00

1

MEL BNE PER ADE HOB CAN DAR E–D

1.0

Class

(c) 0.2

0.5

2.0

5.0

20.0

daily rainfall depth (mm), R

Fig. 8. Comparison of the covariance structure: (a) rank distribution for E–D December, (b) actual values for E-D December and (c) actual values for GRE December for one simulation run using all capital cities’ (8CC) OSDSPs.

for sampling the qw statistic from the lognormal distribution Eq. (6) required to disaggregate the rainfall data for the whole year. A depth resolution of 0.01 mm and 6 min timescale were used in the simulation. For the Gregory case study site, the simulated data were re-sampled with the 0.18 mm depth resolution of the tipping bucket used to record its data. The statistics used to compare the observed and simulated rainfall are the dry probability, variance, lags 1–3 autocorrelations, and the extreme values (intensity–duration– frequency curves). Independently, the capital cities’ and E–D’s OSDSPs were used to create the pool. In addition, all 8 capital cities’ (8CC) OSDSPs were combined to create a pool of over 60,000 OSDSPs to choose from. Results are presented for the disaggregation of the non-continuous E–D daily rainfall data (R P 0:1 mm) used in the parameter estimation and also for the continuous Gregory rainfall data (R P 0). Due to space constraints the results are shown for the simulation (0.1 h) and 1 h aggregation timescales for all months and then for all aggregation timescales for December. At each aggregation timescale, the 2nd and 39th ranks of 8CC are used to define the 95% prediction limits.

Fig. 8 compares the covariance structure of the daily rainfall depth and total wet periods’ duration of one simulation run results and the observed values for December using 8CC OSDSPs. Evidently, the copula model captures the covariance structure very well, and this was the case for all simulation runs. For the evaluation of the daily rainfall disaggregation model, 40 simulation runs were carried out for each stations’ OSDSPs. Only 8CC, SYD and HOB OSDSPs could simulate enough data in the 13th L-quantile class of E–D. Table 5 shows the mean KS-test p-values of the 40 simulation runs comparing E–D observed empirical cumulative distribution function of the qw statistic with the simulation results. Almost all L-quantile class p-values passed the 5% significance level with the exception of ADE which could not generate enough data in the 12th and 13th L-quantile classes. L-quantile class 1 has 4 or fewer wi values for each wet day which is not enough for accurate estimation of the qw statistic. Comparison of the observed and simulated dry probability analysed on the 13 and 4 L-quantile classes for the Emerald–Dingo and the Gregory site is depicted in Fig. 9, with the grey shaded areas

33

Y. Gyasi-Agyei / Advances in Water Resources 45 (2012) 26–36

Table 5 Mean KS-test p-values of the 40 simulation runs by comparing E–D observed empirical cumulative distribution function of the qw statistic with the simulation results for the different L-quantile classes. Class

2

3

4

5

6

7

8

9

10

11

12

13

8CC SYD MEL BNE PER ADE HOB CAN DAR E–D

0.193 0.236 0.259 0.287 0.232 0.176 0.225 0.275 0.225 0.341

0.554 0.440 0.392 0.359 0.174 0.396 0.374 0.347 0.434 0.468

0.486 0.360 0.368 0.419 0.282 0.306 0.481 0.550 0.466 0.628

0.130 0.176 0.202 0.124 0.126 0.126 0.169 0.109 0.170 0.210

0.512 0.313 0.176 0.249 0.028 0.011 0.209 0.068 0.350 0.580

0.259 0.162 0.139 0.049 0.098 0.014 0.092 0.154 0.009 0.173

0.519 0.201 0.285 0.043 0.117 0.000 0.360 0.283 0.164 0.491

0.433 0.175 0.050 0.221 0.255 0.000 0.111 0.429 0.171 0.234

0.227 0.198 0.048 0.072 0.100 0.000 0.036 0.192 0.213 0.441

0.583 0.562 0.560 0.236 0.417 0.008 0.479 0.468 0.484 0.555

0.329 0.110 0.309 0.354 0.148 NA 0.312 0.503 0.529 0.442

0.180 0.155 NA NA NA NA 0.492 NA NA NA

1.0 0.8 0.6

(c)

0.0 1

0.4

(b)

0.0 0.1

8CC

0.0

(a)

dry probability (-)

SYD

obs sim max

0.2

0.8 0.6 0.4

dry probability (-)

8CC

obs sim max

0.2

0.4

0.6

0.8

obs sim max

0.2

dry probability (-)

1.0

1.0

NA – not enough data for KS test; class 1 has 4 or fewer values for each wet day and not enough data for reliable assessment of the KS-test, and all p-values were equal to 0.

10

0.1

aggregation timescale (hours)

1

10

0.1

aggregation timescale (hours)

1

10

aggregation timescale (hours)

4

6

8

12

2

4

month

6

8

10

0.90

obs 8CC SYD MEL BNE PER

0.85

dry probability

0.95

(c)

0.75 0.80

0.98

1.00

10

ADE HOB CAN DAR E-D lim

obs 8CC SYD MEL BNE PER

0.92

0.975

2

0.96

dry probability

0.985

ADE HOB CAN DAR E-D lim

obs 8CC SYD MEL BNE PER

1

(b)

0.94

0.1

(a)

0.965

dry probability

0.995

Fig. 9. Comparison of observed and simulated dry probability curves of classes 1 (top), 4, 6, 8, 9, 10, 11, 12 and 13 (bottom) for Emerald–Dingo using: (a) 8CC (all State capital cities) OSDSPs, (b) SYD (Sydney) OSDSPs and (c) the 4 classes of Gregory using all capital cities OSDSPs, the shaded grey area indicates the 95% prediction limits (2nd and 39th ranks), the max plot is for the case of only one wet period; daily rainfall depth R P 0:1 mm.

12

DEC 0.1

month

ADE HOB CAN DAR E-D lim

1

10

aggregation timescale (hours)

2

4

6

month

8

10

12

(b) 2

4

6

month

8

10

DEC

12

10

2

variance ( mm )

1

obs 8CC SYD MEL BNE PER

1 0.1

8 2

4

6

2

variance ( mm )

0.1

CAN DAR E-D lim

BNE PER ADE HOB

obs 8CC SYD MEL

100

10

(a)

0

0.00

0.10

0.20

CAN DAR E-D lim

BNE PER ADE HOB

obs 8CC SYD MEL

2

variance ( mm )

0.30

Fig. 10. Comparison of observed and simulated dry probability curves: (a) aggregation timescale of 0.1 h, (b) aggregation timescale of 1h and (c) December; the limit (lim) plots indicate the 95% prediction limits (2nd and 39th ranks) for 8CC; daily rainfall depth R P 0 m, Gregory.

(c) 0.1

1

ADE HOB CAN DAR E-D lim

10

aggregation timescale (hours)

Fig. 11. Comparison of observed and simulated variance: (a) aggregation timescale of 0.1 h, (b) aggregation timescale of 1 h and (c) December; the limit (lim) plots indicate the 95% prediction limits (2nd and 39th ranks) for 8CC, daily rainfall depth R P 0:1 mm, Emerald–Dingo.

Y. Gyasi-Agyei / Advances in Water Resources 45 (2012) 26–36 100 6

8

10

ADE HOB CAN DAR E-D lim

12

2

4

6

8

DEC

10

10 1

obs 8CC SYD MEL BNE PER

0.01

0.1

2

variance ( mm )

0.0

(b)

month

(c)

1

0.5

0.01

4

1.0

2

variance ( mm )

0.1

(a) 2

obs 8CC SYD MEL BNE PER

1.5

ADE HOB CAN DAR E-D lim

0.02

0.03

0.04

obs 8CC SYD MEL BNE PER

0.00

2

variance ( mm )

0.05

34

12

0.1

month

1

ADE HOB CAN DAR E-D lim

10

aggregation timescale (hours)

2

6

8

10

0

0.5

12

2

4

6

month

8

10

BNE PER ADE HOB

CAN DAR E-D lim

60

80

obs 8CC SYD MEL

40

intensity (mm/h)

100

(b)

0.1 4

CAN DAR E-D lim

0

(a)

BNE PER ADE HOB

obs 8CC SYD MEL

50

intensity (mm/h)

200

300

CAN DAR E-D lim

100 0

50

intensity (mm/h)

BNE PER ADE HOB

obs 8CC SYD MEL

20

150

100

Fig. 12. Comparison of observed and simulated variance: (a) aggregation timescale of 0.1 h, (b) aggregation timescale of 1 h and (c) December; the limit (lim) plots indicate the 95% prediction limits (2nd and 39th ranks) for 8CC, daily rainfall depth R P 0, Gregory.

12

(c) 2

1 4

month

6

8

10

12

month

2

4

6

8

10

12

80

(b) 2

4

6

8

month

10

12

60

BNE PER ADE HOB

CAN DAR E-D lim

40

0.5

month

obs 8CC SYD MEL

20

intensity (mm/h)

100 60

80

CAN DAR E-D lim

0

0

0.1

BNE PER ADE HOB

0

20

100

(a)

obs 8CC SYD MEL

40

CAN DAR E-D lim

intensity (mm/h)

BNE PER ADE HOB

150

200 250

obs 8CC SYD MEL

50

intensity (mm/h)

300

120

Fig. 13. Comparison of observed and simulated 20-year intensity–frequency–duration curves: (a) aggregation timescale of 0.1 h, (b) aggregation timescale of 0.5 h and (c) aggregation timescale of 1 h; the limit (lim) plots indicate the 95% prediction limits (2nd and 39th ranks) for 8CC, daily rainfall depth R P 0:1 mm, Emerald–Dingo.

(c) 2

1 4

6

8

10

12

month

Fig. 14. Comparison of observed and simulated 10-year intensity–frequency–duration curves: (a) aggregation timescale of 0.1 h, (b) aggregation timescale of 0.5 h and (c) aggregation timescale of 1 h; the limit (lim) plots indicate the 95% prediction limits (2nd and 39th ranks) for 8CC, daily rainfall depth R P 0:1 mm, Gregory.

representing the 95% prediction limits. The uppermost curve (max) defines the maximum dry probability given as 1  hk/24 where hk is the aggregation timescale in hours and it is for the case of only one wet period during the wet day. For the Emerald–Dingo site, the simulated average absolute percentage error was 7.9 and 6.9 for 8CC and SYD respectively, and 4.2 and 5.1 for the Gregory site. DAR recorded the worst error of 13.6% for Emerald–Dingo but only 4.3% for the Gregory site. These small percentages indicate that the predictions of the dry probabilities are considered to be very good. Fig. 10 shows the reproduction of the dry probability for Gregory on a continuous (R P 0) and monthly basis. For all data sets the reproduction of the dry probability is excellent, with the average absolute percentage errors being less than 2%.

Figs. 11–15 demonstrate the reproduction of the variance, the intensity–frequency–duration curves and autocorrelations by the daily rainfall disaggregation model presented. Most of the observed values being within the 95% prediction limits indicate that the predictions are good, particularly for 8CC, SYD, MEL, HOB and CAN. With the exception of DAR OSDSPs, the mean absolute percentage error of the standard deviation is less than 5% for all months other than June and July which is about 20%. DAR provided the worst performance of over 40% mean absolute percentage error of the standard deviation. At the simulation timescale, the observed autocorrelations appear to be seasonally invariant with some noise at 0.7, 0.5 and 0.4 for lags 1–3 respectively (Fig. 15). These averages were cap-

35

8

10

2

4

8

10

12

2

4

8

10

0.8 0.6 0.4 0.0

lag-2 autocorrelation

12

8

10

12

0.5

CAN DAR E-D lim

1

0.2

0.3

0.4

BNE PER ADE HOB

(h) 2

4

month

6

8

month

10

obs 8CC SYD MEL BNE PER

ADE HOB CAN DAR E-D lim

(f) 0.1

1

10

aggregation timescale (hours)

0.1

(g)

obs 8CC SYD MEL

DEC

10

12

DEC

0.4

0.6 6

CAN DAR E-D lim

lag-3 autocorrelation

BNE PER ADE HOB

4

1

0.0

0.6 0.4 0.2 0.0

lag-3 autocorrelation

2

CAN DAR E-D lim

month

0.1

1

aggregation timescale (hours)

BNE PER ADE HOB

6

(c) 0.1

(e)

month

obs 8CC SYD MEL

12

0.0 0.1 0.2 0.3 0.4 0.5 0.6

obs 8CC SYD MEL

DEC

obs 8CC SYD MEL BNE PER

0.3

6

10

CAN DAR E-D lim

ADE HOB CAN DAR E-D lim

0.2

4

(d)

CAN DAR E-D lim

lag-2 autocorrelation

BNE PER ADE HOB

0.0 0.1 0.2 0.3 0.4 0.5 0.6 0.7

0.8 0.6 0.4 0.2

lag-2 autocorrelation

0.0

2

8

month 0.1

obs 8CC SYD MEL

6

BNE PER ADE HOB

0.2

0.8 0.6

12

obs 8CC SYD MEL

0.1

6

month

lag-3 autocorrelation

4

(b)

1

0.0

2

CAN DAR E-D lim

BNE PER ADE HOB

lag-1 autocorrelation

1.0

(a)

0.1

obs 8CC SYD MEL

0.4

CAN DAR E-D lim

0.2

BNE PER ADE HOB

lag-1 autocorrelation

obs 8CC SYD MEL

0.0

lag-1 autocorrelation

0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0

Y. Gyasi-Agyei / Advances in Water Resources 45 (2012) 26–36

(i) 0.1

1

10

aggregation timescale (hours)

Fig. 15. Comparison of observed and simulated autocorrelations: column 1 is for aggregation timescale of 0.1 h, column 2 for aggregation timescale of 1 h and column 3 for December; the limit (lim) plots indicate the 95% prediction limits (2nd and 39th ranks) for 8CC, daily rainfall depth R P 0, Gregory.

tured by the simulations using most stations’ OSDSPs with differences being less than 0.05. This is with the exception of PER and DAR which underestimated the simulation timescale observed autocorrelations by about 0.1. In general the simulation results improved with increasing aggregation timescale. While the simulated mean results are very good, the wide prediction limits for some of the statistics examined need to be commented upon. This could be a result of the short 12 years of data disaggregated and it is hoped that, as the number of years of daily data to be disaggregated increases, say to 100 years, the prediction limits would narrow. Introducing more rules in the selection of OSDSPs could also help to narrow the prediction limits. However, it could also be that the variability observed is just the nature of sub-tropical/tropical rainfall, and the application of the daily disaggregation model to data from other climatic regions may support, or otherwise, this hypothesis. 5. Summary and conclusions This paper contributes to the potential use of the long-term records of daily rainfall data in hydrological and environmental mod-

elling through disaggregation. It combines the power of copula modelling and the use of simple rules to select representative observed scaled daily storm profiles (OSDSPs) from records of rainfall stations across Australia. The properties of interest for a daily storm profile have been outlined as the total depth (R), total wet periods’ duration (L), the distribution of wet periods, and the distribution of R over the wet periods. For a given R, the normal copula model is used to determine the corresponding value of L. It has been observed that the normal copula model parameters exhibit a clear seasonality which has been captured by a first harmonic Fourier series, thus considerably reducing the number of parameters. An important question addressed in this paper was whether a wet day’s storm profile scaled by the daily depth at a given site could be applied to another wet day of the same value, or close, to L at a different climatic site. This question was examined by grouping the observed daily data into classes according to the value of L before analysing the groups at carefully chosen aggregation levels. A major advantage in analysing the data on a daily basis is the use of every good quality wet day data irrespective of whether there are missing data within that same month which would otherwise be discarded. The OSDSPs are classified into four types

36

Y. Gyasi-Agyei / Advances in Water Resources 45 (2012) 26–36

with respect to whether the first and/or the last periods experienced rain. For two consecutive Julian wet days, a transition probability is introduced to determine whether the rain on the first wet day will cross midnight and spread to the next wet day. For each wet day a pool of the selected type of OSDSP having L values equal, or close, to that obtained from the bivariate copula is created. Finally, the OSDSP in the pool that has the closest value of lag-1 autocorrelation of the natural logarithm of the wet period depths to a simulated value conditioned on the total wet periods’ duration is selected, and multiplication of the OSDSP by the daily depth gives the complete storm profile. The developed daily rainfall disaggregation model was tested by applying it to 12 years continuous rainfall data observed at an erosion control experimental site, comparing several observed and simulated statistics on a monthly basis. A good prediction of the total wet periods’ duration was observed, indicating the adequacy of the copula model for the covariance structure. This was also true for the dependence structure of the lag-1 autocorrelation of the natural logarithm of wet period depths. The disaggregation model nearly perfectly reproduced the dry probabilities irrespective of which site’s OSDSPs used. The variance was also very well predicted for all months except for June and July. With respect to the reproduction of the extreme values the performance of most capital cities’ OSDSPs are also very good. For the autocorrelations, nearly all the observed values were within the prediction range, with the results of Perth and Darwin OSDSPs showing a slight departure from the observed. This may be due to their short record length with limited OSDSPs to choose from. Overall the performance of the disaggregation is very good, bearing in mind the short number of years of the case study site data. It is hoped that, as the number of years of daily data to be disaggregated increases, the prediction range would narrow. The modelling philosophy presented will be very useful where fine timescale data are very rare, with the possibility that OSDSPs could be transferred across continents. The full potential of the presented model will be realised after testing on other climatic sites within and outside Australia.

Acknowledgements Continuing financial support by Queensland Rail and QR National to the author through the HEFRAIL Project is gratefully acknowledged. Comments by anonymous reviewers were very helpful in revising the paper.

References [1] Bardossy A, Pegram GGS. Copula based multisite model for daily precipitation simulation. Hydrol Earth Syst Sci 2009;13:2299–314. [2] Bo Z, Islam S, Eltahir EAB. Aggregation–disaggregation properties of a stochastic rainfall model. Water Resour Res 1994;30(12):3423–35.

[3] De Michele C, Salvadori G. A Generalized Pareto intensity–duration model of storm rainfall exploiting 2-copulas. J Geophys Res 2003;108(D2):4067. doi:10.1029/2002JD002534. [4] Econopouly TW, Davis DR, Woolhiser DA. Parameter transferability for a daily rainfall disaggregation model. J Hydrol 1990;118:209–28. [5] Favre AC, Adlouni SE, Perreault L, Thiémonge N, Bobée B. Multivariate hydrological frequency analysis using copulas. Water Resour Res 2004;40:W01101. [6] Genest C, Favre AC. Everything you always wanted to know about copula modelling but were afraid to ask. J Hydrol Eng 2007;12:347–68. [7] Glasbey CA, Cooper G, McGechan MB. Disaggregation of daily rainfall by conditional simulation from a point-process model. J Hydrol 1995;165:1–9. [8] Güntner A, Olsson J, Calver A, Gannon B. Cascade-based disaggregation of continuous rainfall time series: the influence of climate. Hydrol Earth Syst Sc 2001;5(2):145–64. [9] Gyasi-Agyei Y. Identification of regional parameters of a stochastic model for rainfall disaggregation. J Hydrol 1999;223(3–4):148–63. [10] Gyasi-Agyei Y, Sibley J, Ashwath N. Quantitative evaluation of strategies for erosion control on a railway embankment batter. Hydrol Process 2001;15:3249–68. [11] Gyasi-Agyei Y. Stochastic disaggregation of daily rainfall into one-hour time scale. J Hydrol 2005;309:178–90. [12] Gyasi-Agyei Y. Erosion risk assessment of controlled burning of grasses established on steep slopes. J Hydrol 2006;317:276–90. [13] Gyasi-Agyei Y, Mahbub PB. A stochastic model for daily rainfall disaggregation into fine time scale for a large region. J Hydrol 2007;347:358–70. [14] Gyasi-Agyei Y. Copula-based daily rainfall disaggregation model. Water Resour Res 2011;47. doi:10.1029/2011WR010519. W07535, 17 pp. [15] Hershenhorn J, Woolhiser DA. Disaggregation of daily rainfall. J Hydrol 1987;95:299–322. [16] Joe H. Multivariate models and dependence concepts. London: Chapman and Hall; 1997. [17] Kojadinovic I, Yan J. Modeling multivariate distributions with continuous margins using the copula R package. J Stat Softw 2010;34(9):1–20. [18] Koutsoyiannis D, Onof C. Rainfall disaggregation using adjusting procedures on a Poisson cluster model. J Hydrol 2001;246:109–22. [19] Molnar P, Burlando P. Preservation of rainfall properties in stochastic disaggregation by a simple random cascade model. Atmos Res 2005;77:137–51. [20] Nelson RB. An introduction to copulas. New York: Springer-Verlag; 1999. [21] Remillard B, Scaillet O. Testing for equality between two copulas. J Multivariate Anal 2009;100(3):377–86. [22] Salvadori G, De Michele C. Analytical calculation of storm volume statistics involving Pareto-like intensity-duration marginals. Geophys Res Lett 2004;31:L04502. doi:10.1029/2003GL018767. [23] Salvadori G, De Michele C, Kottegoda N, Rosso R. Extremes in nature: an approach using copulas. Springer; 2007. [24] Scaillet O. A Kolmogorov–Smirnov type test for positive quadrant dependence. Can J Stat 2005;33:415–27. [25] Serinaldi F. A multisite daily rainfall generator driven by bivariate copulabased mixed distributions. J Geophys Res – Atmos 2009:114. doi:10.1029/ 2008JD011258. [26] Serinaldi F. Multifractality, imperfect scaling and hydrological properties of rainfall time series simulated by continuous universal multifractal and discrete random cascade models. Nonlin Process Geophys 2010;17:697–714. doi:10.5194/npg-17-697-2010. [27] Vandenberghe S, Verhoest NEC, Buyse E, De Baets B. A stochastic design rainfall generator based on copulas and mass curves. Hydrol Earth Syst Sci 2010;7:3613–48. [28] Vandenberghe S, Verhoest NEC, De Baets B. Fitting bivariate copulas to the dependence structure between storm characteristics: a detailed analysis based on 105 year 10 min rainfall. Water Resour Res 2010;46:W01512. doi:10.1029/ 2009WR007857. [29] Zhang L, Singh VP. Gumbel–Hougaard copula for trivariate rainfall frequency analysis. J Hydrol Eng 2007;12:409–19. [30] Zhang L, Singh VP. Bivariate rainfall frequency distributions using Archimedean copulas. J Hydrol 2007;332(1–2):93–109.