HYDROLOGICAL PROCESSES Hydrol. Process. (2007) Published online in Wiley InterScience (www.interscience.wiley.com) DOI: 10.1002/hyp.6851
Rain-gauge network evaluation and augmentation using geostatistics Ke-Sheng Cheng,* Yun-Ching Lin and Jun-Jih Liou Department of Bioenvironmental Systems Engineering, National Taiwan University, Taipei, Taiwan, ROC
Abstract: Rain-gauge networks are often used to provide estimates of area average rainfall or point rainfalls at ungauged locations. The level of accuracy a network can achieve depends on the total number and locations of gauges in the network. A geostatistical approach for evaluation and augmentation of an existing rain-gauge network is proposed in this study. Through variogram analysis, hourly rainfalls are shown to have higher spatial variability than annual rainfalls, with hourly Mei-Yu rainfalls having the highest spatial variability. A criterion using ordinary kriging variance is proposed to assess the accuracy of rainfall estimation using the acceptance probability defined as the probability that estimation error falls within a desired range. Based on the criterion, the percentage of the total area with acceptable accuracy Ap under certain network configuration can be calculated. A sequential algorithm is also proposed to prioritize rain-gauges of the existing network, identify the base network, and relocate non-base gauges. Percentage of the total area with acceptable accuracy is mostly contributed by the base network. In contrast, non-base gauges provide little contribution to Ap and are subject to removal or relocation. Using a case study in northern Taiwan, the proposed approach demonstrates that the identified base network which comprises of approximately two-thirds of the total rain-gauges can achieve almost the same level of performance (expressed in terms of percentage of the total area with acceptable accuracy) as the complete network for hourly Mei-Yu rainfall estimation. The percentage of area with acceptable accuracy can be raised from 56% to 88% using an augmented network. A threshold value for the percentage of area with acceptable accuracy is also recommended to help determine the number of non-base gauges which need to be relocated. Copyright 2007 John Wiley & Sons, Ltd. KEY WORDS
geostatistics; network evaluation; network augmentation; variogram analysis
Received 14 January 2007; Accepted 12 June 2007
INTRODUCTION Rainfall data are essential in many hydrological analyses and engineering design projects, including water budget analysis, frequency analysis and stormwater drainage design. Direct measurement of rainfall can only be achieved by rain-gauges, and rain-gauge networks are often installed to provide measurements that characterize the temporal and spatial variations of rainfall. However, even though modern rain-gauges are capable of providing rainfall rate in real time and at very fine resolution in time, the spatial variation of rainfall is still difficult to characterize without a rain-gauge network of enough density in space. Although recent advances in satellite remote sensing seem to have the potential to provide full spatial coverage of pixel rainfall estimates and have caused deterioration of rain-gauge networks in some cases (Ali et al., 2005), satellite images alone still cannot provide accurate rainfall estimates at the spatial resolution to match rain-gauge measurements. Satellite rainfall estimation algorithms must be calibrated and validated using rain-gauge networks, and in a few applications the remote-sensing satellite images were coupled with
* Correspondence to: Ke-Sheng Cheng, Department of Bioenvironmental Systems Engineering, National Taiwan University, Taipei, Taiwan, ROC. E-mail:
[email protected] Copyright 2007 John Wiley & Sons, Ltd.
measurements by rain-gauge networks (Krajewski, 1987) to yield better estimates of point or areal rainfalls. In addition, selection of rain-gauge locations is affected by many factors, such as accessibility, easiness of maintenance, topographical aspects, etc. Furthermore, the required minimum density of a rain-gauge network is dependent on the time resolution (or scale) of the desired rainfall measurements. Therefore, a methodology for performance evaluation and augmentation of an existing rain-gauge network is crucial, in that it can help to understand the capability of an existing network and the quality of the data it provides. Many approaches to optimal selection of rainfall gauges take into account the number and location of rainfall gauges to yield greater accuracy of areal rainfall estimation with minimum cost. These approaches are generally known as the variance reduction method (Bras and Rodriguez-Iturbe, 1976; Hughes and Lettenmaier, 1981; Bastin et al., 1984; Bogardi and Bardossy, 1985; Rouhani, 1985), which involves searching for the appropriate number of rainfall gauges and their locations. Geostatistics has been demonstrated to have many potential applications in hydrological research (Delhomme, 1978), in particular the optimal estimation of the average value over a region using the concept of variance reduction. A few methods and applications of rain-gauge
K.-S. CHENG, Y.-C. LIN AND J.-J. LIOU
network evaluation using geostatistics have been presented in the literature, including those by Bastin et al. (1984), Kassim and Kottegoda (1991), Pardo-Ig´uzquiza (1998), Tsintikidis et al. (2002), and St-Hilaire et al. (2003). Other approaches for rain-gauge network evaluation using information entropy were also proposed by Krstanovic and Singh (1992a,b) and Al-Zahrani and Husain (1998). Most approaches and applications of rain-gauge network optimization aim at providing accurate areal rainfall estimation. For such applications, performance evaluation of a network focuses on reducing the estimation variance of the areal rainfall, but not that of point rainfalls across the study area. In many applications, using areal rainfall as input to hydrological models may be enough for flood forecasting since catchments behave as low-pass filters and the necessary intervals for rainfall measurements in time and in space may be determined by variability in the discharge at the watershed outlet (Eagleson, 1967). However, most tributaries of river systems in Taiwan have steep slopes and respond quickly to local rainfalls which often are very intense. As a result, monitoring localized high-intensity rainfalls becomes immensely important, and rain-gauge networks have been installed tailoring to such needs. Therefore, the purpose of this study is to propose a rain-gauge network evaluation and augmentation approach focusing on accuracy assessment of point rainfalls across the whole study area. In the following sections we first describe the topographic and rainfall characteristics of the study area and the data selected for use in this study. We then go on to discuss the spatial variation characteristics of hourly and annual rainfalls using variogram analysis. A detailed account of the proposed network evaluation and augmentation approach is subsequently given. The final section gives the conclusions drawn from this research.
STUDY AREA AND DATA The Danshuei River basin downstream from Shihmen Reservoir in northern Taiwan was selected for this study (Figure 1). It covers an area of approximately 2200 km2 and encompasses three major tributaries and two reservoirs that are used to supply water for paddy irrigation and domestic and industrial usages and to attenuate peak flows during typhoon seasons. The furthest area downstream includes the Taipei metropolitan area and is prone to flood inundation. A network of 27 rain-gauges (see Figure 1) is jointly maintained by the Central Weather Bureau and the Water Resources Agency. Hourly rainfall data from 1990 to 2004 were recorded at all rain-gauge sites and were used for this study. In general, rainfall in Taiwan is brought by four season-specific dominant storm types: winter frontal rainfall (from November to early April of the next year), Mei-Yu (from late April to mid to end of June), convective rainfall and typhoons (both from July to October). The study area has an average annual rainfall of 2800 mm and its terrain elevation Copyright 2007 John Wiley & Sons, Ltd.
2 4 23
1
10 21
12 15
26
11
3 22
5
16
9
8
18
13
14 7
20 19
25
17
6
27 24 Taiwan
Figure 1. The study area and rain-gauge locations
Table I. Properties of the selected hourly rainfall data Storm type
Rainfall depth range (mm)
Total time (h)
Mei-Yu Convective Typhoon Frontal
0Ð49–10Ð27 0Ð46–10Ð30 1Ð39–23Ð62 0Ð38–4Ð36
118 128 274 236
varies from near sea level to about 1000 m above mean sea level. A data selection procedure was adopted for selecting hourly rainfall data of individual storm events to be used for this study. One criterion of the selection procedure was that, for each particular hour, at least twothirds of the 27 rain-gauges have non-zero hourly rainfall. Details of the selected hourly rainfall data are listed in Table I. CHARACTERIZING THE SPATIAL VARIATION OF RAINFALL DEPTHS Based on the criteria recommend by the World Meteorological Organization (WMO, 1974), the minimum raingauge density for the study area is between 100 and 250 km2 /gauge, or equivalently 9 to 22 stations in the study area. However, WMO’s criteria only represent a generally applicable minimum requirement and are not particularly tailored to certain storm types. Rainfall of different storm types may exhibit different levels of spatial variation. For example, hourly rainfall of convective storms and typhoons exhibits high spatial variation, whereas hourly frontal rainfall is less variable in space. In addition, long-duration (e.g. monthly, seasonal and annual) rainfall tends to be more spatially homogeneous than hourly and daily rainfall. For assessment of potential water resources, long-duration rainfall is of major concern, whereas short-duration rainfall is often used for design and modelling of flood mitigation projects. As such, evaluation and design of a rain-gauge network Hydrol. Process. (2007) DOI: 10.1002/hyp
RAIN-GAUGE NETWORK EVALUATION AND AUGMENTATION USING GEOSTATISTICS
should be based on its functional objective, which in turn also determines the duration of rainfall for analysis. In the following sections the spatial variation of rainfall of various durations and storm types is characterized using variogram analysis. Variogram analysis The random nature of spatial variation of many natural phenomena can be described by a random field Zx, where x represents the spatial location and Z is the random phenomenon under investigation, e.g. rainfall depth in our study. The spatial variation structure of a random field Zx with a stationary mean can be represented by its variogram, defined as
where EZ and Var(Z) respectively represent the expected value and the variance of the random variable Z. Equation (1) indicates that the variogram is independent of spatial locations, i.e. it only depends on the distance between the two points under consideration. The influence range of a variogram is the minimum distance jxi xj j beyond which the two random variables Zxi and Zxj become independent. For a second-order stationary random field, as the distance h increases, the variogram will reach an asymptotic value, known as the sill, which is numerically the same as the variance of the random variable ZX. Readers are referred to Journel and Huijbregts (1978) for details of variogram calculation. Storm events of various seasons and intensities may have different spatial variation structures; therefore, it may not always be realistic to adopt a unique variogram for all storm events, irrespective of rainfall intensities, and seasonal and meteorological conditions (Cheng et al., 2003). Bastin et al. (1984) suggested a climatological variogram model of the form m, h D ˛m Ł h
2
where m is an index for storm events, h is the Euclidian distance and ˛m is the scaling factor. With this structure, all the time (or event) nonstationarity is concentrated in the scaling factor ˛m, and the scaled component Ł h is time (or event) invariant and is called the dimensionless variogram. In regions of high meteorological variability, the factor ˛m accounts mainly for the scale effect due to the time/event variation of rainfall intensity. In this study, dimensionless variograms of annual and hourly rainfall were developed and used for evaluation of the rain-gauge network. Spatial variation of annual rainfall Under the second-order stationarity assumption, a variogram approaches its sill, which is numerically the same Copyright 2007 John Wiley & Sons, Ltd.
Rm D 2Ð35H C 2540Ð6 r D 0Ð61
1
3
It worth noting that gauges 10 and 21 are located immediately outside of the watershed boundary and both have about 58% of their average annual rainfalls contributed by frontal rainfall, whereas the frontal rainfall of other gauges accounts for only 30 to 40% of the annual rainfall. If the two gauges are excluded in trend analysis, then this will result in a new regression model with a correlation coefficient of 0Ð76. The scaling factor in Equation (2) is equivalent to the sill ω, or the variance of the rainfall field. In general, random fields (annual rainfall in our case) with higher measurement values tend to have a higher degree of variability, and this is known as the proportional effect in geostatistics. Therefore, years with higher annual rainfall depths correspond to higher sill values and a set of experimental variograms with different sill values will result. In order to remove the spatial and temporal nonstationarity and to construct the dimensionless variogram, annual rainfall data were preprocessed by Ri j Rm,i i D 1, 2, Ð Ð Ð , n; j D 1, 2, Ð Ð Ð , N sj 4 where Ri j and Rm,i respectively represent the jth year’s annual rainfall and the mean annual rainfall of the ith RiŁ j D
Average annual rainfall, Rm (mm)
1 Var[Zxi Zxj ] 2 1 D E[Zxi Zxj ]2 2
jxi xj j D
as the variance of the random field under consideration. However, if a spatial trend or orographic effect of the expected value of Zx exists, then the random field Zx is not stationary and the variogram may increase without a sill. Therefore, it is necessary to check for an orographic effect or trend existence of E[Zx] prior to experimental variogram fitting. Average annual rainfall depths Rm (mm) were calculated for each of the 27 rain-gauges using 15 years of hourly rainfall data. Figure 2 demonstrates that average annual rainfall depth varies with elevation H (m). In the study area, terrain elevation is lowest at the river mouth at the northwest corner and is highest in the southeast mountainous region, resulting in a spatial trend for average annual rainfall. The orographic effect of the average annual rainfall is characterized by a simple linear regression model:
5000
Gauge 10 Gauge 21
4000 3000 2000 Rm = 2.35H + 2540.6 (r = 0.61)
1000 0
0
200
400
600
800
1000
Elevation, H(m)
Figure 2. Relationship between average annual rainfall depth and elevation Hydrol. Process. (2007) DOI: 10.1002/hyp
K.-S. CHENG, Y.-C. LIN AND J.-J. LIOU
rain-gauge, and sj is the standard deviation of the jth year’s detrended annual rainfalls, i.e. n 1 sj D Rd,i j Rd j2 5 n 1 iD1 Rd,i j D Ri j Rm,i
6
1 Rd,i j n iD1 n
Rd j D
7
Hereafter, the preprocessed rainfall RiŁ j, i D 1, 2, Ð Ð Ð , n; j D 1, 2, Ð Ð Ð , N, is referred to as the rescaled rainfall. The rescaled rainfall is dimensionless and has zero mean and unit standard deviation. The experimental and fitted dimensionless variograms of the rescaled annual rainfall are shown in Figure 3. Parameters of the fitted exponential variogram model are given in Table II. The influence range (3a for the exponential variogram model) of the rescaled annual rainfall is about 63 km.
where Rm (mm) represents the average hourly rainfall of typhoons at a rain-gauge site. The proportional effect of hourly rainfalls (exemplified by hourly rainfall of Mei-Yu) is also evident, as can be seen in Figure 5; therefore, hourly rainfall data were also preprocessed in the same manner as the annual rainfall, except that no detrending was needed for Mei-Yu, convective, and winter frontal rainfalls. The parameters of dimensionless variograms for hourly rainfall of different storm types are tabulated in Table II. The influence range of hourly rainfall varies with storm type. In particular, hourly rainfall of Mei-Yu has the smallest influence range of 24 km, suggesting the highest spatial variation among all storm types. Figure 6 shows the experimental and fitted dimensionless variograms of hourly Mei-Yu rainfall. In comparison, the influence range of annual rainfall is significantly larger than that of hourly rainfalls due to the smoothing effect of rainfall accumulation. Such an effect has also been observed by Lebel et al. (1987) and Pardo-Ig´uzquiza (1998).
Spatial variation of hourly rainfall As mentioned previously, the spatial variation structure of hourly rainfall may vary with storm type. Therefore, variogram modelling of hourly rainfall was conducted separately for different storm types. No significant relation between the average hourly rainfall and terrain elevation was found for all storm types except typhoons (Figure 4). Similarly, the orographic effect of the average hourly rainfall of typhoons is characterized by Rm D 0Ð0034H C 5Ð01 r D 0Ð52
8
Dimensionless variogram, r (h)
2
1.5
1
0.5 Experimental 00
10000
Fitted
20000 30000 Distance, h (m)
40000
Figure 3. Experimental and fitted variograms of annual rainfall
Table II. Parameters for dimensionless variogram fitting. Variogram model:a h D ω[1 exph/a] Parameters Annual rainfall
Hourly rainfall Mei-Yu Convective Typhoon Frontal
ω a (m) a Influence
1Ð30 21 000
1Ð05 8 000
1Ð10 9 000
range: 3a.
Copyright 2007 John Wiley & Sons, Ltd.
1Ð20 11 500
1Ð20 10 300 Figure 4. Relationship between average hourly rainfall depth and elevation Hydrol. Process. (2007) DOI: 10.1002/hyp
RAIN-GAUGE NETWORK EVALUATION AND AUGMENTATION USING GEOSTATISTICS
Standard deviation of rainfall depth (mm)
25
20
15
10 N = 118 5
0
0
5
10 15 Average rainfall depth (mm)
20
Figure 5. Proportional effect of Mei-Yu hourly rainfall
Dimensionless variogram, r (h)
1.4 1.2 1 0.8 0.6 0.4 0.2 0
Experimental 0
10000
20000 30000 Distance, h (m)
Fitted 40000
50000
Figure 6. Experimental and fitted variograms of hourly Mei-Yu rainfall
Mei-Yu rainfalls are produced by surface frontal systems which advance southeastward from southern China to Taiwan from mid to late spring through early to mid summer each year. The fronts are usually accompanied by a synoptic-scale cloud band with embedded mesoscale convective systems (MCSs), extending several thousand kilometres from southern Japan to southern China with an approximately east–west orientation (Yeh et al., 2002). During the passage of a Mei-Yu frontal system, a few very active mesoscale convective cells may develop repeatedly, causing heavy and localized rainfall for the area. Although the synoptic-scale frontal system may last for a few days, the MCSs generally have lifetime of a few hours to 1 day only. We suspect that the small influence range of Mei-Yu rainfall may be attributed to redevelopment of MCSs.
RAIN-GAUGE NETWORK EVALUATION AND AUGMENTATION A rain-gauge network encompasses a number of gauges within an area of interest. Two major concerns for the evaluation and design of a network are the density and locations of these rain-gauges. The minimum density of Copyright 2007 John Wiley & Sons, Ltd.
a network depends on its functional objective and characteristics of dominant storm types. For water resources planning, estimates of the area average of long-duration rainfall, such as monthly, seasonal and annual rainfall, are needed. However, for the purpose of flood mitigation, interpolation of hourly rainfall at ungauged locations is often required and measurements of hourly or even smaller duration rainfalls are necessary. Whereas at smaller time scales (e.g. hourly and event rainfall) the rain fields may be treated as a purely stochastic process, rainfall of longer durations tends to exhibit spatial nonstationarity and the presence of a deterministic trend can be observed (Ali et al., 2003, 2005). With the trend presence, long-duration rainfall is stochastically less variable in space and has higher spatial correlation, as indicated by the influence range (3a) of annual rainfall in Table II. In contrast, short-duration rainfall can be characterized by second-order stationary random fields with smaller influence ranges or, equivalently, higher spatial variability. As such, for the purpose of water resources planning, a lowdensity network may be enough to capture the spatial variability of longer duration rainfall, whereas a highdensity rain-gauge network is required for flood forecasting and mitigation. In addition, the influence range of rainfall fields also varies with storm type. Among the four major storm types, hourly rainfall of Mei-Yu has the highest spatial variability and, therefore, is considered as the basis for subsequent network evaluation and design. In evaluating and designing a rain-gauge network, the major concern is how to quantify the performance of an existing or proposed network. The proposed network may be a chosen subset of the existing network, augmentation of a subset (including the whole set) of the existing network, or a new network including the augmented network and newly added gauges. Intuitively, a measure of the network performance is the accuracy associated with estimates of the areal average or point rainfall. The accuracy should reflect not only the absolute error, but also the error variance, since the errors depend heavily on rain field variability, and this variability is far from being consistent from one event to another and from one year to the next. If only unbiased estimators are considered, then the characteristics of estimation accuracy can thus be focused on the variance of estimation error. Using the error variance as a measure of network performance implies that the evaluation is conducted on an ensemble basis. Geostatistical methods, a set of spatial interpolation algorithms which yield unbiased estimators with minimum error variance, have been applied for many monitoring network evaluation practices. Among these methods, ordinary kriging is most widely used owing to its computational simplicity and data availability in most applications. Interpolation algorithm of ordinary kriging is briefly given in the following section. For a full account of the theoretical development of the method, readers are referred to Journel and Huijbregts (1978). Hydrol. Process. (2007) DOI: 10.1002/hyp
K.-S. CHENG, Y.-C. LIN AND J.-J. LIOU
Ordinary kriging Let Zx be a random variable defined at location x and fZx, x 2 g be a second-order stationary random field in a spatial domain . We now attempt to estimate an unknown value of Z at x0 , i.e.zx0 , using the observed values zxi , i D 1, 2, . . . , n, at neighbouring locations and the following linear equation: zO x0 D
n
i zxi
9
site x0 , i.e. Zx0 , t, can then be estimated using measurements Zx, t, i D 1, Ð Ð Ð , n, and Equation (9), and the estimation accuracy is given by the ordinary kriging variance in Equation (11). Since the estimation is made using the same-time measurements only, the time dependence of rainfall Zx, t will be dropped hereafter. Intuitively, an estimate is considered acceptable only if the estimate falls within a given range of the true value, i.e. 12 jQzx0 j D jOzx0 zx0 j < r
iD1
where i are weights assigned to measurements zxi . Weights i are obtained by solving the following ordinary kriging system equation: Ð Ð Ð 1i Ð Ð Ð 1n 1 1 10 11 .. .. .. .. .. .. .. .. . . . . . . . . i1 Ð Ð Ð ii Ð Ð Ð in 1 i i0 10 . .. .. .. .. .. .. .. . . . . . . . . . n n0 n1 Ð Ð Ð ni Ð Ð Ð nn 1 1 ÐÐÐ 1 ÐÐÐ 1 0 1 where is a Lagrange multiplier and ij represents the semivariance between zxi andzxj . Variance of the estimation error, also known as the ordinary kriging variance, is given by
k2 x0 D Ef[Zx0 Zx0 ]2 g D C
n
i i0
11
iD1
Bastin et al. (1984) and Kassim and Kottegoda (1991) used the kriging variance as a tool in rain-gauge network design for optimal estimation of the areal average rainfall. An iterative selection procedure was implemented, and, in each iteration, the rain-gauge which is associated with the minimum kriging variance was selected. All rain-gauges, therefore, are prioritized and such information is used for adding or deleting rain-gauges. Defining the acceptable accuracy In many applications, rain-gauge networks are designed to provide a good estimate of areal avearge rainfall. For such applications, the focus is placed on the estimation accuracy of areal average rainfall, and estimation accuracies of point rainfall at individual ungauged sites are not of major concern. On the contrary, if a rain-gauge network is tailored for flood mitigation, it is often desired to yield high estimation accuracy for point rainfall across the study area. The estimation accuracy of point rainfall varies within the study area and is dependent on the number and geometric layout of rain-gauges. Practically speaking, a good rain-gauge network should yield acceptable accuracy for most points within the study area. In this study, a criterion described in the following section is adopted to define the acceptable accuracy. Suppose that hourly rainfalls Zx, t are measured by a network of n rain-gauges at location xi , i D 1, Ð Ð Ð , n, for a period of time t1 t tp . Rainfall at an ungauged Copyright 2007 John Wiley & Sons, Ltd.
where r > 0. However, at a given location x0 , the estimation accuracy varies from hour to hour and from event to event; thus, it should be evaluated on an ensemble basis. Also, the given range r should be determined by taking into account the variance of the rainfall field Zx, i.e. z2 . Therefore, a revised and more reasonable criterion is given by P[jOzx0 zx0 j < kz ] ½ ˛
13
In the above equation, the acceptable range of the estimation error is expressed in terms of the standard deviation of the random variable Zx, and the multiple k and the minimum probability ˛ are to be chosen according to factors such as available budget for gauge installation and maintenance and desired estimation accuracies. For our study we chose k D 1 and value of ˛ is either 0Ð9 or 0Ð8. For ordinary kriging estimation, variance of the estimation error (i.e. the kriging variance) is given by Equation (11). Since the ordinary kriging estimator is unbiased, the estimation error at x0 has zero mean and variance k2 x0 . If the estimation error is further assumed to be normally distributed, then the probability for the estimation error zQ x0 to fall within the desired range (z , z ) can be readily determined using the cumulative probability of the standard normal distribution: P[jQzx0 j < z ] D P[jOzx0 zx0 j < z ] z jQzx0 j < DP k x0 k x0 z Ł D P jQz x0 j < D pA x0 14 k x0 In Equation (14), zQ Ł x0 is referred to as the standardized estimation error and has a standard normal distribution, i.e. zQ Ł x0 ¾ N0, 1. Also, pA x0 is termed the acceptance probability at x0 , which represents the probability that the estimation error at x0 is less than z . Estimation accuracy at an ungauged point is acceptable only if the associated acceptance probability is no less than ˛; accordingly, estimation at that point is said to have acceptable accuracy. Figure 7 illustrates the probability distribution of estimation error with different kriging vari2 2 ances k1 x0 and k2 x0 . Points associated with higher kriging variances correspond to lower acceptance probabilities, and all points with acceptance probability lower than ˛ are considered as not having acceptable accuracy. It is worth noting that, for our study, z in Equation (14) Hydrol. Process. (2007) DOI: 10.1002/hyp
RAIN-GAUGE NETWORK EVALUATION AND AUGMENTATION USING GEOSTATISTICS
is determined as the sill value of the dimensionless variogram, since rainfall data were preprocessed using Equation (4). Network performance evaluation based on the percentage of area with acceptable accuracy In the previous section we elaborated that the estimation accuracy at every point within the study can be evaluated using the acceptance probability. Therefore, performance evaluation of a rain-gauge network can be made based on the percentage of area with acceptable accuracy (hereafter expressed by Ap ). Also, since the acceptance probability can be calculated at each point within the study area, a contour map of acceptance probability can be developed to assist in evaluation and augmentation of an existing network. In this study, hourly rainfall of Mei-Yu and annual rainfall were used to establish contour maps of acceptance probability. The whole study area was discretized into a set of grid nodes with 1 km interval and the acceptance probability at each node was calculated. Figure 8 shows contour maps of acceptance probability for hourly rainfall of Mei-Yu and annual rainfall. In comparison, for a fixed value of ˛, say 0Ð8 or 0Ð9, annual rainfall is associated with a higher value of Ap due to its higher spatial correlation. For example, at ˛ D 0Ð8, about 56% of the total area has acceptable accuracy for hourly rainfall, whereas 82% of the total area has acceptable accuracy for annual rainfall. It can also be seen that Ap values are very low at ˛ D 0Ð9 (16% and 36% for hourly Mei-Yu and annual rainfall respectively), suggesting that ˛ D 0Ð9 may not be a realistic choice for the study area. If a threshold value of Ap , say 60%, is set as the network evaluation and augmentation criterion, then, at ˛ D 0Ð8, the current network can barely meet the criterion. However, if the threshold is set at a higher level, say 80%, then the current network fails the evaluation and network augmentation is necessitated. It is also worth mentioning that, at a rain-gauge site, the acceptance probability always equals unity, since ordinary kriging is an exact estimator and yields zero estimation error at the measurement points (de Marsily, 1986).
2
N(0, σk1) σz 2
2
σk1 < σk2
2
N(0, σk2)
0
Z(x)
Figure 7. Probability distribution of the estimation error under different kriging variances Copyright 2007 John Wiley & Sons, Ltd.
Augmentation of an existing rain-gauge network An existing network can be augmented in terms of the percentage of area with acceptable accuracy by relocating the existing gauges or adding new gauges. The sequential algorithm described below was devised to prioritize the existing rain-gauges and to evaluate sequentially the joint performance of a subset of rain-gauges. 1. One gauge at location xi , i 2 f1, Ð Ð Ð , mg, hereafter referred to as gauge xi , is chosen from a set of m m n remaining gauges S1 D fxj , j D 1, Ð Ð Ð , mg. Initially, the set of remaining gauges consists of all existing raingauges (m D n). Then, the kriging variance at each grid node is calculated using Equation (11) with only (m 1) gauges involved, i.e. S2 D fxj , j D 1, Ð Ð Ð , m; j 6D ig. The percentage of area with acceptable accuracy, when only gauges S2 D fxj , j D 1, Ð Ð Ð , m; j 6D ig are in place, is then obtained. 2. Return xi to S2 and select another gauge xk k 6D i from S1 . Then, the kriging variance at each grid node and the Ap value are recalculated. This process is repeated until all gauges in S1 have been chosen and a set of m values of Ap are obtained. 3. Remove the gauge associated with the highest value of Ap in step (2) from S1 . Then, reduce the number of remaining gauges by one and repeat steps (1) and (2). This step (3) is repeated until there is only one gauge remaining in S1 . Upon completion of the sequential algorithm, all raingauges are prioritized according to their order of removal in step (3). In addition, at each stage when a gauge is removed in step (3), contour maps of acceptance probability for hourly rainfall of Mei-Yu and annual rainfall and their corresponding Ap values were also established using only the remaining gauges. This information forms the basis of the proposed rain-gauge network augmentation approach. Using the values of Ap corresponding to removal of individual rain-gauges, or equivalently a set of remaining gauges, illustrative figures can be constructed to show the prioritized order of rain-gauges and performance of a subset of rain-gauges. Figure 9 demonstrates the result of rain-gauge prioritization and the corresponding values of Ap for hourly Mei-Yu and annual rainfall. For hourly Mei-Yu rainfall, it can be seen clearly that, at ˛ D 0Ð8, with only 19 gauges (gauges 4, 22, 2, 10, 11, 16, 18, and 15 were sequentially removed) the remaining network can achieve almost the same level of Ap as the complete network of 27 gauges. These remaining gauges form the base network and will not be relocated in the subsequent network augmentation process. Figure 10 shows the contour map of acceptance probability for hourly rainfall using only the base network. The contour maps in Figures 8a and 10 are almost identical, except for minor variations in close vicinity of the removed gauges, validating the removal of these gauges. Similarly, for annual rainfall, the base network is comprised of only Hydrol. Process. (2007) DOI: 10.1002/hyp
K.-S. CHENG, Y.-C. LIN AND J.-J. LIOU
2790000
2780000
TM2 Northing (m)
2770000
2760000
2750000
2740000
(a) 2730000 280000
300000 TM2 Easting (m)
320000
340000
2790000
2780000
TM2 Northing (m)
2770000
2760000
2750000
2740000
(b) 2730000 280000
300000 TM2 Easting (m)
320000
340000
Figure 8. Contour map of acceptance probabilities of (a) hourly rainfall of Mei-Yu and (b) annual rainfall (number of rain-gauges: 27)
15 gauges (gauges 2, 4, 5, 11, 12, 22, 9, 10, 25, 13, 16 and 14 were sequentially removed) and still can achieve almost the same level of Ap as the complete network. Rain-gauges not being included in the base network are redundant and provide little contribution to the network. These rain-gauges can either be abolished to reduce the maintenance cost or be relocated to achieve higher level of Ap . Relocating rain-gauges can again be conducted in a sequential manner. The first relocated point is determined by searching one point among all points with pA < ˛ (grey region in Figure 10) that, together with the base network, yields the highest value of Ap . Such a search generally does not take much effort, since, with the Copyright 2007 John Wiley & Sons, Ltd.
base network in place, most points already have pA ½ ˛. The rest of the relocated points can then be sequentially determined in a similar manner. Figure 11 shows the locations of the eight relocated gauges and the contour map of acceptance probability of hourly rainfall established by the augmented network, which includes the base network and all relocated gauges. Comparing with Figure 10, the area with acceptable accuracy (non-grey area) in Figure 11 increases significantly. The level of Ap that can be achieved by sequentially adding relocated gauges to the base network is also illustrated in Figure 9. For the hourly Mei-Yu rainfall case, Ap increases quickly and reaches 88% when all non-base Hydrol. Process. (2007) DOI: 10.1002/hyp
RAIN-GAUGE NETWORK EVALUATION AND AUGMENTATION USING GEOSTATISTICS
Existing network
Augmented network
Raingauges to be relocated for network augmentation.
(a)
density and configuration of a rain-gauge network that can be practically implemented is dependent on the available budget. A final note about the network augmentation approach presented in this paper is that, even though hourly MeiYu and annual rainfall data were both used for evaluation and augmentation of the existing network, the resultant network using hourly Mei-Yu data should be adopted for practical operation since it has higher rain-gauge density and satisfies the functional objectives of both water resources planning and flood mitigation. Inclusion of network augmentation using annual rainfall data are merely for comparative illustration purposes. CONCLUSIONS
Existing network
Augmented network
Raingauges to be relocated for network augmentation. (b)
Figure 9. Raingauge prioritization and corresponding Ap values. (a) Hourly rainfall of Mei-Yu. (b) Annual rainfall. Note: (1) ˛ D 0.8, (2) the number adjacent to a point represents the gauge number
gauges are relocated (see Figure 9a). The amount of increase in Ap (about 32%) is very significant and demonstrates the applicability of the proposed network augmentation approach. However, readers are also reminded that it may not always be necessary to relocate all non-base gauges, as illustrated by Figure 9b. For the annual rainfall case, the value of Ap almost reaches 100% after only three or four non-base gauges are relocated; therefore, relocating more gauges is unnecessary. Also, as was mentioned in the ‘Network performance evaluation based on the percentage of area with acceptable accuracy’ section, if the threshold value of 80% is set for Ap of hourly Mei-Yu, then after relocation of five non-base gauges the augmented network already achieves higher than 80% of Ap and relocation of the other three non-base gauges is not necessary. Comparing this approach with other approaches that aim to yield greater accuracy of areal rainfall estimation, the approach proposed in this study focuses on accuracy of point rainfall across the whole study area. It is particularly applicable in watersheds with localized highintensity rainfall and a short time of concentration. It also has the advantage of being highly flexible in parameters related to accuracy assessment, such as k and ˛ in Equation (13) and the percentage of area with acceptable accuracy Ap . Such flexibility is necessary because the Copyright 2007 John Wiley & Sons, Ltd.
A geostatistical approach for rain-gauge network evaluation and augmentation is presented. The main concept of the proposed approach is to evaluate a network based on the percentage of the total area with acceptable accuracy. Variances of the estimation error at ungauged sites can be calculated by ordinary kriging and are assumed to be normally distributed. The acceptance probability is defined as the probability that the estimation error falls within a desired range expressed in terms of the standard deviation of rainfall. Estimates with acceptance probability no less than a given value (0Ð8 in our study) are said to have acceptable accuracies. The acceptance probability can be determined according to the cumulative probability of the standard normal distribution. Finally, a threshold value for the percentage of area with acceptable accuracy can be set as the network evaluation and augmentation criterion. A sequential algorithm is also proposed to prioritize the existing rain-gauges, identify the base network, and relocate non-base gauges. A few concluding remarks are drawn as follows: 1. The spatial variation characteristics of hourly rainfall of different storm types and annual rainfall were analysed. Annual rainfall is shown to exhibit a significant orographic effect and less spatial variability, whereas hourly rainfall exhibits higher variability in space and the spatial variation structures vary among different storm types. In particular, hourly Mei-Yu rainfall has the highest spatial variability, which may be attributed to the redevelopment of MCSs embedded in Mei-Yu fronts. 2. Given a threshold value of 0Ð8 for acceptance probability, the existing network is shown to achieve 56% and 82% of the total area with acceptable accuracy for hourly Mei-Yu and annual rainfall estimations respectively. 3. Using the sequential algorithm, a base network for hourly rainfall comprising only about two-thirds of the existing gauges was identified. The base network can achieve almost the same level of the percentage of area with acceptable accuracy as the complete network. The base network for annual rainfall comprises even less gauges. Hydrol. Process. (2007) DOI: 10.1002/hyp
K.-S. CHENG, Y.-C. LIN AND J.-J. LIOU
2790000
2780000
TM2 Northing (m)
2770000
2760000
2750000
2740000
2730000 280000
300000 TM2 Easting (m)
320000
340000
Figure 10. Contour map of the acceptance probability established by the base network of hourly Mei-Yu rainfall. The non-grey region represents area with acceptable accuracy. (Number of raingauges D 19)
2790000
TM2 Northing (m)
2780000
2770000
2760000
2750000
: relocated gauge
2740000
2730000 280000
300000 TM2 Easting (m)
320000
340000
Figure 11. Contour map of the acceptance probability established by the augmented network (including base stations and all relocated gauges) of hourly Mei-Yu rainfall. The non-grey region represents area with acceptable accuracy. (Number of raingauges D 27)
4. Non-base gauges were sequentially relocated to achieve higher levels of percentage of area with acceptable accuracy. For hourly Mei-Yu rainfall, the percentage of area with acceptable accuracy can be raised from 56% for the existing network to 88% for an augmented network, demonstrating the applicability of the proposed network augmentation approach. 5. Given a threshold value for the percentage of area with acceptable accuracy, it may not always be necessary to relocate all non-base gauges. For example, in this Copyright 2007 John Wiley & Sons, Ltd.
study, the threshold value was set at 80% for hourly Mei-Yu rainfall and only five of the eight non-base gauges needed to be relocated to achieve the desired level of accuracy.
ACKNOWLEDGEMENTS
We are grateful for the assistance in data collection offered by the Central Weather Bureau and Water Resources Agency. The first author is also grateful to the Hydrol. Process. (2007) DOI: 10.1002/hyp
RAIN-GAUGE NETWORK EVALUATION AND AUGMENTATION USING GEOSTATISTICS
National Science Council of Taiwan, ROC, and Kyoto University, Japan, for providing financial and facility supports for his sabbatical leave at Kyoto University during which this manuscript was prepared. Our thanks also extend to two anonymous reviewers for their insightful comments.
REFERENCES Ali A, Lebel T, Amani A. 2003. Invariance in the spatial structure of Sahelian rain fields at climatological scales. Journal of Hydrometeorology 4: 996–1011. Ali A, Lebel T, Amani A. 2005. Rainfall estimation in the Sahel. Part I: error function. Journal of Applied Meteorology 44: 1691– 1706. Al-Zahrani M, Husain T. 1998. An algorithm for designing a precipitation network in the south-western region of Saudi Arabia. Journal of Hydrology 205: 205– 216. Bastin G, Lorent B, Duque C, Gevers M. 1984. Optimal estimation of the average areal rainfall and optimal selection of raingauge locations. Water Resources Research 20: 463– 470. Bogardi I, Bardossy A. 1985. Multicriterion network design using geostatistics. Water Resources Research 21: 199– 208. Bras RF, Rodriguez-Iturbe I. 1976. Network design for the estimation of areal mean rainfall events. Water Resources Research 12: 1185– 1195. Cheng KS, Wei C, Cheng YB, Yeh HC. 2003. Effect of spatial variation characteristics on contouring of design storm depth. Hydrological Processes 17: 1755– 1769. DOI: 10Ð1002/hyp.1209. Delhomme JP. 1978. Kriging in the hydrosciences. Advances in Water Resources 1: 251– 266. De Marsily G. 1986. Quantitative Hydrogeology. Academic Press: Orlando, FL. Eagleson PS. 1967. Optimum density of rainfall networks. Water Resources Research 3: 1021– 1033.
Copyright 2007 John Wiley & Sons, Ltd.
Hughes JP, Lettenmaier DP. 1981. Data requirements for kriging: estimation and network design. Water Resources Research 17: 1641– 1650. Journel AG, Huijbregts CJ. 1978. Mining Geostatistics. Academic Press: London. Kassim AHM, Kottegoda NT. 1991. Rainfall network design through comparative kriging methods. Hydrological Sciences Journal 36: 223– 240. Krajewski WF. 1987. Cokriging radar-rainfall and rain gage data. Journal of Geophysical Research 92: 9571– 9580. Krstanovic PF, Singh VP. 1992a. Evaluation of rainfall network using entropy: I. Theoretical development. Water Resources Management 6: 279– 293. Krstanovic PF, Singh VP. 1992b. Evaluation of rainfall network using entropy: II. Application. Water Resources Management 6: 295–314. Lebel T, Bastin G, Obled C, Creutin JD. 1987. On the accuracy of areal rainfall estimation: a case study. Water Resources Research 23: 2123– 2134. Pardo-Ig´uzquiza E. 1998. Optimal selection of number and location of rainfall gauges for areal rainfall estimation using geostatistics and simulated annealing. Journal of Hydrology 210: 206– 220. Rouhani S. 1985. Variance reduction analysis. Water Resources Research 21: 837– 846. St-Hilaire A, Ouarda TBMJ, Lachance M, Bob´ee B, Gaudet J, Gignac C. 2003. Assessment of the impact of meteorological network density on the estimation of basin precipitation and runoff: a case study. Hydrological Processes 17: 3561– 3580. Tsintikidis D, Georgakakos KP, Sperfslage JA, Smith DE, Carpenter TM. 2002. Precipitation uncertainty and raingauge network design within Folsom Lake watershed. Journal of Hydrologic Engineering 7: 175– 184. WMO. 1974. Guide to Hydrometeorological Practices, 3rd edition. World Meteorological Organization: Geneva. Yeh HC, Chen GT-J, Liu WT. 2002. Kinematic characteristics of a Meiyu front detected by the QuikSCAT oceanic winds. Monthly Weather Review 130: 700– 711.
Hydrol. Process. (2007) DOI: 10.1002/hyp
