- Email: [email protected]
A spatial rainfall generator for small spatial scales
Journal of Hydrology 252 (2001) 126±144
www.elsevier.com/locate/jhydrol
A spatial rainfall generator for small spatial scales Patrick Willems* Hydraulics Laboratory, Katholieke Universiteit Leuven, Kasteelpark Arenberg 40, B-3001 Heverlee, Belgium Received 29 March 2000; revised 9 May 2001; accepted 22 May 2001
Abstract A stochastic spatial rainfall generator is developed for use at the small spatial scale of urban and small hydrographic catchments. The generator is based on a spatial rainfall model of the conceptual and hierarchical type. It describes the spatial rainfall ®eld in a macroscopic physically-based way by distinguishing rainfall entities with different scales: rain cells, cell clusters, small and large mesoscale areas (or rain storms). For applications at small spatial scales, the individual rain cells need a detailed description. Data of a dense network of rain gauges at Antwerp, enclosing 5940 rain cells in 807 rain storms are used to derive such description. For separation of the rain cells in the rainfall time series, an algorithm is developed based on the identi®cation of increasing and decreasing rain cell ¯anks. The rain cells observed at different rain gauges are linked together by applying criteria for testing the similarity in rain cell properties. After separating and linking the rain cells and storms, the spatial rainfall model is calibrated to many storms by two methods (e.g. Kalman ®lter). The derived model structure and model parameter distributions apply to the stochastic generation of long-term time series of spatial rainfall. The model is tested by comparing intensity±duration±frequency relationships and temporal scaling properties of the generated and historical rainfall series. q 2001 Elsevier Science B.V. All rights reserved. Keywords: Rain cells; Time correlation; Rain storm
1. Introduction Much research has already been done trying to describe the organization in spatial rainfall ®elds; see for example Bras and Rodriguez-Iturbe (1976), Gupta and Waymire (1979), Waymire and Gupta (1981), Waymire et al. (1984) and Sivapalan and Wood (1987), to name a few. In this research, the `clustering' of high intensity rainfall areas embedded within rainfall areas of lower intensity has been determined as a characteristic pattern. The pattern is observed at many different scales. At the smallest scale, the individual rain cell forms the building* Fax: 132-16-32-18-98. E-mail address: [email protected] (P. Willems).
block of the spatial rainfall structure. The rain cells are embedded in a clustered way within `small mesoscale areas' (10 2 ±10 3 km 2). They move in the same small mesoscale area with nearly identical velocity. At larger scales, small mesoscale areas occur in a clustered way within lower intensity `large mesoscale areas', which in turn are embedded within some synoptic-scale lowest intensity rainfall ®eld. The large mesoscale areas are observed as (sometimes elongated) bands and their spatial extent varies from 10 3 to 10 4 km 2. The precipitation associated with cyclonic and frontal storms are of the level of these large mesoscale areas. Although rain cells and cell clusters most often appear in large and small mesoscale areas, they can also occur isolated outside these regions (e.g. in air mass thunderstorms).
0022-1694/01/$ - see front matter q 2001 Elsevier Science B.V. All rights reserved. PII: S 0022-169 4(01)00446-2
P. Willems / Journal of Hydrology 252 (2001) 126±144
127
Fig. 1. Gaussian-shapped rain cells on the basis of the structure of spatial rainfall (map area of Belgium).
This spatial rainfall structure, which is also represented in the radar view of Fig. 1, has served as the basis for many space-time rainfall modeling applications. Wilson et al. (1979), Bras et al. (1988), Bartolini and ValdeÂs (1994), Bierkens and Puente (1990), Northrop (1998) and Sivapalan and Wood (1987) and many others used such a rainfall model to analyse the in¯uence of rainfall input data on the results of hydrological models. Most of these applications deal with large watersheds. Almost no applications at very small spatial scales exist, such as applications to urban hydrology and small hydrographic catchments. For these applications, the accurate description of individual rain cells and cell clusters is very important, while the description of the spatial rainfall structure at the scale of large mesoscale and synoptic areas is of minor importance. The spatial extent of individual rain cells (10±50 km 2 following Waymire et al., 1984) is indeed of the same order of magnitude as the spatial extent of most urban and small hydrographic catchments. Rain cells, however, have a rather super®cial description in the above mentioned research. There is most often lack of a solid empirical basis for their modeling. This basis should consist of a very dense network of rain gauges. Its availability is however rather scarce. The present study focusses on the accurate spatiotemporal description of rain cells and rain storms at the small spatial scale of urban and small hydrographic catchments. The `solid empirical basis' for
Fig. 2. Rain gauge network of the City of Antwerp, together with the ®ve sub-networks A1,¼,A5.
the study consists of a dense network of 12 tippingbucket rain gauges at the City of Antwerp. The gauges are located at distances of approximately three kilometres and they cover a large urban area of about 100 km 2 (Fig. 2). This area is larger than the maximal spatial extent of rain cells and it is of the same order of magnitude as the spatial extent of cell clusters. It also approximates the area of hydrographic subcatchments and the largest urban catchments. The rain gauges record rainfall intensities every minute with a resolution of 0.1 mm since July 1994. 2. The structure of rain cells On the basis of the data of the Antwerp network, the basic structure of rain cells has been analyzed by Luyckx et al. (1998), Willems and Berlamont (1998) and Willems (1999). The structure includes two aspects: the spatial distribution of rainfall intensities in a rain cell (the shape of the rain cell) and the
128
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 3. Calibration of the average value of the parameters rmax and sy of the rain cell model (the points indicate the twelve rain gauges P1, P2,¼, P12).
movement of the rain cell in time and space (transport of the rainfall intensities). It has been shown by Luyckx et al. (1998) that the spatial distribution of rainfall intensities in a single rain cell can be described well by a bivariate Gaussian distribution: r x; y; t rmax exp 2 2
x 2 x0 2 u t 2 t0 2 2s2x !
y 2 y0 2 2 g t 2 t0 2s2y
1
In this formula, r x; y; t represents the rainfall intensity at spatial co-ordinates (x,y) and time t. The x-direction is taken equal to the direction of rain cell movement u . The other rain cell model parameters consist of: the maximum rainfall intensity rmax (the amplitude of the rain cell), the spatial extent in the x- and y-direction) in terms of the standard deviation of the Gaussian distribution) sx and sy, the mean velocity u and the decay velocity g . x0 ; y0 are the spatial co-ordinates of the rain cell at a reference time t0. To describe the movement of the rain cells, Jinno et al. (1993) have shown that a two-dimensional advection±diffusion model can be used. 2r 2r 22 r 22 r 1u dx 2 1 dy 2 2 gr 2t 2x 2x 2y
2
with dx and dy the diffusion coef®cients in the x- and the y-direction. This advection±diffusion model has already been applied with success by Jinno et al. (1993) and Kawamura et al. (1997) in real-time
control applications. Together with the decay of cell volumes, advection and diffusion are indeed observed as the most important processes in the movement of rainfall intensities. On the other hand, a Gaussian distribution can analytically be determined as a possible solution to the advection±diffusion model r0 p exp r x; y; t 4p dx dy t 2 t0 x 2 x0 2 u t 2 t0 2 4dx t 2 t0 ! y 2 y0 2 2 g t 2 t0 2 4dy t 2 t0 £
2
3
Comparing Eqs. (1) and (3), it can easily be seen that, at time t, the spatial extent of the cell equals p p sx 2dx t 2 t0 ; sy 2dy t 2 t0 4 By diffusion, the cell extent will increase in time and the maximum rainfall intensity of the cell rmax will decrease r0 p rmax 5 4p dx dy t 2 t0 where r0 represents the total rain volume of the cell before decay of the cell volume starts. 3. Calibration of the rain cell model To calibrate the rain cell model, two methods are tested and compared: the semi-automatic (and
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 4. Calibration of the average value of the parameter sx of the rain cell model of optimizing the agreement with a standard Gaussian curve.
statistical) Kalman ®lter technique applied as recursive parameter estimator (Meinhold and Singpurwalla, 1983; Willems and Berlamont, 1998) and a more manual method in which the model parameters are related to physical properties that can be derived from the rainfall time series. In this latter method, the moving-velocity u and the moving-direction u are ®rst calibrated by considering two couples of rain gauges. The moving-times between the two gauges of each couple can be expressed in terms of u and u and equalized to the real moving-times (see also Luyckx et al., 1998). By considering two couples, the two unknowns (u and u ) can be calculated from the two equations. The real moving-times are not known, but they can be estimated from the time between the measured peak intensities. By calculating u and u for each two couples of rain gauges and by averaging the results, the model parameters u and u can be estimated accurately. When the moving-direction is known, the peak intensities at the different rain gauges can be plotted against the distances of the rain gauges in the direction perpendicular to the movingdirection (the y-direction (Fig. 3). From this plot, the average values of the maximum rainfall intensity rmax and the spatial extent in the y-direction sy can be derived. In Fig. 3, y0 is the distance (in the y-direction) of the rain cell peak. The spatial extent in the moving(or x-) direction sx (or the small scale spatial correlation in the direction of storm movement) can be calibrated on the basis of the small scale time correlation as observed in the rainfall times series of the different rain gauges. This time correlation can be transformed to space using the mean storm velocity
129
u. The transformation is based on the hypothesis of ¯uid turbulence by Taylor (1938), which states that the spatial pattern of turbulent motion, as re¯ected in the spatial correlations, is carried past a ®xed point (here the rain gauge) by the mean ¯ow speed without undergoing any essential change. The hypothesis is only valid for uniform mean ¯ow and a low level of turbulence (or small `noise' on the average motion). The usefulness of Taylor's hypothesis for a spacetime rainfall storm was shown earlier by Zawadski (1973), Bras and Rodriguez-Iturbe (1976) and Waymire et al. (1984). The following transformation formula is therefore applied: h
xi 2 x0 2 u t 2 t0 2 y 2 y 2 2 i 20 2sx 2sy
6
In this equations, x0 is the position (in the x-direction) of the rain cell peak at time t0; xi ; yi are the spatial co-ordinates of the rain gauge i. The value x0 can be calibrated by maximizing the agreement between the measured time moments of the peak rainfall intensities and the modelled ones ti t0 1 u xi 2 x0
7
In Fig. 4, the rainfall intensities at the different rain gauges i i 1; ¼; 12 and at different moments t are plotted against the relative distance h to the rain cell centre. For the calibrated value of sx and for Gaussianshaped rain cells, Fig. 4 approximates the standard Gaussian curve r exp 2h 8 rmax By scaling the rainfall intensities to the maximum rain cell intensity rmax ; the maximum of this curve equals 1. After testing and comparing the two calibration techniques (Willems, 1999), a combination of the techniques is recommended: application of the manual method ®rst to derive `a priori' estimates of the model parameters for the Kalman ®lter. The Kalman ®lter technique has major advantages because it takes into consideration rainfall measurement errors and uncertainties in the rain cell model structure. As rainfall measurement errors depend on the rainfall intensity, this consideration of the errors is important for the derivation of unbiased parameter values. In the Kalman ®lter technique, the model parameters are continuously updated in time (as
130
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 6. Results of the linking of rain cells at different rain gauges for the rain storm on 14.11.94 (identi®ed rain cells with signi®cant volumes are represented by bubbles, bubbles with the same ®lled color are identi®ed as the same rain cell, white bubbles represent rain cells that cannot be identi®ed).
4. Rain cell separation Fig. 5. Results of the linking of rain cells at different rain gauges for the rain storm on 14.11.94 (identi®ed rain cells are represented by triangles, triangles with the same ®lled color are identi®ed as the same rain cell).
more measurements become available) and a statistical estimate is made of their uncertainties. In this way, the variability in time of model parameters can be recognized. By the automatic nature of this calibration method, calibration times are very small. An important disadvantage of the method is the dif®cult control of the physical realism of the parameter values during the calibration procedure. By using the results of the manual method as initial values for the model parameters, this problem is however often solved. This consideration can be made since the experience of the modeller is used in the manual method. The manual calibration procedure is however more time-consuming, although it can be partly automated by different operations in a spreadsheet computer programme. To reduce the in¯uence of irregularities in the measurements (e.g. measurement noise) on the calibration results, the original 1 min time series are processed ®rst by moving average over 10 min intervals. The smoothing of 1 min peak intensities, caused by this processing, is of no relevance to hydrology because a much higher smoothing is described in the rainfall-runoff transfer of the smallest hydrological system.
For rain storms with more than one rain cell, the different rain cells have to be identi®ed and separated in the rainfall time series. An algorithm is derived for this purpose. It is described in detail in the appendices and consists of two parts. In the ®rst part (Appendix A), the different rain cells in the point rainfall time series are separated for each rain gauge independently. This separation is based on the identi®cation of peaks, decreasing ¯anks and increasing ¯anks in the time series. In the second part (Appendix B), the identi®ed rain cells for the different rain gauges are linked together. This means that the algorithm identi®es rain cells in the different time series, which correspond with the same spatial cell. As a result of the rain cell separation and identi®cation process, the different spatial rain cells in each rain storm are numbered and some physical rain cell properties are derived. In Fig. 5, the results of the algorithm are shown for the rain storm on 14.11.94. Identi®ed rain cells are represented in this ®gure as triangles. Rain cells with the same color are linked together. The cells with insigni®cant volume (volume smaller than 0.4 mm) are not considered in the linking procedure. They are represented with white color. As can be seen in the ®gure, the linking procedure performs well. In Fig. 6, the same results are represented in a different way. The identi®ed rain cells with signi®cant volume are represented by bubbles for which the magnitude is proportional to the cell volume. The non-colored cells are the ones that could not be
P. Willems / Journal of Hydrology 252 (2001) 126±144
131
Fig. 7. Calibration results for three rain cells in the rain storm on 28.08.96.
linked. They are considered as unidenti®able. The result of the rain cell separation and linking procedure are used for two purposes. They are ®rstly used for the identi®cation of subperiods in the rainfall time series for which the individual rain cell model has to be calibrated. Secondly, they are used for counting the number of rain cells in rain storms, as explained in Section 6.
5. Probability distributions of rain cell parameters According to the procedures of Sections 3 and 4, about 104 rain cells in 58 rain storms are calibrated. They are randomly selected in the period October 1994±February 1997. Fig. 7 shows an example of
132
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 8. Calibration results of the probability distributions for the rain cell model parameters rmax, s, u and u .
the calibration results for three rain cells in the rain storm on 28.08.96. After calibration of the rain cell model for a large number of rain cells, statistical properties of the model parameters are derived in terms of probability distributions and correlations. For the rain cell parameters, clear distribution functions are determined, as shown in Fig. 8 for the parameters rmax ; s, u and u . The derived distributions and correlation coef®cients of the different model parameters are summarized in Tables 2 and 3. As mentioned in Section 2, the spatial extent of a rain cell sx and sy is not constant during its movement over a spatial area. It was explained by the diffusion process. For the detailed description of rain cells, as needed in applications of real time prediction, a detailed consideration of the diffusion and decay processes is important. However, for purposes of stochastic spatial rainfall generation, this consideration is less important. The time the rain cells need to move over an urban or small hydrographic catchment is rather small and the inclusion of the diffusion and decay processes will only change the statistical properties of the generated spatial rainfall in an insignificant way. Therefore, for each rain cell the average
spatial extent during its movement over the rain gauge network is calibrated and the decay velocity is neglected. The spatial extent in both x- and y-directions is not signi®cantly different. Therefore, also a combined distribution is derived for sx and sy : The mean spatial extent of a rain cell in both directions equals 2.5 km (in terms of the standard deviation of the Gaussian shape). If the diameter of a rain cell is considered as six times the standard deviation, the mean diameter equals 15 km. This diameter encircles almost the total rainfall volume (99.7%) of the rain cell. 6. The cell cluster model The spatial organization of the rain cell centers in small mesoscale areas or cell clusters is studied using the numbering results as explained in Section 4. From the temporal density of rain cells in rain storms (or the average number of rain cells per time unit), the spatial density can be calculated using the moving velocity and direction of the rain storm. Indeed, the spatial extent (in the moving direction u ) of the rain storm
P. Willems / Journal of Hydrology 252 (2001) 126±144
133
Fig. 9. Schematic representation of the cell cluster model (rain storm location shown for time moment t t0 ).
ls u can be calculated from the temporal duration ds and the moving velocity u of the rain storm ls uds 2 la u
9
la u is in this formula the spatial extent of the area covered by the rain gaue network in the direction u (see Fig. 9). The rain storm duration ds is de®ned by the time span from the beginning of the ®rst rain cell (the earliest beginning for all rain gauges) to the ending of the last rain cell (the latest ending for all rain gauges). It is clear that the time the rain storm needs to move over the network area la u=u should be subtracted from the storm duration ds to derive the real rain storm length ls : The rain storm width (the
extent of the rain storm in the direction perpendicular to the moving direction) is assumed larger than the extent of the area covered by the rain gauge network in the same direction wa u: To study the organization of the rain cell centres, the number of rain cells is then analysed in relation with the spatial area V s over which the cells as observed by the rain gauge network are spread (see also Fig. 9)
V s ls wa u 1 2w
10
It is assumed in this formula that all rain cells with a cell centre closer than a distance w from the most remote rain gauges are observed by the rain gauge network. The distance w is taken equal to the average
Fig. 10. Analysis of the number of rain cells per unit area in rain storms.
134
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 11. Cumulative distribution function of the rain storm duration ds (comparison of the empirical distribution and the distribution simulated from an exponential rainstorm length distribution).
spatial radius of a single rain cell (7.5 km, see Section 5). For the most simple structure of a spatial rain cell process, being the simple spatial Poisson process, the only parameter of the process (l ) can be calculated as the average number of rain cells per unit area. In Fig. 10, the observed number of rain cells is for all calibrated rain storms plotted against the spatial area of the rain gauge network V s. The latter area is calculated by Eq. (10). The slope in the ®gure equals the parameter l
l
average number of rain cells area V s
11
To derive the ®gure for a wide range of areas, ®ve sub-networks of the Antwerp rain gauge network are considered. They are indicated in Fig. 2 as A1, A2, A3, A4 and A5. Sub-network A1, for instance, only contains rain gauge no. 6, while sub-network A2 also includes the gauges 3, 4, 5, 7 and 8. The ®ve sub-networks have their own values for la u and wa u; which are calculated from rain gauge map coordinates. From Fig. 10, an average value of l 0:002 is derived. There is however a large scattering in the ®gure. This may indicate that other more complicated spatial process descriptions are needed. For the current study, however, a simple Poisson process is considered. This process is assumed very often for the stochastic description of sequences of rain cells in rain storms both in time and space. This is the case, for instance, in the cluster-based rectangular pulses models (e.g. Neyman±Scott model,
Bartlett±Lewis model), which have been calibrated earlier for temporal Belgian rainfall data by Velghe et al. (1993) and Verhoest et al. (1997). It is clear that the larger the area of the rain gauge network, the larger the maximum extent of the spatial organization of rain cells that can be considered in the analysis. For the Antwerp rain gauge network, it can be studied well up to the spatial level of small mesoscale areas. To ®nalize the cell cluster model, it is assumed that all rain cells in the same cell cluster move with the same velocity and direction. It is thus assumed that the cells' velocities relative to the storm velocity are essentially zero. This assumption is based on the early work of Zawadski (1973) and, after calibration of many rain storms to the Antwerp data, it is seen to be a good assumption. 7. Probability distributions of cell cluster parameters Also for the other cell cluster parameters, such as the inter-arrival time of rain storms and the rain storm duration, the information is derived from the rain cell separation and linking results as explained in Section 4. All rain cells and rain storms in the period under study (August 1994±May 1997) can be included in this analysis. They consist of 5940 rain cells in 807 rain storms. For each rain storm, the following information is derived: the inter-arrival time, the storm duration and the number of rain cells.
P. Willems / Journal of Hydrology 252 (2001) 126±144
135
Fig. 12. Exponential quantile plot for the inter-arrival time ta of rain storms.
The length of the cell cluster areas (or small mesoscale areas) can only be calculated for the rain storms with known cell movement direction. Because this information is only available for the 104 calibrated rain cells, a different approach is used for calibration of the distribution of ls : On the basis of an assumed distribution for ls and the known Weibull and normal distributions for u and u , the distribution of the rain storm duration ds is simulated on the basis of Eq. (9). A comparison with the empirical distribution of ds allows to determine the best distribution for ls : An exponential distribution is assumed by many authors for the rain storm lengths. This distribution is therefore tested ®rst. It is seen in Fig. 11 that for an exponentially distributed ls with a mean storm length of 1=l 35:7 km; the simulated distribution of ds matches well the empirical distribution. The range of rain storm lengths described by this exponential distribution corresponds with the range indicated by Waymire et al. (1984): 10±10 3 km 2 for cell clusters or small mesoscale areas. For the inter-arrival time of rain storms, the probability distribution is calibrated directly to the empirical distribution. A two-component exponential distribution is determined (Fig. 12) ta days 2 0:0354 F ta 0:15 1 2 exp 2 5:5 t days 2 0:0354 1 0:85 1 2 exp 2 a 0:5 12
The two-components show the existence of two populations of exponentially distributed inter-arrival times. For each population, the exponential distribution indicates that the rain event arrivals correspond to a temporal Poisson process. The second population with small inter-arrival times (mean inter-arrival time of 0.5 days) corresponds to inter-arrival times between different rain storms (or small meso-scale areas) in the same large meso-scale area (e.g. frontal storm). They vary from 50 min to many hours. The ®rst population with large inter-arrival times (mean interarrival time of 5.5 days) corresponds to inter-arrival times between different large meso-scale or synoptic areas. They may have values up to 30 days (see Fig. 12). The two-component distribution for ta can probably also be derived when the spatial organization of small meso-scale areas is described by a spatial Poisson process (in the same way as the spatial organization of rain cells in small meso-scale areas is described). However, because of the limited spatial extent of the Antwerp rain gauge network, this spatial clustering cannot be studied for spatial scales larger than small meso-scale areas. The in¯uence of the spatial organization of small meso-scale areas is therefore modeled in an indirect way by describing the temporal clustering of the smaller inter-arrival times, which is a result of this spatial organization. This temporal clustering of the second population of inter-arrival times could be described by an exponential distribution for the duration of a cluster of small
136
P. Willems / Journal of Hydrology 252 (2001) 126±144
Table 1 Overview of the probability distribution function for the model parameters Model
Distribution
Distribution parameter values
Mean
u (km/min) u (degrees) sx (km) sy (km) s (sx and sy together) rmax (mm/h) ls (km)
Weibull Normal Lognormal Lognormal Lognormal Generalized Pareto (GPD) Exponential
t Weibull index 0:70; b 2:36 m u 221:75; s u 68:46 m ln s x 0:90; s ln s x 0:41 m ln sy 0:92; s ln s y 0:46 m ln s 0:91; s ln s 0:44 g Extreme value index 0:396; b 1:65; r t 0:6 l 0:028
0.6 2 21.75 2.5 2.5 2.5
and properties are derived for the combined spatiotemporal rain storm process at the small spatial scale. They consist of the rain cell and cell cluster model structures (Fig. 9, Eq. (1)), the probability distributions of the rain cell parameters and their correlations (Tables 1 and 2, Fig. 8), the l parameter of the spatial Poisson process of rain cells in rain storms and the statistical models for the rain storm duration and the inter-arrival time of rain storms (Eq. (12)). On the basis of these properties, a spatial rainfall generator can be easily constructed. Using a random number generator, a large number of spatial rain storms can be randomly simulated. The long-term time series of spatial rainfall that are generated in this way have the same statistical properties compared with those derived from the rain gauge network data. The simulated long-term time series of spatiotemporal rainfall can be used for many applications. Two examples are given hereafter. The time series can be used as input to hydrological models to study the in¯uence of the spatial variability of rainfall on the model results. They also can be used to simulate spatially correlated rainfall simultaneously for different subsystems in an integrated catchment model, such as the sewer system and the hydrographic catchment. Due to the limitation of the spatial scale, the
meso-scale areas (mean duration of 3.0 days), combined with a Poisson process for the description of the different small inter-arrival times inside such a cluster. In Fig. 12, a comparison is shown for the inter-arrival times simulated by this model and the ones derived from the Antwerp time series by the rain cell identi®cation procedure. The model shows a small underestimation for the intermediate inter-arrival times. Marien and Vandewiele (1986) also determined a compound distribution for the inter-arrival times of rain storms in a point rainfall series. They considered three components, one of them corresponding with the inter-arival times between different rain cells in a rain storm (modeled here in an indirect way by the spatial Poisson process). For the two other components, an exponential and a Weibull distribution were considered. Exponential and Weibull distributions were also used by Eagleson (1978) and Amorocho and Wu (1977) for the inter-arrival times for large mesoscale areas. 8. Spatial rainfall generator In the previous sections, basic conceptual features Table 2 Correlation matrix for the model parameters
u (km/min) u (degrees) sx (km) sy (km) rmax (m/h)
35.7
u (km/min)
u (degrees)
sx (km)
sy (km)
rmax (mm/h)
1 ± ± ± ±
0.15 1 ± ± ±
0.86 0.15 1 ± ±
0.13 0.006 0.35 1 ±
2 0.25 0.12 2 0.22 ±0.13 1
P. Willems / Journal of Hydrology 252 (2001) 126±144
137
Fig. 13. Model validation by comparison of IDF-relationships. (a) for the ®nal rainfall generator. (b) for a rectangular rain cell model with and without clustering of small inter-arrival times.
generator is, of course, only applicable to sewer system catchments and hydrological cachments of limited spatial extent. 9. Model testing by IDF-relationships The accuracy of the spatial rainfall generator is tested by comparing the statistical properties of the stochastically generated long-term time series of spatial rainfall and these of historical rainfall series.
The comparison is made for different statistical properties: extreme value distributions, IDF-relationships, autocorrelations and scale properties. Apart from the Antwerp data on the basis of which the rainfall generator is built and calibrated, also the Ukkel/ Uccle series for 27 years (1967±1993) is used as historical series. Through the greater time-span of this latter series, the accuracy of model extrapolations to rare events (which are not available in the spatial rainfall data at Antwerp) can be tested. The extreme value analysis of the generated and historical time
138
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 14. Model validation by comparison of the extreme value distribution (example of an exponential quantile plot for an aggregation-level of 30 min).
series and the calculation of IDF-relationships is performed according to the methodology described in Willems (2000). In Fig. 13(a), the comparison is shown for intensity±duration±frequency-relationships up to 27 years, which were derived for the Uccle series in Willems (2000). The con®dence limits are indicating the differences in IDF values between the 12 individual rain gauges of the Antwerp network. Good results are shown, both for rainfall frequency and for a large range of temporal scales. Even the temporal scaling properties, indicated by the linear shape of the IDFcurves in Fig. 13, are well represented. Apart from this good general performace of the generator, two minor differences can be observed. Asymptotically to the smallest aggregation-levels (smaller than 20 min), the IDF-relationships show a higher bending down for the rainfall generator. This is explained by the full deterministic nature of the rain cell model. By adding a random component in terms of high frequency noise to this model, the bending down can be corrected. Another difference is the small underestimation of the rainfall volumes at large aggregation-levels (larger than one day). A sensitivity analysis has shown that this underestimation is very sensitive to the model for the smaller inter-arrival
times. If the short-term dependency in inter-arrival times is totally ignored in the model, the underestimation will increase with a factor three. This is shown in Fig. 13(b) by randomly changing the time order of the different rain cells that were identi®ed in the original Antwerp rainfall time series. Also by ignoring the smaller rain cells in the model, the underestimation of rainfall volumes for large aggregation-levels increases. For instance, by ignoring all rain cells with a volume less than 0.4 mm, the decrease in rainfall intensities (Dr) given in relation to the aggregation-level da by the following formula: Dr mm=h 0:02 1
0:3 da days
13
In addition to the IDF-relationships, also the individual distributions of rainfall intensities at the different aggregation-levels are compared. An example is shown in Fig. 14 for an aggregation-level of 30 min. The shape of the distribution's tail (which is linear for an exponential distribution in the exponential quantile plot of Fig. 14) is well reproduced by the generator. This is of primary importance for an accurate extrapolation to extreme events.
P. Willems / Journal of Hydrology 252 (2001) 126±144
139
10. Conclusions A spatial rainfall model is derived for use in hydrological applications at small spatial scales. The model is of the conceptual and hierarchical type and its structure is similar to many existing models for larger meso-scale areas. It differs from these models by the more detailed description of the smallest buildingblocks of a spatial rainfall structure: the individual rain cells and the small scale rain cell clusters (or small meso-scale area). The description is based on a detailed analysis of the observed cell cluster patterns by a dense network of rain gauges at the city of Antwerp. Another feature is the development of an algorithm to identify and separate different rain cells in the rain gauge time series and to link together identi®ed rain cells for different rain gauges. The rain cell separation is based on the identi®cation of increasing and decreasing rain cell ¯anks, while the rain cell linking is based on criteria for testing the similarity in rain cell properties. By extending the deterministic model with a stochastic description of model parameters, the model is very useful for the long-term stochastic generation of small scale spatial rainfall. It is shown that the most important statistical properties of rainfall, like temporal scaling properties and extreme value statistics, are well reproduced by the generator. Other properties, such as spatial coverage statistics and spatial scaling properties still have to be validated. It is however expected that also these properties are described reasonably well due to the physically-based nature of the approach that was followed for the structure identi®cation and the calibration. This approach differs from other approaches where the parameters of an assumed model-structure are calibrated by matching global statistical properties such as spatial correlations and spatial coverage statistics. Because of its macroscopic physically-based structure and the applicability of the recursive Kalman ®lter calibration procedure, the model described in this paper may also be useful for short-term rainfall prediction in real-time control applications. Acknowledgements This research was supported by a scholarship of the
Fig. 15. Criteria for associating states to the rainfall intensities at the different 1 min time steps in the rain cell separation procedure.
Flemish Institute for Scienti®c and Technological Research (IWT) and a project for the Programme `Beleidsgericht Onderzoek' of the Ministry of the Flemish Community. The Antwerp rain gauge network data was made available by the `City of Antwerp ± Department voor Werken', while the long-term historical rainfall data at Uccle was processed and delivered by the Royal Meteorological Institute of Belgium (KMI/IRM). The author is also grateful to two reviewers which have made important suggestions for improvement of this paper. Appendix A. Algorithm for rain cell separation and linking In this appendix, a description can be found of the algorithm developed and applied in section 4 for the identi®cation and separation of different rain cells in the rainfall time series. The procedure is described ®rst for a single point rainfall series (Appendix A.1). In Appendix A.2, the more advanced procedure of linking of rain cells observed at different rain gauges is explained.
140
P. Willems / Journal of Hydrology 252 (2001) 126±144
² The rainfall intensity is a minimum (state 0) ² The rainfall intensity is part of an increasing ¯ank (state 1) ² The rainfall intensity is part of a decreasing ¯ank (state-1).
Fig. 16. Criteria of Marien and Vandewiele (1986) for de®ning minima in the 10 min rainfall time series (to de®ne the end of partial storms).
A.1. Rain cell separation in a point rainfall time series The algorithm for the separation of different rain cells in a point rainfall time series is based on the identi®cation of peaks, decreasing ¯anks and increasing ¯anks in the time series. The identi®cation is performed by comparing the rainfall intensity at each time step in the rainfall series with the rainfall intensities at the neighbouring time steps. In this way, the rainfall intensity at each time step can be associated with one of the three following possible conditions:
These possibilities are hereafter called the `states' of the rainfall intensities and are given the number 0,1 and 21.When the identi®cation of a rainfall intensity is not clear, it is given `state 3'. The criteria for discriminating among the four state conditions are summarized in both a graphical and numerical way in Fig. 15 for minima and increasing ¯anks. For decreasing ¯anks, the criteria are analogous to the ones of increasing ¯anks. A rainfall intensity of zero is always given a state number 0. The criteria in Fig. 15 are based on a comparison of the rainfall intensities at two preceding time steps r 21 and r 22 and two succeeding time steps (r(1) and r(2)). In the ®rst two criteria for increasing ¯anks, local minima in increasing ¯anks are identi®ed. They are eliminated from the time series by replacement of the rainfall intensity at the local minimum by the average rainfall intensity for the previous and the next time step. The criteria are based on the work of Marien and Vandewiele (1986). They identi®ed minima in a 10 min point rainfall time series to separate different rain cells (they call it `partial storms'). The criteria that they used to select these minima are summarized in Fig.16 in the form of Fig. 15. All minima in the time series that do not meet these criteria were considered
Fig. 17. Results of the rain cell separation procedure for the rain storm on 14.11.94 (rainfall time series of rain gauge 12).
P. Willems / Journal of Hydrology 252 (2001) 126±144
as local minima in increasing or decreasing ¯anks. Different criteria are used in this study because of the different time step of the time series (1 min instead of 10 min). From the list of state numbers at all time steps, the different rain cells are identi®ed together with the time steps indicating the beginning, the ending and the peak-moment of the rain cells. This is also performed by an algorithm. The following rules are used: ² the beginning of a rain cell is de®ned by the ®rst state 1 after a list of states 0 or 3 (e.g. `¼033331¼' or `¼01¼'; the bold position here indicates the starting position of a rain cell). ² the peak-moment of a rain cell is de®ned by a state between an increasing and a decreasing ¯ank (e.g. `¼1133311-1-1¼'). An increasing ¯ank is de®ned here by a list of states 1 and 3, while a decreasing ¯ank corresponds to a list of states 21 and 3. Whenever the increasing and decreasing ¯anks are separated by a list of states 3, the peakmoment is not de®ned that clear (e.g. `¼113333113333-1-1¼'). The peak-moment is then selected as one of the states 3 more speci®cally the one that is closest to the middle between the starting and ending moments of the rain cell. ² the ending of a rain cell is de®ned by a state 0 after a decreasing ¯ank (e.g. `¼3-133-1-10¼'). Whenever states 3 are listened in between the decreasing ¯ank and the state 0, the ending is considered at the ®rst state 3 (e.g. `¼3-133-1-133330¼'). Such list of states 3 indicates that the ending of the rain cell is not that clear. However, if this list would be considered part of the decreasing ¯ank, the slope of the ¯ank would be too `weak'. For the illustration of such a weak decreasing ¯ank, the reader can look at the third identi®ed rain cell in Fig. 17. The weakness in this ®gure is however not caused by the above-mentioned identi®cation problem. It is just caused by the length of the ¯ank. Sometimes, when two succeeding rain cell are close together, no state 0 exist between the two rain cells. The ending of the ®rst rain cell and the beginning of the second rain cell are then both de®ned by the last state 21 in the decreasing ¯ank (or the ®rst state 3 if both rain cells are separated only by a list of states 3, e.g. `¼3-133-1-1333311131¼').
141
² To avoid that a short and insigni®cant increasing ¯ank in a longer decreasing ¯ank (or local noise in a decreasing ¯ank) is identi®ed as the increasing ¯ank of an additional rain cell, a minimum value (e.g. 3 min) is required for the length of an increasing ¯ank. The following list of states is thus considered as a continuously decreasing ¯ank: `¼3-133-1-11-13-1¼' (the single state 1 (in bold) is ignored). ² Sometimes, there is an overlap in time between two succeeding rain cells. As the rainfall intensity measurements are representing the sum of the rainfall intensities of two rain cells during the overlap period, the starting time of the ®rst rain cell and ending times are not de®ned correctly by the above listed rule. Indeed, the starting time of the ®rst rain cell and the ending time of the second rain cell are assumed equal to the time at which the total rainfall intensity approximately reaches a minimum. To asses the real starting and ending times, the ending time of the ®rst rain cell is repeatedly increased with one time step and the starting time of the second rain cell decreased in the same way until the minimum measured rainfall intensity in the overlap period is reached by the model. Instead of using the Gaussian model for this purpose, a simpli®ed triangular model (see later on) is used. The procedure is applied for each pair of rain cells with less than 10 min between the starting and ending times. To achieve a better separation of rain cell, the 1 min rainfall time series (already smoothed by moving average over intervals of 10 min) is ®rst additionally smoothed by moving average over intervals of 5 min. This is done because the `noise' on the rainfall intensities was still too high, leading up to many wrongidenti®ed rainfall intensities (many local maxima not corresponding with rain cell peaks). As can be seen in Fig. 17, the peak intensities are ¯attened a little by the additional moving-average operation. However, the identi®cation of minima, increasing and decreasing ¯anks is easier and is not disturbed by this operation. In Fig. 17 the results of the state identi®cation and the rain cell separation are shown for the rain storm of 14.11.94. Only for graphical reasons, the rain cells are represented by triangular shapes. The time steps at the beginning, ending and peak-moments of the triangles
142
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 18. Different possibilities for the selection of the most likely rain cell with a certain rain cell number (indicated as known) in a list of three succeeding rain cells with the same rain cell number. The rain cells indicated in `bold' are at each step subject to the criteria for selection of the most likely rain cell from two rain cells with the same number.
are derived from the identi®ed ones of the separated rain cells. The height of the triangle-peak is calculated by equalizing the total rainfall volume in the triangle and the original rainfall measurements during the rain cell period. As can also be seen in the ®gure, there is often an overlap between succeeding rain cells. A.2. Linking of rain cells identi®ed at different rain gauges By the procedure of Appendix A.1, the different rain cells can be separated in the time series of all available rain gauges. Because one spatial rain cell can be measured by many rain gauges, the rain cells in the different time series that belong together have to be identi®ed. The algorithm for such identi®cation is explained in this section. In the algorithm, numbers are given to the different spatial rain cells in each rain storm. This is done in four steps as explained below. Step 1: A priori numbering of rain cells. A ®rst numbering of rain cells is performed for the rain gauge with the largest number of identi®ed rain cells. If more rain gauges have the same maximum number of rain cells, the ®rst numbering is done for the rain gauge with the earliest ®rst rain cell. The motivation for this choice is the higher simplicity to add rain cell numbers at the end of a storm that an the beginning. Afterwards, the rain cells of the other rain gauges are numbered, starting from the rain gauge with the largest number of identi®ed rain cells. For each rain cell, the elapsed times to the peak-moments of all rain cells of the other gauges that have been numbered already are calculated. The number of the
closest rain cell is chosen as the most likely rain cell number. Step 2: Adaptation of double rain cell numbers. When certain rain cell numbers occur twice (or more times) for some rain gauges, the rain cell with the maximum likelihood for having the rain cell number is selected. For this rain cell, the rain cell, the rain cell number is called `known' (indicated here by `K'). For the other rain cells, the number is changed to `unknown' (`U'). For two succeeding rain cells, three criteria are considered to select the most likely rain cell for a certain rain cell number: ² the average time distance between the rain cell center and that of the rain cells with the same number ² the average difference in maximum rainfall intensity between the rain cells and that of the rain cells with the same number ² the average difference in elapsed times between the rain cell center and the centers of the former and next rain cells, on the one hand, and the elapsed times between the rain cell center and the centers of the rain cells with the same number, on the other hand. When looking for rain cells with the same rain cell number in these criteria, only rain cells of gauges that have been adapted already for double rain cell numbers (or with `known' numbers) are considered. Weights are given to the three criteria (e.g. 1/3 min, 2 mm 21 h 21 and 0.2 min 21) and the most likely rain cell for a certain rain cell number is derived by minimization of a weighed average. The largest weight is given to the ®rst criteria, which maximizes the proximity in time of rain cells with the same rain cell number in a rain storm. Whenever more than two succeeding rain cells have the same number, the procedure is repeated by a moving operation for succeeding sets of two rain cells (®rst the two ®rst rain cells, than the second and third rain cell, ¼.). For these succeeding rain cells, the different possibilities are shown in Fig. 18. For the possibility in which the number of the second rain cell is considered unknown for both pairs of two succeeding rain cells (the possibility `¼KUK¼'), the procedure is repeated for the ®rst and third rain cell.
P. Willems / Journal of Hydrology 252 (2001) 126±144
This leads to two possibilities: `¼KUU¼' or `¼UUK¼'. For four succeeding rain cells, the procedure for three succeeding rain cells is repeated twice. For the possibility `¼KUUK¼', the procedure is repeated for the ®rst and the fourth rain cell. It is however a rare possibility as the large elapsed time between the ®rst and the fourth rain cell con¯icts with the minimization of the three criteria. In this way, unique rain cell number are obtained for each rain gauge. Step 3: Check for the minimum proximity of rain cells with the same rain cell number. After the unique assignment of rain cells numbers, the numbers are additionally checked for their proximity in time to cells with the same cell number. For each rain cell, the average elapsed time with rain cells having the same rain cell number is calculated and compared with the time distance to the preceding rain cell (if the average distance is negative) or succeeding rain cell (if the average distance is positive). The average distance with rain cells having the same number should be smaller than the distance with other rain cells. The average distance should also be limited (e.g. smaller than 30 min). For rain cells that do not pass this check, the rain cell number is changed to `unknown'. When certain rain cells numbers do not occur any more after the check for all rain cells, the rain cell numbers are shifted and the number of rain cells in the rain storm are decreased with one. Step 4: Further rain cell number identi®cation for rain cells with unknown number. The possibility of adding rain cell numbers to rain cells with provisionally `unknown' number is additionally tested. Therefore, for each rain cell with unknown number, the `number gap' is calculated as the difference in rain cell number between the two surroundings rain cells with known number. If the number gap equals 0, a number is inserted when the number of rain gauges with gap 0 between the same numbers is suf®ciently high (e.g. 4). If only one rain cell with unknown number exists in the gap, the additional number is connected to this rain cell. If more than one rain cell exist, the additional number is connected to the rain cell with the greatest proximity (see step 3) to the other rain cells with the same `known' number. Therefore, the rain gauges with the lowest number of rain cells in the gap are treated ®rst. More additional
143
numbers are inserted in the same way. The number of added rain cell numbers equals the number of rain cells in the gap, which occurs for a minimum number of rain gauges (e.g. 4). If the gap is 1, the remaining number is simply connected to the single rain cell in the gap. Afterwards, it is checked for its proximity to rain cells with the same number (see step 3). If the gap is higher than 1, the same procedure is repeated. From the different numbers in the gap, the number is chosen on the basis of the three criteria for selecting the most likely rain cell for a certain rain cell number (see step 2). References Amorocho, J., Wu, B., 1977. Mathematical models for the simulation of cyclonic storm sequences and precipitation ®elds. Journal of Hydrology 32, 329±345. Bartolini, P., ValdeÂs, J.B., 1994. Representation of spatial variability of rainfall in aggregated rainfall-runoff models. ASCE Journal of Hydraulic Engineering 120 (10), 1199±1219. Bierkens, M.F.P., Puente, C.E., 1990. Analytically derived models based on rainfall point processes. Water Resources Research 26 (11), 2653±2659. Bras, R., Rodriguez-Iturbe, I., 1976. Rainfall generation: a nonstationary time-varying multidimensional model. Water Resources Research 12 (3), 450±454. Bras, R.L., Tarboton, D.G., Puente, C., 1988. Hydrological sampling: a characterisation in terms of rainfall and basin properties. Journal of Hydrology 102, 113±135. Eagleson, P.S., 1978. Climate, soil, and vegetation, 2. The distribution of annual precipitation derived from observed storm sequences, Water Resources Research 14, 713±721. Gupta, V., Waymire, E., 1979. A stochastic kinematic study of subsynopic space-time rainfall. Water Resources Research 15 (3), 637±644. Jinno, K., Kawamura, A., Berndtsson, R., Larson, M., Niemczynowicz, J., 1993. Real-time rainfall prediction at small space-time scales using a two-dimensional stochastic advection- diffusion model. Water Resources Research 29 (5), 1489±1504. Kawamura, A., Jinno, K., Berndtsson, R., Furukawa, T., 1997. Realtime tracking of convective rainfall properties using a twodimensional advection±diffusion model. Journal of Hydrology 203, 109±118. Luyckx, G., Willems, P., Berlamont, J., 1998. In¯uence of the spatial variability of rainfall on sewer system design. In: Wheater, H., Kirby, C. (Eds.). Hydrology in a Changing Environment, Wheather, vol. 3. Wiley, Chichester, pp. 339±349. Marien, J.L., Vandewiele, G.L., 1986. A point rainfall generator with internal storm structure. Water Resources Research 22 (4), 475±482. Meinhold, R., Singpurwalla, N., 1983. Understanding the Kalman ®lter. The American Statistician 37 (2), 123±127. Northrop, P., 1998. A clustered spatial-temporal model of rainfall.
144
P. Willems / Journal of Hydrology 252 (2001) 126±144
Proceedings of the Royal Society of London, Series A 454, 1875±1888. Sivapalan, M., Wood, E.F., 1986. A multidimensional model of nonstationary space-time rainfall at the catchment scale. Water Resources Research 22 (7), 1289±1299. Taylor, G.I., 1938. The spectrum of turbulence. Proceedings of the Royal Society of London, Series A 164, 476. Velghe, N., Troch, P.A., De Troch, F.P., Van de Velde, J., 1993. Evaluation of cluster-based rectangular pulses point process models for rainfall. Water Resources Research 30 (10), 2847± 2857. Verhoest, N., Troch, P.A., De Troch, F.P., 1997. On the applicability of Bartlett-Lewis rectangular pulses models for calculating design storms at a point. Journal of Hydrology 202 (1-4), 108±120. Waymire, E., Gupta, V., 1981. The mathematical structure of rainfall representations; 1.A review of the stochastic rainfall models. Water Resources Research 17 (5), 1261±1272.
Waymire, E., Gupta, V., Rodriguez-Iturbe, I., 1984. A spectral theory of rainfall intensity at the meso-b scale. Water Resources Research 20 (10), 1453±1465. Willems, P., 1999. Stochastic generation of spatial rainfall for urban areas. Water Science and Technology 39 (9), 23±30. Willems, P., Berlamont, J., 1998. Stochastic modeling of spatial rain cells. In: Wheater, H., Kirby, C. (Eds.). Hydrology in a Changing Environment, vol. 3. Wiley, Chichester, pp. 307± 318. Willems, P., 2000. Compound intensity/duration/frequency-relationships of extreme precipitation for two seasons and two storm types. Journal of Hydrology 233, 189±205. Wilson, C., ValdeÁs, J.B., Rodriguez-Iturbe, I., 1979. On the in¯uence of the spatial distribution of rainfall on storm runoff. Water Resources Research 15 (2), 321±328. Zawadski, I., 1973. Statistical properties of precipitation patterns. Journal of Applied Meteorology 12, 459±472.
www.elsevier.com/locate/jhydrol
A spatial rainfall generator for small spatial scales Patrick Willems* Hydraulics Laboratory, Katholieke Universiteit Leuven, Kasteelpark Arenberg 40, B-3001 Heverlee, Belgium Received 29 March 2000; revised 9 May 2001; accepted 22 May 2001
Abstract A stochastic spatial rainfall generator is developed for use at the small spatial scale of urban and small hydrographic catchments. The generator is based on a spatial rainfall model of the conceptual and hierarchical type. It describes the spatial rainfall ®eld in a macroscopic physically-based way by distinguishing rainfall entities with different scales: rain cells, cell clusters, small and large mesoscale areas (or rain storms). For applications at small spatial scales, the individual rain cells need a detailed description. Data of a dense network of rain gauges at Antwerp, enclosing 5940 rain cells in 807 rain storms are used to derive such description. For separation of the rain cells in the rainfall time series, an algorithm is developed based on the identi®cation of increasing and decreasing rain cell ¯anks. The rain cells observed at different rain gauges are linked together by applying criteria for testing the similarity in rain cell properties. After separating and linking the rain cells and storms, the spatial rainfall model is calibrated to many storms by two methods (e.g. Kalman ®lter). The derived model structure and model parameter distributions apply to the stochastic generation of long-term time series of spatial rainfall. The model is tested by comparing intensity±duration±frequency relationships and temporal scaling properties of the generated and historical rainfall series. q 2001 Elsevier Science B.V. All rights reserved. Keywords: Rain cells; Time correlation; Rain storm
1. Introduction Much research has already been done trying to describe the organization in spatial rainfall ®elds; see for example Bras and Rodriguez-Iturbe (1976), Gupta and Waymire (1979), Waymire and Gupta (1981), Waymire et al. (1984) and Sivapalan and Wood (1987), to name a few. In this research, the `clustering' of high intensity rainfall areas embedded within rainfall areas of lower intensity has been determined as a characteristic pattern. The pattern is observed at many different scales. At the smallest scale, the individual rain cell forms the building* Fax: 132-16-32-18-98. E-mail address: [email protected] (P. Willems).
block of the spatial rainfall structure. The rain cells are embedded in a clustered way within `small mesoscale areas' (10 2 ±10 3 km 2). They move in the same small mesoscale area with nearly identical velocity. At larger scales, small mesoscale areas occur in a clustered way within lower intensity `large mesoscale areas', which in turn are embedded within some synoptic-scale lowest intensity rainfall ®eld. The large mesoscale areas are observed as (sometimes elongated) bands and their spatial extent varies from 10 3 to 10 4 km 2. The precipitation associated with cyclonic and frontal storms are of the level of these large mesoscale areas. Although rain cells and cell clusters most often appear in large and small mesoscale areas, they can also occur isolated outside these regions (e.g. in air mass thunderstorms).
0022-1694/01/$ - see front matter q 2001 Elsevier Science B.V. All rights reserved. PII: S 0022-169 4(01)00446-2
P. Willems / Journal of Hydrology 252 (2001) 126±144
127
Fig. 1. Gaussian-shapped rain cells on the basis of the structure of spatial rainfall (map area of Belgium).
This spatial rainfall structure, which is also represented in the radar view of Fig. 1, has served as the basis for many space-time rainfall modeling applications. Wilson et al. (1979), Bras et al. (1988), Bartolini and ValdeÂs (1994), Bierkens and Puente (1990), Northrop (1998) and Sivapalan and Wood (1987) and many others used such a rainfall model to analyse the in¯uence of rainfall input data on the results of hydrological models. Most of these applications deal with large watersheds. Almost no applications at very small spatial scales exist, such as applications to urban hydrology and small hydrographic catchments. For these applications, the accurate description of individual rain cells and cell clusters is very important, while the description of the spatial rainfall structure at the scale of large mesoscale and synoptic areas is of minor importance. The spatial extent of individual rain cells (10±50 km 2 following Waymire et al., 1984) is indeed of the same order of magnitude as the spatial extent of most urban and small hydrographic catchments. Rain cells, however, have a rather super®cial description in the above mentioned research. There is most often lack of a solid empirical basis for their modeling. This basis should consist of a very dense network of rain gauges. Its availability is however rather scarce. The present study focusses on the accurate spatiotemporal description of rain cells and rain storms at the small spatial scale of urban and small hydrographic catchments. The `solid empirical basis' for
Fig. 2. Rain gauge network of the City of Antwerp, together with the ®ve sub-networks A1,¼,A5.
the study consists of a dense network of 12 tippingbucket rain gauges at the City of Antwerp. The gauges are located at distances of approximately three kilometres and they cover a large urban area of about 100 km 2 (Fig. 2). This area is larger than the maximal spatial extent of rain cells and it is of the same order of magnitude as the spatial extent of cell clusters. It also approximates the area of hydrographic subcatchments and the largest urban catchments. The rain gauges record rainfall intensities every minute with a resolution of 0.1 mm since July 1994. 2. The structure of rain cells On the basis of the data of the Antwerp network, the basic structure of rain cells has been analyzed by Luyckx et al. (1998), Willems and Berlamont (1998) and Willems (1999). The structure includes two aspects: the spatial distribution of rainfall intensities in a rain cell (the shape of the rain cell) and the
128
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 3. Calibration of the average value of the parameters rmax and sy of the rain cell model (the points indicate the twelve rain gauges P1, P2,¼, P12).
movement of the rain cell in time and space (transport of the rainfall intensities). It has been shown by Luyckx et al. (1998) that the spatial distribution of rainfall intensities in a single rain cell can be described well by a bivariate Gaussian distribution: r x; y; t rmax exp 2 2
x 2 x0 2 u t 2 t0 2 2s2x !
y 2 y0 2 2 g t 2 t0 2s2y
1
In this formula, r x; y; t represents the rainfall intensity at spatial co-ordinates (x,y) and time t. The x-direction is taken equal to the direction of rain cell movement u . The other rain cell model parameters consist of: the maximum rainfall intensity rmax (the amplitude of the rain cell), the spatial extent in the x- and y-direction) in terms of the standard deviation of the Gaussian distribution) sx and sy, the mean velocity u and the decay velocity g . x0 ; y0 are the spatial co-ordinates of the rain cell at a reference time t0. To describe the movement of the rain cells, Jinno et al. (1993) have shown that a two-dimensional advection±diffusion model can be used. 2r 2r 22 r 22 r 1u dx 2 1 dy 2 2 gr 2t 2x 2x 2y
2
with dx and dy the diffusion coef®cients in the x- and the y-direction. This advection±diffusion model has already been applied with success by Jinno et al. (1993) and Kawamura et al. (1997) in real-time
control applications. Together with the decay of cell volumes, advection and diffusion are indeed observed as the most important processes in the movement of rainfall intensities. On the other hand, a Gaussian distribution can analytically be determined as a possible solution to the advection±diffusion model r0 p exp r x; y; t 4p dx dy t 2 t0 x 2 x0 2 u t 2 t0 2 4dx t 2 t0 ! y 2 y0 2 2 g t 2 t0 2 4dy t 2 t0 £
2
3
Comparing Eqs. (1) and (3), it can easily be seen that, at time t, the spatial extent of the cell equals p p sx 2dx t 2 t0 ; sy 2dy t 2 t0 4 By diffusion, the cell extent will increase in time and the maximum rainfall intensity of the cell rmax will decrease r0 p rmax 5 4p dx dy t 2 t0 where r0 represents the total rain volume of the cell before decay of the cell volume starts. 3. Calibration of the rain cell model To calibrate the rain cell model, two methods are tested and compared: the semi-automatic (and
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 4. Calibration of the average value of the parameter sx of the rain cell model of optimizing the agreement with a standard Gaussian curve.
statistical) Kalman ®lter technique applied as recursive parameter estimator (Meinhold and Singpurwalla, 1983; Willems and Berlamont, 1998) and a more manual method in which the model parameters are related to physical properties that can be derived from the rainfall time series. In this latter method, the moving-velocity u and the moving-direction u are ®rst calibrated by considering two couples of rain gauges. The moving-times between the two gauges of each couple can be expressed in terms of u and u and equalized to the real moving-times (see also Luyckx et al., 1998). By considering two couples, the two unknowns (u and u ) can be calculated from the two equations. The real moving-times are not known, but they can be estimated from the time between the measured peak intensities. By calculating u and u for each two couples of rain gauges and by averaging the results, the model parameters u and u can be estimated accurately. When the moving-direction is known, the peak intensities at the different rain gauges can be plotted against the distances of the rain gauges in the direction perpendicular to the movingdirection (the y-direction (Fig. 3). From this plot, the average values of the maximum rainfall intensity rmax and the spatial extent in the y-direction sy can be derived. In Fig. 3, y0 is the distance (in the y-direction) of the rain cell peak. The spatial extent in the moving(or x-) direction sx (or the small scale spatial correlation in the direction of storm movement) can be calibrated on the basis of the small scale time correlation as observed in the rainfall times series of the different rain gauges. This time correlation can be transformed to space using the mean storm velocity
129
u. The transformation is based on the hypothesis of ¯uid turbulence by Taylor (1938), which states that the spatial pattern of turbulent motion, as re¯ected in the spatial correlations, is carried past a ®xed point (here the rain gauge) by the mean ¯ow speed without undergoing any essential change. The hypothesis is only valid for uniform mean ¯ow and a low level of turbulence (or small `noise' on the average motion). The usefulness of Taylor's hypothesis for a spacetime rainfall storm was shown earlier by Zawadski (1973), Bras and Rodriguez-Iturbe (1976) and Waymire et al. (1984). The following transformation formula is therefore applied: h
xi 2 x0 2 u t 2 t0 2 y 2 y 2 2 i 20 2sx 2sy
6
In this equations, x0 is the position (in the x-direction) of the rain cell peak at time t0; xi ; yi are the spatial co-ordinates of the rain gauge i. The value x0 can be calibrated by maximizing the agreement between the measured time moments of the peak rainfall intensities and the modelled ones ti t0 1 u xi 2 x0
7
In Fig. 4, the rainfall intensities at the different rain gauges i i 1; ¼; 12 and at different moments t are plotted against the relative distance h to the rain cell centre. For the calibrated value of sx and for Gaussianshaped rain cells, Fig. 4 approximates the standard Gaussian curve r exp 2h 8 rmax By scaling the rainfall intensities to the maximum rain cell intensity rmax ; the maximum of this curve equals 1. After testing and comparing the two calibration techniques (Willems, 1999), a combination of the techniques is recommended: application of the manual method ®rst to derive `a priori' estimates of the model parameters for the Kalman ®lter. The Kalman ®lter technique has major advantages because it takes into consideration rainfall measurement errors and uncertainties in the rain cell model structure. As rainfall measurement errors depend on the rainfall intensity, this consideration of the errors is important for the derivation of unbiased parameter values. In the Kalman ®lter technique, the model parameters are continuously updated in time (as
130
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 6. Results of the linking of rain cells at different rain gauges for the rain storm on 14.11.94 (identi®ed rain cells with signi®cant volumes are represented by bubbles, bubbles with the same ®lled color are identi®ed as the same rain cell, white bubbles represent rain cells that cannot be identi®ed).
4. Rain cell separation Fig. 5. Results of the linking of rain cells at different rain gauges for the rain storm on 14.11.94 (identi®ed rain cells are represented by triangles, triangles with the same ®lled color are identi®ed as the same rain cell).
more measurements become available) and a statistical estimate is made of their uncertainties. In this way, the variability in time of model parameters can be recognized. By the automatic nature of this calibration method, calibration times are very small. An important disadvantage of the method is the dif®cult control of the physical realism of the parameter values during the calibration procedure. By using the results of the manual method as initial values for the model parameters, this problem is however often solved. This consideration can be made since the experience of the modeller is used in the manual method. The manual calibration procedure is however more time-consuming, although it can be partly automated by different operations in a spreadsheet computer programme. To reduce the in¯uence of irregularities in the measurements (e.g. measurement noise) on the calibration results, the original 1 min time series are processed ®rst by moving average over 10 min intervals. The smoothing of 1 min peak intensities, caused by this processing, is of no relevance to hydrology because a much higher smoothing is described in the rainfall-runoff transfer of the smallest hydrological system.
For rain storms with more than one rain cell, the different rain cells have to be identi®ed and separated in the rainfall time series. An algorithm is derived for this purpose. It is described in detail in the appendices and consists of two parts. In the ®rst part (Appendix A), the different rain cells in the point rainfall time series are separated for each rain gauge independently. This separation is based on the identi®cation of peaks, decreasing ¯anks and increasing ¯anks in the time series. In the second part (Appendix B), the identi®ed rain cells for the different rain gauges are linked together. This means that the algorithm identi®es rain cells in the different time series, which correspond with the same spatial cell. As a result of the rain cell separation and identi®cation process, the different spatial rain cells in each rain storm are numbered and some physical rain cell properties are derived. In Fig. 5, the results of the algorithm are shown for the rain storm on 14.11.94. Identi®ed rain cells are represented in this ®gure as triangles. Rain cells with the same color are linked together. The cells with insigni®cant volume (volume smaller than 0.4 mm) are not considered in the linking procedure. They are represented with white color. As can be seen in the ®gure, the linking procedure performs well. In Fig. 6, the same results are represented in a different way. The identi®ed rain cells with signi®cant volume are represented by bubbles for which the magnitude is proportional to the cell volume. The non-colored cells are the ones that could not be
P. Willems / Journal of Hydrology 252 (2001) 126±144
131
Fig. 7. Calibration results for three rain cells in the rain storm on 28.08.96.
linked. They are considered as unidenti®able. The result of the rain cell separation and linking procedure are used for two purposes. They are ®rstly used for the identi®cation of subperiods in the rainfall time series for which the individual rain cell model has to be calibrated. Secondly, they are used for counting the number of rain cells in rain storms, as explained in Section 6.
5. Probability distributions of rain cell parameters According to the procedures of Sections 3 and 4, about 104 rain cells in 58 rain storms are calibrated. They are randomly selected in the period October 1994±February 1997. Fig. 7 shows an example of
132
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 8. Calibration results of the probability distributions for the rain cell model parameters rmax, s, u and u .
the calibration results for three rain cells in the rain storm on 28.08.96. After calibration of the rain cell model for a large number of rain cells, statistical properties of the model parameters are derived in terms of probability distributions and correlations. For the rain cell parameters, clear distribution functions are determined, as shown in Fig. 8 for the parameters rmax ; s, u and u . The derived distributions and correlation coef®cients of the different model parameters are summarized in Tables 2 and 3. As mentioned in Section 2, the spatial extent of a rain cell sx and sy is not constant during its movement over a spatial area. It was explained by the diffusion process. For the detailed description of rain cells, as needed in applications of real time prediction, a detailed consideration of the diffusion and decay processes is important. However, for purposes of stochastic spatial rainfall generation, this consideration is less important. The time the rain cells need to move over an urban or small hydrographic catchment is rather small and the inclusion of the diffusion and decay processes will only change the statistical properties of the generated spatial rainfall in an insignificant way. Therefore, for each rain cell the average
spatial extent during its movement over the rain gauge network is calibrated and the decay velocity is neglected. The spatial extent in both x- and y-directions is not signi®cantly different. Therefore, also a combined distribution is derived for sx and sy : The mean spatial extent of a rain cell in both directions equals 2.5 km (in terms of the standard deviation of the Gaussian shape). If the diameter of a rain cell is considered as six times the standard deviation, the mean diameter equals 15 km. This diameter encircles almost the total rainfall volume (99.7%) of the rain cell. 6. The cell cluster model The spatial organization of the rain cell centers in small mesoscale areas or cell clusters is studied using the numbering results as explained in Section 4. From the temporal density of rain cells in rain storms (or the average number of rain cells per time unit), the spatial density can be calculated using the moving velocity and direction of the rain storm. Indeed, the spatial extent (in the moving direction u ) of the rain storm
P. Willems / Journal of Hydrology 252 (2001) 126±144
133
Fig. 9. Schematic representation of the cell cluster model (rain storm location shown for time moment t t0 ).
ls u can be calculated from the temporal duration ds and the moving velocity u of the rain storm ls uds 2 la u
9
la u is in this formula the spatial extent of the area covered by the rain gaue network in the direction u (see Fig. 9). The rain storm duration ds is de®ned by the time span from the beginning of the ®rst rain cell (the earliest beginning for all rain gauges) to the ending of the last rain cell (the latest ending for all rain gauges). It is clear that the time the rain storm needs to move over the network area la u=u should be subtracted from the storm duration ds to derive the real rain storm length ls : The rain storm width (the
extent of the rain storm in the direction perpendicular to the moving direction) is assumed larger than the extent of the area covered by the rain gauge network in the same direction wa u: To study the organization of the rain cell centres, the number of rain cells is then analysed in relation with the spatial area V s over which the cells as observed by the rain gauge network are spread (see also Fig. 9)
V s ls wa u 1 2w
10
It is assumed in this formula that all rain cells with a cell centre closer than a distance w from the most remote rain gauges are observed by the rain gauge network. The distance w is taken equal to the average
Fig. 10. Analysis of the number of rain cells per unit area in rain storms.
134
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 11. Cumulative distribution function of the rain storm duration ds (comparison of the empirical distribution and the distribution simulated from an exponential rainstorm length distribution).
spatial radius of a single rain cell (7.5 km, see Section 5). For the most simple structure of a spatial rain cell process, being the simple spatial Poisson process, the only parameter of the process (l ) can be calculated as the average number of rain cells per unit area. In Fig. 10, the observed number of rain cells is for all calibrated rain storms plotted against the spatial area of the rain gauge network V s. The latter area is calculated by Eq. (10). The slope in the ®gure equals the parameter l
l
average number of rain cells area V s
11
To derive the ®gure for a wide range of areas, ®ve sub-networks of the Antwerp rain gauge network are considered. They are indicated in Fig. 2 as A1, A2, A3, A4 and A5. Sub-network A1, for instance, only contains rain gauge no. 6, while sub-network A2 also includes the gauges 3, 4, 5, 7 and 8. The ®ve sub-networks have their own values for la u and wa u; which are calculated from rain gauge map coordinates. From Fig. 10, an average value of l 0:002 is derived. There is however a large scattering in the ®gure. This may indicate that other more complicated spatial process descriptions are needed. For the current study, however, a simple Poisson process is considered. This process is assumed very often for the stochastic description of sequences of rain cells in rain storms both in time and space. This is the case, for instance, in the cluster-based rectangular pulses models (e.g. Neyman±Scott model,
Bartlett±Lewis model), which have been calibrated earlier for temporal Belgian rainfall data by Velghe et al. (1993) and Verhoest et al. (1997). It is clear that the larger the area of the rain gauge network, the larger the maximum extent of the spatial organization of rain cells that can be considered in the analysis. For the Antwerp rain gauge network, it can be studied well up to the spatial level of small mesoscale areas. To ®nalize the cell cluster model, it is assumed that all rain cells in the same cell cluster move with the same velocity and direction. It is thus assumed that the cells' velocities relative to the storm velocity are essentially zero. This assumption is based on the early work of Zawadski (1973) and, after calibration of many rain storms to the Antwerp data, it is seen to be a good assumption. 7. Probability distributions of cell cluster parameters Also for the other cell cluster parameters, such as the inter-arrival time of rain storms and the rain storm duration, the information is derived from the rain cell separation and linking results as explained in Section 4. All rain cells and rain storms in the period under study (August 1994±May 1997) can be included in this analysis. They consist of 5940 rain cells in 807 rain storms. For each rain storm, the following information is derived: the inter-arrival time, the storm duration and the number of rain cells.
P. Willems / Journal of Hydrology 252 (2001) 126±144
135
Fig. 12. Exponential quantile plot for the inter-arrival time ta of rain storms.
The length of the cell cluster areas (or small mesoscale areas) can only be calculated for the rain storms with known cell movement direction. Because this information is only available for the 104 calibrated rain cells, a different approach is used for calibration of the distribution of ls : On the basis of an assumed distribution for ls and the known Weibull and normal distributions for u and u , the distribution of the rain storm duration ds is simulated on the basis of Eq. (9). A comparison with the empirical distribution of ds allows to determine the best distribution for ls : An exponential distribution is assumed by many authors for the rain storm lengths. This distribution is therefore tested ®rst. It is seen in Fig. 11 that for an exponentially distributed ls with a mean storm length of 1=l 35:7 km; the simulated distribution of ds matches well the empirical distribution. The range of rain storm lengths described by this exponential distribution corresponds with the range indicated by Waymire et al. (1984): 10±10 3 km 2 for cell clusters or small mesoscale areas. For the inter-arrival time of rain storms, the probability distribution is calibrated directly to the empirical distribution. A two-component exponential distribution is determined (Fig. 12) ta days 2 0:0354 F ta 0:15 1 2 exp 2 5:5 t days 2 0:0354 1 0:85 1 2 exp 2 a 0:5 12
The two-components show the existence of two populations of exponentially distributed inter-arrival times. For each population, the exponential distribution indicates that the rain event arrivals correspond to a temporal Poisson process. The second population with small inter-arrival times (mean inter-arrival time of 0.5 days) corresponds to inter-arrival times between different rain storms (or small meso-scale areas) in the same large meso-scale area (e.g. frontal storm). They vary from 50 min to many hours. The ®rst population with large inter-arrival times (mean interarrival time of 5.5 days) corresponds to inter-arrival times between different large meso-scale or synoptic areas. They may have values up to 30 days (see Fig. 12). The two-component distribution for ta can probably also be derived when the spatial organization of small meso-scale areas is described by a spatial Poisson process (in the same way as the spatial organization of rain cells in small meso-scale areas is described). However, because of the limited spatial extent of the Antwerp rain gauge network, this spatial clustering cannot be studied for spatial scales larger than small meso-scale areas. The in¯uence of the spatial organization of small meso-scale areas is therefore modeled in an indirect way by describing the temporal clustering of the smaller inter-arrival times, which is a result of this spatial organization. This temporal clustering of the second population of inter-arrival times could be described by an exponential distribution for the duration of a cluster of small
136
P. Willems / Journal of Hydrology 252 (2001) 126±144
Table 1 Overview of the probability distribution function for the model parameters Model
Distribution
Distribution parameter values
Mean
u (km/min) u (degrees) sx (km) sy (km) s (sx and sy together) rmax (mm/h) ls (km)
Weibull Normal Lognormal Lognormal Lognormal Generalized Pareto (GPD) Exponential
t Weibull index 0:70; b 2:36 m u 221:75; s u 68:46 m ln s x 0:90; s ln s x 0:41 m ln sy 0:92; s ln s y 0:46 m ln s 0:91; s ln s 0:44 g Extreme value index 0:396; b 1:65; r t 0:6 l 0:028
0.6 2 21.75 2.5 2.5 2.5
and properties are derived for the combined spatiotemporal rain storm process at the small spatial scale. They consist of the rain cell and cell cluster model structures (Fig. 9, Eq. (1)), the probability distributions of the rain cell parameters and their correlations (Tables 1 and 2, Fig. 8), the l parameter of the spatial Poisson process of rain cells in rain storms and the statistical models for the rain storm duration and the inter-arrival time of rain storms (Eq. (12)). On the basis of these properties, a spatial rainfall generator can be easily constructed. Using a random number generator, a large number of spatial rain storms can be randomly simulated. The long-term time series of spatial rainfall that are generated in this way have the same statistical properties compared with those derived from the rain gauge network data. The simulated long-term time series of spatiotemporal rainfall can be used for many applications. Two examples are given hereafter. The time series can be used as input to hydrological models to study the in¯uence of the spatial variability of rainfall on the model results. They also can be used to simulate spatially correlated rainfall simultaneously for different subsystems in an integrated catchment model, such as the sewer system and the hydrographic catchment. Due to the limitation of the spatial scale, the
meso-scale areas (mean duration of 3.0 days), combined with a Poisson process for the description of the different small inter-arrival times inside such a cluster. In Fig. 12, a comparison is shown for the inter-arrival times simulated by this model and the ones derived from the Antwerp time series by the rain cell identi®cation procedure. The model shows a small underestimation for the intermediate inter-arrival times. Marien and Vandewiele (1986) also determined a compound distribution for the inter-arrival times of rain storms in a point rainfall series. They considered three components, one of them corresponding with the inter-arival times between different rain cells in a rain storm (modeled here in an indirect way by the spatial Poisson process). For the two other components, an exponential and a Weibull distribution were considered. Exponential and Weibull distributions were also used by Eagleson (1978) and Amorocho and Wu (1977) for the inter-arrival times for large mesoscale areas. 8. Spatial rainfall generator In the previous sections, basic conceptual features Table 2 Correlation matrix for the model parameters
u (km/min) u (degrees) sx (km) sy (km) rmax (m/h)
35.7
u (km/min)
u (degrees)
sx (km)
sy (km)
rmax (mm/h)
1 ± ± ± ±
0.15 1 ± ± ±
0.86 0.15 1 ± ±
0.13 0.006 0.35 1 ±
2 0.25 0.12 2 0.22 ±0.13 1
P. Willems / Journal of Hydrology 252 (2001) 126±144
137
Fig. 13. Model validation by comparison of IDF-relationships. (a) for the ®nal rainfall generator. (b) for a rectangular rain cell model with and without clustering of small inter-arrival times.
generator is, of course, only applicable to sewer system catchments and hydrological cachments of limited spatial extent. 9. Model testing by IDF-relationships The accuracy of the spatial rainfall generator is tested by comparing the statistical properties of the stochastically generated long-term time series of spatial rainfall and these of historical rainfall series.
The comparison is made for different statistical properties: extreme value distributions, IDF-relationships, autocorrelations and scale properties. Apart from the Antwerp data on the basis of which the rainfall generator is built and calibrated, also the Ukkel/ Uccle series for 27 years (1967±1993) is used as historical series. Through the greater time-span of this latter series, the accuracy of model extrapolations to rare events (which are not available in the spatial rainfall data at Antwerp) can be tested. The extreme value analysis of the generated and historical time
138
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 14. Model validation by comparison of the extreme value distribution (example of an exponential quantile plot for an aggregation-level of 30 min).
series and the calculation of IDF-relationships is performed according to the methodology described in Willems (2000). In Fig. 13(a), the comparison is shown for intensity±duration±frequency-relationships up to 27 years, which were derived for the Uccle series in Willems (2000). The con®dence limits are indicating the differences in IDF values between the 12 individual rain gauges of the Antwerp network. Good results are shown, both for rainfall frequency and for a large range of temporal scales. Even the temporal scaling properties, indicated by the linear shape of the IDFcurves in Fig. 13, are well represented. Apart from this good general performace of the generator, two minor differences can be observed. Asymptotically to the smallest aggregation-levels (smaller than 20 min), the IDF-relationships show a higher bending down for the rainfall generator. This is explained by the full deterministic nature of the rain cell model. By adding a random component in terms of high frequency noise to this model, the bending down can be corrected. Another difference is the small underestimation of the rainfall volumes at large aggregation-levels (larger than one day). A sensitivity analysis has shown that this underestimation is very sensitive to the model for the smaller inter-arrival
times. If the short-term dependency in inter-arrival times is totally ignored in the model, the underestimation will increase with a factor three. This is shown in Fig. 13(b) by randomly changing the time order of the different rain cells that were identi®ed in the original Antwerp rainfall time series. Also by ignoring the smaller rain cells in the model, the underestimation of rainfall volumes for large aggregation-levels increases. For instance, by ignoring all rain cells with a volume less than 0.4 mm, the decrease in rainfall intensities (Dr) given in relation to the aggregation-level da by the following formula: Dr mm=h 0:02 1
0:3 da days
13
In addition to the IDF-relationships, also the individual distributions of rainfall intensities at the different aggregation-levels are compared. An example is shown in Fig. 14 for an aggregation-level of 30 min. The shape of the distribution's tail (which is linear for an exponential distribution in the exponential quantile plot of Fig. 14) is well reproduced by the generator. This is of primary importance for an accurate extrapolation to extreme events.
P. Willems / Journal of Hydrology 252 (2001) 126±144
139
10. Conclusions A spatial rainfall model is derived for use in hydrological applications at small spatial scales. The model is of the conceptual and hierarchical type and its structure is similar to many existing models for larger meso-scale areas. It differs from these models by the more detailed description of the smallest buildingblocks of a spatial rainfall structure: the individual rain cells and the small scale rain cell clusters (or small meso-scale area). The description is based on a detailed analysis of the observed cell cluster patterns by a dense network of rain gauges at the city of Antwerp. Another feature is the development of an algorithm to identify and separate different rain cells in the rain gauge time series and to link together identi®ed rain cells for different rain gauges. The rain cell separation is based on the identi®cation of increasing and decreasing rain cell ¯anks, while the rain cell linking is based on criteria for testing the similarity in rain cell properties. By extending the deterministic model with a stochastic description of model parameters, the model is very useful for the long-term stochastic generation of small scale spatial rainfall. It is shown that the most important statistical properties of rainfall, like temporal scaling properties and extreme value statistics, are well reproduced by the generator. Other properties, such as spatial coverage statistics and spatial scaling properties still have to be validated. It is however expected that also these properties are described reasonably well due to the physically-based nature of the approach that was followed for the structure identi®cation and the calibration. This approach differs from other approaches where the parameters of an assumed model-structure are calibrated by matching global statistical properties such as spatial correlations and spatial coverage statistics. Because of its macroscopic physically-based structure and the applicability of the recursive Kalman ®lter calibration procedure, the model described in this paper may also be useful for short-term rainfall prediction in real-time control applications. Acknowledgements This research was supported by a scholarship of the
Fig. 15. Criteria for associating states to the rainfall intensities at the different 1 min time steps in the rain cell separation procedure.
Flemish Institute for Scienti®c and Technological Research (IWT) and a project for the Programme `Beleidsgericht Onderzoek' of the Ministry of the Flemish Community. The Antwerp rain gauge network data was made available by the `City of Antwerp ± Department voor Werken', while the long-term historical rainfall data at Uccle was processed and delivered by the Royal Meteorological Institute of Belgium (KMI/IRM). The author is also grateful to two reviewers which have made important suggestions for improvement of this paper. Appendix A. Algorithm for rain cell separation and linking In this appendix, a description can be found of the algorithm developed and applied in section 4 for the identi®cation and separation of different rain cells in the rainfall time series. The procedure is described ®rst for a single point rainfall series (Appendix A.1). In Appendix A.2, the more advanced procedure of linking of rain cells observed at different rain gauges is explained.
140
P. Willems / Journal of Hydrology 252 (2001) 126±144
² The rainfall intensity is a minimum (state 0) ² The rainfall intensity is part of an increasing ¯ank (state 1) ² The rainfall intensity is part of a decreasing ¯ank (state-1).
Fig. 16. Criteria of Marien and Vandewiele (1986) for de®ning minima in the 10 min rainfall time series (to de®ne the end of partial storms).
A.1. Rain cell separation in a point rainfall time series The algorithm for the separation of different rain cells in a point rainfall time series is based on the identi®cation of peaks, decreasing ¯anks and increasing ¯anks in the time series. The identi®cation is performed by comparing the rainfall intensity at each time step in the rainfall series with the rainfall intensities at the neighbouring time steps. In this way, the rainfall intensity at each time step can be associated with one of the three following possible conditions:
These possibilities are hereafter called the `states' of the rainfall intensities and are given the number 0,1 and 21.When the identi®cation of a rainfall intensity is not clear, it is given `state 3'. The criteria for discriminating among the four state conditions are summarized in both a graphical and numerical way in Fig. 15 for minima and increasing ¯anks. For decreasing ¯anks, the criteria are analogous to the ones of increasing ¯anks. A rainfall intensity of zero is always given a state number 0. The criteria in Fig. 15 are based on a comparison of the rainfall intensities at two preceding time steps r 21 and r 22 and two succeeding time steps (r(1) and r(2)). In the ®rst two criteria for increasing ¯anks, local minima in increasing ¯anks are identi®ed. They are eliminated from the time series by replacement of the rainfall intensity at the local minimum by the average rainfall intensity for the previous and the next time step. The criteria are based on the work of Marien and Vandewiele (1986). They identi®ed minima in a 10 min point rainfall time series to separate different rain cells (they call it `partial storms'). The criteria that they used to select these minima are summarized in Fig.16 in the form of Fig. 15. All minima in the time series that do not meet these criteria were considered
Fig. 17. Results of the rain cell separation procedure for the rain storm on 14.11.94 (rainfall time series of rain gauge 12).
P. Willems / Journal of Hydrology 252 (2001) 126±144
as local minima in increasing or decreasing ¯anks. Different criteria are used in this study because of the different time step of the time series (1 min instead of 10 min). From the list of state numbers at all time steps, the different rain cells are identi®ed together with the time steps indicating the beginning, the ending and the peak-moment of the rain cells. This is also performed by an algorithm. The following rules are used: ² the beginning of a rain cell is de®ned by the ®rst state 1 after a list of states 0 or 3 (e.g. `¼033331¼' or `¼01¼'; the bold position here indicates the starting position of a rain cell). ² the peak-moment of a rain cell is de®ned by a state between an increasing and a decreasing ¯ank (e.g. `¼1133311-1-1¼'). An increasing ¯ank is de®ned here by a list of states 1 and 3, while a decreasing ¯ank corresponds to a list of states 21 and 3. Whenever the increasing and decreasing ¯anks are separated by a list of states 3, the peakmoment is not de®ned that clear (e.g. `¼113333113333-1-1¼'). The peak-moment is then selected as one of the states 3 more speci®cally the one that is closest to the middle between the starting and ending moments of the rain cell. ² the ending of a rain cell is de®ned by a state 0 after a decreasing ¯ank (e.g. `¼3-133-1-10¼'). Whenever states 3 are listened in between the decreasing ¯ank and the state 0, the ending is considered at the ®rst state 3 (e.g. `¼3-133-1-133330¼'). Such list of states 3 indicates that the ending of the rain cell is not that clear. However, if this list would be considered part of the decreasing ¯ank, the slope of the ¯ank would be too `weak'. For the illustration of such a weak decreasing ¯ank, the reader can look at the third identi®ed rain cell in Fig. 17. The weakness in this ®gure is however not caused by the above-mentioned identi®cation problem. It is just caused by the length of the ¯ank. Sometimes, when two succeeding rain cell are close together, no state 0 exist between the two rain cells. The ending of the ®rst rain cell and the beginning of the second rain cell are then both de®ned by the last state 21 in the decreasing ¯ank (or the ®rst state 3 if both rain cells are separated only by a list of states 3, e.g. `¼3-133-1-1333311131¼').
141
² To avoid that a short and insigni®cant increasing ¯ank in a longer decreasing ¯ank (or local noise in a decreasing ¯ank) is identi®ed as the increasing ¯ank of an additional rain cell, a minimum value (e.g. 3 min) is required for the length of an increasing ¯ank. The following list of states is thus considered as a continuously decreasing ¯ank: `¼3-133-1-11-13-1¼' (the single state 1 (in bold) is ignored). ² Sometimes, there is an overlap in time between two succeeding rain cells. As the rainfall intensity measurements are representing the sum of the rainfall intensities of two rain cells during the overlap period, the starting time of the ®rst rain cell and ending times are not de®ned correctly by the above listed rule. Indeed, the starting time of the ®rst rain cell and the ending time of the second rain cell are assumed equal to the time at which the total rainfall intensity approximately reaches a minimum. To asses the real starting and ending times, the ending time of the ®rst rain cell is repeatedly increased with one time step and the starting time of the second rain cell decreased in the same way until the minimum measured rainfall intensity in the overlap period is reached by the model. Instead of using the Gaussian model for this purpose, a simpli®ed triangular model (see later on) is used. The procedure is applied for each pair of rain cells with less than 10 min between the starting and ending times. To achieve a better separation of rain cell, the 1 min rainfall time series (already smoothed by moving average over intervals of 10 min) is ®rst additionally smoothed by moving average over intervals of 5 min. This is done because the `noise' on the rainfall intensities was still too high, leading up to many wrongidenti®ed rainfall intensities (many local maxima not corresponding with rain cell peaks). As can be seen in Fig. 17, the peak intensities are ¯attened a little by the additional moving-average operation. However, the identi®cation of minima, increasing and decreasing ¯anks is easier and is not disturbed by this operation. In Fig. 17 the results of the state identi®cation and the rain cell separation are shown for the rain storm of 14.11.94. Only for graphical reasons, the rain cells are represented by triangular shapes. The time steps at the beginning, ending and peak-moments of the triangles
142
P. Willems / Journal of Hydrology 252 (2001) 126±144
Fig. 18. Different possibilities for the selection of the most likely rain cell with a certain rain cell number (indicated as known) in a list of three succeeding rain cells with the same rain cell number. The rain cells indicated in `bold' are at each step subject to the criteria for selection of the most likely rain cell from two rain cells with the same number.
are derived from the identi®ed ones of the separated rain cells. The height of the triangle-peak is calculated by equalizing the total rainfall volume in the triangle and the original rainfall measurements during the rain cell period. As can also be seen in the ®gure, there is often an overlap between succeeding rain cells. A.2. Linking of rain cells identi®ed at different rain gauges By the procedure of Appendix A.1, the different rain cells can be separated in the time series of all available rain gauges. Because one spatial rain cell can be measured by many rain gauges, the rain cells in the different time series that belong together have to be identi®ed. The algorithm for such identi®cation is explained in this section. In the algorithm, numbers are given to the different spatial rain cells in each rain storm. This is done in four steps as explained below. Step 1: A priori numbering of rain cells. A ®rst numbering of rain cells is performed for the rain gauge with the largest number of identi®ed rain cells. If more rain gauges have the same maximum number of rain cells, the ®rst numbering is done for the rain gauge with the earliest ®rst rain cell. The motivation for this choice is the higher simplicity to add rain cell numbers at the end of a storm that an the beginning. Afterwards, the rain cells of the other rain gauges are numbered, starting from the rain gauge with the largest number of identi®ed rain cells. For each rain cell, the elapsed times to the peak-moments of all rain cells of the other gauges that have been numbered already are calculated. The number of the
closest rain cell is chosen as the most likely rain cell number. Step 2: Adaptation of double rain cell numbers. When certain rain cell numbers occur twice (or more times) for some rain gauges, the rain cell with the maximum likelihood for having the rain cell number is selected. For this rain cell, the rain cell, the rain cell number is called `known' (indicated here by `K'). For the other rain cells, the number is changed to `unknown' (`U'). For two succeeding rain cells, three criteria are considered to select the most likely rain cell for a certain rain cell number: ² the average time distance between the rain cell center and that of the rain cells with the same number ² the average difference in maximum rainfall intensity between the rain cells and that of the rain cells with the same number ² the average difference in elapsed times between the rain cell center and the centers of the former and next rain cells, on the one hand, and the elapsed times between the rain cell center and the centers of the rain cells with the same number, on the other hand. When looking for rain cells with the same rain cell number in these criteria, only rain cells of gauges that have been adapted already for double rain cell numbers (or with `known' numbers) are considered. Weights are given to the three criteria (e.g. 1/3 min, 2 mm 21 h 21 and 0.2 min 21) and the most likely rain cell for a certain rain cell number is derived by minimization of a weighed average. The largest weight is given to the ®rst criteria, which maximizes the proximity in time of rain cells with the same rain cell number in a rain storm. Whenever more than two succeeding rain cells have the same number, the procedure is repeated by a moving operation for succeeding sets of two rain cells (®rst the two ®rst rain cells, than the second and third rain cell, ¼.). For these succeeding rain cells, the different possibilities are shown in Fig. 18. For the possibility in which the number of the second rain cell is considered unknown for both pairs of two succeeding rain cells (the possibility `¼KUK¼'), the procedure is repeated for the ®rst and third rain cell.
P. Willems / Journal of Hydrology 252 (2001) 126±144
This leads to two possibilities: `¼KUU¼' or `¼UUK¼'. For four succeeding rain cells, the procedure for three succeeding rain cells is repeated twice. For the possibility `¼KUUK¼', the procedure is repeated for the ®rst and the fourth rain cell. It is however a rare possibility as the large elapsed time between the ®rst and the fourth rain cell con¯icts with the minimization of the three criteria. In this way, unique rain cell number are obtained for each rain gauge. Step 3: Check for the minimum proximity of rain cells with the same rain cell number. After the unique assignment of rain cells numbers, the numbers are additionally checked for their proximity in time to cells with the same cell number. For each rain cell, the average elapsed time with rain cells having the same rain cell number is calculated and compared with the time distance to the preceding rain cell (if the average distance is negative) or succeeding rain cell (if the average distance is positive). The average distance with rain cells having the same number should be smaller than the distance with other rain cells. The average distance should also be limited (e.g. smaller than 30 min). For rain cells that do not pass this check, the rain cell number is changed to `unknown'. When certain rain cells numbers do not occur any more after the check for all rain cells, the rain cell numbers are shifted and the number of rain cells in the rain storm are decreased with one. Step 4: Further rain cell number identi®cation for rain cells with unknown number. The possibility of adding rain cell numbers to rain cells with provisionally `unknown' number is additionally tested. Therefore, for each rain cell with unknown number, the `number gap' is calculated as the difference in rain cell number between the two surroundings rain cells with known number. If the number gap equals 0, a number is inserted when the number of rain gauges with gap 0 between the same numbers is suf®ciently high (e.g. 4). If only one rain cell with unknown number exists in the gap, the additional number is connected to this rain cell. If more than one rain cell exist, the additional number is connected to the rain cell with the greatest proximity (see step 3) to the other rain cells with the same `known' number. Therefore, the rain gauges with the lowest number of rain cells in the gap are treated ®rst. More additional
143
numbers are inserted in the same way. The number of added rain cell numbers equals the number of rain cells in the gap, which occurs for a minimum number of rain gauges (e.g. 4). If the gap is 1, the remaining number is simply connected to the single rain cell in the gap. Afterwards, it is checked for its proximity to rain cells with the same number (see step 3). If the gap is higher than 1, the same procedure is repeated. From the different numbers in the gap, the number is chosen on the basis of the three criteria for selecting the most likely rain cell for a certain rain cell number (see step 2). References Amorocho, J., Wu, B., 1977. Mathematical models for the simulation of cyclonic storm sequences and precipitation ®elds. Journal of Hydrology 32, 329±345. Bartolini, P., ValdeÂs, J.B., 1994. Representation of spatial variability of rainfall in aggregated rainfall-runoff models. ASCE Journal of Hydraulic Engineering 120 (10), 1199±1219. Bierkens, M.F.P., Puente, C.E., 1990. Analytically derived models based on rainfall point processes. Water Resources Research 26 (11), 2653±2659. Bras, R., Rodriguez-Iturbe, I., 1976. Rainfall generation: a nonstationary time-varying multidimensional model. Water Resources Research 12 (3), 450±454. Bras, R.L., Tarboton, D.G., Puente, C., 1988. Hydrological sampling: a characterisation in terms of rainfall and basin properties. Journal of Hydrology 102, 113±135. Eagleson, P.S., 1978. Climate, soil, and vegetation, 2. The distribution of annual precipitation derived from observed storm sequences, Water Resources Research 14, 713±721. Gupta, V., Waymire, E., 1979. A stochastic kinematic study of subsynopic space-time rainfall. Water Resources Research 15 (3), 637±644. Jinno, K., Kawamura, A., Berndtsson, R., Larson, M., Niemczynowicz, J., 1993. Real-time rainfall prediction at small space-time scales using a two-dimensional stochastic advection- diffusion model. Water Resources Research 29 (5), 1489±1504. Kawamura, A., Jinno, K., Berndtsson, R., Furukawa, T., 1997. Realtime tracking of convective rainfall properties using a twodimensional advection±diffusion model. Journal of Hydrology 203, 109±118. Luyckx, G., Willems, P., Berlamont, J., 1998. In¯uence of the spatial variability of rainfall on sewer system design. In: Wheater, H., Kirby, C. (Eds.). Hydrology in a Changing Environment, Wheather, vol. 3. Wiley, Chichester, pp. 339±349. Marien, J.L., Vandewiele, G.L., 1986. A point rainfall generator with internal storm structure. Water Resources Research 22 (4), 475±482. Meinhold, R., Singpurwalla, N., 1983. Understanding the Kalman ®lter. The American Statistician 37 (2), 123±127. Northrop, P., 1998. A clustered spatial-temporal model of rainfall.
144
P. Willems / Journal of Hydrology 252 (2001) 126±144
Proceedings of the Royal Society of London, Series A 454, 1875±1888. Sivapalan, M., Wood, E.F., 1986. A multidimensional model of nonstationary space-time rainfall at the catchment scale. Water Resources Research 22 (7), 1289±1299. Taylor, G.I., 1938. The spectrum of turbulence. Proceedings of the Royal Society of London, Series A 164, 476. Velghe, N., Troch, P.A., De Troch, F.P., Van de Velde, J., 1993. Evaluation of cluster-based rectangular pulses point process models for rainfall. Water Resources Research 30 (10), 2847± 2857. Verhoest, N., Troch, P.A., De Troch, F.P., 1997. On the applicability of Bartlett-Lewis rectangular pulses models for calculating design storms at a point. Journal of Hydrology 202 (1-4), 108±120. Waymire, E., Gupta, V., 1981. The mathematical structure of rainfall representations; 1.A review of the stochastic rainfall models. Water Resources Research 17 (5), 1261±1272.
Waymire, E., Gupta, V., Rodriguez-Iturbe, I., 1984. A spectral theory of rainfall intensity at the meso-b scale. Water Resources Research 20 (10), 1453±1465. Willems, P., 1999. Stochastic generation of spatial rainfall for urban areas. Water Science and Technology 39 (9), 23±30. Willems, P., Berlamont, J., 1998. Stochastic modeling of spatial rain cells. In: Wheater, H., Kirby, C. (Eds.). Hydrology in a Changing Environment, vol. 3. Wiley, Chichester, pp. 307± 318. Willems, P., 2000. Compound intensity/duration/frequency-relationships of extreme precipitation for two seasons and two storm types. Journal of Hydrology 233, 189±205. Wilson, C., ValdeÁs, J.B., Rodriguez-Iturbe, I., 1979. On the in¯uence of the spatial distribution of rainfall on storm runoff. Water Resources Research 15 (2), 321±328. Zawadski, I., 1973. Statistical properties of precipitation patterns. Journal of Applied Meteorology 12, 459±472.