- Email: [email protected]
Identification of regional parameters of a stochastic model for rainfall disaggregation
Journal of Hydrology 223 (1999) 148–163 www.elsevier.com/locate/jhydrol
Identification of regional parameters of a stochastic model for rainfall disaggregation Y. Gyasi-Agyei Centre for Railway Engineering, James Goldston Faculty of Engineering and Physical Systems, Central Queensland University, Rockhampton QLD4702, Australia Received 22 June 1998; accepted 14 July 1999
Abstract This paper demonstrates how the Gyasi-Agyei–Willgoose hybrid model for point processes could be regionalised for daily rainfall disaggregation using limited high resolution data within a region of interest. Their model is a product of the binary nonrandomised Bartlett–Lewis rectangular pulse model and a lognormal autoregressive model used as a jitter. The computationally efficient multinormal approximation to parameter uncertainty is used to group the monthly parameter values of the binary model. For central Queensland, Australia, it has been established that the parameters of cell origins and the duration of the rectangular pulse of the binary model could have constant values for all months. Second harmonic Fourier series is used to represent the seasonal variation of the parameter governing the storm lifetime. The storm arrival rate is a function of the daily dry probability and the other parameters. Additive properties of random variables with finite variances were used to scale down the daily mean and variance of the historical data to the simulation timescale, values required by the jitter model. The results of using observed daily rainfall statistics to capture sub-daily statistics by the regionalised model are very encouraging. The model is therefore a valuable tool for disaggregating daily rainfall data. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Rainfall; Disaggregation; Stochastic point process; Regionalisation; Parameter uncertainty; Hybrid model
1. Introduction The design of strategies to control runoff and erosion often requires information on individual rain storms at a fine timescale (1 min for example). In particular, for regions where few rain storms of short duration within a year account for a large proportion of runoff and erosion damage. The need for minute-by-minute information on short duration hydrological processes increases where slopes are steep, as on railway or highway batters. Unfortunately, fine timescale rainfall data are limited worldwide due to the high cost of obtaining and storing such E-mail address: [email protected] (Y. Gyasi-Agyei)
data. Usually, available long records of rainfall data are on a daily timescale. The need for a model to disaggregate daily rainfall into a sequence of individual storms of finer timescale cannot be overemphasised. Nowadays most Australian farmers and mining companies monitor and store daily rainfall data, increasing the density of daily rain gauges without a significant cost. Bo et al. (1994) used the modified Bartlett–Lewis model (Rodriguez-Iturbe et al., 1988) to capture the statistics of finer timescale rainfall from the observed daily rainfall statistics. Their approach involved using daily statistics to estimate the model parameters and then simulate the sequence of rainfall events at any desired timescale. Glasbey et al. (1995) also used the
0022-1694/99/$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(99)00114-6
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
modified Bartlett–Lewis model to disaggregate daily rainfall by conditional simulation. They calibrated the model using statistics of hourly data, without taking into account the seasonality in the data. The model was used to generate an archive of a large number of years of rainfall data. Then, given daily rainfall data, the best matches to sequences of daily totals in the archive, and the hourly data from those past events were used to provide disaggregated sequences. As they pointed out, the disaggregated data are constrained to follow only the same trend in the mean rainfall. Therefore, the model could not necessarily predict the second-order characteristics and the dry probability. For the rainfall data within central Queensland, the Bartlett–Lewis models failed to adequately reproduce statistics of the finer timescale, even where these statistics were included in the model parameter calibration (Gyasi-Agyei and Willgoose, 1997). Cowpertwait et al. (1996a) presented a fitting procedure for using Newman–Scott model (Rodriguez-Iturbe et al., 1987 and Cowpertwait, 1991) to disaggregate daily rainfall data. When only daily aggregated properties were used to fit the model, they found that the 1-h aggregation level statistics were poorly reproduced. They underlined that since the parameters of the model relate to rain cells that generally have lifetimes less than 1 h, higher aggregation levels are unlikely to contain enough information from which the properties of the cells can be found. Cowpertwait et al. (1996a) then estimated the sample hourly variance from the sample daily variance using a derived empirical relationship. These sub-daily variance values together with the daily second-order moments and the wet and dry transition probabilities were used to estimate the five parameters of the model for each calendar month. The objectives of this paper are: 1. to use a stochastic point process model to capture the statistics of finer timescale rainfall from readily available observed daily rainfall statistics; 2. to minimise the number of parameters that need to be jointly calibrated by incorporating more realistic seasonal variation than monthly grouping of data; 3. to determine a parameter set for the region of interest.
149
A hybrid point process model presented in GyasiAgyei and Willgoose (1997, 1999) is used to achieve the above objectives. The numerics of the hybrid model are summarised in Section 2 and the data used are presented in Section 3. This is followed by the parameter identification procedure of the binary model in Section 4. Additive properties of random variables with finite variances used to scale down the daily mean and variance of the historical data are presented in Section 5. Section 6 evaluates the disaggregation methodology. The paper is concluded in Section 7.
2. Gyasi-Agyei–Willgoose hybrid model The hybrid model {Y t} is a product of a binary chain model {Y t} and an autocorrelated jitter model {A t} expressed as {Y t} {A t}{Y t}
1
where A t expZ t and {Z t} is a stationary Gaussian process of mean mZ h; variance s2Z h and lag-t autocovariance function cZ h; t; h is the timescale and t is time. Further, it is assumed that {Z t} is independent of {Y t}: The exponential form of the jitter process {A t} is necessary to ensure positivity 0 only if Y t 0: This ensures that and that Y t the probability of dry periods remain the same for both {Y t} and {Y t}: Choosing mA h 1; that is making mY h mY h; 2mZ h 2s2Z h; and the second order properties of {Z t} have been derived as (Gyasi-Agyei and Willgoose, 1997) " 2 # sY h 1 m2Y h 2 sZ h ln 2 2 sY h 1 m2Y h} and "
c h; t 1 m2Y h cZ h; t ln Y cY h; t 1 m2Y h
# 3
The random process {Z t} is modelled by a firstorder autoregressive model given as Zt h mZ h 1 rZ h; 1 Zt21 h 2 mZ h 1 Et 4 where rZ h; 1 is the lag-1 autocorrelation and {Et } is normally distributed, with a mean value of zero and
150
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Fig. 1. The locations of pulviograph and daily rain gauge stations within central Queensland, Australia: each of the daily rain gauge stations has over 50 years of data.
variance s2E h given by
s2E h 1 2 rZ h; 12 s2Z h
(Gyasi-Agyei and Willgoose, 1999) 5
The condition urZ h; 1u , 1 must be satisfied. The binary chain model generates a string of two numbers: Yi 0 for a dry period and a constant value Yi w h for a wet period, where Yi is the cumulative rainfall depth over time interval i. The moments of the binary chain have been derived analytically as
w h
mY h 1 2 P h
6
the variance s2Y h as
s2Y h mY h2
P h 1 2 P h
7
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
151
Table 1 Number of years of data Station
Month
ID
Name
No.
1
2
3
4
5
6
7
8
9
10
11
12
Min
Max
33087 33099 33119 34002 35025 35069 35098 35104 39006 39083 39090 39104 39123 Min Max
Barganal Teemburra Mackay Charters Dingo Tambo Emerald Kilmacolm Biloela Rockhampton Theodore Monto Gladstone
1 2 3 4 5 6 7 8 9 10 11 12 13
20 21 24 20 14 19 15 15 44 29 24 22 14 14 44
27 21 17 11 15 20 15 14 41 23 22 17 11 11 41
19 20 18 24 15 23 16 16 37 31 21 15 10 10 37
21 20 16 20 16 22 18 10 43 39 26 21 23 10 43
19 18 14 19 13 21 20 12 45 38 24 14 25 12 45
24 12 20 23 16 21 14 10 44 42 23 21 26 10 44
22 14 13 17 13 24 19 10 39 35 21 17 19 10 39
25 13 16 21 14 24 20 13 41 43 26 18 21 13 43
21 12 19 14 15 26 14 11 40 46 23 24 21 11 46
23 16 23 22 16 23 23 15 41 43 24 24 18 15 43
25 17 28 26 17 25 18 16 48 45 26 20 17 16 48
24 19 26 20 15 23 19 16 45 41 22 17 16 15 45
19 12 13 11 13 19 14 10 37 23 21 14 10
27 21 28 26 17 26 23 16 48 46 26 24 26
where P h is the probability that an interval i of duration h hours is dry, and lag-1 joint probabilities are defined as
g (h 21). It is assumed that one cell arrives at the storm origin. The duration of the rectangular pulse associated with each is exponentially distributed with rate h (h 21). Rectangular pulses of cells are allowed to overlap with cells of the same storm and across cells of different storms. The sequences of the random variables defining the process are assumed mutually independent. For convenience, dimensionless parameters k b=h and f g=h are introduced. The probability that an interval i of duration h hours is dry, P h; is given by (Rodriguez-Iturbe et al., 1987)
P00 h P Yi11 0; Yi 0 P 2h
P h exp2 h 1 uT 2 Cl
P01 h P Yi11 . 0; Yi 0 P h 2 P 2h
where
P11 h P Yi11 . 0; Yi . 0 1 2 2P h 1 P 2h
C
and the lag-1 autocovariance cY h; 1 as " ( cY h; 1 mY h
2
P h P00 h; 1 2 2P01 h; 1 1 2 P h "
P h 1 P11 h; 1 1 2 P h
#
#2 ) 8
9 In Eq. (9), P 2h is the probability that two consecutive time intervals are dry. In this paper, the non-randomised Bartlett–Lewis model (Rodriguez-Iturbe et al., 1987) is used to generate the binary chain. For the non-randomised Bartlett–Lewis model, storms arrive in a Poisson Process of rate l (h 21). Followed by each storm origin is the arrival of cell origins that are also governed by a Poisson Process but of rate b (h 21). Cell arrivals cease after a time, exponentially distributed with rate
10
GP f 1 k e2 k1fhh k1f
11
The expected duration of a storm, mT, is given by Z1 Z1 dn dtn21 tf21 uT f21 h21 1 fh21 0
0
1 2 1 2 nt e2kn 12t and the function GP is defined as Z1 GP h21 e2k tf21 1 2 t ekt dt 0
12
13
Eqs. (12) and (13) are solved numerically. The time
152
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Fig. 2. Stage 1 calibrated parameter values with a band of twice the standard deviation, Rockhampton station.
series generated by the non-randomised Bartlett– Lewis model is represented by a string of two numbers: zero for a dry period and a constant value, w h; for a wet period.
3. Data Six-minute continuous rainfall data from 13 stations within central Queensland, Australia, were used in this paper. The location of the pluviograph stations together with daily rainfall stations having over 50 years of data is shown in Fig. 1. The altitude of these stations varies between 10 and 400 m above
mean sea level. The rainfall data, some dating back to 1939, were digitised from syphon (float) type rain gauge electronic tracing charts by the Bureau of Meteorology, Melbourne, Australia. For a given year, months with missing data were excluded. This resulted in rainfall data of different months for the same station having a different number of years of record. Table 1 shows the number of years of data of the various months of a given rainfall station used. Only stations with at least 10 years of good quality data for each month were used. The data were analysed at 20 aggregation levels: 6, 12, 18, 24, 30, 36, 42, 48 and 54 min and 1, 1.5, 2, 3, 4, 6, 8, 10, 12, 18 and 24 h.
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
153
Table 2 Correlation coefficient between parameter pairs–stage 1 Month
l –b
l –g
l –h
b –g
b –h
g –h
January February March April May June July August September October November December
0.59790 0.48320 0.65910 0.60670 0.58730 0.58070 0.58540 0.56940 0.56940 0.57170 0.57700 0.60150
0.93890 0.94380 0.95970 0.93790 0.94250 0.95330 0.94350 0.93050 0.90140 0.90030 0.88270 0.92170
0.39570 0.39300 0.46270 0.40820 0.39780 0.39400 0.39600 0.38440 0.37980 0.37620 0.38870 0.41170
0.71820 0.65190 0.77340 0.73520 0.66170 0.62710 0.65600 0.61680 0.71630 0.73600 0.78400 0.74560
0.92940 0.93780 0.93630 0.93290 0.93950 0.94140 0.93950 0.94120 0.93670 0.93170 0.93730 0.93900
0.46980 0.40640 0.55340 0.49570 0.41970 0.38920 0.41380 0.36740 0.46890 0.48530 0.55290 0.51250
4. Identification of regional parameters of the binary model The binary Bartlett–Lewis model has four parameters (l , b , g and h ) to be calibrated for each month, resulting in 48 parameters for the whole year for a given rainfall station. The parameters of the model are examined for uniformity over the months and the possibility of fitting a curve in order to reduce the number of parameters. This may provide a smooth representation of the parameters with a view to capturing the seasonal variability of the data. Moreover, regionalisation is easier with a reduced number of parameters. In modelling rainfall over Britain using the random Bartlett–Lewis model, Onof and Wheater (1993) fitted polynomial curves of degree varying between 5 and 12 to the individual parameters. Their objective was to represent the parameters as daily variables in order to avoid the artificiality of the monthly subdivision of the year and allowing for possible generalisation of parameters through regional analysis. The parameters of the binary model were identified through four stages using rainfall data of Rockhampton. At each stage the parameters were calibrated by minimising the objective function, J, J
X
Po h 2 Ps h2
14
where Po(h) and Ps(h) are the historical and analytical dry probabilities of aggregation level h, respectively. The shuffled complex evolution global probabilistic
search strategy (Duan et al., 1992) option of the Bayesian non-linear regression software NLFIT (Kuczera, 1994) was used to estimate the parameters. The number of complexes was set equal to the number of fitted parameters. The global optimum parameter set was then used as starting values for a local gradient search. NLFIT uses the computationally efficient classical approximations to parameter uncertainty based on the multinormal approximation of the posterior distribution to estimate the standard deviations and correlations of the parameter sets. The covariance matrix of the multinormal approximation can help identify the poorly defined parameters and their interactions (Kuczera and Mroczkowski, 1998). In stage one, the four parameters were calibrated for each month. Fig. 2 depicts the values of the calibrated parameters, with a band of twice the standard deviation. It is observed that only parameter l is well identified with a clear seasonal pattern, high values for the wet season (October to March) and low values for the dry season. Parameter g also shows a seasonal trend. However, the envelopes of twice the standard deviation of the parameters b and h suggest the possibility of adopting constant values for these parameters over the months without undue loss of accuracy. The very high correlations between the parameter pairs b and h (Table 2) suggest nearly indistinguishable simulation results from parameters spanning a wide range of parameter space. Thus only one is assumed to have a constant value for all months in stage two of the parameter identification process. During stage two of the parameter identification
154
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163 Table 3 Calibrated parameter values of the binary model
Fig. 3. Stage 2 calibrated parameter h values with a band of twice the standard deviation, Rockhampton station.
process, parameter b was considered constant for all months, reducing the number of parameters calibrated to 37. Since the correlation between the parameter pairs l – b and g – b in stage one was not very high, there were slight changes in the values of the parameters l and g . The stage two calibrated values of h indicated narrow fluctuations with envelopes of two standard deviations (Fig. 3), yet a constant value could be assumed. Therefore in the third stage of the parameter identification process, parameters b and h were assumed to have constant values for all months,
Parameter (h 21)
Rockhampton
Regional
b h a0 a1 b1 a2 b2
0.7095 1.3721 0.1407 0.0432 2 0.0428 2 0.0094 2 0.0369
0.4808 1.3175 0.1180 0.0316 2 0.0335 2 0.0026 2 0.0226
reducing the number of calibrated parameters to 26. While there were changes in the values of the parameters l and g , the trend of their variation with the month of the year remained largely unchanged. What emerged at this stage was a strong logarithmic relationship between parameter l and the dry probabilities. This is as a result of Eq. (10) for given values of b , g and h . Moreover, the correlation between the parameter pair l – g was in excess of 83% for all months, an indication of multiple optima within the feasible parameter space. In stage 4, parameters b and h were still considered constant for all months and second harmonic Fourier series
g m a0 1
2 X ak cos k2pm=12 1 bk sin k2pm=12 k1
15
Fig. 4. Fitted second harmonic Fourier series of parameter g .
was used to represent the variation of parameter g over the months, m. In Eq. (15) a0, ak and bk are coefficients to be determined. For given parameter values of b , g and h , and the daily dry probability, P(24), Eq. (10) was used to estimate parameter l . Fig. 4 shows second harmonic Fourier series fitted to parameter g and Table 3 gives the calibrated parameter values at stage 4. Essentially, the number of calibrated parameters in stage 4 was 7. Eq. (10) and the second harmonic Fourier series for parameter g capture the seasonality effects in the sequence of dry and wet spells. The binary model error in reproducing the dry probabilities through the four stages was examined using two error statistics: the average absolute percentage error E1 and the maximum absolute percentage error
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
155
Fig. 5. Comparison of error statistics E1 and E2 during stages 1–4, Rockhampton station.
E2, defined as E1 %
n 1 X Po k 2 Ps k 100 n k1 Po k
P k 2 Ps k 100 E2 % max o k1;n Po k
16
17
where Po(k) and Ps(k) are as defined in Eq. (14) and n is the number of values. If the error statistics are estimated for a given month then n is 20 the number of aggregation levels, but if they are for a given aggregation level then n is 12 the number months. The error statistics are presented in Fig. 5. Apparently, there were no significant differences between the error statistics from stages 1–3. Stage 4 error statistics
were slightly higher than the values of the other stages. However, for all intents and purposes the stage 4 fit was very good since the error statistics were very small. The representativeness of Rockhampton’s calibrated parameter values (Table 3) for central Queensland region was investigated. These calibrated values and the P(24) values of the 12 other stations were used to generate the analytical dry probabilities of the 20 aggregation levels. Fig. 6 shows the error statistics of the 12 stations for the individual months and also for the aggregation levels, error statistic E1 being averaged over the 12 stations and E2 the maximum of the 12 stations. Clearly, the error statistics were small (all less than 1.4% for E1 and less than 7% for E2), the dry months having the least values. In Fig. 7 is shown the
156
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Fig. 6. Variation of error statistics E1 and E2 over the months and aggregation levels.
error statistics for the individual stations, estimated over the aggregation levels and months. Next, all 13 stations were jointly calibrated for a representative regional parameter set (see Table 3 and Fig. 4). The Rockhampton and the regional parameter sets were different, but the average error statistics were not significantly different (Figs. 6 and 7). It is therefore concluded that a set of parameter values could be used to disaggregate daily dry probabilities to finer aggregation levels for all stations within the region. The values of the first three stations were highest, suggesting the possibility of clustering since these stations are grouped geographically. To investigate this issue the stations were divided into two groups separated by latitude 238 (Fig. 1). Joint calibration was carried out for each group. Results not shown here
indicated a slight reduction in the error statistics in the average sense, but too small to warrant grouping the stations into clusters.
5. Estimation of sub-daily historical moments In order to apply the jitter model, the historical moments at the simulation timescale are required. The contained analysis is a practical solution for estimating rainfall statistics for different timescales. This analysis is not a consequence of the hybrid model discussed in the previous sections and an alternative approach could be adopted. The rainfall depths of timescale h, Y k h; are
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
157
Fig. 7. Variation of error statistics E1 and E2 over the stations.
connected with the daily value Y 24 by Y 24 Y 1 h 1 Y 2 h 1 … 1 Y n h
18
Fig. 8. Variation of 1-h timescale lag-1 autocorrelation over the months (top) and of 1-h timescale lag-2 autocorrelation with 1-h timescale lag-1 autocorrelation (bottom).
where n 24=h: From the additive properties of random variables with finite variances
mY 24
n X
mYk h nmY h
19
k1
and
s2Y 24
n X k1
s2Yk h 1 2
X
h i cov Y j h; Y k h
20
j,k
Fig. 8 (top) shows the 1-h lag-1 autocorrelation, rY 1; 1; values of all station-months. In order to reduce the number of parameters, the monthly mean
Table 4 Parameter values of the second harmonic Fourier series fitted to 1-h lag-1 autocorrelation, rY 1; 1 Parameter
Value
f0 f1 g1 f2 g2
0.4189 2 0.1038 0.0312 0.0375 0.0476
158
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
be simplified to
mY 24 24mY 1
22
and Eq. (20) to
s2Y 24 24 1 2 24 2 1rY 1; 1 1 erY 1; 12 s2Y 1 23
Fig. 9. Observed and predicted 1-h variance and lag-1 autocovariance. Solid line is a line of perfect match.
values (open circle) were fitted with second harmonic Fourier series
rY 1; 1 m f0 1
2 X
fk cos k2pm=12
where the third term on the right is an error term accounting for the higher lag autocovariances. Parameter e varied over the station-months without apparent regional trend. However, a constant regional value of parameter e was sought in addition to the rY 1; 1 regional function to fit s2Y 1: The s2Y 24 values of 14 station-months, mostly from the summer months of stations 33087, 33099 and 33119, exceeded 200 mm 2 and required a separate parameter e value. The parameter e was estimated as e 201:5; s2Y 24 # 200 and e 670:5; s2Y 24 . 200 mm2 : Fig. 9 compares the observed and predicted values of 1-h variance and lag-1 autocovariance. It is observed in Fig. 9 that the predicted second-order moments matched the observed values very well. The 14 station-months have exceptionally high variances and are the wettest of all stationmonths, explaining the huge range between the two constant values of parameter e. A higher value of parameter e indicates the importance of the higher lag autocorrelations. Sub-hourly timescales can be similarly handled. For the 6-min-timescale, it was observed that rY 1; 1 was nearly constant with an average value of 0.8. However, the variation of parameter e over the station-months was very high. In such circumstances it is recommended to work with hourly values and distribute the hourly rainfall depths over the wet sub-hourly intervals as explained in the next section.
k1
1gk sin k2pm=12
21
where f0 , fk and gk are coefficients given in Table 4. The higher lag autocorrelations of the 1-h timescale were observed to increase with the square of the lag-1 autocorrelation, as shown in Fig. 8 (bottom) for lag-2. Therefore, for 1-h simulation timescale Eq. (19) can
6. Evaluation of the regionalised hybrid model The procedure for estimating the parameters of the hybrid model when only daily data are available is summarised as follows: 1. calculate mY 24, s2Y 24 and P(24) for each month;
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
159
Fig. 10. Comparison of historical (shaded circle) and simulated (dashed line) dry probability.
2. use Eq. (21) to estimate rY 1; 1; Eq. (22) to estimate mY 1; Eq. (23) to estimate s2Y 1 and then calculate cY 1; 1 rY 1; 1s2Y 1; 3. estimate parameter l in Eq. (10) using P(24) and the regional parameters given in Table 3; 4. select the simulation timescale, h, and use Eq. (10) to estimate P h and P 2h; 5. estimate the moments of the binary chain model using Eqs. (6)–(9);
6. estimate the parameters of the jitter model using Eqs. (2), (3) and (5). The binary process can be simulated at any desired timescale but 6-min (the historical data timescale) is used in this paper. However, the jitter model was simulated at 1-h timescale. Analysis of the distribution of the historical hourly rainfall depth over the 6-min wet intervals indicated that the fractional
160
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Fig. 11. Comparison of historical (shaded circle) and simulated (dashed line) mean.
weights vary randomly and could be sampled from a uniform distribution, U 0:1; 1: Hence the fractional weights within a wet hour were drawn from U 0:1; 1 and each divided by their sum in order to preserve a sum of unity. The scaled fractional weights were multiplied by the hourly rainfall depth to obtain the 6-min interval’s rainfall depth. For example, if the rainfall depth within an hour was 50 mm, and there were five 6-min wet periods within that hour, then five
fractional weights 0.405, 0.775, 0.271, 0.152 and 0.410 drawn from U 0:1; 1 were divided by 2.013, their sum. Each of the scaled fractional weights 0.201, 0.385, 0.135, 0.075 and 0.204 was then multiplied by 50 mm to obtain the rainfall depths 10.05, 19.25, 6.75, 3.75 and 10.20 mm of the five 6-min wet periods. This approach preserved the dry probabilities. Figs. 10–13 compare the historical statistics
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
161
Fig. 12. Comparison of historical (shaded circle) and simulated (dashed line) variance.
(shaded circle) with 1000 years of simulation results (dashed line) of 6 station-months, used as examples. These station-months are from December and January, the wettest months, and across the region. The dry probabilities were very well reproduced. It must be underlined that the simulated and analytical dry probabilities were very close. In all cases the historical means were very well matched by the simulated values. The variances and lag-1 autocovariances
reproduced by the hybrid model were in good agreement with the historical ones. Although for some station-months the absolute percentage error of variance and lag-1 autocovariance were in excess of 50% for a given aggregation level, lag-1 autocovariance in particular. Station-month Dec-33099, showing the worst scenario, is one of the station-months with s2Y 24 in excess of 200 mm 2. As indicated by the examples shown, the method of distributing hourly
162
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Fig. 13. Comparison of historical (shaded circle) and simulated (dashed line) lag-1 autocovariance.
rainfall depths over 6-min wet periods during the simulation worked very well. Data size could play a significant role in the errors in the reproduction of the variances and lag-1 autocovariances. For example, Cowpertwait et al. (1996b) established that for a 10-year record of British rainfall data the coefficient of variation of mY 24; s2Y 24 and P(24) was of the order of 16, 31 and 12%, respectively. These values were reduced by about 30%
when the record length was doubled. Hence, with the availability of more data it is hoped that the reproduction of the variance and lag-1 autocovariance by the hybrid model would be improved.
7. Conclusions In this paper it is demonstrated how the
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Gyasi-Agyei–Willgoose hybrid model for point processes could be regionalised for daily rainfall disaggregation using limited high resolution rainfall data within a region of interest. Their model is a product of the binary non-randomised Bartlett– Lewis rectangular pulse model and a lognormal autoregressive model used as a jitter. The time series generated by the non-randomised Bartlett–Lewis model is replaced by a string of two numbers: zero for a dry period and a constant value for a wet period chosen to preserve the mean of the historical rainfall data. The four stages of the parameter identification process of the binary model involve joint calibration of the parameters of all months. The computationally efficient multinormal approximation to parameter uncertainty is used to group the monthly parameter values of the binary model. For central Queensland, Australia, it has been established that the parameters b and h could have constant values for all months. Second harmonic Fourier series was used to represent the seasonal variation of the parameter g . The parameter l is a function of the daily dry probability and the parameters b , h and g . Thus only a minimal number of parameters need to be calibrated making the model more parsimonious. The approach removes the artificiality of grouping the yearly data on a monthly basis. Also the additive properties of random variables with finite variances were used to scale down the daily mean and variance of the historical data to the simulation timescale, values required by the jitter model. The regionalised model is structured so that it requires only the daily values of dry probability, mean and variance, the regional 1-h lag-1 autocorrelation function and parameter e to simulate rainfall at a fine timescale. The results on using observed daily rainfall statistics to capture subdaily statistics by the regionalised model are very encouraging. The model is therefore a valuable tool for disaggregating daily rainfall, data of which abound worldwide. Further research is, however, required to improve the reproduction of the second-order moments.
163
Acknowledgements The reviews by E.N. Anagnostou, T. Lebel and D. Koutsoyiannis were very helpful and are gratefully acknowledged. The author wishes to thank G. Kuczera for providing NLFIT software. Financial support from the Coal and Minerals Infrastructure Division of Queensland Rail, Australia, is very much appreciated. References Bo, Z., Islam, S., Eltahir, E.A.B., 1994. Aggregation–disaggregation properties of a stochastic rainfall model. Water Resour. Res. 30 (12), 3423–3435. Cowpertwait, P.S.P., 1991. Further developments of the Newman– Scott clustered point process for modelling rainfall. Water Resour. Res. 27 (7), 1431–1438. Cowpertwait, P.S.P., O’Connell, P.E., Metcalfe, A.V., Mawdsley, J.A., 1996a. Stochastic point process modelling of rainfall. I. Single-site fitting and validation. J. Hydrol. 175, 17–46. Cowpertwait, P.S.P., O’Connell, P.E., Metcalfe, A.V., Mawdsley, J.A., 1996b. Stochastic point process modelling of rainfall. II. Regionalisation and disaggregation. J. Hydrol. 175, 47–65. Duan, Q., Sorooshian, S., Gupta, V., 1992. Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resour. Res. 28 (4), 1015–1031. Glasbey, C.A., Cooper, G., McGehan, M.B., 1995. Disaggregation of daily rainfall by conditional simulation from a point-success model. J. Hydrol. 165, 1–9. Gyasi-Agyei, Y., Willgoose, G.R., 1997. A hybrid model for point rainfall modelling. Water Resour. Res. 33 (7), 1699–1706. Gyasi-Agyei, Y., Willgoose, G.R., 1999. Generalisation of a hybrid model for point rainfall. J. Hydrol. 219, 218–224. Kuczera, G., 1994. NLFIT: a Bayesian nonlinear regression program suite, Department of Civil Engineering and Survey, the University of Newcastle, Callaghan, NSW, Australia, 162pp. Kuczera, G., Mroczkowski, M., 1998. Assessment of hydrologic parameter uncertainty and the worth of multiresponse data. Water Resour. Res. 34 (6), 1481–1489. Onof, C., Wheater, H.S., 1993. Modelling of British rainfall using a random parameter Barlett–Lewis rectangular pulse model. J. Hydrol. 149, 67–95. Rodriguez-Iturbe, I., Cox, D.R., Isham, V., 1987. Some models for rainfall based on stochastic point processes. Proc. R. Soc. London A 410, 269–288. Rodriguez-Iturbe, I., Cox, D.R., Isham, V., 1988. A point success model for rainfall: further developments. Proc. R. Soc. London A 417, 283–298.
Identification of regional parameters of a stochastic model for rainfall disaggregation Y. Gyasi-Agyei Centre for Railway Engineering, James Goldston Faculty of Engineering and Physical Systems, Central Queensland University, Rockhampton QLD4702, Australia Received 22 June 1998; accepted 14 July 1999
Abstract This paper demonstrates how the Gyasi-Agyei–Willgoose hybrid model for point processes could be regionalised for daily rainfall disaggregation using limited high resolution data within a region of interest. Their model is a product of the binary nonrandomised Bartlett–Lewis rectangular pulse model and a lognormal autoregressive model used as a jitter. The computationally efficient multinormal approximation to parameter uncertainty is used to group the monthly parameter values of the binary model. For central Queensland, Australia, it has been established that the parameters of cell origins and the duration of the rectangular pulse of the binary model could have constant values for all months. Second harmonic Fourier series is used to represent the seasonal variation of the parameter governing the storm lifetime. The storm arrival rate is a function of the daily dry probability and the other parameters. Additive properties of random variables with finite variances were used to scale down the daily mean and variance of the historical data to the simulation timescale, values required by the jitter model. The results of using observed daily rainfall statistics to capture sub-daily statistics by the regionalised model are very encouraging. The model is therefore a valuable tool for disaggregating daily rainfall data. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Rainfall; Disaggregation; Stochastic point process; Regionalisation; Parameter uncertainty; Hybrid model
1. Introduction The design of strategies to control runoff and erosion often requires information on individual rain storms at a fine timescale (1 min for example). In particular, for regions where few rain storms of short duration within a year account for a large proportion of runoff and erosion damage. The need for minute-by-minute information on short duration hydrological processes increases where slopes are steep, as on railway or highway batters. Unfortunately, fine timescale rainfall data are limited worldwide due to the high cost of obtaining and storing such E-mail address: [email protected] (Y. Gyasi-Agyei)
data. Usually, available long records of rainfall data are on a daily timescale. The need for a model to disaggregate daily rainfall into a sequence of individual storms of finer timescale cannot be overemphasised. Nowadays most Australian farmers and mining companies monitor and store daily rainfall data, increasing the density of daily rain gauges without a significant cost. Bo et al. (1994) used the modified Bartlett–Lewis model (Rodriguez-Iturbe et al., 1988) to capture the statistics of finer timescale rainfall from the observed daily rainfall statistics. Their approach involved using daily statistics to estimate the model parameters and then simulate the sequence of rainfall events at any desired timescale. Glasbey et al. (1995) also used the
0022-1694/99/$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S0022-169 4(99)00114-6
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
modified Bartlett–Lewis model to disaggregate daily rainfall by conditional simulation. They calibrated the model using statistics of hourly data, without taking into account the seasonality in the data. The model was used to generate an archive of a large number of years of rainfall data. Then, given daily rainfall data, the best matches to sequences of daily totals in the archive, and the hourly data from those past events were used to provide disaggregated sequences. As they pointed out, the disaggregated data are constrained to follow only the same trend in the mean rainfall. Therefore, the model could not necessarily predict the second-order characteristics and the dry probability. For the rainfall data within central Queensland, the Bartlett–Lewis models failed to adequately reproduce statistics of the finer timescale, even where these statistics were included in the model parameter calibration (Gyasi-Agyei and Willgoose, 1997). Cowpertwait et al. (1996a) presented a fitting procedure for using Newman–Scott model (Rodriguez-Iturbe et al., 1987 and Cowpertwait, 1991) to disaggregate daily rainfall data. When only daily aggregated properties were used to fit the model, they found that the 1-h aggregation level statistics were poorly reproduced. They underlined that since the parameters of the model relate to rain cells that generally have lifetimes less than 1 h, higher aggregation levels are unlikely to contain enough information from which the properties of the cells can be found. Cowpertwait et al. (1996a) then estimated the sample hourly variance from the sample daily variance using a derived empirical relationship. These sub-daily variance values together with the daily second-order moments and the wet and dry transition probabilities were used to estimate the five parameters of the model for each calendar month. The objectives of this paper are: 1. to use a stochastic point process model to capture the statistics of finer timescale rainfall from readily available observed daily rainfall statistics; 2. to minimise the number of parameters that need to be jointly calibrated by incorporating more realistic seasonal variation than monthly grouping of data; 3. to determine a parameter set for the region of interest.
149
A hybrid point process model presented in GyasiAgyei and Willgoose (1997, 1999) is used to achieve the above objectives. The numerics of the hybrid model are summarised in Section 2 and the data used are presented in Section 3. This is followed by the parameter identification procedure of the binary model in Section 4. Additive properties of random variables with finite variances used to scale down the daily mean and variance of the historical data are presented in Section 5. Section 6 evaluates the disaggregation methodology. The paper is concluded in Section 7.
2. Gyasi-Agyei–Willgoose hybrid model The hybrid model {Y t} is a product of a binary chain model {Y t} and an autocorrelated jitter model {A t} expressed as {Y t} {A t}{Y t}
1
where A t expZ t and {Z t} is a stationary Gaussian process of mean mZ h; variance s2Z h and lag-t autocovariance function cZ h; t; h is the timescale and t is time. Further, it is assumed that {Z t} is independent of {Y t}: The exponential form of the jitter process {A t} is necessary to ensure positivity 0 only if Y t 0: This ensures that and that Y t the probability of dry periods remain the same for both {Y t} and {Y t}: Choosing mA h 1; that is making mY h mY h; 2mZ h 2s2Z h; and the second order properties of {Z t} have been derived as (Gyasi-Agyei and Willgoose, 1997) " 2 # sY h 1 m2Y h 2 sZ h ln 2 2 sY h 1 m2Y h} and "
c h; t 1 m2Y h cZ h; t ln Y cY h; t 1 m2Y h
# 3
The random process {Z t} is modelled by a firstorder autoregressive model given as Zt h mZ h 1 rZ h; 1 Zt21 h 2 mZ h 1 Et 4 where rZ h; 1 is the lag-1 autocorrelation and {Et } is normally distributed, with a mean value of zero and
150
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Fig. 1. The locations of pulviograph and daily rain gauge stations within central Queensland, Australia: each of the daily rain gauge stations has over 50 years of data.
variance s2E h given by
s2E h 1 2 rZ h; 12 s2Z h
(Gyasi-Agyei and Willgoose, 1999) 5
The condition urZ h; 1u , 1 must be satisfied. The binary chain model generates a string of two numbers: Yi 0 for a dry period and a constant value Yi w h for a wet period, where Yi is the cumulative rainfall depth over time interval i. The moments of the binary chain have been derived analytically as
w h
mY h 1 2 P h
6
the variance s2Y h as
s2Y h mY h2
P h 1 2 P h
7
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
151
Table 1 Number of years of data Station
Month
ID
Name
No.
1
2
3
4
5
6
7
8
9
10
11
12
Min
Max
33087 33099 33119 34002 35025 35069 35098 35104 39006 39083 39090 39104 39123 Min Max
Barganal Teemburra Mackay Charters Dingo Tambo Emerald Kilmacolm Biloela Rockhampton Theodore Monto Gladstone
1 2 3 4 5 6 7 8 9 10 11 12 13
20 21 24 20 14 19 15 15 44 29 24 22 14 14 44
27 21 17 11 15 20 15 14 41 23 22 17 11 11 41
19 20 18 24 15 23 16 16 37 31 21 15 10 10 37
21 20 16 20 16 22 18 10 43 39 26 21 23 10 43
19 18 14 19 13 21 20 12 45 38 24 14 25 12 45
24 12 20 23 16 21 14 10 44 42 23 21 26 10 44
22 14 13 17 13 24 19 10 39 35 21 17 19 10 39
25 13 16 21 14 24 20 13 41 43 26 18 21 13 43
21 12 19 14 15 26 14 11 40 46 23 24 21 11 46
23 16 23 22 16 23 23 15 41 43 24 24 18 15 43
25 17 28 26 17 25 18 16 48 45 26 20 17 16 48
24 19 26 20 15 23 19 16 45 41 22 17 16 15 45
19 12 13 11 13 19 14 10 37 23 21 14 10
27 21 28 26 17 26 23 16 48 46 26 24 26
where P h is the probability that an interval i of duration h hours is dry, and lag-1 joint probabilities are defined as
g (h 21). It is assumed that one cell arrives at the storm origin. The duration of the rectangular pulse associated with each is exponentially distributed with rate h (h 21). Rectangular pulses of cells are allowed to overlap with cells of the same storm and across cells of different storms. The sequences of the random variables defining the process are assumed mutually independent. For convenience, dimensionless parameters k b=h and f g=h are introduced. The probability that an interval i of duration h hours is dry, P h; is given by (Rodriguez-Iturbe et al., 1987)
P00 h P Yi11 0; Yi 0 P 2h
P h exp2 h 1 uT 2 Cl
P01 h P Yi11 . 0; Yi 0 P h 2 P 2h
where
P11 h P Yi11 . 0; Yi . 0 1 2 2P h 1 P 2h
C
and the lag-1 autocovariance cY h; 1 as " ( cY h; 1 mY h
2
P h P00 h; 1 2 2P01 h; 1 1 2 P h "
P h 1 P11 h; 1 1 2 P h
#
#2 ) 8
9 In Eq. (9), P 2h is the probability that two consecutive time intervals are dry. In this paper, the non-randomised Bartlett–Lewis model (Rodriguez-Iturbe et al., 1987) is used to generate the binary chain. For the non-randomised Bartlett–Lewis model, storms arrive in a Poisson Process of rate l (h 21). Followed by each storm origin is the arrival of cell origins that are also governed by a Poisson Process but of rate b (h 21). Cell arrivals cease after a time, exponentially distributed with rate
10
GP f 1 k e2 k1fhh k1f
11
The expected duration of a storm, mT, is given by Z1 Z1 dn dtn21 tf21 uT f21 h21 1 fh21 0
0
1 2 1 2 nt e2kn 12t and the function GP is defined as Z1 GP h21 e2k tf21 1 2 t ekt dt 0
12
13
Eqs. (12) and (13) are solved numerically. The time
152
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Fig. 2. Stage 1 calibrated parameter values with a band of twice the standard deviation, Rockhampton station.
series generated by the non-randomised Bartlett– Lewis model is represented by a string of two numbers: zero for a dry period and a constant value, w h; for a wet period.
3. Data Six-minute continuous rainfall data from 13 stations within central Queensland, Australia, were used in this paper. The location of the pluviograph stations together with daily rainfall stations having over 50 years of data is shown in Fig. 1. The altitude of these stations varies between 10 and 400 m above
mean sea level. The rainfall data, some dating back to 1939, were digitised from syphon (float) type rain gauge electronic tracing charts by the Bureau of Meteorology, Melbourne, Australia. For a given year, months with missing data were excluded. This resulted in rainfall data of different months for the same station having a different number of years of record. Table 1 shows the number of years of data of the various months of a given rainfall station used. Only stations with at least 10 years of good quality data for each month were used. The data were analysed at 20 aggregation levels: 6, 12, 18, 24, 30, 36, 42, 48 and 54 min and 1, 1.5, 2, 3, 4, 6, 8, 10, 12, 18 and 24 h.
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
153
Table 2 Correlation coefficient between parameter pairs–stage 1 Month
l –b
l –g
l –h
b –g
b –h
g –h
January February March April May June July August September October November December
0.59790 0.48320 0.65910 0.60670 0.58730 0.58070 0.58540 0.56940 0.56940 0.57170 0.57700 0.60150
0.93890 0.94380 0.95970 0.93790 0.94250 0.95330 0.94350 0.93050 0.90140 0.90030 0.88270 0.92170
0.39570 0.39300 0.46270 0.40820 0.39780 0.39400 0.39600 0.38440 0.37980 0.37620 0.38870 0.41170
0.71820 0.65190 0.77340 0.73520 0.66170 0.62710 0.65600 0.61680 0.71630 0.73600 0.78400 0.74560
0.92940 0.93780 0.93630 0.93290 0.93950 0.94140 0.93950 0.94120 0.93670 0.93170 0.93730 0.93900
0.46980 0.40640 0.55340 0.49570 0.41970 0.38920 0.41380 0.36740 0.46890 0.48530 0.55290 0.51250
4. Identification of regional parameters of the binary model The binary Bartlett–Lewis model has four parameters (l , b , g and h ) to be calibrated for each month, resulting in 48 parameters for the whole year for a given rainfall station. The parameters of the model are examined for uniformity over the months and the possibility of fitting a curve in order to reduce the number of parameters. This may provide a smooth representation of the parameters with a view to capturing the seasonal variability of the data. Moreover, regionalisation is easier with a reduced number of parameters. In modelling rainfall over Britain using the random Bartlett–Lewis model, Onof and Wheater (1993) fitted polynomial curves of degree varying between 5 and 12 to the individual parameters. Their objective was to represent the parameters as daily variables in order to avoid the artificiality of the monthly subdivision of the year and allowing for possible generalisation of parameters through regional analysis. The parameters of the binary model were identified through four stages using rainfall data of Rockhampton. At each stage the parameters were calibrated by minimising the objective function, J, J
X
Po h 2 Ps h2
14
where Po(h) and Ps(h) are the historical and analytical dry probabilities of aggregation level h, respectively. The shuffled complex evolution global probabilistic
search strategy (Duan et al., 1992) option of the Bayesian non-linear regression software NLFIT (Kuczera, 1994) was used to estimate the parameters. The number of complexes was set equal to the number of fitted parameters. The global optimum parameter set was then used as starting values for a local gradient search. NLFIT uses the computationally efficient classical approximations to parameter uncertainty based on the multinormal approximation of the posterior distribution to estimate the standard deviations and correlations of the parameter sets. The covariance matrix of the multinormal approximation can help identify the poorly defined parameters and their interactions (Kuczera and Mroczkowski, 1998). In stage one, the four parameters were calibrated for each month. Fig. 2 depicts the values of the calibrated parameters, with a band of twice the standard deviation. It is observed that only parameter l is well identified with a clear seasonal pattern, high values for the wet season (October to March) and low values for the dry season. Parameter g also shows a seasonal trend. However, the envelopes of twice the standard deviation of the parameters b and h suggest the possibility of adopting constant values for these parameters over the months without undue loss of accuracy. The very high correlations between the parameter pairs b and h (Table 2) suggest nearly indistinguishable simulation results from parameters spanning a wide range of parameter space. Thus only one is assumed to have a constant value for all months in stage two of the parameter identification process. During stage two of the parameter identification
154
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163 Table 3 Calibrated parameter values of the binary model
Fig. 3. Stage 2 calibrated parameter h values with a band of twice the standard deviation, Rockhampton station.
process, parameter b was considered constant for all months, reducing the number of parameters calibrated to 37. Since the correlation between the parameter pairs l – b and g – b in stage one was not very high, there were slight changes in the values of the parameters l and g . The stage two calibrated values of h indicated narrow fluctuations with envelopes of two standard deviations (Fig. 3), yet a constant value could be assumed. Therefore in the third stage of the parameter identification process, parameters b and h were assumed to have constant values for all months,
Parameter (h 21)
Rockhampton
Regional
b h a0 a1 b1 a2 b2
0.7095 1.3721 0.1407 0.0432 2 0.0428 2 0.0094 2 0.0369
0.4808 1.3175 0.1180 0.0316 2 0.0335 2 0.0026 2 0.0226
reducing the number of calibrated parameters to 26. While there were changes in the values of the parameters l and g , the trend of their variation with the month of the year remained largely unchanged. What emerged at this stage was a strong logarithmic relationship between parameter l and the dry probabilities. This is as a result of Eq. (10) for given values of b , g and h . Moreover, the correlation between the parameter pair l – g was in excess of 83% for all months, an indication of multiple optima within the feasible parameter space. In stage 4, parameters b and h were still considered constant for all months and second harmonic Fourier series
g m a0 1
2 X ak cos k2pm=12 1 bk sin k2pm=12 k1
15
Fig. 4. Fitted second harmonic Fourier series of parameter g .
was used to represent the variation of parameter g over the months, m. In Eq. (15) a0, ak and bk are coefficients to be determined. For given parameter values of b , g and h , and the daily dry probability, P(24), Eq. (10) was used to estimate parameter l . Fig. 4 shows second harmonic Fourier series fitted to parameter g and Table 3 gives the calibrated parameter values at stage 4. Essentially, the number of calibrated parameters in stage 4 was 7. Eq. (10) and the second harmonic Fourier series for parameter g capture the seasonality effects in the sequence of dry and wet spells. The binary model error in reproducing the dry probabilities through the four stages was examined using two error statistics: the average absolute percentage error E1 and the maximum absolute percentage error
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
155
Fig. 5. Comparison of error statistics E1 and E2 during stages 1–4, Rockhampton station.
E2, defined as E1 %
n 1 X Po k 2 Ps k 100 n k1 Po k
P k 2 Ps k 100 E2 % max o k1;n Po k
16
17
where Po(k) and Ps(k) are as defined in Eq. (14) and n is the number of values. If the error statistics are estimated for a given month then n is 20 the number of aggregation levels, but if they are for a given aggregation level then n is 12 the number months. The error statistics are presented in Fig. 5. Apparently, there were no significant differences between the error statistics from stages 1–3. Stage 4 error statistics
were slightly higher than the values of the other stages. However, for all intents and purposes the stage 4 fit was very good since the error statistics were very small. The representativeness of Rockhampton’s calibrated parameter values (Table 3) for central Queensland region was investigated. These calibrated values and the P(24) values of the 12 other stations were used to generate the analytical dry probabilities of the 20 aggregation levels. Fig. 6 shows the error statistics of the 12 stations for the individual months and also for the aggregation levels, error statistic E1 being averaged over the 12 stations and E2 the maximum of the 12 stations. Clearly, the error statistics were small (all less than 1.4% for E1 and less than 7% for E2), the dry months having the least values. In Fig. 7 is shown the
156
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Fig. 6. Variation of error statistics E1 and E2 over the months and aggregation levels.
error statistics for the individual stations, estimated over the aggregation levels and months. Next, all 13 stations were jointly calibrated for a representative regional parameter set (see Table 3 and Fig. 4). The Rockhampton and the regional parameter sets were different, but the average error statistics were not significantly different (Figs. 6 and 7). It is therefore concluded that a set of parameter values could be used to disaggregate daily dry probabilities to finer aggregation levels for all stations within the region. The values of the first three stations were highest, suggesting the possibility of clustering since these stations are grouped geographically. To investigate this issue the stations were divided into two groups separated by latitude 238 (Fig. 1). Joint calibration was carried out for each group. Results not shown here
indicated a slight reduction in the error statistics in the average sense, but too small to warrant grouping the stations into clusters.
5. Estimation of sub-daily historical moments In order to apply the jitter model, the historical moments at the simulation timescale are required. The contained analysis is a practical solution for estimating rainfall statistics for different timescales. This analysis is not a consequence of the hybrid model discussed in the previous sections and an alternative approach could be adopted. The rainfall depths of timescale h, Y k h; are
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
157
Fig. 7. Variation of error statistics E1 and E2 over the stations.
connected with the daily value Y 24 by Y 24 Y 1 h 1 Y 2 h 1 … 1 Y n h
18
Fig. 8. Variation of 1-h timescale lag-1 autocorrelation over the months (top) and of 1-h timescale lag-2 autocorrelation with 1-h timescale lag-1 autocorrelation (bottom).
where n 24=h: From the additive properties of random variables with finite variances
mY 24
n X
mYk h nmY h
19
k1
and
s2Y 24
n X k1
s2Yk h 1 2
X
h i cov Y j h; Y k h
20
j,k
Fig. 8 (top) shows the 1-h lag-1 autocorrelation, rY 1; 1; values of all station-months. In order to reduce the number of parameters, the monthly mean
Table 4 Parameter values of the second harmonic Fourier series fitted to 1-h lag-1 autocorrelation, rY 1; 1 Parameter
Value
f0 f1 g1 f2 g2
0.4189 2 0.1038 0.0312 0.0375 0.0476
158
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
be simplified to
mY 24 24mY 1
22
and Eq. (20) to
s2Y 24 24 1 2 24 2 1rY 1; 1 1 erY 1; 12 s2Y 1 23
Fig. 9. Observed and predicted 1-h variance and lag-1 autocovariance. Solid line is a line of perfect match.
values (open circle) were fitted with second harmonic Fourier series
rY 1; 1 m f0 1
2 X
fk cos k2pm=12
where the third term on the right is an error term accounting for the higher lag autocovariances. Parameter e varied over the station-months without apparent regional trend. However, a constant regional value of parameter e was sought in addition to the rY 1; 1 regional function to fit s2Y 1: The s2Y 24 values of 14 station-months, mostly from the summer months of stations 33087, 33099 and 33119, exceeded 200 mm 2 and required a separate parameter e value. The parameter e was estimated as e 201:5; s2Y 24 # 200 and e 670:5; s2Y 24 . 200 mm2 : Fig. 9 compares the observed and predicted values of 1-h variance and lag-1 autocovariance. It is observed in Fig. 9 that the predicted second-order moments matched the observed values very well. The 14 station-months have exceptionally high variances and are the wettest of all stationmonths, explaining the huge range between the two constant values of parameter e. A higher value of parameter e indicates the importance of the higher lag autocorrelations. Sub-hourly timescales can be similarly handled. For the 6-min-timescale, it was observed that rY 1; 1 was nearly constant with an average value of 0.8. However, the variation of parameter e over the station-months was very high. In such circumstances it is recommended to work with hourly values and distribute the hourly rainfall depths over the wet sub-hourly intervals as explained in the next section.
k1
1gk sin k2pm=12
21
where f0 , fk and gk are coefficients given in Table 4. The higher lag autocorrelations of the 1-h timescale were observed to increase with the square of the lag-1 autocorrelation, as shown in Fig. 8 (bottom) for lag-2. Therefore, for 1-h simulation timescale Eq. (19) can
6. Evaluation of the regionalised hybrid model The procedure for estimating the parameters of the hybrid model when only daily data are available is summarised as follows: 1. calculate mY 24, s2Y 24 and P(24) for each month;
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
159
Fig. 10. Comparison of historical (shaded circle) and simulated (dashed line) dry probability.
2. use Eq. (21) to estimate rY 1; 1; Eq. (22) to estimate mY 1; Eq. (23) to estimate s2Y 1 and then calculate cY 1; 1 rY 1; 1s2Y 1; 3. estimate parameter l in Eq. (10) using P(24) and the regional parameters given in Table 3; 4. select the simulation timescale, h, and use Eq. (10) to estimate P h and P 2h; 5. estimate the moments of the binary chain model using Eqs. (6)–(9);
6. estimate the parameters of the jitter model using Eqs. (2), (3) and (5). The binary process can be simulated at any desired timescale but 6-min (the historical data timescale) is used in this paper. However, the jitter model was simulated at 1-h timescale. Analysis of the distribution of the historical hourly rainfall depth over the 6-min wet intervals indicated that the fractional
160
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Fig. 11. Comparison of historical (shaded circle) and simulated (dashed line) mean.
weights vary randomly and could be sampled from a uniform distribution, U 0:1; 1: Hence the fractional weights within a wet hour were drawn from U 0:1; 1 and each divided by their sum in order to preserve a sum of unity. The scaled fractional weights were multiplied by the hourly rainfall depth to obtain the 6-min interval’s rainfall depth. For example, if the rainfall depth within an hour was 50 mm, and there were five 6-min wet periods within that hour, then five
fractional weights 0.405, 0.775, 0.271, 0.152 and 0.410 drawn from U 0:1; 1 were divided by 2.013, their sum. Each of the scaled fractional weights 0.201, 0.385, 0.135, 0.075 and 0.204 was then multiplied by 50 mm to obtain the rainfall depths 10.05, 19.25, 6.75, 3.75 and 10.20 mm of the five 6-min wet periods. This approach preserved the dry probabilities. Figs. 10–13 compare the historical statistics
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
161
Fig. 12. Comparison of historical (shaded circle) and simulated (dashed line) variance.
(shaded circle) with 1000 years of simulation results (dashed line) of 6 station-months, used as examples. These station-months are from December and January, the wettest months, and across the region. The dry probabilities were very well reproduced. It must be underlined that the simulated and analytical dry probabilities were very close. In all cases the historical means were very well matched by the simulated values. The variances and lag-1 autocovariances
reproduced by the hybrid model were in good agreement with the historical ones. Although for some station-months the absolute percentage error of variance and lag-1 autocovariance were in excess of 50% for a given aggregation level, lag-1 autocovariance in particular. Station-month Dec-33099, showing the worst scenario, is one of the station-months with s2Y 24 in excess of 200 mm 2. As indicated by the examples shown, the method of distributing hourly
162
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Fig. 13. Comparison of historical (shaded circle) and simulated (dashed line) lag-1 autocovariance.
rainfall depths over 6-min wet periods during the simulation worked very well. Data size could play a significant role in the errors in the reproduction of the variances and lag-1 autocovariances. For example, Cowpertwait et al. (1996b) established that for a 10-year record of British rainfall data the coefficient of variation of mY 24; s2Y 24 and P(24) was of the order of 16, 31 and 12%, respectively. These values were reduced by about 30%
when the record length was doubled. Hence, with the availability of more data it is hoped that the reproduction of the variance and lag-1 autocovariance by the hybrid model would be improved.
7. Conclusions In this paper it is demonstrated how the
Y. Gyasi-Agyei / Journal of Hydrology 223 (1999) 148–163
Gyasi-Agyei–Willgoose hybrid model for point processes could be regionalised for daily rainfall disaggregation using limited high resolution rainfall data within a region of interest. Their model is a product of the binary non-randomised Bartlett– Lewis rectangular pulse model and a lognormal autoregressive model used as a jitter. The time series generated by the non-randomised Bartlett–Lewis model is replaced by a string of two numbers: zero for a dry period and a constant value for a wet period chosen to preserve the mean of the historical rainfall data. The four stages of the parameter identification process of the binary model involve joint calibration of the parameters of all months. The computationally efficient multinormal approximation to parameter uncertainty is used to group the monthly parameter values of the binary model. For central Queensland, Australia, it has been established that the parameters b and h could have constant values for all months. Second harmonic Fourier series was used to represent the seasonal variation of the parameter g . The parameter l is a function of the daily dry probability and the parameters b , h and g . Thus only a minimal number of parameters need to be calibrated making the model more parsimonious. The approach removes the artificiality of grouping the yearly data on a monthly basis. Also the additive properties of random variables with finite variances were used to scale down the daily mean and variance of the historical data to the simulation timescale, values required by the jitter model. The regionalised model is structured so that it requires only the daily values of dry probability, mean and variance, the regional 1-h lag-1 autocorrelation function and parameter e to simulate rainfall at a fine timescale. The results on using observed daily rainfall statistics to capture subdaily statistics by the regionalised model are very encouraging. The model is therefore a valuable tool for disaggregating daily rainfall, data of which abound worldwide. Further research is, however, required to improve the reproduction of the second-order moments.
163
Acknowledgements The reviews by E.N. Anagnostou, T. Lebel and D. Koutsoyiannis were very helpful and are gratefully acknowledged. The author wishes to thank G. Kuczera for providing NLFIT software. Financial support from the Coal and Minerals Infrastructure Division of Queensland Rail, Australia, is very much appreciated. References Bo, Z., Islam, S., Eltahir, E.A.B., 1994. Aggregation–disaggregation properties of a stochastic rainfall model. Water Resour. Res. 30 (12), 3423–3435. Cowpertwait, P.S.P., 1991. Further developments of the Newman– Scott clustered point process for modelling rainfall. Water Resour. Res. 27 (7), 1431–1438. Cowpertwait, P.S.P., O’Connell, P.E., Metcalfe, A.V., Mawdsley, J.A., 1996a. Stochastic point process modelling of rainfall. I. Single-site fitting and validation. J. Hydrol. 175, 17–46. Cowpertwait, P.S.P., O’Connell, P.E., Metcalfe, A.V., Mawdsley, J.A., 1996b. Stochastic point process modelling of rainfall. II. Regionalisation and disaggregation. J. Hydrol. 175, 47–65. Duan, Q., Sorooshian, S., Gupta, V., 1992. Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resour. Res. 28 (4), 1015–1031. Glasbey, C.A., Cooper, G., McGehan, M.B., 1995. Disaggregation of daily rainfall by conditional simulation from a point-success model. J. Hydrol. 165, 1–9. Gyasi-Agyei, Y., Willgoose, G.R., 1997. A hybrid model for point rainfall modelling. Water Resour. Res. 33 (7), 1699–1706. Gyasi-Agyei, Y., Willgoose, G.R., 1999. Generalisation of a hybrid model for point rainfall. J. Hydrol. 219, 218–224. Kuczera, G., 1994. NLFIT: a Bayesian nonlinear regression program suite, Department of Civil Engineering and Survey, the University of Newcastle, Callaghan, NSW, Australia, 162pp. Kuczera, G., Mroczkowski, M., 1998. Assessment of hydrologic parameter uncertainty and the worth of multiresponse data. Water Resour. Res. 34 (6), 1481–1489. Onof, C., Wheater, H.S., 1993. Modelling of British rainfall using a random parameter Barlett–Lewis rectangular pulse model. J. Hydrol. 149, 67–95. Rodriguez-Iturbe, I., Cox, D.R., Isham, V., 1987. Some models for rainfall based on stochastic point processes. Proc. R. Soc. London A 410, 269–288. Rodriguez-Iturbe, I., Cox, D.R., Isham, V., 1988. A point success model for rainfall: further developments. Proc. R. Soc. London A 417, 283–298.