Generalisation of a hybrid model for point rainfall

Journal of Hydrology 219 (1999) 218–224

Technical note

Generalisation of a hybrid model for point rainfall Y. Gyasi-Agyei 1,*, G.R. Willgoose Department of Civil, Surveying and Environmental Engineering, The University of Newcastle, Callaghan, NSW 2308, Australia Received 4 September 1998; accepted 8 April 1999

Abstract A generalised hybrid model to generate point rainfall for a wide range of aggregation levels is presented in this paper. The rainfall process is expressed as a product of a binary chain model, which preserves the dry and wet sequences as well as the mean, and a correlated jitter used to improve the deficiencies in the second-order properties of the binary chain. Analytical derivations of the moments of a binary chain are presented. As the jitter model the exponential of a second-order autoregressive Gaussian process is selected. Two possible binary chain models are analysed, a non-randomised Bartlett–Lewis model and a Markov chain. Although both binary chain models perform equally well, the Bartlett–Lewis model is preferred for reasons of parameter parsimony. q 1999 Elsevier Science B.V. All rights reserved. Keywords: Rainfall; Aggregation; Random processes; Binary chain; Jitter model; Hybrid model

1. Introduction A rainfall time series is generally represented as a combination of an occurrence process (i.e. the sequence of wet and dry periods) and an intensity process, which assigns values to wet periods. Models taking into account these two processes may either model them separately or together. Chain-dependent process models have been used for daily (e.g. Chin, 1977) and hourly (Katz and Parlange, 1995) rainfall but without much attention to finer time scales and aggregation properties. Katz and Parlange (1995), in fitting variants of first-order Markov chain models to hourly rainfall data, observed that the models failed to reproduce the statistics of 12 * Corresponding author. Fax: 1 61-7493-09-382. E-mail address: [email protected] (Y. Gyasi-Agyei) 1 Now at James Goldston Faculty of Engineering and Physical Systems, Central Queensland University, Rockhampton, Queensland 4702, Australia.

and 24 h aggregation levels. They recommended fitting higher-order Markov chain models to the hourly occurrence process. However, lack of parsimony has been noted as a drawback of the Markov chain models as the order is increased. In general, Markov chain models are attractive because of their largely non-parametric nature (i.e. the parameters are derived directly from data), their ease of application and interpretability, and well-developed literature (Lall et al., 1996). Lall et al. (1996) presented a non-parametric wet/dry spell model for resampling daily rainfall. They acknowledged the need to evaluate the concurrent representation of the rainfall process at different time scales. Cluster-based point rainfall models (Bartlett– Lewis and Neyman–Scott) pioneered by RodriguezIturbe et al. (1987, 1988) have had a wide appeal in the hydrological modelling community. A common approach to the estimation of the parameters of these models is by using the method of moments. Some selected statistics of the historical data,

0022-1694/99/$ - see front matter q 1999 Elsevier Science B.V. All rights reserved. PII: S0022-1694(99 )00 054-2

Y. Gyasi-Agyei, G.R. Willgoose / Journal of Hydrology 219 (1999) 218–224

covering a range of aggregation levels, are matched to those of the models. The number of statistics used in the calibration is usually set to the number of parameters of the model. One problem with the method of moments is that different combinations of the statistics may give different results of the statistics not used in calibration particularly the probability that an interval is dry (e.g. Velghe et al., 1994; Gyasi-Agyei and Willgoose, 1997). The hybrid model presented by Gyasi-Agyei and Willgoose (1997) avoids the need to select a set of moments for the parameter calibration. Their model is a product of two random process models, being the non-randomised Bartlett–Lewis and an autoregressive model, referred to in the remainder of this paper as BLAR. For the BLAR model, the five parameters of the Bartlett–Lewis model were estimated by the method of moments using the mean of one aggregation level and the dry probabilities of all aggregation levels considered. The parameters of the autoregressive model were then derived from the second-order properties of the historical data and those given by the Bartlett–Lewis model, without additional parameter calibration being necessary. For the 15 min point rainfall data from Capella, Central Queensland, Australia, the BLAR model reproduced the second-order moments of the historical data very well at all aggregation levels. The purpose of this paper is to generalise the hybrid model by replacing the Bartlett–Lewis model with a binary chain model while maintaining the autoregressive model used as a jitter. The binary chain consists of a string of two numbers, zero for a dry period and a constant value, w, for a wet period. A Markov chain and the Bartlett– Lewis models are used as examples. Where the Bartlett–Lewis model is used to generate the binary chain, the distribution of the rectangular pulse (or cell) intensity within the rainfall event is not required, reducing the number of parameters to be calibrated by one. Moreover, regionalisation is easier with a reduced number of parameters. The formulation of the hybrid model and the autoregressive model parameter estimation given in Gyasi-Agyei and Willgoose (1997) are summarised in Section 2. This is followed by the analytical derivation of the moments of a binary chain in Section 3. In Section 4, the models for generating a binary chain are

219

presented. Section 5 compares and evaluates the hybrid models using a rainfall data set. The same data set, 15 min point rainfall data from Capella, Central Queensland, Australia (Carroll et al., 1997), given in Gyasi-Agyei and Willgoose (1997) and nine aggregation levels (0.25, 0.5, 1, 2, 4, 6, 8, 12 and 24 h) are used. 2. The hybrid model and jitter process formulation Following Rodriguez-Iturbe et al. (1987), Onof and Wheater (1994) and Gyasi-Agyei and Willgoose (1997), the hybrid model, {Y…t†}; is a product of two random processes, a jitter model {A(t)} and a binary chain model {Y(t)}, expressed as {Y…t†} ˆ {A…t†}{Y…t†}

…1†

The jitter model {A(t)} is the exponential of an autoregressive Gaussian process {Z(t)} which has the mean, m Z, variance s2Z and lag-t autocovariance function cZ(t ), and t is time. 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 and that Y…t† ˆ 0 only if Y(t) ˆ 0. This ensures that the probability of dry periods, the number and the duration of rainfall events remain the same for both {Y(t)} and {Y…t†}: Choosing m A ˆ 1, that is making mY ˆ mY ;and with the properties of the log-normal distribution, it can be shown that 2mZ ˆ 2s2Z : With this result and the moment generating function of a bivariate Gaussian distribution, it can be shown that E‰A…t†·A…t 1 t†Š ˆ exp‰cz …t†Š. These assumptions lead to the second-order properties of {Z(t)} as (Gyasi-Agyei and Willgoose, 1997) ! s2Y 1 m2Y 2 sz ˆ ln 2 …2† sY 1 m2Y and c …t† 1 m2Y cz …t† ˆ ln Y cY …t† 1 m2Y

! …3†

Eqs. (2) and (3) show that the second-order properties of the stationary Gaussian process can be derived from those of the historical data and the random process {Y(t)}. The random process {Z(t)} is modelled by a

220

Y. Gyasi-Agyei, G.R. Willgoose / Journal of Hydrology 219 (1999) 218–224

second-order autoregressive model, given as (e.g. Bras and Rodriguez-Iturbe, 1985, p. 19) Zt ˆ mz 1 v1 …Zt21 2 mz † 1 v2 …Zt22 2 mz † 1 Et

…4†

where {Et} is normally distributed, with a mean value of zero and variance s2E given by

s2E ˆ ‰1 2 v1 rz …1† 2 v2 rz …2†Šs2z

…5†

In Eq. (5), r Z(i) is the lag-i autocorrelation and v i a constant. The parameters v i are the solutions of the Yule–Walker equations expressed as (e.g. Bras and Rodriguez-Iturbe, 1985, p. 20)

v1 ˆ rz …1†

‰1 2 rz …2†Š ‰1 2 rz …1†2 Š

v2 ˆ

‰rz …2† 2 rz …1†2 Š ‰1 2 rz …1†2 Š …6†

In all cases the conditions urz …i†u , 1; 2‰rz …1†Š , rz …2† 1 1; v1 1 v2 , 1; v2 2 v1 , 1 and 21 , v2 , 1 must be satisfied. Note that parameters of the jitter process are derived without additional calibration. They are given by the second-order properties of the finest aggregation level of the historical data and those of the random process {Y(t)}.

s2Y of a binary chain is derived as s2Y ˆ E‰…Yi 2 mYi †2 Š ˆ P…i†…0 2 mYi †2 1 ‰1 2 P…i†Š…w 2 mYi †2 ( ˆ …mYi †

2

ˆ …mYi †2

mYi wˆ 1 2 P…i†

…7†

where P(i) is the probability that an interval i is dry, also written as P(Yi ˆ 0). Assuming second order stationarity, the variance

P…i† ‰1 2 P…i†Š

T00 ˆ P…Yi11 ˆ 0uYi ˆ 0† ˆ ˆ

…8†

P…Yi11 ˆ 0; Yi ˆ 0† P…Yi ˆ 0†

P…2i† P…i†

…9†

T01 ˆ P…Yi11 . 0uYi ˆ 0† ˆ 1 2 P…Yi11 ˆ 0uYi ˆ 0† ˆ

P…i† 2 P…2i† P…i†

T10 ˆ P…Yi11 ˆ 0uYi . 0† ˆ ˆ

A binary chain model generates a string of two numbers, Yi ˆ 0 for a dry period and a constant value Yi ˆ w for a wet period, where Yi is the cumulative rainfall depth over time interval i. The moments of a binary chain are derived analytically as functions of the dry probabilities. For the historical data and a binary chain of the same time scale to have the same mean, mYi , the cumulative rainfall depth over a wet period w must be given by

!2 )

Before the autocovariance of a binary chain is derived, the following conditional or transition probabilities are first defined:

2

3. Derivation of the moments of a binary chain

w P…i† 1 ‰1 2 P…i†Š 21 mYi

P…Yi11 ˆ 0; Yi . 0† P…Yi . 0†

P…i† 2 P…2i† 1 2 P…i†

T11 ˆ P…Yi11 . 0uYi . 0† ˆ 1 2 P…Yi11 ˆ 0uYi . 0† ˆ

1 2 2P…i† 1 P…2i† 1 2 P…i†

In Eq. (9), P(2i) is the probability that two consecutive time intervals of total length 2i are dry, i.e. P(Yi ˆ 0, Yi11 ˆ 0). The variable T01 means a transition from a dry (0) interval to a wet (1) one. The lag-l joint probabilities are also defined as follows: P00 …l† ˆ P…Yi11 ˆ 0; Yi ˆ 0† P01 …l† ˆ P…Yi11 . 0; Yi ˆ 0† P10 …l† ˆ P…Yi11 ˆ 0; Yi . 0† P11 …l† ˆ P…Yi11 . 0; Yi . 0†

…10†

Y. Gyasi-Agyei, G.R. Willgoose / Journal of Hydrology 219 (1999) 218–224

Noting that the joint probabilities are the proportions of the occurrences of two rainfall values separated by the lag in the chain, the lag-l autocovariance Cov…Yi1l ; Yi † is given by Cov…Yi11 ; Yi † ˆ E‰…Yi 2 mYi †…Yi1l 2 mYi †Š

221

two branches at a node, higher lag joint probabilities can be derived. It is easily shown that for all lags, P01 ˆ P10 . As the moments of a binary chain are functions of mYi , P(i) and P(2i), two binary chains having the same values of these statistics will give the same parameter set for the jitter process.

ˆ P00 …l†…0 2 mYi †…0 2 mYi † 1P01 …l†…0 2 mYi †…w 2 mYi † 1 P10 …l†…w 2 mYi †…0 2 mYi † 1P11 …l†…w 2 mYi †…w 2 mYi † ( ˆ …mYi †2

"

P…i† P00 …l† 2 ‰P01 …l† 1 P10 …l†Š 1 2 P…i† "

P…i† 1 P11 …l† 1 2 P…i†

#

4. Models for the random process {Y(t)} Any model which preserves the dry probabilities of the aggregation levels of interest can be used for the random process {Y(t)}. The non-randomised Bartlett– Lewis and Markov Chain models used as examples in this paper are briefly described.

#2 ) …11†

4.1. The non-randomised Bartlett–Lewis model

From Eq. (10), lag-1 joint probabilities are derived as P00 …1† ˆ P…Yi11 ˆ 0uYi ˆ 0†P…Yi ˆ 0† ˆ T00 P…i† P01 …1† ˆ P…Yi11 . 0uYi ˆ 0†P…Yi ˆ 0† ˆ T01 P…i† P10 …1† ˆ P…Yi11 ˆ 0uYi . 0†P…Yi . 0† ˆ T10 ‰1 2 P…i†Š P11 …1† ˆ P…Yi11 . 0uYi . 0†P…Yi . 0† ˆ T11 ‰1 2 P…i†Š

…12†

Similarly, lag-2 joint probabilities are P00 …2† ˆ P…i†‰P…Yi12 ˆ 0uYi11 ˆ 0† × P…Yi11 ˆ 0uYi ˆ 0†Š 1P…i†‰P…Yi12 ˆ 0uYi11 . 0†P…Yi11 . 0uYi ˆ 0†Š ˆ P…i†‰T00 T00 1 T01 T10 Š P01 …2† ˆ P…i†‰T00 T01 1 T01 T11 Š

The non-randomised Bartlett–Lewis rectangular pulse model (Rodriguez-Iturbe et al., 1987) is the first model considered. This model is the same as reported in Gyasi-Agyei and Willgoose (1997) except that the distribution of the rectangular pulse (or cell) intensity is not required, reducing the number of parameters by one. As a result, the number of parameters to be calibrated is reduced by twelve, if the data are grouped on a monthly basis. In this case the objective function (20) of that paper reduces to a least square fit of the dry probabilities. Using the same optimisation procedure (Duan et al., 1992; Kuczera, 1997) the calibrated values of the four remaining parameters are unchanged. Rather, the time 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, for a wet period. The moments of the binary chain are estimated analytically from the equations presented in Section 3. In what follows, the hybrid of this binary chain and the second-order autoregressive models is referred to as BBLAR.

P10 …2† ˆ ‰1 2 P…i†Š‰T10 T00 1 T11 T10 Š 4.2. Markov chain model

P11 …2† ˆ ‰1 2 P…i†Š‰T10 T01 1 T11 T11 Š: …13† By means of a probability decision tree diagram with

The random process {Y(t)} is said to be an n-state Markov chain of order q if (Foufoula-Georgiou and

222

Y. Gyasi-Agyei, G.R. Willgoose / Journal of Hydrology 219 (1999) 218–224

Table 1 Parameters of the hybrid models (l ˆ 0.0083748 h 21; k ˆ 0.48638; f ˆ 0.09943; h ˆ 1.49326 h 21) Model

s2Z

r Z(1)

r Z(2)

s2E

v1

v2

w

BLAR BBLAR

1.00981 1.83457

0.57717 0.77505

2 0.03424 0.40012

0.46906 0.54768

0.89511 1.16440

2 0.55087 2 0.50230

– 0.59225

Georgiou, 1987; Lloyd, 1974). P…Yi ˆ xi uYi21 ˆ xi21 ; Yi22 ˆ xi22 ; …; Yi2q ˆ xi2q ; Yi2q21 ˆ xi2q21 ; …; Y0 ˆ x0 † ˆ P…Yi ˆ xi uYi21 ˆ xi21 ; Yi22 ˆ xi22 ; …; Yi2q ˆ xi2q †

…14†

where xi takes on values from a set of n states, and P…Yi ˆ xi uC† denotes the conditional probability that Yi takes the value xi, given the condition C. That is P…Yi ˆ xi uC† is the relative frequency of the event Yi ˆ xi in the particular subset of combinations of previous values specified by C. Eq. (14) implies that the probability distribution of the next value Yi is determined by the q most recent values. It can be shown (Lloyd, 1974) that the random sequence of C formed out of q tuples of the sequence {xi} is an N-state Markov chain of order 1, where N ˆ n q. The N × N transition probability is a sparse matrix, with at most n non-zero

Fig. 1. Variation of dry probability with aggregation level.

entries for each row or column. An N-state Markov chain of order 1 is used in the hybrid model. In this paper, there are only two states, dry and wet periods, so order q Markov chain is replaced by 2 qstate Markov chain of order 1. The transition probabilities (i.e. the parameters of the model) are derived directly from the historical data. Each column or row has only two entries. It is proposed that the lowest order to give satisfactory fit of the aggregated dry probabilities be chosen. A hybrid model of a Markov chain of order q and the second-order autoregressive model is referred to as MCARq.

5. Comparison and evaluation of the hybrid models In this section the hybrid models are compared and evaluated using 1000 years of simulated data. Table 1 gives the estimated parameters of the jitter process of the hybrid model BBLAR, as well as those of BLAR (Gyasi-Agyei and Willgoose, 1997), for comparison. The calibrated parameters of the non-randomised Bartlett–Lewis model are given in the table. Also given in this table is the constant depth for wet periods, w. Note that autocorrelation values of the jitter component of the BBLAR model are all positive and are larger than the corresponding values of the BLAR model. The same is true of the variance values. These observations are due to the smoothing effects of rainfall depths of wet periods. Fig. 1 depicts the variation of dry probability with aggregation levels of the hybrid models. It is not surprising that the BBLAR model simulated values match the historical ones, because the dry probabilities were used in calibrating its parameters. For the MCAR model, the dry probabilities are not satisfactorily reproduced until the order reaches 12. Note that the hybrid models have identical jitter model parameters as they reproduced very well the dry probabilities at 15 and 30 min aggregation levels, statistics controlling Eqs. (7)–(13). The BBLAR and

Y. Gyasi-Agyei, G.R. Willgoose / Journal of Hydrology 219 (1999) 218–224

223

Fig. 2. Aggregation statistics of the hybrid models-mean and variance.

MCAR12 models are examined in the remainder of the paper. The MCAR12 model has 4096 states of Markov chain of order 1. From the historical data only 451 states have marginal probability different from zero and greater that 1.E-5. Out of these states, 112 have marginal probability greater than 1.E-4, 10 greater than 1.E-3 but less than 1.E-2. The state with only zero (dry) entries has the highest marginal probability of 0.927. The states with zero marginal probability indicated by the historical data may occur during the simulation. In such cases the dry probability of the

Fig. 3. Aggregation statistics of the hybrid models-autocorrelations.

224

Y. Gyasi-Agyei, G.R. Willgoose / Journal of Hydrology 219 (1999) 218–224

simulation time step may be used. This situation occurred only in 0.002% of the 1000 years of rainfall data generated by the MCAR12 model. Hence the number of states could be trimmed down to those with marginal probability greater than zero plus one to cater for the remaining states should they occur during the simulation. The order appears to be high but it is only saying that the probability distribution of the next interval’s rainfall depth is dependent on the values within the previous 3 h, on a 15 min time scale. This implies that order 3 Markov chain model may suffice on a 1 h time scale. In Fig. 2 it is observed that the aggregated mean and variance are very well reproduced by the hybrid models. Fig. 3 demonstrates that autocorrelations are also well reproduced by the models, the BBLAR model results being slightly better than those of the MCAR12 model. It can be seen that the two models have a higher potential to mimic the cycles observed in the aggregated autocorrelations compared with the BLAR model (Gyasi-Agyei and Willgoose, 1997).

6. Concluding remarks A generalised hybrid model for point rainfall, consisting of a binary chain and an autoregressive model, has been introduced in this paper. The binary chain consists of a string of two states, zero for a dry period and a constant value for a wet period in the case of rainfall. Analytical derivations of the moments of the binary chain have been presented. It is believed that any model, which preserves the sequence of the two states, could be used in the hybrid model. The autoregressive model is used as a jitter to fix deficiencies in the second-order properties of the binary chain. 15 min point rainfall data from Capella, Central Queensland, Australia, were used to compare and evaluate two hybrid models. While these two models reproduced the aggregated statistics very well, the BBLAR model is superior to the Markov chain model because it is parsimonious in terms of the number of model parameters.

Acknowledgements This work was carried out as part of the Post-

Mining Landscapes Project sponsored by the Australian Coal Association Research Program (Grant C1629), BHP Australia Coal Ltd., Callide Coalfields Pvt. Ltd., Capricorn Coal Management Pvt. Ltd., Curragh Queensland Mining Ltd., MIM Holdings Ltd, and Pacific Coal Pvt. Ltd.. This work was also partly sponsored by a TUNRA postdoctoral fellowship to the first author. We thank Chris Carroll of Queensland Department of Natural Resources (Emerald) for providing the Capella rainfall data. Comments by George Kuczera were helpful.

References Bras, R.L., Rodriguez-Iturbe, I., 1985. Random Functions and Hydrology. Addison-Wesley, Reading, MA 559 pp. Carroll, C., Halpin, M., Burger, P., Bell, K., Sallaway, M.M., Yule, D.F., 1997. The effect of crop type, crop rotation and tillage practice, on runoff and soil loss on a vertisol in Central Queensland. Aust. J. Soil Res. 35, 925–939. Chin, E.H., 1977. Modelling daily precipitation occurrence process with Markov Chain. Water Resour. Res. 13, 949–956. Duan, Q., Sorooshian, S., Gupta, V., 1992. Effective and efficient global optimization for conceptual rainfall-runoff models. Water Resour. Res. 28 (4), 1015–1031. Foufoula-Georgiou, E., Georgiou, T.T., 1987. Interpolation of binary series based on discrete-time Markov chain models. Water Resour. Res. 23 (3), 515–518. Gyasi-Agyei, Y., Willgoose, G.R., 1997. A hybrid model for point rainfall modelling. Water Resour. Res. 33 (7), 1699–1706. Katz, R.W., Parlange, M.B., 1995. Generalization of chain-dependent processes. Application to hourly precipitation. Water Resour. Res. 31 (5), 1331–1341. Kuczera, G., 1997. Efficient subspace probabilistic parameter optimisation for catchment models. Water Resour. Res. 33 (1), 177–185. Lall, U., Rajagopalan, B., Tarboton, D.G., 1996. A nonparametric wet/dry spell model for resampling daily precipitation. Water Resour. Res. 32 (9), 2803–2823. Lloyd, E.H., 1974. What is, and what is not a Markov chain?. J. Hydrology 22, 1–28. Onof, C., Wheater, H.S., 1994. Improvements to the modelling of British rainfall using a modified random parameter Bartlett– Lewis rectangular pulse model. J. Hydrol. 157, 177–195. Rodriguez-Iturbe, I., Cox, D.R., Isham, V., 1987. Some models for rainfall based on stochastic point processes. Proc. R. Soc. Lond. A 410, 269–288. Rodriguez-Iturbe, I., Cox, D.R., Isham, V., 1988. A point process model for rainfall: further developments. Proc. R. Soc. Lond. A 417, 283–298. Velghe, T., Troch, P.A., De Troch, F.P., Van de Velde, J., 1994. Water Resour. Res. 30 (10), 2847–2857.