- Email: [email protected]
Modelling of British rainfall using a random parameter Bartlett-Lewis Rectangular Pulse Model
Journal of Hydrology, 149 (1993) 67 95 0022-1694/93/$06.00 © 1993 - Elsevier Science Publishers B.V. All rights reserved
67
[4]
Modelling of British rainfall using a random parameter Bartlett-Lewis Rectangular Pulse Model Christian Onof*, Howard S. Wheater Dept. of Civil Engineering, Imperial College of Science, Technology and Medicine, London SW7 2BU, UK (Received 26 June 1992; revision accepted 27 January 1993)
Abstract The performance of a stochastic model, in which Poisson cluster processes represent the arrival of rain cells within a rainstorm, is assessed for hourly point rainfall data from Birmingham, UK. The BartlettLewis Rectangular Pulse Model is extended to include random cell duration, and its performance is compared with that of the basic model analysed previously. Both models perform well in reproducing statistics of the rainfall depth distribution at different time-scales for individual months, but the randomised model also enables a good reproduction of the proportions of dry periods of different durations. Performance is improved by an optimisation method which addresses issues of parameter identification. However, some deficiencies remain with respect to extreme values. Finally, the seasonal variability of parameters is represented on a daily, rather than monthly basis, with comparable results.
Introduction Background M o d e l s o f p o i n t rainfall time-series h a v e p o t e n t i a l a p p l i c a t i o n to a r a n g e o f h y d r o l o g i c a l p r o b l e m s , for example, the g e n e r a t i o n o f rainfall across a r a n g e o f time-scales for h y d r o l o g i c a l design a n d the d i s a g g r e g a t i o n o f large timeinterval d a t a for s h o r t d u r a t i o n a p p l i c a t i o n . Because o f the c o m p l e x i t y a n d s t r o n g d e p e n d e n c e u p o n initial c o n d i t i o n s o f the p r e c i p i t a t i o n process, a s t o c h a s t i c a p p r o a c h is likely to be preferable to a p u r e l y physical m o d e l (Smith, 1981). A time-series a p p r o a c h is n o t a p p r o p r i a t e for small time-scales such as the h o u r l y time-scale because o f the large n u m b e r o f p a r a m e t e r s r e q u i r e d (Pattison, 1965). A l t e r n a t i v e l y , rainfall events c a n be m o d e l l e d s e p a r a t e l y * Corresponding author.
68
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67 95
from dry spells. This can be achieved by using a point process to generate the arrival of rainfall cells while modelling their duration and depth with independent distributions, assuming, for example, a simple representation of a rainfall cell as a rectangular pulse. Recently, the use of Poisson cluster processes in stochastic modelling of rainfall has been investigated by a number of authors. The cell arrivals are modelled by a Poisson cluster process, i.e. storm arrivals form a Poisson process and a cell arrival distribution is assigned to each storm; the depth and duration of the cell are modelled by exponential distributions. The cluster process may be a Neyman-Scott process for which the number of cells in a storm is randomly distributed (e.g. Poisson or geometrically) and the cell arrival times are exponentially distributed; this model was examined by Rodriguez-Iturbe et al. (1987a,b) for US data and Cowpertwait (1991) for British data. Alternatively the Bartlett-Lewis Rectangular Pulse Model (BLRPM), for which the storm duration is exponentially distributed and cell arrivals are Poisson distributed, was applied to a single month of Denver data by Rodriguez-Iturbe et al. (1987a,b). Although this model proved efficient in reproducing second-order properties of the depth distribution at all time intervals considered (1 to 24h), a major deficiency was its inability to reproduce the proportion of dry periods. Rodriguez-Iturbe et al. (1988) and Entekhabi et al. (1989) overcame this problem by randomising a model parameter: the exponential parameter of the cell duration distribution was chosen as randomly varying from storm to storm, with a gamma distribution. The present study is motivated by the need to examine the applicability of a randomised Bartlett-Lewis model to hourly British data and the ability of the model to reproduce seasonality. Moreover, problems of parameter identifiability, observed by Rodriguez-Iturbe et al. (1987a), are dealt with by an optimisation procedure.
Bartlett-Lew& Rectangular Pulse Model (BLRPM) In a previous paper (Onof and Wheater, 1993), the Bartlett-Lewis model was examined for the modelling of British rainfall using 38.5 years of hourly data from the Elmdon raingauge (Birmingham, UK). In this model, storms arrive according to a Poisson process of parameter A; within each storm, cells arrive according to a Poisson process of parameter /3; with each cell arrival is associated a rectangular pulse which has an exponentially distributed duration (parameter r/) and an exponentially distributed depth of mean #x; the storm has an exponentially distributed duration (parameter 7).
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
69
The parameters n and ~ were introduced by Rodriguez-Iturbe et al. (1987a), where
and
The B L R P M is thus a model with five parameters (A, n, ~, 77, #x). The fits of this and the randomised model discussed below are assessed by comparing both analytical and simulated values of certain variables characteristic of the rainfall process with their historical values. The notation from Rodriguez-Iturbe et al, (1987a,b, 1988) is used: y~h) represents the ith value in the time-series of h hourly rainfall depths. The characteristic variables are (for h = l, 6, 12 and 24h):
(1) the average rainfall depth at time-scale h hours, E[Y{h)]; (2) the standard deviation of the same, S[y~h)]; (3) the autocorrelation of lag k of rainfall depths for time-scale h hours, (h) (h)
A(k)
__cov[Yi
,Yi+k]
{S[y~h)]2}
(4) the average duration of an event in the hourly data, row, where an event is a wet period preceded and followed by at least one dry time-interval; (5) the mean inter-event time in the hourly data, md; (6) the standard deviation of the same, Sd; (7) the average number of rainfall events per month, rnn; (8) the proportion of dry h hour time intervals, p(h). It should be noted that these variables are all features of the discrete-time process for a given time-scale. Using analytical expressions which are functions of the model parameters, we can derive features of the aggregated process which approximate the discrete-time characteristics; these variables ! ! are given primed nomenclatures: rn~, s~, row, mn and p(h)'. Thus the distribution of lET', which is the time separating the end of an event from the beginning of the next one, is defined by Prob(IET' = h) = Prob[(t, t + h) dry, (t + h, t + h + k) wet I (t - e, t) wet for e > 0] and m~ and s~ are the mean and standard deviation of this distribution. The Bartlett-Lewis model performed well at all levels of aggregation for the reproduction of the first- and second-order statistics of rainfall depth, but not
70
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
for the temporal properties of rainfall events. An optimisation procedure was suggested, allowing both for a better identification of the model parameters and improving the reproduction of the time distribution of rainfall events. However, the main remaining shortcoming observed was the fact that it produced too few wet periods at all time-scales. This overestimation of p(1), p(6), p(12) and p(24) had already been observed by Rodriguez-Iturbe et al. (1987a) for Denver summer rainfall.
Randomisation of the Bartlett-Lewis Rectangular Pulse Model (BLRPM)
The randomised model One solution to the overestimation, as adopted by Rodriguez-Iturbe et al. (1988), is to randomise the parameters characteristic of the cellular structure of the storm. The parameter rI of the exponential distribution of the cell duration is now chosen to be gamma distributed with shape parameter c~ and scale parameter l/u: /,,a r/c~ 1
fr(r/)-
r(a)
e-~
The parameters ~ and ~bare chosen as fixed; therefore,/3 = ~- 7/and 7 = ~b- ,/ are respectively U(c~, ~/u) and P(~, ~b/u) distributed. Thus, the whole of the cellular structure varies from one storm to the other, since for each new storm, a new value of 7/(which determines/3 and 7) is picked from the P(~, l/u) distribution. This should give the model greater flexibility to represent the proportion of dry periods at all time-scales. Finally, A and #x are defined as for the Bartlett-Lewis model, as the rate of arrival of storms and the average cell depth. The model therefore has six parameters (A, #x, c~, u, ~, ~b) and is referred to as the randomised model.
Model equations The second-order properties of the continuous-time randomised process can easily be deduced from those for the Bartlett-Lewis model by replacing the average period of activity of a storm (#T), by its expected value in the randomised process. This is 1
E(#T)
= c)
U
1
1
-(1-vt)e-XU(l-t)]+~ •/
d
0
0
-1
c~--u1
71
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67 95
which may be approximated by (for n and 05 small) .05 1 [1 ÷ 05(t~ ÷ 05/2) - 105(505t; + t~2 ÷ 205 2) E ( # T ) ~ ct -- 1
+ 1 05(4t;3 ÷ 31t;205 + 99t~052 + 36053)]
(1)
The first- and second-order properties of the aggregated process are given in Rodriguez-Iturbe et al. (1988). The equations are E[Y} h)] =
Ahu#x 1 + ~;/05
(2)
a-1
var[Y~ h)] = 2 A l [ ( a - 3)h. 2-~ - u 3-`~ + ( . + h) 3-~]
(3)
- 2 AZ[05(a - 3)hu 2-~ - u 3-~ + ( . + Oh) 3-~]
cov[r} h) v(h)l --i+kJ
=
Al{[u+(k+
1)h]3-~-2(u+kh)3-o+[u+(k
1)h] 3-~}
- A2{[u + (k + l)05h] 3-~ - 2(. + k05h) 3-~ + [. + (k - 1)05h]3 ~} where A1 =
(a-
A#c"~ 1)(a-2)(a-
3)
]
E(:) +
and A2=
A#c~#2u ~ 052(052 __ 1 ) ( ~ -- 1 ) ( a - -
2)(a - 3)
In these equations, X is the cell depth. The proportion dry is also given by Rodriguez-Iturbe et al. (1988):
p(h)' = exp { - ) ~ h - , ~ T ÷ / ~ G ; ( 0 , 0)
[
1°-I} 05+~
(4)
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
72
where @(0,0) is defined in that paper. This is approximately (for n and 0 small)
p(h)' ..~ exp
-
Aq~-l// 1 + 0 n +
Ah
-- 0(50t~-k-
-1- 262)
A~_I v + ~1 0(4n3 + 31ecz~b+ 99n02 + 3663 ) + - 72 c~- 1 x ( 1 - n - 0 + S n 0 +3
02 + ~ 2 ){~__~___£n+ 0 - -n- ~ Iu + ( g +u0 ) / ~ 1~-1})
(5) Equations (2)-(5) are used for the fitting of the model.
Other character&tic variables Inter-event properties The inter-event time properties may be derived in a similar fashion to those of the non-randomised Bartlett-Lewis model, by first calculating the probability that the end of a wet spell is followed by a dry period of duration at least h: PWD'(h) = Prob [(0, h) dry [ t = 0 is the last instant of a wet period] The details of the calculation are shown in the appendix. The result is I (~-7)1 u _ ~ e-Ah @" Ik PNCB(h) k-1 PWD'(h) = 1 - E (u+ u) ~ ~ ~-~I2~33 3
+ E(~)e-Ah k~=olkPNCB(h)k
(6)
where C', I0,11, I2 and 13 are defined in the appendix. The probability distribution of the inter-event time is then given by Ji PWD'(t + Prob(IET' _> n) =
n)f~(t)dt
i'l
]0 PWD'(t +
1)f~(t)dt
where f~ is the density function the distribution of the time t(0 < t < 1) from the beginning of the dry period to the first clock hour; 1 - t is
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67 95
73
therefore the time from the last clock hour in the storm to the end of the storm. This distribution is complex: if there were only one cell alive at the beginning of the last hour, then for a given storm, this time would have the distribution 77e-,~(1-t) -1 ~ - g
fT,,,(t) -
However, there can be more than one cell and new cells can be born; we shall therefore use as approximation for this expression: Prob(IET' > n) ~
Ii PWD'(t + n) dt PWD'(1)
The mean inter-event time m~ is obtained by summing Prob(IET'_> n) from n=lto~: m~ = Z P r ° b ( I E T ' - > n) n=l or
,
J' PWD'(t + n)dt
md ~ n=l 0 PWDI(1)
(7)
and similarly s 'a ~ ~/ 1 + Z~( 2 n -
1) I1 PWD,(t + n) d t - m~t2 0 PWD'(1)
n=2
(8)
Number o f events and average event duration
As for the non-randomised model, these may be derived from the expression for m~: ,
mn =
p(1)' x 24 × N(i) ,
(%)
md
where N(i) is the number of days in month i and ,
1 -p(1)'
,
mw-
p(1)'
× md
(lOa)
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
74
It must, however, be noted that these two variables may also be derived without using m~, as follows: !
m, = [ p ( l ) ' - p ( 2 ) ' ] x 24 ×
NM
(9b)
and mtw = [ 1 - p ( 1 ) ' ] / [ p ( 1 ) ' - p ( 2 ) ' ]
(10b)
Formulae (9b) and (10b) will be preferred to the others when the probability PSTOV of overlap of storms becomes too large: PSTOV _> 20%. Parameters and results for the randomised model
Parameter estimation As for the Bartlett-Lewis model, parameter estimation consists of the identification of certain variables characteristic of the precipitation process, for which Eqs. (2)-(5) above are solved using a modified version of the Powell hybrid method. Several alternative sets of variables could be selected. They must include features both of the depth process and of the proportion of dry periods; thus, two sets of variables were examined: (1) the R1 set
{E[Yi(U], var[Yi(1)],cov[Y~ ') ~v(1)],p(1),cov[Yi (6), yi(6+l],p(24)} i+1 (2) the R2 set
{E[Yi (l)1, var[Yi (1)], cov[Y} 1), V~i+,],p(1), (1) var[Yi(6)],p(24)} The two methods were applied to individual months of the data record and yield similar values of A and #x, but the other parameters are very different according to the method. It is moreover easy to see that parameters such as u and ~ are not well identified, since for very different values of u and c~, it is possible to obtain parameter sets yielding characteristic values within 5% of the historical values. This problem of identification is examined below. For the m o n t h of July, no parameters were obtained for the R1 set because the model was incapable of reproducing both the relatively high value of p(24) and the relatively low value of the six-hourly lag-1 autocorrelation A(1). Parameters for the R2 set are shown in Table 1.
Analytical results for R2 variables The second-order properties of the rainfall depth are well reproduced by the
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67 95
w i
o
y_ .=_
0
0
i
e. o
75
76
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67--95 0.50
--
0.40
--
6-hourly
-- - - -
A(1)
analytical values
BLRRPM Historical
data
- ~ ' ~ -~ / ~ ~
r -
/~--
\\//
0.30--
0.20
--
0.10---
I
I
I
2
3
,~
[ ~
l
6
I
7
I
s
I
9
I
~e
I
~
[ ~2
Months
Fig. 1. Analytical autocorrelations of lag-1 of six-hourly rainfall for the randomised model.
Bartlett Lewis Randomised Rectangular Pulse Model (BLRRPM) (Figs. 1 and 2), with a tendency to overestimate some of the autocorrelations at timescales larger than 12 h. The chief improvement brought about by the randomisation is an adequate representation of the proportion dry even at those levels of aggregation not used in the fitting (Fig. 3(a),(b)). The inter-event time properties are somewhat improved on the two methods used for the BLRPM (Fig. 4), which were based on the identification of the sets {E[y~I)],
var[Yi(l)],c°v[Yi 0), Yi~l], var[yi(6)],c°v[Yi (6) , ~i+lJJY (6)]~
and {E[ ri(l)], var[ yi(l)], cov[ yi(1) yfill] 'cov[ yi(6) yi(6+l],cov[ ri (12>, yi(l~)] } denoted as H = {1;6} and H = {1;6; 12} respectively. Nevertheless, clear discrepancies remain in the reproduction of rod, mw and ran.
77
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
~.se-
24-hourly
- -- -
0.40--
B L R R P M analytical Historical data
A(1 )
values
0.30--
.¢ 0.20
/
--
0.10--
0.~
I 2
r 3
I
f
l
I
1
I
4
~
6
7
e
9
I le
I ~1
1 1~
Months Fig. 2. A n a l y t i c a l a u t o c o r r e l a t i o n s o f lag-1 o f 2 4 - h o u r l y rainfall f o r the r a n d o m i s e d m o d e l .
Simulated results for R2 variables These results were confirmed by simulation of 144 years of hourly data. The discrepancies between analytical and simulated results do not exceed 10% apart from Im~w- mwl and ]mln- m,I for which the discrepancies sometimes reach 15% because of the approximations involved in their analytical expressions. A small discrepancy (smaller than 3% for p(1) and smaller than 10% for p(24)) remains between the analytical and simulated values of the proportions dry. This is analysed in Onof (1992).
Parameter optimisation Parameter identification with simulated results As was observed above, only A and #x are properly identified by the R1 or R2 methods. This identification problem was also observed by RodriguezIturbe et al. (1988). Moreover, the features characteristic of the event distribution are not adequately reproduced. Therefore, as in the case of the
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67 95
78
6-hourly proportion dry BLRRPM
++I //
0,
iI
I III /
75
~k"l, xx
/
;
]II/
x
1~i%x ;I
]
x,, ]
.... - -
0.70 I
R2 method mlto~l.cL1,
.....
data
I
I
I
1
[
I
I
I
I
I
2
3
B
6
7
8
9
10
11
12
Months
(a)
tg-hourly proportion dry
0.7B --
BLRRPM
0.70
/
I
x
l
x
tillI
tlttl
0.60
.... - -
0.m5 (b)
I~ 'method llimt,o r i c a l data
I
I
I
I
I
I
I
I
t
I
I
2
3
~1
B
6
7
8
9
10
11
12
Months
Fig. 3. (a) Analytical six-hourly proportion dry for the randomiscd modcl; (b) analytical 12-hourly proportion dry for the randomiscd modcl
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
Mean Inter-Event
79
Time
BLRPM/BLRRPM m o d e l s \
/\ /
ES--
\
/
\
f 60-k.
\
/
\
\ / \/
/ /
.o
\
/ \
/
m
/
\ \ k~
/ f
~
/
....
/
-
J
analyt/cal results Historical data B L E P M H = ~ I ; 6 ~ analytical results ~
--
xxx~ ~
25-
2e 1
2
3
4
~
6
7
8
9
10
11
12
Months Fig. 4. Analytical mean inter-event time for the non-randomised and the randomised models.
BLRPM (Onof and Wheater, 1993) which also gave rise to a parameter identification problem, an optimisation was carried out with the following objective function:
/[md(H ) -- md(i)] 2
[mn(H) - mn(i)] 2
where md(H), m,(H) and rnd(i), m,(i) are respectively the historical and simulated (iteration i) values of the two variables md and m, (mw is not used since it is a function of the other two). To investigate and improve parameter identifiability, a poorly identified parameter was selected and held constant while the five remaining parameters were optimised as before (by solving a set of five equations). The value of di was then computed; this procedure was repeated for different values of the fixed parameter, until an optimal value of the fixed parameter had been found. The five variables selected for the parameter estimation and the parameter which is to be fixed remain to be chosen. The following investigations were first conducted on the month of April.
C. Onof, H.S. Wheater/Journal of Hydrology 149 (1993) 67 95
80
Using the set of variables
{E[Yi(1)],c°v[Yi (1) , Y(1)I,p(1),c°v[Yi k
, Ji+l],p(24)
and fixing a or u, no optimum was obtained. However, if instead of the autocovariances, the autocorrelations are used, i.e if the set
{ E[Yi(I)],c°v[Yi(')'YiO+l]{S[Yi(1)]} 2
,p(1), c°v[Yi(6)'yi(6+l]}~,__ ~)]~5 ,p(24)
is used, then: (1) if a is fixed, the convergence of the Powell hybrid method is difficult to obtain; (2) if u is fixed, then we obtain an optimal value//opt (as verified by simulation of 200 years of data) which corresponds approximately to a minimum for the number of cells per storms #c. This corresponds to a small value of E(3) as had been found for the BLRPM optimisation but unlike the BLRPM case, EQ/) also reaches a minimum in the neighbourhood of the optimum. Applying this procedure to all months and choosing for initialisation values of the parameter sets those obtained by an examination of the parameters from grouped months for method R2, we obtain the parameters of Table 2. (The grouped month parameters correspond to larger quantities of data (three or four months) and although of limited significance, are useful in providing a value with which to initialise iterations for any given month in that group.) Table 2 Optimal parameters for the randomised model Month
A (h I)
#x (mm h-l)
a
u (h)
n
0
#~
January February March April May June July August September October November December
0.0261 0.0221 0.0212 0.0219 0.0226 0.0194 0.0111 0.0231 0.0215 0.0206 0.0282 0.0317
1.0329 1.1851 1.2122 1.4182 2.3427 2.0769 4.1869 3.3444 1.5922 1.2486 1.1771 0.7938
2.5032 2.6175 2.1800 2.3390 3.0652 3.3702 2.2894 2.7016 6.6513 4.0159 5.1272 3.6535
1.0600 0.9000 0.6000 0.4700 0.5000 0.9000 0.2476 0.4331 3.8000 2.0000 3.6000 1.1328
0.2069 0.2431 0.1662 0.2496 0.1816 0.1845 0.0421 0.1593 0.2476 0.2728 0.2584 0.8483
0.0669 0.0639 0.0386 0.0518 0.0366 0.0467 0.0064 0.0408 0.0986 0.0834 0.1188 0.1297
4.0941 4.8076 5.3047 5.8216 5.9638 4.9512 7.5515 4.9009 3.5099 4.2711 3.1740 7.5408
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
81
Analytical determination of the optima The same optimisation process was carried out using the analytical values of the characteristic variables, i.e. with the objective function
/[md(H ) - m~t(i)]2 di : v
md ~I--I7
[mn(H) - mrn(i)]2 -~
mn ( H ) 2
(12)
where ma(t), , . mn(~) , . are the analytical (iteration i) values of the two variables / ! md, mn. The optimal values of u which are obtained are very close to those obtained with the simulated results above: //opt ~'~ /Jopta
is the analytical optimum. This provides us with an analytical method for the approximated determination of Vop t. Figure 5(a)-(c) (where RO is the optimised version of the randomised model) shows the improvement brought about by this optimisation in the reproduction of the variables md, mw and mn. Moreover, the optimal distribution of the simulated inter-event time population was compared with its historical counterpart using a Student test for the average, a Fischer test for the standard deviations and a M a n n - W h i t n e y test, and showed no significant discrepancy for all months except July and December. For these two months, the above method did not provide an optimum. Table 3 shows the analytical values of the characteristic variables for the parameters in Table 2. where/Jopta
Extreme values of the model
Hourly data For most months of the year, the extreme-values of the hourly rainfall distribution are overestimated by the R1 parameters for all return periods and by the R2 parameters for return periods longer than the data set (38.5 years). Figure 6(a),(b) shows the extreme value distributions for the historical and the simulated data for the months of January and June. The results for the optimised parameters are similar to these.
Daily data For daily rainfall, for all months apart from January and November,
82
C.
Onof, H.S. Wheaterl Journal of Hydrology 149 (1993) 67 95
Average E v e n t D u r a t i o n
4--
," - - . -
BLRRPM
/ ~xxx
,
/
i#
/
3--
It
i
/ // \',
/
/ I
2--
\,
I
I
I
', I
,'
v_---" .... - -
-
'
I
I
2
I
3
I
4
I
5
I
6
(a)
7
I
--
R2 s/mulated H/storlcal data RO simulated
I
8
I
9
1
10
I
11
12
Months
Mean I n t e r - E v e n t
40-
Time
BLRRPM
Ixxx\ !
x
; ii m
',
I ," l
..--. x s I
\
~
",
/
\
I
t
It
x\
A
/
ix ii
~
(b)
I \
I 2
~
\/
I
I 3
~
/
I 4
I El
# / / z"I
'~
,
/
\\1
'~ \ 1
I 6
I. 7
I IB
I ?
I le
I 11
I 12
Months
Fig. 5. (a) S i m u ] a t e d o p t i m a l a v e r a g e e v e n t d u r a t i o n ; (b) s i m u l a t e d o p t i m a l m e a n i n t e r - e v e n t time; (c) s i m u l a t e d o p t i m a l n u m b e r o f events.
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
Number
83
of E v e n t s
BLRRPM
2K-v
x
¢
20--
.... - -
1B
[ 2
(c)
1~ ~ul.at, lllrl~rlcs1
RO simulated
--
ed
data
I
I
I
I
I
I
I
I
1
3
4
B
6
7
5
~, ?
10
11
12
Months
F i g . 5. C o n t i n u e d .
method R2 reproduces the extreme values of short return period but diverges from the historical data for larger return periods; method R1 often overestimates the historical values for all return periods (Fig. 6(c),(d)). Again, the results for the optimised parameters show no improvement.
Daily representation of the randomised parameters From the analysis of monthly data, twelve sets of parameters have been obtained for each of the estimation methods. Two problems arise: (1) is it possible to provide a smoother representation (e.g. polynomial) of the variation of the parameters over the year, eliminating the arbitrary monthly grouping of data?; (2) can this variation be represented in such a form (e.g. Fourier analysis) that the coefficients of this representation might be used for a regionalisation of the parameters?
84
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
H o u r l y E x t r e m e Values
12--
J~uauary; BLRRPM
11--
//
10--
//
/
/
/
•1
9--
/• B--
/
// / 6--
/
/
i"
/•
/
/
/
/ /
/~ F { / /
/ J // /~
5~ 4~ a-
¢ S . /
....
R2
SI----
I
-1.0
I
0.0
I
1.0
I
2.0
Return
period
1 2
(yrs)
I
3,0
GumbeI
reduced 1
1 10
5
]
I
4.0
5.0
6.0
variate
1 ZO
1
I
100
200
(a)
~-
H o u r l y E x t r e m e Values June;
BLRRPM
2~-
//
/"
17
18-
i 16 -
/
-
F¢
14-
./ /
~f 1 2 - /
10 -8--
,14
/
./
-/
/-
/
~
11 / FI
/~
i
/I
/
/
/
6--
2 0 -1.0
I
l
0.0
1.0
I
Return
period
(yrs)
[ 2
I
I
reduced
variate
2.0
GumbeI i 5
3.0
[ S.O
4.0
I
r
10
20
I 6+0
I
T
I~
200
(b) F i g . 6. (a) H o u r l y e x t r e m e v a l u e s f o r J a n u a r y ; (b) h o u r l y e x t r e m e v a l u e s f o r J u n e ; (c) d a i l y e x t r e m e v a l u e s f o r J a n u a r y ; (d) d a i l y e x t r e m e v a l u e s f o r J u n e .
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67 95
~0-
Daily Extreme January;
Values
BLRRPM 7
6o-
7
50-
f
j"
7
I~1
.f
~0--
J" j-"
./" 7
/ I SO--
7
J" "''"
~
lO-
o
I 0.0
-1.0
I 1.0
I 8.0
I 3.0
Cumbel
Return
period
(yrs)
I 4.0
reduced
I 5.0
1 6.0
variate
1
1
1
1
2
5
10
20
t
l
100
200
(c)
Daily Extreme
100-
June;
Values
BLRRPM /-
90-
J J J
00-
J J 70-
J f J
i
J
Et
J J j j
i ~ ~"
10-0
I 0.0
-I .0
Return
period
(d) Fig. 6. Continued.
I 1.0
(yrs)
I 2
I I I 2.0 3.0 4.0 GumbeI reduced variate I 5
I 10
I 20
i 5,o
T
1~
i 6.0
T
200
85
86
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67--95
~
.
.
.
.
~
.
.
.
.
C. Onof, H, S. Wheater Journal of Hydrology 149 (1993) 67-95
¢,.q
t"q
~D
H
[..,
O
[-.,
az
~
87
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
88
1.0
I 31
I 61
I 91
I 121
I 1 I IE;I 181 211 Days
I 241
I 271
I 301
l I 331:361
I 241
I 271
I 301
I 331
(a)
BLRRPM daily
0.7 0. a0. ¢L,
0. 0. 0. 0.1
0"0
I 1
I 31
1 61
I 91
I 121
1 151
I I 181 211 Days
I 361
(b) F i g 7 (a) D a i l y variation o f A for the randomised mode] (R2); (b) daily variation o f e~for the randomised mode] (R2); (c) daily variation o f ~ for randomised mode] (R2).
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
89
1:R.O-12.0-
~'~
BLRRPM daily
11 . 0 -
10.09.@o,m
8.0-
~
7.e6.05.04°0:3.02.01,0
I 31
I 61
I 91
I 121
(c)
I 151
I I 181 211 Days
I 241
[ 271
I 3~1
I 331
I 361
Fig. 7, Continued.
Polynomial curve fitting F o r each parameter ah[h E ( 1 , . . . , 6 ) ] , polynomial coefficients for an expression of the parameter in a Chebyshev series form are calculated: n=k
~'~h,k = ½ao(ah)To(t) + Z an(ah)Tn(t) n=l
where
T~(t) = cos(n arcos t) = ½[(t + iv/1 - t2)" + ( t - ix,/1 - t2) "] The approximation of a h by Y]h,k is carried out by a least-squares method. On the whole, for parameters A, n and 05, a representation of degree 5 is sufficient (Fig. 7a-c), whereas for u, a polynomial of degree 7 appears necessary. F o r the other two parameters, #x and a, it is preferable to choose a degree k = 12. The parameters thus s m o o t h e d produce very satisfactory results for simulations of the m o n t h l y characteristic variables. They have the additional advantage of yielding more realistic results for periods which are not clusters of the calendar months.
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67 95
90
Table 4 Fourier coefficient moduli for the randomised model, N = 12 (R2)
Degree
A (h-1)
#x (mmh-l)
a
u (h)
~
0 1 2 3 4 5 6
0.0803 0.0126 0.0064 0.0027 0.0040 0.0018 0.0028
5.9714 2.4249 0.9357 0.3420 0.5920 0.5912 0.4482
10.2632 0.8729 0.3700 0.0489 0.5950 0.5148 0.4496
2.4343 0.6602 0.1886 0.0950 0.4232 0.4818 0.0726
1.6591 0.6273 0.0909 0.0379 0.1700 0.1314 0.0674
0.2412 0.0885 0.0232 0.0074 0.0079 0.0415 0.0151
Fourier analysis of the continuous-time representation A Fourier analysis of the parameters may also be carried out to isolate the predominant periodicities of the monthly parameter values. The Fourier coefficients for parameter crh are calculated in the complex form: l
N-I
{
.27rjk'~
F~h,k -- vT~ j~=oCrh( j ) exp ~-' --~ -) Choosing N = 12 enables us to represent the monthly components which dominate the picture (the coefficients decrease w i t h j for A and #x) and longer trends (two to three monthly), as can be seen from Table 4. This provides us with a representation of the parameters which could be used for a spatial generalisation of the model; the coefficients of lower frequency are expected to vary less than those of higher frequency and thus represent underlying seasonal trends which may not be as location specific as variations of higher frequency. Conclusion
The extension of the Bartlett-Lewis Rectangular Pulse Model (BLRPM) to include a random cell duration is shown to improve the simulation of British rainfall data. Good reproduction of dry periods is achieved at all time-scales (hourly to daily), in addition to the previous reproduction of daily rainfall depth distributions. This corroborates results obtained by Rodriguez-Iturbe et al. (1988) for data from Denver, CO, USA. Problems of parameter identification are apparent, and a two-stage optimisation is proposed which allows for optimal reproduction of features of the
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
91
event distribution and which can be carried out analytically. Representation of the parameters as daily variables avoids the artificiality of the monthly subdivision of the year and allows for the possible generalisation of parameters through regional analysis. Two problems remain which require further examination: (1) the excessive values of some of the autocorrelations for lags greater than 12 h, in particular for the optimised parameters; (2) the inadequate reproduction of extreme values for return periods greater than the length of the data set. These problems can be addressed through an extension of this model, e.g. by the introduction of a jitter and by using a gamma distribution for the cell depths; although additional parameters are expected to give rise to further problems of identifiability. These issues will be examined in a subsequent paper.
Acknowledgments The helpful advice and comments from Sir David Cox and Valerie Isham during the progress of this work are gratefully acknowledged.
References Cowpertwait, P.S.P., 1991. Further developments of the Neyman-Scott clustered point process for modelling rainfall. Water Resour. Res., 27(7): 1431-1438. Entekhabi, D., Rodriguez-Iturbe, I. and Eagleson, P.S., 1989. Probabilistic representation of the temporal rainfall process by a modified Neyman-Scott Rectangular Pulses Model: parameter estimation and validation. Water Resour. Res., 25(2): 295-302. Onof, C., 1992. Stochastic modelling of British rainfall using Poisson processes. Ph.D. Thesis, Imperial College, London. Onof, C. and Wheater, H.S., 1993. Optimisation of the Bartlett-Lewis Rectangular Pulse Model for British Rainfall (submitted). Pattison, A,, 1965. Synthesis of hourly rainfall data. Water Resour. Res., 1(4): 489-498. Rodriguez-Iturbe, I., Cox, D.R. and Isham, V, 1987a. Some models for rainfall based on stochastic point processes. Proc. R. Soc. London, Ser. A, 410: 269-288. Rodriguez-lturbe, I, Febres de Power, B. and Valdes, J.B., 1987b. Rectangular pulses point process models for rainfall: analysis of empirical data. J. Geophys. Res., 92D8: 96459656. Rodriguez-lturbe, I., Cox, D.R. and Isham, V., 1988. A point process model for rainfall: further developments. Proc. R. Soc. London, Ser. A, 417: 283-298. Smith, J.A., 1981. Point process models of rainfall. Ph.D. Thesis, The Johns Hopkins University, Baltimore, MD.
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
92
Appendix To calculate PWD(h) ', let us first suppose there is only one storm alive at t = O:
One storm alive at t = 0 If we suppose this, then if P1 is the probability that an interval o f length h is dry conditional upon the first instant wet P1 = Prob [(0, h) dryIt = 0 is the last instant of a wet period]
E
c
I C e -'~h × e -~oh + E ~-s [ C e -'~h /
where C is the number o f cells in the storm and C' the number of those whose cell end times are followed by a dry period. This obtains since t = 0 m a y be any cell end time within the storm with probability [E(C'[C) - 1]/E(C'[C) or the last cell end time with probability 1/E(C~[ C)', where C' is the number of cell end times which are followed by a dry period. To calculate E [1/C'[C] as a function o f C, the number o f cells in the storm, we need to compute the probability P S U P that two consecutive cells overlap. :3C
P S U P = / nr/e -~u e -~u e -*~u du d
0
because for each time u, we multiply the probabilities that the end of the cell is not reached, that the cell is still alive and that a cell arrives at u. We obtain P S U P = n / ( n + ~ + 1). The probability that a cell end time is dry is approximately that o f two consecutive cells not overlapping if we neglect the probability PSUP3 of a first cell overlapping with the two others arriving later which are mutually non-overlapping. This is 0(3
X
X--U
X
PSUP3=Irle-'TXe-~'TXdxlnrle-'~'TUdulrle-'WdvlmTe-'~'mdw 0
0
0
u+v
which computes to ~+ 1 PSUP3 =
2(1 - ~;)(t~ + 1)(2~ + 1 + ~b)
1 ~+1+05
+-
2 (1 - t~)(~ + 1)(t~ + 2 + ~b)
1 2(~ + 1)(1 + ~h)
For the parameter values obtained, we have PSUP3 < 0.05.
C. OnoJ; H. S. Wheater/ Journal of Hydrology 149 (1993) 67 95
93
Therefore, the distribution of C' conditional upon C is binomial:
P(C' -- C - q l C ) where Cab =
PSUp)c-I-qPSUP q
= cc~-~-q(1 -
a!/[b!(a - b)!] and we therefore obtain
E(1"-1TI C) ~}--~CPc-1 (1 - PSUP)PPSUPP p=0
p+l
This yields the probability P1 for one particular storm with C cells; to obtain the probability P2 that an interval of length h is 'dry' within one storm knowing that it starts at the end of a wet spell, we must decondition PI upon C and take expectations over r]: E ( )~-5 l
"~
ZZCcLl°CC-'
(I-PSUP)PPSUPP ( I3 ~"7
c=lp=0 P 2 = E(P1) ~ E,j
and this yields p2~ [I_E(1)]
p+ 1 1- E
\ c )1
u~
• ~(~ + u)
\ ~ J
e -ah x e -~vh +
e -~h +
E(1).e-,h
where E ( 1 / C ' ) has to be evaluated numerically.
Multi-storm case P2 above is the sum of two terms: P21 =
1-
. ( u +/Ju ) ~ e
-Ah
and
If there are k storms active at t = 0, then conditional upon this P21 must be multiplied by [PNCB(h)] k-l, where PNCB(h) is the probability that no cells are born during [0, h] for any o f the storms and similarity P22 must be multiplied by PNCB(h) k. Removing the condition, we obtain PWD'(h) ~ P21 x A21 + P22 x A22
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67--95
94
where ~-, Prob(exactly k storms alive) PNCB(h) k-1 A21 = ~ ~ e or more storms alive) k=l
and O<3
A22 = Z
Prob(exactly k storms alive)PNCB(h) k
k=0
Remark: this equation for PWD~(h) supposes that the distribution of C' is the same as that for the single-storm case. This approximation is acceptable if PSTOV (defined below) is not too large, e.g. if PSTOV_< 20%. Then I k = Prob(number of storms alive = k) = E~ e -A/~
k!
or
1 Ik -/~!r(~) e a/~tk~t/ t ~
lu~e-~t at
0
Practically, the sums in A21 and A22 will be approximated by sums up to k = 3; indeed, the probability of more than three storms superposed at any given instant is 3
PSTOV3~ [E~(iAe-Ate-~Vtdt)] or 3
PSTOV 3 ~
~
dt
which can be evaluated numerically. PSTOV 3 _< 0.01 for the parameter sets obtained below. Thus
,~ ~ IkPNCB(h)k-1 A21
~
I1 + / 2 + / 3 3
A22 ~ Z k=0
IkPNCB(h)~'
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
95
where PNCB(h) is the probability that no cells occur in a given storm during the time interval [0, h]. The calculation of the probability PNCB(h) now follows. The calculation is as in the BLRPM case, but the integral which is obtained must be averaged over all possible values of r/: the storm may end at u (0 _< u < h) and no cells occur during [0,u] (with probability O~le-°~Ue-~'lu),or may end after t = h and no cells occur during [0,hi (with probability e e,lh e ~h). Thus PNCB(h) = E+l
Or/e -4'~u e '+~' du + e -0~h e -'+'Th
We obtain PNCB(h) = E,, ~ O
+ ~ n
e_(~+0)@]]
which may be simplified by using the expression of the Laplace transform of the gamma distribution: PNCB(h)_
0
+
~
[
u
]~
We obtain finally [ PWD'(h)~
(~5)] 1-E
u'~
3IkPNCB(h)k' Z I1+I2+I~
-.Xh
(u+u) ~e
k=l
+ E(l e h , PNCS/h/
67
[4]
Modelling of British rainfall using a random parameter Bartlett-Lewis Rectangular Pulse Model Christian Onof*, Howard S. Wheater Dept. of Civil Engineering, Imperial College of Science, Technology and Medicine, London SW7 2BU, UK (Received 26 June 1992; revision accepted 27 January 1993)
Abstract The performance of a stochastic model, in which Poisson cluster processes represent the arrival of rain cells within a rainstorm, is assessed for hourly point rainfall data from Birmingham, UK. The BartlettLewis Rectangular Pulse Model is extended to include random cell duration, and its performance is compared with that of the basic model analysed previously. Both models perform well in reproducing statistics of the rainfall depth distribution at different time-scales for individual months, but the randomised model also enables a good reproduction of the proportions of dry periods of different durations. Performance is improved by an optimisation method which addresses issues of parameter identification. However, some deficiencies remain with respect to extreme values. Finally, the seasonal variability of parameters is represented on a daily, rather than monthly basis, with comparable results.
Introduction Background M o d e l s o f p o i n t rainfall time-series h a v e p o t e n t i a l a p p l i c a t i o n to a r a n g e o f h y d r o l o g i c a l p r o b l e m s , for example, the g e n e r a t i o n o f rainfall across a r a n g e o f time-scales for h y d r o l o g i c a l design a n d the d i s a g g r e g a t i o n o f large timeinterval d a t a for s h o r t d u r a t i o n a p p l i c a t i o n . Because o f the c o m p l e x i t y a n d s t r o n g d e p e n d e n c e u p o n initial c o n d i t i o n s o f the p r e c i p i t a t i o n process, a s t o c h a s t i c a p p r o a c h is likely to be preferable to a p u r e l y physical m o d e l (Smith, 1981). A time-series a p p r o a c h is n o t a p p r o p r i a t e for small time-scales such as the h o u r l y time-scale because o f the large n u m b e r o f p a r a m e t e r s r e q u i r e d (Pattison, 1965). A l t e r n a t i v e l y , rainfall events c a n be m o d e l l e d s e p a r a t e l y * Corresponding author.
68
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67 95
from dry spells. This can be achieved by using a point process to generate the arrival of rainfall cells while modelling their duration and depth with independent distributions, assuming, for example, a simple representation of a rainfall cell as a rectangular pulse. Recently, the use of Poisson cluster processes in stochastic modelling of rainfall has been investigated by a number of authors. The cell arrivals are modelled by a Poisson cluster process, i.e. storm arrivals form a Poisson process and a cell arrival distribution is assigned to each storm; the depth and duration of the cell are modelled by exponential distributions. The cluster process may be a Neyman-Scott process for which the number of cells in a storm is randomly distributed (e.g. Poisson or geometrically) and the cell arrival times are exponentially distributed; this model was examined by Rodriguez-Iturbe et al. (1987a,b) for US data and Cowpertwait (1991) for British data. Alternatively the Bartlett-Lewis Rectangular Pulse Model (BLRPM), for which the storm duration is exponentially distributed and cell arrivals are Poisson distributed, was applied to a single month of Denver data by Rodriguez-Iturbe et al. (1987a,b). Although this model proved efficient in reproducing second-order properties of the depth distribution at all time intervals considered (1 to 24h), a major deficiency was its inability to reproduce the proportion of dry periods. Rodriguez-Iturbe et al. (1988) and Entekhabi et al. (1989) overcame this problem by randomising a model parameter: the exponential parameter of the cell duration distribution was chosen as randomly varying from storm to storm, with a gamma distribution. The present study is motivated by the need to examine the applicability of a randomised Bartlett-Lewis model to hourly British data and the ability of the model to reproduce seasonality. Moreover, problems of parameter identifiability, observed by Rodriguez-Iturbe et al. (1987a), are dealt with by an optimisation procedure.
Bartlett-Lew& Rectangular Pulse Model (BLRPM) In a previous paper (Onof and Wheater, 1993), the Bartlett-Lewis model was examined for the modelling of British rainfall using 38.5 years of hourly data from the Elmdon raingauge (Birmingham, UK). In this model, storms arrive according to a Poisson process of parameter A; within each storm, cells arrive according to a Poisson process of parameter /3; with each cell arrival is associated a rectangular pulse which has an exponentially distributed duration (parameter r/) and an exponentially distributed depth of mean #x; the storm has an exponentially distributed duration (parameter 7).
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
69
The parameters n and ~ were introduced by Rodriguez-Iturbe et al. (1987a), where
and
The B L R P M is thus a model with five parameters (A, n, ~, 77, #x). The fits of this and the randomised model discussed below are assessed by comparing both analytical and simulated values of certain variables characteristic of the rainfall process with their historical values. The notation from Rodriguez-Iturbe et al, (1987a,b, 1988) is used: y~h) represents the ith value in the time-series of h hourly rainfall depths. The characteristic variables are (for h = l, 6, 12 and 24h):
(1) the average rainfall depth at time-scale h hours, E[Y{h)]; (2) the standard deviation of the same, S[y~h)]; (3) the autocorrelation of lag k of rainfall depths for time-scale h hours, (h) (h)
A(k)
__cov[Yi
,Yi+k]
{S[y~h)]2}
(4) the average duration of an event in the hourly data, row, where an event is a wet period preceded and followed by at least one dry time-interval; (5) the mean inter-event time in the hourly data, md; (6) the standard deviation of the same, Sd; (7) the average number of rainfall events per month, rnn; (8) the proportion of dry h hour time intervals, p(h). It should be noted that these variables are all features of the discrete-time process for a given time-scale. Using analytical expressions which are functions of the model parameters, we can derive features of the aggregated process which approximate the discrete-time characteristics; these variables ! ! are given primed nomenclatures: rn~, s~, row, mn and p(h)'. Thus the distribution of lET', which is the time separating the end of an event from the beginning of the next one, is defined by Prob(IET' = h) = Prob[(t, t + h) dry, (t + h, t + h + k) wet I (t - e, t) wet for e > 0] and m~ and s~ are the mean and standard deviation of this distribution. The Bartlett-Lewis model performed well at all levels of aggregation for the reproduction of the first- and second-order statistics of rainfall depth, but not
70
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
for the temporal properties of rainfall events. An optimisation procedure was suggested, allowing both for a better identification of the model parameters and improving the reproduction of the time distribution of rainfall events. However, the main remaining shortcoming observed was the fact that it produced too few wet periods at all time-scales. This overestimation of p(1), p(6), p(12) and p(24) had already been observed by Rodriguez-Iturbe et al. (1987a) for Denver summer rainfall.
Randomisation of the Bartlett-Lewis Rectangular Pulse Model (BLRPM)
The randomised model One solution to the overestimation, as adopted by Rodriguez-Iturbe et al. (1988), is to randomise the parameters characteristic of the cellular structure of the storm. The parameter rI of the exponential distribution of the cell duration is now chosen to be gamma distributed with shape parameter c~ and scale parameter l/u: /,,a r/c~ 1
fr(r/)-
r(a)
e-~
The parameters ~ and ~bare chosen as fixed; therefore,/3 = ~- 7/and 7 = ~b- ,/ are respectively U(c~, ~/u) and P(~, ~b/u) distributed. Thus, the whole of the cellular structure varies from one storm to the other, since for each new storm, a new value of 7/(which determines/3 and 7) is picked from the P(~, l/u) distribution. This should give the model greater flexibility to represent the proportion of dry periods at all time-scales. Finally, A and #x are defined as for the Bartlett-Lewis model, as the rate of arrival of storms and the average cell depth. The model therefore has six parameters (A, #x, c~, u, ~, ~b) and is referred to as the randomised model.
Model equations The second-order properties of the continuous-time randomised process can easily be deduced from those for the Bartlett-Lewis model by replacing the average period of activity of a storm (#T), by its expected value in the randomised process. This is 1
E(#T)
= c)
U
1
1
-(1-vt)e-XU(l-t)]+~ •/
d
0
0
-1
c~--u1
71
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67 95
which may be approximated by (for n and 05 small) .05 1 [1 ÷ 05(t~ ÷ 05/2) - 105(505t; + t~2 ÷ 205 2) E ( # T ) ~ ct -- 1
+ 1 05(4t;3 ÷ 31t;205 + 99t~052 + 36053)]
(1)
The first- and second-order properties of the aggregated process are given in Rodriguez-Iturbe et al. (1988). The equations are E[Y} h)] =
Ahu#x 1 + ~;/05
(2)
a-1
var[Y~ h)] = 2 A l [ ( a - 3)h. 2-~ - u 3-`~ + ( . + h) 3-~]
(3)
- 2 AZ[05(a - 3)hu 2-~ - u 3-~ + ( . + Oh) 3-~]
cov[r} h) v(h)l --i+kJ
=
Al{[u+(k+
1)h]3-~-2(u+kh)3-o+[u+(k
1)h] 3-~}
- A2{[u + (k + l)05h] 3-~ - 2(. + k05h) 3-~ + [. + (k - 1)05h]3 ~} where A1 =
(a-
A#c"~ 1)(a-2)(a-
3)
]
E(:) +
and A2=
A#c~#2u ~ 052(052 __ 1 ) ( ~ -- 1 ) ( a - -
2)(a - 3)
In these equations, X is the cell depth. The proportion dry is also given by Rodriguez-Iturbe et al. (1988):
p(h)' = exp { - ) ~ h - , ~ T ÷ / ~ G ; ( 0 , 0)
[
1°-I} 05+~
(4)
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
72
where @(0,0) is defined in that paper. This is approximately (for n and 0 small)
p(h)' ..~ exp
-
Aq~-l// 1 + 0 n +
Ah
-- 0(50t~-k-
-1- 262)
A~_I v + ~1 0(4n3 + 31ecz~b+ 99n02 + 3663 ) + - 72 c~- 1 x ( 1 - n - 0 + S n 0 +3
02 + ~ 2 ){~__~___£n+ 0 - -n- ~ Iu + ( g +u0 ) / ~ 1~-1})
(5) Equations (2)-(5) are used for the fitting of the model.
Other character&tic variables Inter-event properties The inter-event time properties may be derived in a similar fashion to those of the non-randomised Bartlett-Lewis model, by first calculating the probability that the end of a wet spell is followed by a dry period of duration at least h: PWD'(h) = Prob [(0, h) dry [ t = 0 is the last instant of a wet period] The details of the calculation are shown in the appendix. The result is I (~-7)1 u _ ~ e-Ah @" Ik PNCB(h) k-1 PWD'(h) = 1 - E (u+ u) ~ ~ ~-~I2~33 3
+ E(~)e-Ah k~=olkPNCB(h)k
(6)
where C', I0,11, I2 and 13 are defined in the appendix. The probability distribution of the inter-event time is then given by Ji PWD'(t + Prob(IET' _> n) =
n)f~(t)dt
i'l
]0 PWD'(t +
1)f~(t)dt
where f~ is the density function the distribution of the time t(0 < t < 1) from the beginning of the dry period to the first clock hour; 1 - t is
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67 95
73
therefore the time from the last clock hour in the storm to the end of the storm. This distribution is complex: if there were only one cell alive at the beginning of the last hour, then for a given storm, this time would have the distribution 77e-,~(1-t) -1 ~ - g
fT,,,(t) -
However, there can be more than one cell and new cells can be born; we shall therefore use as approximation for this expression: Prob(IET' > n) ~
Ii PWD'(t + n) dt PWD'(1)
The mean inter-event time m~ is obtained by summing Prob(IET'_> n) from n=lto~: m~ = Z P r ° b ( I E T ' - > n) n=l or
,
J' PWD'(t + n)dt
md ~ n=l 0 PWDI(1)
(7)
and similarly s 'a ~ ~/ 1 + Z~( 2 n -
1) I1 PWD,(t + n) d t - m~t2 0 PWD'(1)
n=2
(8)
Number o f events and average event duration
As for the non-randomised model, these may be derived from the expression for m~: ,
mn =
p(1)' x 24 × N(i) ,
(%)
md
where N(i) is the number of days in month i and ,
1 -p(1)'
,
mw-
p(1)'
× md
(lOa)
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
74
It must, however, be noted that these two variables may also be derived without using m~, as follows: !
m, = [ p ( l ) ' - p ( 2 ) ' ] x 24 ×
NM
(9b)
and mtw = [ 1 - p ( 1 ) ' ] / [ p ( 1 ) ' - p ( 2 ) ' ]
(10b)
Formulae (9b) and (10b) will be preferred to the others when the probability PSTOV of overlap of storms becomes too large: PSTOV _> 20%. Parameters and results for the randomised model
Parameter estimation As for the Bartlett-Lewis model, parameter estimation consists of the identification of certain variables characteristic of the precipitation process, for which Eqs. (2)-(5) above are solved using a modified version of the Powell hybrid method. Several alternative sets of variables could be selected. They must include features both of the depth process and of the proportion of dry periods; thus, two sets of variables were examined: (1) the R1 set
{E[Yi(U], var[Yi(1)],cov[Y~ ') ~v(1)],p(1),cov[Yi (6), yi(6+l],p(24)} i+1 (2) the R2 set
{E[Yi (l)1, var[Yi (1)], cov[Y} 1), V~i+,],p(1), (1) var[Yi(6)],p(24)} The two methods were applied to individual months of the data record and yield similar values of A and #x, but the other parameters are very different according to the method. It is moreover easy to see that parameters such as u and ~ are not well identified, since for very different values of u and c~, it is possible to obtain parameter sets yielding characteristic values within 5% of the historical values. This problem of identification is examined below. For the m o n t h of July, no parameters were obtained for the R1 set because the model was incapable of reproducing both the relatively high value of p(24) and the relatively low value of the six-hourly lag-1 autocorrelation A(1). Parameters for the R2 set are shown in Table 1.
Analytical results for R2 variables The second-order properties of the rainfall depth are well reproduced by the
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67 95
w i
o
y_ .=_
0
0
i
e. o
75
76
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67--95 0.50
--
0.40
--
6-hourly
-- - - -
A(1)
analytical values
BLRRPM Historical
data
- ~ ' ~ -~ / ~ ~
r -
/~--
\\//
0.30--
0.20
--
0.10---
I
I
I
2
3
,~
[ ~
l
6
I
7
I
s
I
9
I
~e
I
~
[ ~2
Months
Fig. 1. Analytical autocorrelations of lag-1 of six-hourly rainfall for the randomised model.
Bartlett Lewis Randomised Rectangular Pulse Model (BLRRPM) (Figs. 1 and 2), with a tendency to overestimate some of the autocorrelations at timescales larger than 12 h. The chief improvement brought about by the randomisation is an adequate representation of the proportion dry even at those levels of aggregation not used in the fitting (Fig. 3(a),(b)). The inter-event time properties are somewhat improved on the two methods used for the BLRPM (Fig. 4), which were based on the identification of the sets {E[y~I)],
var[Yi(l)],c°v[Yi 0), Yi~l], var[yi(6)],c°v[Yi (6) , ~i+lJJY (6)]~
and {E[ ri(l)], var[ yi(l)], cov[ yi(1) yfill] 'cov[ yi(6) yi(6+l],cov[ ri (12>, yi(l~)] } denoted as H = {1;6} and H = {1;6; 12} respectively. Nevertheless, clear discrepancies remain in the reproduction of rod, mw and ran.
77
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
~.se-
24-hourly
- -- -
0.40--
B L R R P M analytical Historical data
A(1 )
values
0.30--
.¢ 0.20
/
--
0.10--
0.~
I 2
r 3
I
f
l
I
1
I
4
~
6
7
e
9
I le
I ~1
1 1~
Months Fig. 2. A n a l y t i c a l a u t o c o r r e l a t i o n s o f lag-1 o f 2 4 - h o u r l y rainfall f o r the r a n d o m i s e d m o d e l .
Simulated results for R2 variables These results were confirmed by simulation of 144 years of hourly data. The discrepancies between analytical and simulated results do not exceed 10% apart from Im~w- mwl and ]mln- m,I for which the discrepancies sometimes reach 15% because of the approximations involved in their analytical expressions. A small discrepancy (smaller than 3% for p(1) and smaller than 10% for p(24)) remains between the analytical and simulated values of the proportions dry. This is analysed in Onof (1992).
Parameter optimisation Parameter identification with simulated results As was observed above, only A and #x are properly identified by the R1 or R2 methods. This identification problem was also observed by RodriguezIturbe et al. (1988). Moreover, the features characteristic of the event distribution are not adequately reproduced. Therefore, as in the case of the
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67 95
78
6-hourly proportion dry BLRRPM
++I //
0,
iI
I III /
75
~k"l, xx
/
;
]II/
x
1~i%x ;I
]
x,, ]
.... - -
0.70 I
R2 method mlto~l.cL1,
.....
data
I
I
I
1
[
I
I
I
I
I
2
3
B
6
7
8
9
10
11
12
Months
(a)
tg-hourly proportion dry
0.7B --
BLRRPM
0.70
/
I
x
l
x
tillI
tlttl
0.60
.... - -
0.m5 (b)
I~ 'method llimt,o r i c a l data
I
I
I
I
I
I
I
I
t
I
I
2
3
~1
B
6
7
8
9
10
11
12
Months
Fig. 3. (a) Analytical six-hourly proportion dry for the randomiscd modcl; (b) analytical 12-hourly proportion dry for the randomiscd modcl
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
Mean Inter-Event
79
Time
BLRPM/BLRRPM m o d e l s \
/\ /
ES--
\
/
\
f 60-k.
\
/
\
\ / \/
/ /
.o
\
/ \
/
m
/
\ \ k~
/ f
~
/
....
/
-
J
analyt/cal results Historical data B L E P M H = ~ I ; 6 ~ analytical results ~
--
xxx~ ~
25-
2e 1
2
3
4
~
6
7
8
9
10
11
12
Months Fig. 4. Analytical mean inter-event time for the non-randomised and the randomised models.
BLRPM (Onof and Wheater, 1993) which also gave rise to a parameter identification problem, an optimisation was carried out with the following objective function:
/[md(H ) -- md(i)] 2
[mn(H) - mn(i)] 2
where md(H), m,(H) and rnd(i), m,(i) are respectively the historical and simulated (iteration i) values of the two variables md and m, (mw is not used since it is a function of the other two). To investigate and improve parameter identifiability, a poorly identified parameter was selected and held constant while the five remaining parameters were optimised as before (by solving a set of five equations). The value of di was then computed; this procedure was repeated for different values of the fixed parameter, until an optimal value of the fixed parameter had been found. The five variables selected for the parameter estimation and the parameter which is to be fixed remain to be chosen. The following investigations were first conducted on the month of April.
C. Onof, H.S. Wheater/Journal of Hydrology 149 (1993) 67 95
80
Using the set of variables
{E[Yi(1)],c°v[Yi (1) , Y(1)I,p(1),c°v[Yi k
, Ji+l],p(24)
and fixing a or u, no optimum was obtained. However, if instead of the autocovariances, the autocorrelations are used, i.e if the set
{ E[Yi(I)],c°v[Yi(')'YiO+l]{S[Yi(1)]} 2
,p(1), c°v[Yi(6)'yi(6+l]}~,__ ~)]~5 ,p(24)
is used, then: (1) if a is fixed, the convergence of the Powell hybrid method is difficult to obtain; (2) if u is fixed, then we obtain an optimal value//opt (as verified by simulation of 200 years of data) which corresponds approximately to a minimum for the number of cells per storms #c. This corresponds to a small value of E(3) as had been found for the BLRPM optimisation but unlike the BLRPM case, EQ/) also reaches a minimum in the neighbourhood of the optimum. Applying this procedure to all months and choosing for initialisation values of the parameter sets those obtained by an examination of the parameters from grouped months for method R2, we obtain the parameters of Table 2. (The grouped month parameters correspond to larger quantities of data (three or four months) and although of limited significance, are useful in providing a value with which to initialise iterations for any given month in that group.) Table 2 Optimal parameters for the randomised model Month
A (h I)
#x (mm h-l)
a
u (h)
n
0
#~
January February March April May June July August September October November December
0.0261 0.0221 0.0212 0.0219 0.0226 0.0194 0.0111 0.0231 0.0215 0.0206 0.0282 0.0317
1.0329 1.1851 1.2122 1.4182 2.3427 2.0769 4.1869 3.3444 1.5922 1.2486 1.1771 0.7938
2.5032 2.6175 2.1800 2.3390 3.0652 3.3702 2.2894 2.7016 6.6513 4.0159 5.1272 3.6535
1.0600 0.9000 0.6000 0.4700 0.5000 0.9000 0.2476 0.4331 3.8000 2.0000 3.6000 1.1328
0.2069 0.2431 0.1662 0.2496 0.1816 0.1845 0.0421 0.1593 0.2476 0.2728 0.2584 0.8483
0.0669 0.0639 0.0386 0.0518 0.0366 0.0467 0.0064 0.0408 0.0986 0.0834 0.1188 0.1297
4.0941 4.8076 5.3047 5.8216 5.9638 4.9512 7.5515 4.9009 3.5099 4.2711 3.1740 7.5408
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
81
Analytical determination of the optima The same optimisation process was carried out using the analytical values of the characteristic variables, i.e. with the objective function
/[md(H ) - m~t(i)]2 di : v
md ~I--I7
[mn(H) - mrn(i)]2 -~
mn ( H ) 2
(12)
where ma(t), , . mn(~) , . are the analytical (iteration i) values of the two variables / ! md, mn. The optimal values of u which are obtained are very close to those obtained with the simulated results above: //opt ~'~ /Jopta
is the analytical optimum. This provides us with an analytical method for the approximated determination of Vop t. Figure 5(a)-(c) (where RO is the optimised version of the randomised model) shows the improvement brought about by this optimisation in the reproduction of the variables md, mw and mn. Moreover, the optimal distribution of the simulated inter-event time population was compared with its historical counterpart using a Student test for the average, a Fischer test for the standard deviations and a M a n n - W h i t n e y test, and showed no significant discrepancy for all months except July and December. For these two months, the above method did not provide an optimum. Table 3 shows the analytical values of the characteristic variables for the parameters in Table 2. where/Jopta
Extreme values of the model
Hourly data For most months of the year, the extreme-values of the hourly rainfall distribution are overestimated by the R1 parameters for all return periods and by the R2 parameters for return periods longer than the data set (38.5 years). Figure 6(a),(b) shows the extreme value distributions for the historical and the simulated data for the months of January and June. The results for the optimised parameters are similar to these.
Daily data For daily rainfall, for all months apart from January and November,
82
C.
Onof, H.S. Wheaterl Journal of Hydrology 149 (1993) 67 95
Average E v e n t D u r a t i o n
4--
," - - . -
BLRRPM
/ ~xxx
,
/
i#
/
3--
It
i
/ // \',
/
/ I
2--
\,
I
I
I
', I
,'
v_---" .... - -
-
'
I
I
2
I
3
I
4
I
5
I
6
(a)
7
I
--
R2 s/mulated H/storlcal data RO simulated
I
8
I
9
1
10
I
11
12
Months
Mean I n t e r - E v e n t
40-
Time
BLRRPM
Ixxx\ !
x
; ii m
',
I ," l
..--. x s I
\
~
",
/
\
I
t
It
x\
A
/
ix ii
~
(b)
I \
I 2
~
\/
I
I 3
~
/
I 4
I El
# / / z"I
'~
,
/
\\1
'~ \ 1
I 6
I. 7
I IB
I ?
I le
I 11
I 12
Months
Fig. 5. (a) S i m u ] a t e d o p t i m a l a v e r a g e e v e n t d u r a t i o n ; (b) s i m u l a t e d o p t i m a l m e a n i n t e r - e v e n t time; (c) s i m u l a t e d o p t i m a l n u m b e r o f events.
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
Number
83
of E v e n t s
BLRRPM
2K-v
x
¢
20--
.... - -
1B
[ 2
(c)
1~ ~ul.at, lllrl~rlcs1
RO simulated
--
ed
data
I
I
I
I
I
I
I
I
1
3
4
B
6
7
5
~, ?
10
11
12
Months
F i g . 5. C o n t i n u e d .
method R2 reproduces the extreme values of short return period but diverges from the historical data for larger return periods; method R1 often overestimates the historical values for all return periods (Fig. 6(c),(d)). Again, the results for the optimised parameters show no improvement.
Daily representation of the randomised parameters From the analysis of monthly data, twelve sets of parameters have been obtained for each of the estimation methods. Two problems arise: (1) is it possible to provide a smoother representation (e.g. polynomial) of the variation of the parameters over the year, eliminating the arbitrary monthly grouping of data?; (2) can this variation be represented in such a form (e.g. Fourier analysis) that the coefficients of this representation might be used for a regionalisation of the parameters?
84
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
H o u r l y E x t r e m e Values
12--
J~uauary; BLRRPM
11--
//
10--
//
/
/
/
•1
9--
/• B--
/
// / 6--
/
/
i"
/•
/
/
/
/ /
/~ F { / /
/ J // /~
5~ 4~ a-
¢ S . /
....
R2
SI----
I
-1.0
I
0.0
I
1.0
I
2.0
Return
period
1 2
(yrs)
I
3,0
GumbeI
reduced 1
1 10
5
]
I
4.0
5.0
6.0
variate
1 ZO
1
I
100
200
(a)
~-
H o u r l y E x t r e m e Values June;
BLRRPM
2~-
//
/"
17
18-
i 16 -
/
-
F¢
14-
./ /
~f 1 2 - /
10 -8--
,14
/
./
-/
/-
/
~
11 / FI
/~
i
/I
/
/
/
6--
2 0 -1.0
I
l
0.0
1.0
I
Return
period
(yrs)
[ 2
I
I
reduced
variate
2.0
GumbeI i 5
3.0
[ S.O
4.0
I
r
10
20
I 6+0
I
T
I~
200
(b) F i g . 6. (a) H o u r l y e x t r e m e v a l u e s f o r J a n u a r y ; (b) h o u r l y e x t r e m e v a l u e s f o r J u n e ; (c) d a i l y e x t r e m e v a l u e s f o r J a n u a r y ; (d) d a i l y e x t r e m e v a l u e s f o r J u n e .
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67 95
~0-
Daily Extreme January;
Values
BLRRPM 7
6o-
7
50-
f
j"
7
I~1
.f
~0--
J" j-"
./" 7
/ I SO--
7
J" "''"
~
lO-
o
I 0.0
-1.0
I 1.0
I 8.0
I 3.0
Cumbel
Return
period
(yrs)
I 4.0
reduced
I 5.0
1 6.0
variate
1
1
1
1
2
5
10
20
t
l
100
200
(c)
Daily Extreme
100-
June;
Values
BLRRPM /-
90-
J J J
00-
J J 70-
J f J
i
J
Et
J J j j
i ~ ~"
10-0
I 0.0
-I .0
Return
period
(d) Fig. 6. Continued.
I 1.0
(yrs)
I 2
I I I 2.0 3.0 4.0 GumbeI reduced variate I 5
I 10
I 20
i 5,o
T
1~
i 6.0
T
200
85
86
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67--95
~
.
.
.
.
~
.
.
.
.
C. Onof, H, S. Wheater Journal of Hydrology 149 (1993) 67-95
¢,.q
t"q
~D
H
[..,
O
[-.,
az
~
87
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
88
1.0
I 31
I 61
I 91
I 121
I 1 I IE;I 181 211 Days
I 241
I 271
I 301
l I 331:361
I 241
I 271
I 301
I 331
(a)
BLRRPM daily
0.7 0. a0. ¢L,
0. 0. 0. 0.1
0"0
I 1
I 31
1 61
I 91
I 121
1 151
I I 181 211 Days
I 361
(b) F i g 7 (a) D a i l y variation o f A for the randomised mode] (R2); (b) daily variation o f e~for the randomised mode] (R2); (c) daily variation o f ~ for randomised mode] (R2).
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
89
1:R.O-12.0-
~'~
BLRRPM daily
11 . 0 -
10.09.@o,m
8.0-
~
7.e6.05.04°0:3.02.01,0
I 31
I 61
I 91
I 121
(c)
I 151
I I 181 211 Days
I 241
[ 271
I 3~1
I 331
I 361
Fig. 7, Continued.
Polynomial curve fitting F o r each parameter ah[h E ( 1 , . . . , 6 ) ] , polynomial coefficients for an expression of the parameter in a Chebyshev series form are calculated: n=k
~'~h,k = ½ao(ah)To(t) + Z an(ah)Tn(t) n=l
where
T~(t) = cos(n arcos t) = ½[(t + iv/1 - t2)" + ( t - ix,/1 - t2) "] The approximation of a h by Y]h,k is carried out by a least-squares method. On the whole, for parameters A, n and 05, a representation of degree 5 is sufficient (Fig. 7a-c), whereas for u, a polynomial of degree 7 appears necessary. F o r the other two parameters, #x and a, it is preferable to choose a degree k = 12. The parameters thus s m o o t h e d produce very satisfactory results for simulations of the m o n t h l y characteristic variables. They have the additional advantage of yielding more realistic results for periods which are not clusters of the calendar months.
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67 95
90
Table 4 Fourier coefficient moduli for the randomised model, N = 12 (R2)
Degree
A (h-1)
#x (mmh-l)
a
u (h)
~
0 1 2 3 4 5 6
0.0803 0.0126 0.0064 0.0027 0.0040 0.0018 0.0028
5.9714 2.4249 0.9357 0.3420 0.5920 0.5912 0.4482
10.2632 0.8729 0.3700 0.0489 0.5950 0.5148 0.4496
2.4343 0.6602 0.1886 0.0950 0.4232 0.4818 0.0726
1.6591 0.6273 0.0909 0.0379 0.1700 0.1314 0.0674
0.2412 0.0885 0.0232 0.0074 0.0079 0.0415 0.0151
Fourier analysis of the continuous-time representation A Fourier analysis of the parameters may also be carried out to isolate the predominant periodicities of the monthly parameter values. The Fourier coefficients for parameter crh are calculated in the complex form: l
N-I
{
.27rjk'~
F~h,k -- vT~ j~=oCrh( j ) exp ~-' --~ -) Choosing N = 12 enables us to represent the monthly components which dominate the picture (the coefficients decrease w i t h j for A and #x) and longer trends (two to three monthly), as can be seen from Table 4. This provides us with a representation of the parameters which could be used for a spatial generalisation of the model; the coefficients of lower frequency are expected to vary less than those of higher frequency and thus represent underlying seasonal trends which may not be as location specific as variations of higher frequency. Conclusion
The extension of the Bartlett-Lewis Rectangular Pulse Model (BLRPM) to include a random cell duration is shown to improve the simulation of British rainfall data. Good reproduction of dry periods is achieved at all time-scales (hourly to daily), in addition to the previous reproduction of daily rainfall depth distributions. This corroborates results obtained by Rodriguez-Iturbe et al. (1988) for data from Denver, CO, USA. Problems of parameter identification are apparent, and a two-stage optimisation is proposed which allows for optimal reproduction of features of the
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
91
event distribution and which can be carried out analytically. Representation of the parameters as daily variables avoids the artificiality of the monthly subdivision of the year and allows for the possible generalisation of parameters through regional analysis. Two problems remain which require further examination: (1) the excessive values of some of the autocorrelations for lags greater than 12 h, in particular for the optimised parameters; (2) the inadequate reproduction of extreme values for return periods greater than the length of the data set. These problems can be addressed through an extension of this model, e.g. by the introduction of a jitter and by using a gamma distribution for the cell depths; although additional parameters are expected to give rise to further problems of identifiability. These issues will be examined in a subsequent paper.
Acknowledgments The helpful advice and comments from Sir David Cox and Valerie Isham during the progress of this work are gratefully acknowledged.
References Cowpertwait, P.S.P., 1991. Further developments of the Neyman-Scott clustered point process for modelling rainfall. Water Resour. Res., 27(7): 1431-1438. Entekhabi, D., Rodriguez-Iturbe, I. and Eagleson, P.S., 1989. Probabilistic representation of the temporal rainfall process by a modified Neyman-Scott Rectangular Pulses Model: parameter estimation and validation. Water Resour. Res., 25(2): 295-302. Onof, C., 1992. Stochastic modelling of British rainfall using Poisson processes. Ph.D. Thesis, Imperial College, London. Onof, C. and Wheater, H.S., 1993. Optimisation of the Bartlett-Lewis Rectangular Pulse Model for British Rainfall (submitted). Pattison, A,, 1965. Synthesis of hourly rainfall data. Water Resour. Res., 1(4): 489-498. Rodriguez-Iturbe, I., Cox, D.R. and Isham, V, 1987a. Some models for rainfall based on stochastic point processes. Proc. R. Soc. London, Ser. A, 410: 269-288. Rodriguez-lturbe, I, Febres de Power, B. and Valdes, J.B., 1987b. Rectangular pulses point process models for rainfall: analysis of empirical data. J. Geophys. Res., 92D8: 96459656. Rodriguez-lturbe, I., Cox, D.R. and Isham, V., 1988. A point process model for rainfall: further developments. Proc. R. Soc. London, Ser. A, 417: 283-298. Smith, J.A., 1981. Point process models of rainfall. Ph.D. Thesis, The Johns Hopkins University, Baltimore, MD.
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67-95
92
Appendix To calculate PWD(h) ', let us first suppose there is only one storm alive at t = O:
One storm alive at t = 0 If we suppose this, then if P1 is the probability that an interval o f length h is dry conditional upon the first instant wet P1 = Prob [(0, h) dryIt = 0 is the last instant of a wet period]
E
c
I C e -'~h × e -~oh + E ~-s [ C e -'~h /
where C is the number o f cells in the storm and C' the number of those whose cell end times are followed by a dry period. This obtains since t = 0 m a y be any cell end time within the storm with probability [E(C'[C) - 1]/E(C'[C) or the last cell end time with probability 1/E(C~[ C)', where C' is the number of cell end times which are followed by a dry period. To calculate E [1/C'[C] as a function o f C, the number o f cells in the storm, we need to compute the probability P S U P that two consecutive cells overlap. :3C
P S U P = / nr/e -~u e -~u e -*~u du d
0
because for each time u, we multiply the probabilities that the end of the cell is not reached, that the cell is still alive and that a cell arrives at u. We obtain P S U P = n / ( n + ~ + 1). The probability that a cell end time is dry is approximately that o f two consecutive cells not overlapping if we neglect the probability PSUP3 of a first cell overlapping with the two others arriving later which are mutually non-overlapping. This is 0(3
X
X--U
X
PSUP3=Irle-'TXe-~'TXdxlnrle-'~'TUdulrle-'WdvlmTe-'~'mdw 0
0
0
u+v
which computes to ~+ 1 PSUP3 =
2(1 - ~;)(t~ + 1)(2~ + 1 + ~b)
1 ~+1+05
+-
2 (1 - t~)(~ + 1)(t~ + 2 + ~b)
1 2(~ + 1)(1 + ~h)
For the parameter values obtained, we have PSUP3 < 0.05.
C. OnoJ; H. S. Wheater/ Journal of Hydrology 149 (1993) 67 95
93
Therefore, the distribution of C' conditional upon C is binomial:
P(C' -- C - q l C ) where Cab =
PSUp)c-I-qPSUP q
= cc~-~-q(1 -
a!/[b!(a - b)!] and we therefore obtain
E(1"-1TI C) ~}--~CPc-1 (1 - PSUP)PPSUPP p=0
p+l
This yields the probability P1 for one particular storm with C cells; to obtain the probability P2 that an interval of length h is 'dry' within one storm knowing that it starts at the end of a wet spell, we must decondition PI upon C and take expectations over r]: E ( )~-5 l
"~
ZZCcLl°CC-'
(I-PSUP)PPSUPP ( I3 ~"7
c=lp=0 P 2 = E(P1) ~ E,j
and this yields p2~ [I_E(1)]
p+ 1 1- E
\ c )1
u~
• ~(~ + u)
\ ~ J
e -ah x e -~vh +
e -~h +
E(1).e-,h
where E ( 1 / C ' ) has to be evaluated numerically.
Multi-storm case P2 above is the sum of two terms: P21 =
1-
. ( u +/Ju ) ~ e
-Ah
and
If there are k storms active at t = 0, then conditional upon this P21 must be multiplied by [PNCB(h)] k-l, where PNCB(h) is the probability that no cells are born during [0, h] for any o f the storms and similarity P22 must be multiplied by PNCB(h) k. Removing the condition, we obtain PWD'(h) ~ P21 x A21 + P22 x A22
C. Onof, H.S. Wheater/ Journal of Hydrology 149 (1993) 67--95
94
where ~-, Prob(exactly k storms alive) PNCB(h) k-1 A21 = ~ ~ e or more storms alive) k=l
and O<3
A22 = Z
Prob(exactly k storms alive)PNCB(h) k
k=0
Remark: this equation for PWD~(h) supposes that the distribution of C' is the same as that for the single-storm case. This approximation is acceptable if PSTOV (defined below) is not too large, e.g. if PSTOV_< 20%. Then I k = Prob(number of storms alive = k) = E~ e -A/~
k!
or
1 Ik -/~!r(~) e a/~tk~t/ t ~
lu~e-~t at
0
Practically, the sums in A21 and A22 will be approximated by sums up to k = 3; indeed, the probability of more than three storms superposed at any given instant is 3
PSTOV3~ [E~(iAe-Ate-~Vtdt)] or 3
PSTOV 3 ~
~
dt
which can be evaluated numerically. PSTOV 3 _< 0.01 for the parameter sets obtained below. Thus
,~ ~ IkPNCB(h)k-1 A21
~
I1 + / 2 + / 3 3
A22 ~ Z k=0
IkPNCB(h)~'
C. Onof, H. S. Wheater/ Journal of Hydrology 149 (1993) 67-95
95
where PNCB(h) is the probability that no cells occur in a given storm during the time interval [0, h]. The calculation of the probability PNCB(h) now follows. The calculation is as in the BLRPM case, but the integral which is obtained must be averaged over all possible values of r/: the storm may end at u (0 _< u < h) and no cells occur during [0,u] (with probability O~le-°~Ue-~'lu),or may end after t = h and no cells occur during [0,hi (with probability e e,lh e ~h). Thus PNCB(h) = E+l
Or/e -4'~u e '+~' du + e -0~h e -'+'Th
We obtain PNCB(h) = E,, ~ O
+ ~ n
e_(~+0)@]]
which may be simplified by using the expression of the Laplace transform of the gamma distribution: PNCB(h)_
0
+
~
[
u
]~
We obtain finally [ PWD'(h)~
(~5)] 1-E
u'~
3IkPNCB(h)k' Z I1+I2+I~
-.Xh
(u+u) ~e
k=l
+ E(l e h , PNCS/h/