Environ Monit Assess (2008) 145:375–386 DOI 10.1007/s10661-007-0046-z

Optimal extension of the rain gauge monitoring network of the Apulian Regional Consortium for Crop Protection E. Barca & G. Passarella & V. Uricchio

Received: 10 May 2007 / Accepted: 30 October 2007 / Published online: 23 November 2007 # Springer Science + Business Media B.V. 2007

Abstract The goal of this paper is to provide a methodology for assessing the optimal localization of new monitoring stations within an existing rain gauge monitoring network. The methodology presented, which uses geostatistics and probabilistic techniques (simulated annealing) combined with GIS instruments, could be extremely useful in any area where an extension of whatever existing environmental monitoring network is planned. The methodology has been applied to the design of an extension to a rainfall monitoring network in the Apulia region (South Italy). The considered monitoring network is managed by the Apulian Regional Consortium for Crop Protection (ARCCP), and, currently consists of 45 gauging stations distributed over the regional territory, mainly located on the basis of administrative needs. Fifty new stations have been added to the existing monitoring network, split in two groups: 15 fixed and 35 mobile stations. Two different methods were applied and tested: the Minimization of the Mean of Shortest Distances method (MMSD) and Ordinary Kriging (OK) whose related objective function is estimation variance. The MMSD, being a E. Barca : G. Passarella (*) : V. Uricchio Water Research Institute, CNR, V.le De Blasio, 5, 70123 Bari, Italy e-mail: [email protected]

purely geometric method, produced a spatially uniform configuration of the gauging stations. On the contrary, the approach based on the minimization of the average of the kriging estimation variances, produced a less regular configuration, though a more reliable one from a spatial standpoint. Nevertheless, the MMSD approach was chosen, since the ARCCP’s intention was to create a new monitoring network characterized by uniform spatial distribution throughout the regional territory. This was the most important constraint given to the project by the ARCCP, whose main objective was to accomplish a territorial network capable of detecting hazardous events quickly. A seasonal aggregation of the available rainfall data was considered. The choice of the temporal aggregation in quarterly averages allowed four different optimal configurations to be determined per season. The overlapping of the four configurations allowed a number of new station locations, which tended to remain fixed season after season, to be identified. Other stations, instead, changed their coordinates considerably over the four seasons. Constraints were defined in order to avoid placing new monitoring locations either near existing stations, belonging to other Agencies, or near the coast line. Keywords Monitoring . Rain gauges . Computational statistics . Simulated annealing . Geostatistics . GIS

376

Introduction With the growth of public environmental awareness and the contemporary improvement in national and EU legislation regarding the environment, monitoring has assumed great importance in the frame of all those managerial activities related to monitoring and safeguarding the environment. In particular, over the last decade, a number of public agencies whose purpose is to monitor meteorological, hydrological and hydrogeological parameters etc., have invested great economic, technical and human resources in planning and operating improvements on existing monitoring networks within their catchment areas. The problem of extending an environmental monitoring network (EMN) has consequently increased its importance in scientific literature because of the need to produce reliable managerial tools (Arbia and Espa 2001; Bogardi et al. 1985; Carrera and Szidarovzsky 1985; Cox and Cox 1994; Fedorov and Hackl 1994; Harmancioglu et al. 1999; Knopman and Voss 1989; Meyer et al. 1994; Nunes et al. 2002; Van Groenigen and Stein 1998; Wu 2004). Meteorological monitoring networks and particularly those devoted to rainfall monitoring, are among those which have received most attention from researchers, with a consequent abundance in the production of scientific papers, undoubtedly due to the importance of this resource (Al-Zahrani and Husain 1998; Bastin et al. 1984; Bras and Rodríguez-Iturbe 1975; Bras and Rodríguez-Iturbe 1976; Goovaerts 2000; Lebel et al. 1987; Papamichail and Metaxa 1996; Rodríguez-Iturbe and Mejía 1974). In particular, in the scientific community, the problem of the extension of rainfall monitoring networks has been tackled by searching for optimal criteria for the positioning of new measuring gauges. However, the sparse spatial coverage of regional territories, and/or the technological obsolescence of the gauges already installed, often provides scarce information on which to base reliable decision processes. Recent scientific literature has provided various approaches, characterized by different levels of complexity according to the level of detail required, capable of supporting both the design and realization of such networks (Ashraf et al. 1996). The methodology proposed in this paper integrates processes of stochastic and geostatistical theory with

Environ Monit Assess (2008) 145:375–386

optimisation methods based on simulation tools (Pardo-Igùzquiza 1998). The methodology was applied to the design of an extension to a rainfall monitoring network doubling the number of the gauging stations. The considered monitoring network is managed by the Apulian Regional Consortium for Crop Protection (ARCCP) and it currently consists of 45 gauging stations distributed randomly over the regional territory of Apulia (South Italy) mainly as a result of administrative needs. The Apulian territory is also covered by a second and denser network (about 150 stations), managed by the Hydrographic Regional Office (HRO). Obviously, the institutional aims of the two networks differ, nevertheless, from a general managerial and economic perspective, it is desirable that the monitoring locations designed to widen the existing ARCCP network should not overlap those belonging to the concurrent network. The methodology needs, therefore, to be sufficiently flexible to exclude a new proposed position, if the location is already covered by a gauge belonging to the second network. In general, the methodology allows buffer zones having a different amplitude to be defined, where new monitoring sites (e.g.: the coastal area) are not needed.

Methodology The proposed methodology consists of two steps; in the first it is necessary to define an objective function to be minimized. There are two possible choices: the Minimization of the Mean of Shortest Distances method (MMSD) and Ordinary Kriging (OK) whose related objective function is the estimation variance. The MMSD criterion was defined by Van Groenigen and Stein (1998) and modified, to take into account secondary information (weights), (Van Groenigen et al. 2000). This criterion, as inferred by its definition, is independent of the measured values and entirely based on the relative position of the considered points; therefore it is a geometric criterion and it provides extremely regular final configurations. In the modified version (Van Groenigen et al. 2000), it is possible to introduce weights, conditioning the related objective function, in order to define areas with a greater or lesser need for monitoring sites.

Environ Monit Assess (2008) 145:375–386

The OK criterion, derived from the theory of regionalized variables (Matheron 1970), allows the value of the estimation variance to be calculated in every location of the new configuration (Journel and Huijbrechts 1978; Isaaks and Srivastava 1989; Goovaerts 1997). In this case it is possible to define as the objective function, the average or the maximum estimation variance. Also in this case, the estimation variances depend uniquely on the sampling configuration; nevertheless a variogram model needs to be defined (Journel and Huijbrechts 1978; Isaaks and Srivastava 1989; Goovaerts 1997) that implicitly models the spatial behaviour of the considered variable. This phase is usually defined as structural analysis. The choice between the MMSD and OK optimisation approach cannot be based on a stringent theoretical and quantitative criterion. MMSD, being a geometrical driven method, produces spatially even distributions of monitoring point locations, while OK, based on estimation variance minimization, produces a monitoring network capable of providing better estimations in non-sampled points. In short, the first method appears to be more useful for designing “alert monitoring networks”, while the second is necessary when a reliable statistic description of a spatial phenomenon is required. The second step of the proposed methodology consists in the application of so-called “simulated annealing”, which provides a number of random configurations “driven” by the objective function. This method, implemented by Deutsch and Journel (1992), is used for finding the optimum in combinatorial problems, when the optimal solution of a given problem needs to be selected among a large number of possible available solutions without exploring them all. The theory of “simulated annealing” is based on the analogy with the organization of the atom network of a metal when it undergoes a process of temperature change (abrupt heating and slow cooling). Following this process, the atoms of the metal change their arrangement to a configuration of low energetic maintenance cost. In the analogy, the configuration of the atoms corresponds to that of the sampling points while the objective function corresponds to the energy of the system (Pardo-Igùzquiza 1998, Deutsch and Cockerham 1994). In algorithmic terms, with reference to the described metallurgical analogy (Metropolis et al.

377

1953), we assign an initial value to the temperature of the system, then we randomly choose a starting configuration from all the possible configurations, and we determine the corresponding value of the objective function, which is called energy. The temperature drives the duration of the process and, at every following step it decreases down to zero, which is the final temperature; the slower the cooling the higher the probability of finding the optimal configuration is, while the greater the initial temperature, the higher is the probability that the final configuration matches the absolute optimum that is the absolute minimum for the objective function. The starting configuration is perturbed in a randomised way, varying the position of only one sampling point of the monitoring network at a time, and the corresponding value of the objective function is computed again. If the perturbed configuration is better than the previous one (i.e.: the value of the objective function decreases) it is assumed as a transitory excellent solution; otherwise, the new configuration is not automatically discharged, as would happen with a classical method of optimisation, but it is submitted to a probabilistic criterion of acceptance which compares it again with the transitory optimal configuration. If this probabilistic criterion establishes that the configuration is acceptable, it is accepted as a transitory optimal solution. In detail, this happens by verifying that the following expression:   ΔE exp
Ti

ð1Þ

where ΔE represents the variation of the objective function, and Ti the current value of the temperature parameter, is smaller than a randomly generated number. This test allows the method to avoid the process of converging to a local optimum rather than the global one. Independently from the criterion chosen, another specific requirement was considered for improving the monitoring network. In fact, once the results had been obtained from one of the two methods, the option of determining a number of “mobile” stations, among the new monitoring locations was investigated. This option would allow the stations to be moved within a given distance, during the seasons of the year. The main reason why this option was investigated is that a “fixed” monitoring network may

378

sometimes be considered too rigid by agency managers, who ask for a certain flexibility in the gauge positioning. Therefore, the methodology was repeated, considering seasonal rainfall means for the OK, and the four best realizations for the MMSD criterion, thus four different configurations were determined. Successively, a possible maximum tolerance of 10 km was introduced. In practice, a regular mesh of 10 km side cells was overlaid over the study area map and the locations of the new gauging stations were associated to the correspondent mesh cell by means of GIS software (ESRI 1996). As a result, it was possible to distinguish two types of stations: those whose position remained fixed in the same mesh cell, throughout the four seasons and those whose position changed. The new gauging stations belonging to the former group, i.e., those which did not move from their original mesh cell even when a reduction of the position tolerance to 5 km was considered, were defined “fixed stations”. On the contrary, those stations which moved from one cell to another during the different monitoring seasons were labelled as “mobile stations”. All the remaining stations were defined as “potentially mobile stations” which means that, even considering these stations as fixed stations, it would be possible to choose some mobile stations among them, to be moved to some other location in particular seasons and conditions. Fig. 1 Monitoring network of the Apulian Regional Consortium for Crop Protection and the Regional Hydrographic Office. Stations lying outside the Apulian boundaries belong to interregional catchments

Environ Monit Assess (2008) 145:375–386

Study case The methodology was applied to the design of an extension to the rainfall monitoring network located in the Apulian region. The monitoring network considered is managed by the Apulian Regional Consortium for Crop Protection (ARCCP) and currently consists of 45 stations irregularly spread over the regional territory. A second meteorological monitoring network exists covering the area considered, managed by the Hydrographic Regional Office and consisting of about 150 stations. Figure 1 shows the location of all the existing stations belonging to the two networks, including also some stations outside the regional boundaries but belonging to interregional hydrographic basins. As stated above, even though the two networks have different institutional goals, a design constraint given by the ARCCP was that the monitoring locations should not overlap those belonging to the concurrent network. Nevertheless, in the present study, measurements from the HRO stations were used to improve our knowledge of the spatial behaviour of the mean seasonal rainfall, but their locations were constrained so as not to allow points in the optimisation algorithm. Other constraints were defined related to the distance of new monitoring points from both the existing ARCCP stations and the coast line.

Environ Monit Assess (2008) 145:375–386

379

The simulations were performed by using the software SANOS (Van Groenigen and Stein 2000; Van Groenigen 2000). It is an established fact that any spatial analysis is strongly dependent on the available data. In the study case, they were taken from the electronic files of the two regional agencies. As already mentioned previously, the data were preliminarily submitted to a statistical exploratory analysis. In particular, for every station, rainfall data were aggregated per season. In the following table the percentages of stations belonging to each of the five Apulian provinces are reported. As Table 1 clearly shows, an equal number of stations has been allocated to each province, providing an almost uniform distribution from an administrative standpoint. However, this distribution does not guarantee spatial uniformity, since the provinces vary in size. It was therefore decided to plan design simulations in order to re-equilibrate the coverage percentages throughout the regional territory, favouring those provinces having a worse gauge/km2 ratio. ARCCP provided a series of daily rainfall data related to the period 2000–2003. Obviously, a 3 year temporal series is not enough to get representative seasonal average values. Consequently, in order to obtain the best possible characterisation of the spatial behaviour of rainfall over the regional territory, a decision was taken to use the historical series provided by the HRO. This choice, however, did not affect the correctness of the application of the methodology; in fact, these data were used only to determine the spatial law characterizing the mean seasonal rainfall rate throughout the Apulian territory. The precision in determining the spatial law depends,

Table 1 Density of rain gauging stations, in each province, of the Apulian Regional Consortium for Crop Protection Province

Area (km2)

No. of gauging stations

Density (stations/km2)

Density (%)

Bari Brindisi Foggia Lecce Taranto Apulia

5,127,609 1,843,752 7,157,063 2,770,077 2,438,445 19,336,946

9 9 9 9 9 45

0.00176 0.00488 0.00126 0.00325 0.00369

11.83 32.91 8.48 21.90 24.88

obviously, on the abundance of the available data throughout the territory. Thus, the historical series published by the HRO were used only to appraise the spatial behaviour of the mean seasonal rainfall, but were ignored during the actual optimisation phase, when, instead, only the positions of the existing ARCCP gauging stations were considered. The historical series provided by the HRO cover about a 50 year period, approximately from the 1950s to today. This interval is long enough to define the quarterly mean behaviour of rainfall, filtering possible distortions due to intense phenomena and, in particular, rainy or dry periods. As stated above, there are about 150 monitoring stations belonging to the HRO network, but 27 of them are located in the provinces of Potenza and Avellino (Fig. 1), outside the Apulian borders, for monitoring the inter-regional basins of the rivers Ofanto, Candelaro and Carapelle. Unfortunately, for various reasons, only 93 gauging stations were actually usable for the simulations, instead of 150, corresponding to a spatial density of about 0.005 stations per squared kilometre. Using the aggregated values of these stations, some preliminary, descriptive statistics were computed in order to evaluate the related PDFs. In fact, it is preferable that these distributions should be normal to respect the ordinary kriging hypotheses. Table 2 reports the main descriptive statistics for each period of 3 months.

Table 2 Descriptive statistics of precipitation (mm) in the four quarters of the year

No. of cases Minimum Maximum Range Median Mean 95% CI. sup. 95% CI. inf. Std. error Standard dev. Variance C.V. Skewness Kurtosis

Quart 1

Quart 2

Quart 3

Quart 4

93 38.0 94.7 56.7 63.3 63.8 66.1 61.5 1.2 11.2 126.0 0.2 0.4 −0.1

93 23.4 66.4 43.0 36.8 39.7 41.6 37.7 1.0 9.6 92.7 0.2 0.7 −0.2

93 22.6 59.2 36.6 33.8 35.1 36.4 33.8 0.7 6.3 39.3 0.2 0.9 1.7

93 51.2 119.0 67.8 73.8 77.9 81.2 74.6 1.7 16.2 262.2 0.2 0.7 −0.2

380

Table 2 shows that the values of the mean and the median related to each quarter are almost the same, pointing out a tendency to symmetry of the sample distributions (Ott 1995). In fact, since the values used are quarterly means of daily values, unless phenomena of casual or systematic distortion affect the starting values, the distributions are expected to be normal, because of the Central Limit Theorem. The normal distribution of data is at the basis of the geostatistical approach (ordinary kriging); consequently, the statistical analysis described below was functional to the verification of this hypothesis, with the purpose of guaranteeing non-distorted results when applying the ordinary kriging, rather than the MMSD, method. The non-parametric Kolmogorov–Smirnov (K–S) test, with a 99% level of significance (Massey 1951; Lilliefors 1969), was applied to all the seasonal data of all the gauging stations, outlining the following results: in all the seasons nothing suggests the sample distributions are non-normal, or better, no meaningful differences were shown among the Gaussian distribution and the four sample distributions. The box and whiskers and the stem and leaf diagrams confirmed, even though at a qualitative level, that the frequency distributions for each of the considered seasons are approximately symmetrical. Figure 2 shows an overview of the four box and whiskers diagrams of every season and allows the

Fig. 2 Box and whiskers diagrams of the average rainfall values recorded at the 93 considered gauging stations of the Hydrographic Regional Office

Environ Monit Assess (2008) 145:375–386

shape and the position of the four distributions to be compared. These diagrams represent schematically the main characteristics of the distributions. In particular, the box represents the first and third quartiles and the median, while the whiskers represent the range between the first and the 99th percentile; outliers, outside this range, are also marked. Observing the diagrams in Fig. 2, it appears that, during the summer, the median is smaller than in autumn and winter, confirming that it rains less in warm seasons. A wide dispersion of rainfall values is, finally, evident around the median during the rainiest seasons, which is symptomatic of a great nonhomogeneity of the phenomenon over the territory. All this information, jointly with the results of the normality test, confirms the hypotheses made about the average behaviour of the considered phenomenon, which was, partially, already known. Following the preliminary phase of statistical investigation, the experimental variograms, representing the spatial behaviour of the mean seasonal rainfall rate for each season, were calculated and the theoretical models were determined. In the following Table 3, the parameters of the four theoretical variograms are reported; all of them were “spherical” and anisotropic. The reported ranges are those related to the principal axis of anisotropy.

Environ Monit Assess (2008) 145:375–386

381

Table 3 Characteristic parameters of the theoretical variograms for the four seasons

Quart Quart Quart Quart

1 2 3 4

Nugget (mm2)

Sill (mm2)

Range (m)

Anisotropy angle (deg)

Anisotropy ratio

40 10 5 40

130 80 60 200

65,000 65,000 70,000 65,000

150 150 150 150

2.5 3.3 5.0 2.0

From a comparison of the values in Table 3 it can be seen that: –

–

–

–

During the cold seasons (I and IV quarter) the scale of the considered phenomenon is wider; this indicates that the differences in rainfall values for different areas of the region are greater; this is most evident observing the dispersion of values around the median in the box and whiskers diagrams; Likewise, during these seasons the discontinuity at the origin (nugget) increases; this indicates that the rainfall values tend to differ even over short distances. This can be explained by spatial discontinuity, which is a specific peculiarity of rainfall phenomena. Similar values of the ranges point out, instead, that the distance of the spatial correlation remains nearly constant, around 65 km, throughout the four seasons; Throughout all the seasons there is a strong anisotropy (ratio from 2 to 5) of the phenomenon with a principal direction more or less parallel to the coast line.

The proposed methodology was applied and four different configurations were determined both for the OK and the MMSD criteria. As expected, the configurations produced by the two approaches are notably different. In fact, while for OK, the optimisation process uses seasonal variograms, the MMSD criterion involves only the locations of the existing 45 monitoring stations. Consequently, while the four configurations obtained by OK are really seasonal configurations, the four obtained by MMSD are simply the best four among several realizations. However, for a matter of clarity, in both the cases the four configurations have been labelled as “seasonal”. Figure 3 shows, as an example, the simulated configurations for the first season, using as objective function the average of the OK estimation variances

(Fig. 3a) and the average of the distances between an arbitrary point and its nearest neighbour (Fig. 3b); white circles indicate the new monitoring station locations. The grey area along the coast and the inner borders represents the part of the regional territory where the algorithm was constrained to avoid the placement of new monitoring points. Obviously, the second approach, being purely geometric, produced a new configuration, with a very regular distribution of the gauging stations. On the contrary, the approach based on the minimization of the average of the kriging estimation variances, produced a less regular configuration, but, more reliable from a spatial standpoint, in terms of estimation variance. The managers of the ARCCP preferred the approach based on the minimization of the average distance among the points since it allowed the spatial density of the gauging stations to be made consistent at a provincial level. Thus, the MMSD approach was chosen as the working criterion, simply on the basis of managerial requirements. The ARCCP also asked for the 50 new stations to be set up, divided into two groups: 15 fixed and 35 mobile stations. This request can be explained, in managerial terms, by the necessity of getting more detailed information from different parts of the Apulian territory according to the current season, with the purpose of defining and circulating reliable forecasts of rainfall availability, among the pooled consumers. The overlapping of the four seasonal configurations allowed a number of new station locations to be located automatically, which tended to remain unchanged season after season and certain others that, on the contrary, sometimes changed their coordinates considerably. The methodology, as described above, yielded the required number of locations where the new gauging stations could be placed. Nevertheless, this result was considered too rigid by the managers of the ARCCP, and they asked for a certain flexibility in the gauge positioning, since some locations were not achievable. Consequently, a possible maximum tolerance of 10 km was introduced between the determined and actual gauge position. In practice, a regular mesh of 10 km side cells was overlaid on the Apulia map and the locations of the 50 new gauging stations were associated to the correspondent mesh cell by means of GIS software (ESRI 1996). Doing so, it was possible

382

Environ Monit Assess (2008) 145:375–386

Fig. 3 Configurations of the new monitoring network of the Apulian Regional Consortium for Crop Protection resulting for the first season using a the average of estimated variances of ordinary kriging; b the average of the points distance from those at the border line

to distinguish two types of stations: those whose position remained fixed in the same mesh cell, throughout the four seasons and those whose position changed. Forty-two of the 50 new gauging stations belonged to the first group, and at least 17 of them were kept strictly to their original position, allowing a reduction of the position tolerance to 5 km. Only eight stations moved from one cell to another during the different monitoring seasons. These eight gauging stations were labelled as mobile stations, while the previous 17 were, obviously, considered fixed stations. The remaining 25 stations were defined as potentially mobile stations which means that even considering these stations as fixed stations, it would be possible to choose some mobile stations among them, to be moved to some other location in particular seasons and conditions. Figures 4 and 5 show graphically what was said above. In particular, Fig. 4 summarizes the four seasonal configurations; the squared boxes represent those cases where the gauging stations remain almost fixed throughout the year. Figure 5 shows the four final configurations of the new gauging stations per season, labelled according to type, achieved by means of the minimization of the mean of shortest distances method (MMSD). Finally, Table 4 reports the number of stations of the upgraded monitoring network per Province. The last column of Table 4 shows that the percent density of stations per Province has been rebalanced as required by the ARCCP.

Conclusion One of the main institutional assignments of the Regional Consortium for Crop Protection is the elaboration of data gained from the meteorological monitoring network with the purpose of obtaining information about the state of crops, hazards related to the actual and predicted meteorological conditions, and advising on agricultural practices to safeguard agricultural production and the environment. The precision of the predicted information and the climatological characterization of the territory are strongly conditioned by the optimal spatial arrangement of the monitoring stations. The present study holds particular importance also because it aims to improve the efficiency and effectiveness of the whole agro-meteorological monitoring system. The broadening of the agro-meteorological monitoring network is aimed at achieving more and more precise and reliable predictive information, able to satisfy the increasing need of knowledge regarding meteorological and climatic phenomena. An optimal location of the monitoring stations was established through the application of geostatistical methodologies, which followed a climatic characterization of the region based on a time series analysis of available data. In fact, the characterization of any spatial or time–space phenomena represents the first and most important step of any geostatistical study. Geostatistics methods, including kriging and cokrig-

Environ Monit Assess (2008) 145:375–386

383

Fig. 4 Final configuration from the elaborations in the four quarters of a year

ing techniques, were used to finalize the estimation of spatial variables that is intrinsically linked to the territory. In particular, kriging and its modifications besides providing an estimation of the considered spatial variable, also gives a measure of the precision of the estimation in terms of “estimation variance”. The proposed study highlighted, qualitatively and quantitatively, the variability of the considered phenomenon, specifying its typology, with regard to the presence of possible anisotropies and to the existence of different space or time scales of variability. A first phase of statistic analysis, on a quarterly level data, was followed by the computation of the experimental variograms and the variogram model fitting; the analysis of these variogram models (spherical), clearly highlighted some peculiar charac-

teristics of the Apulian climate, consisting in a great spatial variability of rainfall during the Winter, even over relatively small distances with a constant spatial correlation distance of about 65 km. This characteristic of seasonal variability was also confirmed by other approaches, including a qualitative analysis carried out by means of GIS instruments. Comparing the results obtained, it was possible to define the co-ordinates of the optimal locations where the 50 new monitoring stations should be placed. Subsequently, a distinction was made between fixed and mobile stations: those, among the 50 new stations, characterized by a strong convergence of the seasonal optimal locations were classified as fixed. The results of the present study were put into practice with the actual setting of the 50 monitoring

384

Environ Monit Assess (2008) 145:375–386

Fig. 5 Results from the elaborations for the 4 year seasons and groupings

stations in the suggested locations. This practical evolution of the methodology gives it an added value related to the possibility of continuously checking the efficiency of the proposed solution. Moreover, it allows further experiments to be made aimed at improving the methodology itself. One of the critical aspects of the proposed methodology is the choice between the MMSD and OK optimisation approach. Nevertheless, in our opinion it cannot be based on a stringent theoretical and quantitative criterion. In fact, MMSD is a geometrically driven method, and consequently produces spatially even distributions of monitoring point locations. OK, instead, is based on estimation variance minimization and produces a monitoring network capable of providing better estimations in non-sampled points. In short, the former method is more useful for designing “alert monitoring networks”, while the second is necessary when a reliable statistical description of a spatial phenomenon is required. A further planned improvement of the methodology consists in adopting a combined or mixed two- or multiple-step approach, which integrates the two approaches, so that the drawbacks of one method are partly compensated by the advantages of the other one.

A brand new feature of the proposed methodology consists in the possibility of designing a flexible monitoring network. Providing a criterion for distinguishing between “mobile” and “fixed” gauging stations allows the monitoring administrator to change the network configuration over a period of time, taking into account the seasonal behaviour of the considered natural phenomenon. Finally, a simplification was made with regard to installation costs. A total of 50 new stations was adopted in this paper, neglecting the trade-off between cost and accuracy of results, since the given budget

Table 4 New density of rain gauging stations, in each province, of the Apulian Regional Consortium for Crop Protection Province

Area (km2)

N° stations

Density (stations/km2)

Density (%)

Bari Brindisi Foggia Lecce Taranto Puglia

5,127,609 1,843,752 7,157,063 2,770,077 2,438,445 19,336,946

25 10 34 14 12 95

0.00488 0.00542 0.00475 0.00505 0.00492

19.52 21.51 19.12 20.32 19.52

Environ Monit Assess (2008) 145:375–386

available for improving the monitoring network was already defined by the ARCCP. Nevertheless, a further development of the methodology could be the possibility of introducing a satisfactory balance between costs and results, allowing a reliability threshold to be defined in terms of either monitoring station density in the MMSD case, or estimation variance in the OK case.

Acknowledgements The authors wish to acknowledge the courtesy of the Apulian Regional Consortium for Crop Protection (ARCCP) and Hydrographic Regional Office (HRO) in providing data used throughout the paper.

References 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. Arbia, G., & Espa, G. (2001). Optimal spatial sampling strategies for agricultural and environmental data. (Paper presented at the Conference on Agricultural and Environmental Statistical Applications, Rome – CAESAR), June. Ashraf, M., Loftis, J. C., & Hubbard, K. G. (1996). Application of geostatistics to evaluate partial weather station networks. Agricultural and forest meteorology, 84, 255–271. Bastin, G., Lorent, B., Duqué, C., & Gervers, M. (1984). Optimal estimation of the average areal rainfall and optimal selection of rain gauge locations. Water Resources Research, 20, 463–470. Bogardi, I., Bardossy, A., & Duckstein, L. (1985). Multicriterion network design using geostatistics. Water Resources Research, 21, 199–208. Bras, R. L., & Rodríguez-Iturbe, I. (1975). Rainfall-runoff as spatial stochastic processes. Data collection and synthesis. (Cambridge, USA: Ralph M. Parsons Lab. for Water Resources and Hydrodynamics, dep. of Civil Eng., Mass. Inst. of Technology). Bras, R. L., & Rodríguez-Iturbe, I. (1976). Network design for the estimation of areal mean of rainfall events. Water Resources Research, 12, 1185–1195. Carrera, J., & Szidarovzsky, F. (1985). Numerical comparison of network design algorithms for regionalized variables. Applied Mathematics and Computation, 16, 189–202. Cox, D. D., & Cox, L. H. (1994). Spatial sampling and the environment: report on a workshop. (Paper presented at NISS/USEPA workshop, Chapel Hill, NC), September. Deutsch, C. V., & Cockerham, P. W. (1994). Practical considerations in the application of simulated annealing to stochastic simulation. Mathematical Geology, 26, 67– 82. Deutsch, C., & Journel, A. (1992). GSLIB: geostatistical software library and user’s guide. New York: Oxford University Press.

385 ESRI. (1996). ArcView GIS. Redlands, CA: Environmental Systems Research Institute Inc. Fedorov, V. V., & Hackl, P. (1994). Optimal experimental design: spatial sampling. Calcutta Statistical Association Bulletin, 44, 73–174. Goovaerts, P. (1997). Geostatistics for natural resources evaluation. New York: Oxford University Press. Goovaerts, P. (2000). Geostatistical approaches for incorporating elevation into the spatial interpolation of rainfall. Journal of Hydrology, 228, 113–129. Harmancioglu, N. B., Alpaslan, M. N., Singh, V. P., Fistikoglu, O., & Ozkul, S. D. (1999). Water quality monitoring network design. Boston: Kluwer Academic Publishers. Isaaks, E. H., & Srivastava, R. M. (1989). An introduction to applied geostatistics. New York: Oxford Unversity Press. Journel, A. G., & Huijbrechts, C. J. (1978). Mining geostatistics. London: Academic Press. Knopman, D. S., & Voss, C. I. (1989). Multiobjective sampling design for parameter estimation and model discrimination in groundwater solute transport. Water Resources Research, 25, 2245–2258. Lebel, T., Bastin, G., Obled, C., & Creutin, J. D. (1987). On the accuracy of areal rainfall estimation. A case study. Water Resources Research, 23, 2123–2134. Lilliefors, H. W. (1969). On the Kolmogorov–Smirnov test for the exponential distribution with mean unknown. Journal of the American Statistical Association, 64, 387–389. Massey Jr, F. J. (1951). The Kolmogorov–Smirnov test of goodness of fit. Journal of the American Statistical Association, 46, 68–78. Matheron, G. (1970). The theory of regionalized variables and its applications. Paris: Les Cahiers de Morphologie Mathematique. Fascicule I, Ecole des Mines de Paris. Metropolis, N., Rosenbluth, A., Rosenbluth, M., Teller, A., & Teller, E. (1953). Equation of state calculations by fast computing machines. Journal of Chemical Physics, 21, 1087–1092. Meyer, D., Valocchi, A. J., & Eheart, J. W. (1994). Monitoring network design to provide initial detection of groundwater contamination. Water Resources Research, 30, 2647– 2659. Nunes, L. M., Cunha, M. C., & Ribeiro L. (2002). Monitoring network optimisation using both spatial and temporal information. (Presented at the 3rd International Conference on Decision Making in Urban & Civil Engineering, London), November. Ott, W. R. (1995). Environmental Statistics and data analysis. Boca Raton, FL: Lewis Publisher. Papamichail, D. M., & Metaxa, I. G. (1996). Geostatistical analysis of spatial variability of rainfall and optimal design of a rain gauge network. Water Resources Management, 10, 107–127. Pardo-Igùzquiza, 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. Rodríguez-Iturbe, I., & Mejía, J. M. (1974). The design of rainfall networks in time and space. Water Resources Research, 10, 713–728. Van Groenigen, J. W. (2000). SANOS 0.1 for Windows 95, preliminary version.

386 Van Groenigen, J. W., Pieters, G., & Stein, A. (2000). Optimizing spatial sampling for multivariate contamination in urban areas. Environmetrics, 11, 227–244. Van Groenigen, J. W., & Stein, A. (1998). Constrained optimization of spatial sampling using continuous simulated annealing. Journal of Environmental Quality, 27, 1078–1086. Van Groenigen, J. W., & Stein, A. (2000). Constrained optimization of spatial sampling in a model-based setting

Environ Monit Assess (2008) 145:375–386 using SANOS software. In G. B. M. Heuvelink, & M. J. P. M. Lemmens (Eds.) Accuracy 2000 – Proceedings of the 4th international symposium on spatial accuracy assessment in nature resources and environmental sciences (pp. 679–687). Delft: Delft Univ. Press. Wu, Y. (2004). Optimal design of a groundwater monitoring network in Daqing, China. Environmental Geology, 45, 527–535.