Annual Journal of Hydraulic Engineering, JSCE, Vol.52, 2008, February

POT NORMALIZED VARIANCE PARAMETERSEARCH OF THE TEMPORAL NEYMAN-SCOTTRAINFALL MODEL
Carlo Arturo MONDONEDO 1, Yasuto TACHIKAWA2, and Kaoru TAKARA3
1Student Member, Graduate student, Dept. of Urban and Environmental Eng., Kyoto University(Kyoto 615-8540, Japan) E-mail:mondonedo@flood.dpri.kyoto-u.ac.jp
2Member of JSCE, Dr. Eng., Associate Professor, Dept. of Urban and Environmental Eng., Kyoto University(Kyoto 615-8540, Japan)
3Fellow of JSCE, Dr. Eng., Professor, DPRI, Kyoto University (Gokasho, Uji 611-0011, Japan)

In previous investigations on the Neyman-Scott clustered Poisson rectangular pulse rainfall model(NSM), model parameters were often estimated using rainfall moments, covariances, correlations, and dry
moments while model performance was often assessed through the accuracy of extreme values. No
former contribution included extreme value information in improving the parameter search. This was
the motivation of the authors in finding a link between the scaling apparent in the historical rainfall
maxima and the NSM parameters. The normalized variance or Fano factor of the Peaks Over Threshold(POT) rainfall process was adopted to define this link. In general, the results obtained from thisapproach showed acceptable performance with minor limitations.
Key Words: point processes, Neyman-Scott rainfall model, Fano factor, fractals

as Telesca et al.8) demonstrated scaling in rainfall
occurrence vs. counting time in historical data.
This study focuses on linking the scaling of
historical extremes to the NSM parameters.
Historical rainfall from the AutomatedMeteorological
Data
Acquisition
System(AMeDAS) station at Kamishiiba in thesouthwestern island of Kyushu, Japan, wasgathered to include mechanisms such as fronts,
typhoons, etc. A brief description of the NSM
appears in Section 2 followed by a discussion on
the Peaks Over Threshold (POT) rainfall Fanofactor and its eventual inclusion in the NSM in
Section 3. Model application and verification
appear in Section 4 followed by concluding
remarks in Section 5.

1. INTRODUCTION
Stochastic rainfall modeling under point
processes1) is an essential tool used byhydrologists for numerous applications.
Athorough review of this subject was presented byWaymire and Gupta2, 3, 4). A variant of thisapproach, the Neyman-Scott clustered Poisson
point process, NSM here for brevity, was further
developed by Rodriguez-Iturbe et al. 5)
For flood control applications, a synthetic
rainfall record generated from a point process
such as the NSM should replicate the historicalextreme values. It would thus be advantageous
to consider the rainfall extreme values in
estimating the parameters of the NSM.
However, previous investigations did not employ
this information in the parameter search, adopting
mainly a method of moments approach which
incorporated sample moments and correlations
instead5). A means to include the maximal
information in the method of moments parametersearch is proposed here.
Often, the scaling of statistics in point
processes leads mathematically to power-law
dependencies. This tendency for scaling has
been a key focus in the investigations of fractal
stochastic point process6, 7). In fact, studies such

2. MODEL DESCRIPTIONIn the NSM, rainfall is formulated as atemporal cluster of rain cells which occur basedon a Poisson process with occurrence rate λ.
Cluster size is based on a geometric distribution
with mean μc. Rain cells lag cluster occurrences
randomly following an exponential distribution
with mean 1/β. Rain cell durations follow an
exponential distribution as well with mean 1/η.
Rain cell intensities are modeled by a gamma- 97 -

A1 (h ) = ηh − 1 + e −ηh
B1 (h ) = βh − 1 + e − βh()
B (h, k ) = 0.5(1 − e ) e ( )A2 (h, k ) = 0.5 1 − e −ηh e −ηh (k −1)
2

Synthetic
RainfallRainfallintensity− βh 2− βh k −12where:
n = time interval counterh = integer specifying time step interval ofdata (1 for 1 hour, 24 for 1 day, etc.)h
Yn = rainfall depth in the n-th time of intervalh
E Yn(h )= mean rainfall depth record at

Time
Cell
Intensity

Cell
Duration

Time

Cell
ArrivalTime

var Yn(h )(h )2

3. THE POT RAINFALL AND NSM
FANO FACTOR EXPONENT

Time

3

(1) POT Rainfall count process Ni(t)Figure 2 shows the construction of the special
count process Ni(t) to be used in determining theFano factor FF(T) at an arbitrary window T. With
rainfall pooled for each month i (i=1: January,i=2: February, etc.), we define the process Ui(t)(Fig. 2a), as the rainfall magnitude per fixed
duration t (i.e.: t=1 hour in this study). For the
M months in each Ui(t), we determine the M
monthly or block maxima Zi(t) and thresholdvalue zi, the minimum of Zi(t). We define the
POT rainfall point process Qi(t) based on Eq. (4)(Fig. 2b).

Cluster
ArrivalTime

Fig. 1. Schematic drawing of the Neyman-Scott Rainfall Model.

distribution with scale parameter θ and shapeparameter α. The resulting rainfall is thesuperposition of these processes as shown in Fig.
1.
The moment and correlation expressions forthe aggregated rainfall point process in terms of
the NSM parameters were derived byRodriguez-Iturbe et al.9)
The succeeding
Equations (1) – (3) list these expressions in theaggregated form, specific to the distribution
parameters mentioned: β, λ, μc, η, α, and θ.

= λμ c (αθ )h / η()[β A (h) − η B (h)] +3314λμ c (αθ ) A1 (h )2

cor Yn(h ) , Yn(+hk) =

(2)

η3

1

[β A (h, k ) − η B (h, k )]
3

3

2

2

)

(3)

⎧1 ; Qi (t ) > 0
Bi (t ) = ⎨⎩0 ; Qi (t ) = 0(5)(2) Fano factor and governing relationships
The Fano factor FF(T) is defined as the ratio of
variance of count Ni(T) and mean of Ni(T), or:

4λμ c (αθ ) A2 (h, k )
2

+

(4)

Using adjacent windows of arbitrary value T,
we define the count process Ni(T), the sum of Bi(t)within each segment T (Fig. 2d).

2λμ cαθ 2 (1 + α )(μ c − 1)×var Yi (h ) βη 3 β 2 − η 2(⎧U (t ) ; U i (t ) ≥ zi
Qi (t ) = ⎨ i⎩ 0 ; U i (t ) < zi

For each rainfall occurrence in Qi(t), we define
the binary process Bi(t) based on Eq.(5) (Fig. 2c).
This process counts all durations of rainfall
greater than or equal to the threshold zi and serves
as the basis for defining the Fano factor used here:(1)

2λμ cαθ 2 (1 + α )(μ c − 1)
var Yn(h ) =×βη 3 β 2 − η 2= autocorrelation of rainfall
record at h-hours at lag k

Time

3

E Yn

h-hours(h )cor Yn , Yn+k

1

Clustersize (cells)(h )h-hours= variance of rainfall record at

η3- 98 -

(b)(a)
SampleUi(t)Sample

Qi(t)

19

25

rainfall (mm)rainfall (mm)
30

30

25

25

20

20

15

15
10

10
z

i

5
0 1

7

13

19

25

5
31

37

43

49

1

0

7

13

time (hr)

(c)

37

43

49

37

43

49

(d)
Sample Ni(T=20 hours)

Sample Bi(t)Ni(T1)=11

1

31

time (hr)

7

13

19

25

31

37

43

49

1

7

13

Ni(T2)=2

19

25

31

time (hr)time (hr)Fig. 2. Determining the Peaks Over Threshold (POT) rainfall point process and count process Ni(t). (a) sample rainfall
Ui(t) with previously determined threshold zi, (b) POT rainfall Ui(t) values greater than or equal to zi (c) unit counts
Bi(t) assigned for each rainfall occurrence, (d) count process Ni(t).

E [N i (T )] − E N i (T )2

FFi (T ) =are used to explicitly determine FFHi(TMi), theFano factor of the historical Qi(t) at window T =mean storm duration TMi such that:2

E N i (T )(6)where:
E = operation to obtain expected value.

FFHi (TMi ) ≈⎛T ⎞FFSi (TMi ) = 1 + ⎜⎜ Mi ⎟⎟⎝ T0 ⎠

ξin which T0 is the basic data duration (i.e.: T0 = 1
hour here). Strictly speaking, the rainfall data
used by Telesca et al. 8) were pooled yearly (zi = 0)
instead of Qi(t). It was assumed here that Eq. (7)is valid for Qi(t) throughout the set of windows Tœ Ω = {2, 10, 20, …, 100 hours}, a conservative
range for possible mean storm durations (Fig. 3b).
Alternatively, within Ω, the approximations forvariance and mean of Ni(T) are proposed here as
Eqs. (8) and (9) such that the Fano factor can be
independently estimated as Eq. (10).
(8)

E N i (T ) = CiT(9)AFFi (T ) ≈ i T Bi −1
Ci

ξ Si

(12)

For our purposes, an ideal simulation should
yield synthetic rainfall with Qi(t) such that
historical Fano factor FFHi(TMi) and synthetic
Fano factor FFSi(TMi) are equal, or based on Eqs.
(10) and (12):(7)

Var N i (T ) = AiT Bi(11)

The relationship for the synthetic equivalent
FFSi(TMi) can be written explicitly using Eq. (7)such that:

Based on data from independent studies,
Lowen and Teich6,7) and Telesca et al.8) proposed
the power law relationship of Eq. (7) to describe
the scaling that occurs over several decadal valuesof T (Fig. 3):⎛T ⎞FF (T ) = 1 + ⎜⎜ ⎟⎟⎝ T0 ⎠

Ai
B −1
TMi i
Ci

⎛T ⎞
Ai
B −1
TMi i = 1 + ⎜⎜ Mi ⎟⎟Ci⎝ T0 ⎠ξ Si(13)Isolating the synthetic Fano factor exponent ξSileads to an expression relating properties of thehistorical Qi(t) to the unknown mean duration TMi.

(10)(3) NSM Fano factor exponent
Curve fitting operations for Eqs. (8), (9), andEq. (10) (i.e.: least squares fit for Ai, Bi, and Ci)- 99 -⎡A⎤( B −1)ln ⎢ i TMi i − 1 ⎥Ci⎦ξ Si = ⎣ln[TMi ] − ln[T0 ]

(14)

⎛
⎡
η ⎤⎞
TMi = 1 ⎜⎜ γ + ln ⎢(μ c − 1)
⎥ ⎟⎟
η
η
β
−
⎦⎠
⎣
⎝

(15)

Actual mean storm duration TMi was estimatedhere based on Cowpertwait’s expression10)
derived from the NSM parameters shown here as
Eq. (15).

(a)

(b)

⎛τ ⎞
FF (τ ) = 1 + ⎜⎜ ⎟⎟
⎝τ0 ⎠

ξ

⎛T⎞
FF(T) = 1 + ⎜⎜ ⎟⎟⎝ T0 ⎠

10

ξ

FF(T)

1

100

10
T (hr)

1

Fig. 3. (a) Scaling observed between count of spikes from a nerve fiber FF(τ) vs. counting time τ (reproduced fromLowen and Teich6)). (b) Scaling in the count process Ni(t) obtained from Kamishiiba POT series Qi(t) in June.

Table 1. Estimated NSM Parameters for Kamishiiba Region.

Jan
Feb
Mar
Apr
May
Jun
Jul
Aug
Sep
Oct
Nov
Dec

α

β (1/hr)η (1/hr)λ (1/hr)μcθ (mm/hr)

0.029520.0000
20.0000
20.0000
20.0000
20.0000
20.0000
20.0000
20.0000
20.0000
20.0000
20.0000

0.2096
0.27410.31410.35480.23020.19250.0771
0.06870.07620.11630.14320.18000.500060.0000
60.0000
60.0000
24.5687
60.0000
60.0000
26.6501
3.452258.2798
9.992460.0000

0.01550.01160.02050.02080.01130.02020.00630.00520.00510.00410.00530.004950.0000
50.0000
21.6258
28.1992
50.0000
37.3840
44.691950.0000
33.2218
38.174224.2467
31.0726

2.20660.73791.73631.35350.67742.9474
6.32612.49560.45493.17180.3746
1.4719TMi (hr)23.9182
16.3187
11.4900
10.9542
19.4561
21.689956.479065.0530
53.4124
36.071026.0940
22.1366

previously derived NSM Fano factor exponent ξSi
is proposed here. The objective function in theparameter search is given by Eq. (17).

By substituting Eq. (15) for TMi in Eq. (14), a
direct link between historical Qi(t) and the NSM
parameters is established as Eq. (16).
( Bi −1)
⎤
⎡ ⎧ ⎛
⎞⎫⎪
⎡
⎤
A
1
η
⎪
i
⎥
⎟
−
ln ⎢ ⎨ ⎜⎜ γ + ln ⎢(μ c − 1)1⎬⎥
⎢ Ci ⎪⎩η ⎝η − β ⎥⎦ ⎟⎠⎪⎭⎣⎦
ξ Si = ⎣
⎡1 ⎛
⎤
⎞
⎡
η ⎤
⎟⎥ − ln[T0 ]ln ⎢ ⎜⎜ γ + ln ⎢(μ c − 1)
− β ⎥⎦ ⎟⎠⎦⎥
η
η
⎣
⎣⎢ ⎝

6 ⎛
NSM m ⎞⎟
O = ∑ ⎜⎜1 −HIS m ⎟⎠
m =1 ⎝

θ=

E Yn

(1)

λμ cα

2

(17)

η
(18)

where:
NSMm = any of equations (1), (2), (3), and (16)
HISm = equivalent historical value of moment.
To determine β, λ, μc, η, α, and θ, six
expressions were required in the search: (1)hourly variance, (2) hourly correlation at lag 1, (3)12-hourly correlation at lag 1, (4) 24-hourly
variance, (5) 24-hourly correlation at lag 1, and(6) hourly Fano exponent ξSi. Eq. (18), the scaleparameter θ in terms of the mean hourly rainfall,
was used to limit the computational difficulty innumerically solving Eq. (17). Historical hourly
rainfall data was obtained from the AMeDAS
station in Kamishiiba, Kyushu, Japan from 1988
to 2002. Estimated parameters appear in Table
1, indicating TMi values within the range of Ω.

(16)

Moreover, with the historical exponent ξHiestimable through curve fitting Eq. (7) to thehistorical Qi(t), it is now possible to include ξSi,
the NSM Fano factor exponent, in the NSM
parameter estimation. Thus, in a limited manner
the proposed parameter search contains
information about historical extreme values thatshould lead to improved synthetic maxima.

4. MODEL APPLICATION(1) Parameter Estimation
A basic scheme consisting of the existing NSM
moments and correlation equations and the- 100 -

average historical monthly total rainfall and its
synthetic counterpart and the estimate of thestandard error. The values were observed to lie
within the 5% significance value of 2.048 (for the
sample size NHis + NSyn - 2 = 28), which indicates
acceptable model performance. The same is true
for other statistics such as those shown in Figs.
4b-c and 5. Monthly maxima of lower rank
were also within the 5% significance levels andwere omitted here.

(a)Monthly Total Rainfall
3
5% Level(His-Syn)/SE

2
1
0
-1
-2

5% Level

-3
1

2

3

4

5

6

7

8

9

10

11

12

Month

(a)

(b)

24-Hourly Variances
3

Hourly Variances(His-Syn)/SE

5% Level

2(His-Syn)/SE

5% Level

2

3

1
0

1
0
-1
-2

-1

-3

-2

1

2

3

4

5

6

-3
1

2

3

4

5

6

7

8

9

10

11

12

Month

9

10

11

12

24-Hourly Lag 1 Autocovariances
3

Hourly Lag 1 Autovariances

5% Level2

3(His-Syn)/SE

5% Level

2(His-Syn)/SE

8

(b)

(c)

1
0

1
0
-1
-2

-1

-3

-2

1

2

3

4

5

-3

6

7

8

9

10

11

12

Month
1

2

3

4

5

6

7

8

9

10

11

12

Month

(c)

(d)

24-Hourly Maxima
3

Hourly Maxima

5% Level2

3(His-Syn)/SE

5% Level

2(His-Syn)/SE

7

Month

1
0

1
0
-1
-2

-1

-3

-2

1

-3
1

2

3

4

5

6

7

8

9

10

11

2

3

4

5

6

7

8

9

10

11

12

Month

12

Fig. 5. Student’s T tests for the hourly rainfall point processfor (a) 24 hourly variances, (b) lag 1 autocovariances, and(c) 24-hourly maxima.

Month

Fig. 4. Student’s T tests for the hourly rainfall point processfor (a) monthly totals, (b) hourly variances, (c) lag 1
autocovariances, and (d) hourly maxima.

(3) Model Limitation
Figure 6 shows a sample Kolmogorov-Smirnov test (KS test) applied to the block
maxima of September. The test was formulated
such that 2 data sets that are more consistent inprobabilistic distribution receive a KS probability
approaching unity11). Based on Table 2, the(2) Model Performance
Fifteen synthetic records were generated permonth corresponding to the 15 monthly historical
records. Each student’s t ordinate tS in Fig. 4acorrespond to the ratio of the difference between- 101 -

synthetic block maxima were significantly
different from the historical counterparts forseveral months (i.e.: values below 0.95). During
these months, rain cell information may by more
relevant than storm cluster information.
Therefore, we propose replacing the 12-hourly lag
1 correlation with the 6-hourly lag 1 correlation to
improve the results for these months.

Table 2. KS Test results for Kamishiiba rainfall.
Block Maxima KS ProbabilityMonth
1-Hour24-Hour
1 (Jan)
0.83
0.83
2
0.98
0.98
3
0.83
0.55
4
0.98
0.28
5
0.99
0.99
6
0.98
0.99
7
0.98
0.99
8
0.55
0.98
9
0.99
0.83
10
0.99
0.99
11
0.98
0.99
12 (Dec)
0.98
0.99

cummulative frequency

September Kamishiiba Rainfall Block Maxima
1.000
0.800
0.600
0.400

Prof. Y. Tachikawa, Kyoto Univ.).

0.200
0.000
0.000

10.00020.00030.00040.00050.000

rainfall (mm)historicalsynthetic

REFERENCES

KS significance prob. = 0.99

1) Cox, D. and Isham, V.: Point processes, Chapman andHall, London, 1980.
2) Waymire, E. and Gupta, V.: The mathematical structureof rainfall representations 1: A review of stochastic
rainfall models, Water Resources Research, Vol. 17, No.
5., pp. 1261-1272, 1981a.
3) Waymire, E. and Gupta, V.: The mathematical structureof rainfall representations 2: A review of the theory ofpoint processes, Water Resources Research, Vol. 17, No.
5., pp. 1273-1285, 1981b.
4) Waymire, E. and Gupta, V.: The mathematical structureof rainfall representations 3: Some applications of thepoint process theory to rainfall processes, Water
Resources Research, Vol. 17, No. 5., pp. 1287-1294,
1981c.
5) Rodriguez-Iturbe, I., Gupta, V. and Waymire, E.: Scaleconsiderations in the modeling of temporal rainfall,
Water Resources Research, Vol. 20, No. 11., pp.
1611-1619, 1984.
6) Lowen, S. and Teich, M.: Estimation and simulation of
fractal stochastic processes, Fractals, Vol. 3., No. 1, pp.
183-210, 1995.
7) Lowen, S. and Teich, M.:
Fractal-Based Point
Processes, Wiley Series in Probability and Statistics,
New Jersey, 2005.
8) Telesca, L., Lapenna, V., Scalcione, E. and Summa, D.:
Searching for time-scaling features in rainfall sequences,
Chaos, Solitons, and Fractals, Vol. 32, pp. 35-41, 2007.
9) Rodriguez-Iturbe, I., Cox, D. and Isham, V. : Somemodels for rainfall based on stochastic point processes,
Proc. R. Soc. London A, Vol. 410, pp. 269-288, 1987.
10) Cowpertwait, P.: Further developments of the
Neyman-Scott clustered Poisson process for modelingrainfall, Water Resources Research, Vol. 27, No. 7, pp.
1431-1438, 1991.
11) Press, W., Teukolsky, S., Vetterling, W. and Flannery,B.: Numerical Recipes in FORTRAN, 2nd ed. Cambridge
University Press, London, 1992.
(Received September 30, 2007)

Fig. 6. Historical vs. synthetic rainfall Kolmogorov-Smirnov
test using 15 monthly maximum values per set.

5. CONCLUSION
A new expression referred to here as the Fanofactor exponent was determined for the
Neyman-Scott Clustered Poisson Rectangular
Rainfall Model, or NSM.
Unlike previous
parameter search methods for the NSM, thisversion linked the model parameters to the Fanofactor of a portion of the historical data based on
the POT rainfall series. Including this Fanofactor exponent expression in the search estimatedNSM parameters that were satisfactory for
generating synthetic rainfall with significant mean
total rainfall, variances, autocovariances andmonthly maxima.
However, in several cases the same parameters
generated synthetic block maxima that wereinconsistent with the historical block maxima in
empirical distribution. To rectify this limitation
during these months, rain cell based information
should be prioritized in the parameter search.
For this purpose, the 6-hourly lag 1
autocorrelation should be used instead of the
12-hourly lag 1 autocorrelation.
ACKNOWLEDGMENTS
This study was supported by Grant-in-Aid for
Scientific Reserch (C) 18560497 (PI: Assoc.

- 102 -