Factorial and indicator kriging
methods using a geographic
information system to delineate
spatial variation and pollution
sources of soil heavy metals
Yu-Pin LinÆ Tsun-Kuo Chang Æ Chiung-Wen Shih Æ Chen-Hui Tseng Abstract This study uses factorial and indicatorkriging methods with a geographic information system (GIS), to overlay and delineate the variation and pollution sources of soil heavy metals in fields in the north of Changhua County, Taiwan. Accord-ingly, the spatial variation and pollution probability patterns of soil cadmium, chromium and copper are estimated to support further environmental moni-toring, remediation and planning. The overlay maps of industrial plants, irrigation channels and kriging estimates, constitute a GIS. This reveals that the short- and long-range variations and soil pollution by Cd, Cr and Cu strongly corresponded to the locations of irrigation systems and industrial plants in the area of study. The maps also verify most locations of suspected pollution sources from the GIS database. Moreover, the probability spatial patterns estimated by indicator kriging, based on the pollution thresholds of Cd, Cr and Cu, were also highly correlated with irrigation channels and industrial plants. The GIS, with factorial and indi-cator kriging under different combinations of vari-ations and pollution probabilities of soil Cd, Cr and
Cu, provide alternative scenarios for delineating the areas that contain high variation and high pollution probability of these heavy metals, to enable decision makers to identify pollution sources, and to further investigate and remedy soil pollution.
Keywords Factorial kriging Æ Indicator kriging Æ Spatial variation Æ Soil
pollution Æ GIS Æ Taiwan
Introduction
Nature and human activity are both important in deter-mining the complex spatial variation of heavy metal con-centrations in soil, via soil parent material mineralogy and waste from industrial plants. Anthropogenic activities cause the most soil pollution. Fully understanding, esti-mating and mapping spatial variations and pollution of soil heavy metals by efficient techniques enables accurate monitoring and remediation. In Taiwan, some soil in ag-ricultural fields has been polluted by waste from industrial plants through irrigation systems. Some industrial plants were suspected as pollution sources. In 1983, the
Environmental Protection Administration (EPA) of the Republic of China initiated a collaborative research programme to determine the presence of As, Cd, Cu, Cr, Hg, Ni, Pb and Zn trace elements in the soil of Taiwan (Lin and Chang 2000). The programme also sought to deter-mine soil properties, such as cation-exchange capacity and pH value. Moreover, in 1998, Taiwan’s EPA also
announced a soil and groundwater pollution act to over-come soil and groundwater pollution in Taiwan. Suspected pollution sources of soil heavy metals were investigated and identified by Taiwan EPA.
Kriging techniques have been widely applied in soil science to analyse the spatial patterns and variability of concen-trations of pollutants. Kriging estimates are weighted sums of the adjacent sample concentrations and weights depend on the correlation structure. The weights are chosen to minimize the variance in estimates of the variable of in-terest. In this framework, kriging estimates can be regarded as the most accurate linear estimator (Best Linear Unbiased
Received: 8 February 2002 / Accepted: 9 April 2002 Published online: 4 June 2002
ª Springer-Verlag 2002
Y.-P. Lin (&)
Department of Landscape Architecture, Chinese Culture University, 55 Hwa-Ken Rd, Yangming Shan, Taipei 11114,
Taiwan, Republic of China E-mail: [email protected] Tel.: +886-2-28626433 Fax: +886-2-28617507 T.-K. Chang
Graduate Institute of Bioenvironmental Systems Engineering, National Taiwan University,
Taipei 10617, Taiwan, Republic of China C.-W. ShihÆ C.-H. Tseng
Graduate Institute of Geography, Chinese Culture University, 55 Hwa-Ken Rd. Yangming Shan,
Estimator). At an unsampled location and for a given variogram, a kriging estimate can simply be considered to be an optimally weighted average of the surrounding sampled data (Cressie 1990). Recent examples of kriging in soil pollution include those of Couto and others (1997), White and others (1997), Juang and Lee (1998a), Meuli and others (1998), Wang (1998), Chang and others (1999), Wang and Zhang (1999), Zhang and others (1999), Lin and Chang (2000), Carlon and others (2001), Facchinelli and others (2001) and Lin and others (2001).
Factorial kriging is a variant of kriging that aims to estimate and map the different sources of spatial variability deter-mined from the experimental variogram (Goovaerts 1992, 1998). Factorial kriging allows the various sources of spatial variation in multivariate data to be distinguished according to the scale on which they operate (Goovaerts 1998). Therefore, variograms decompose the spatial components of interesting variables into structural components. The spatial variations caused by different processes can be dis-tinguished by factorial kriging. The most recent applica-tions of factorial kriging include those of Goovaerts (1994), Webster and others (1994), Dobermann and others (1995), Einax and Soldt (1998), Jime´nez-Espinosa and Chica-Olmo (1999), Bocchi and others (2000), Castrignano` and others (2000a), Batista and others (2001) and Lin (2002). Soil pollution data may contain highly skewed and ex-treme values. Applying indicator kriging, a nonparametric, distribution-free estimation method, may solve these problems. In indicator kriging, no assumptions are made about the underlying invariant distribution, and 0–1 indicator transformations of the data make the predictor robust to outliers (Cressie 1993). The numbers estimated by indicator kriging represent the probability that the unknown values are greater or less than a specified threshold. That is, the expected value at an unsampled location derived from indicator data is equivalent to the cumulative distribution of the variable (Smith and others 1993). Recent works on indicator kriging for soil pollution included those of Smith and Williams (1996), Juang and Lee (1998b), Jime´nez-Espinosa and Chica-Olmo (1999), Castrignano` and others (2000b) and Van Meirvenne and Goovaerts (2001).
GIS, a computer system, can capture, store, check, inte-grate, manipulate, analyse and display geographically ref-erenced information, i.e. data identified according to their locations (Zhu and others 2001). With environmental models, GIS can be an efficient decision system for envi-ronmental management and planning. This study uses factorial and indicator kriging methods with a GIS inte-grating land use, irrigation system and suspected pollution sources data, to overlay, delineate and verify the variation and pollution sources of soil heavy metals in fields in the northern part of Changhua County, Taiwan.
Materials and methods
The county from which the samples were selected, Changhua County, is one of the most important
agricul-tural counties of Taiwan. Figure 1 shows the 194 sampling sites in the northern part of Changhua County and selected from the 878 sites across Taiwan.
The GIS database for this study area included soil sample, land use, irrigation system, suspected pollution sources, spatial variation components of soil heavy metals and probabilities exceeding pollution thresholds. Figure 2 il-lustrates a linkage of this database and integration meth-odologies. Soil heavy metals data in this study, 194 topsoil samples of Cd, Cr and Cu, taken at depths of 0–15 cm and collected between 1981 and 1986, were selected from those of the EPA data programme described above. These sam-ples were taken from geographically distributed sites at target intervals of 100 m and 50 m in a network formation. Table 1 shows descriptive statistics of these 194 topsoil samples of Cd, Cr and Cu. Taiwan’s EPA also defined five classes of soil heavy metals (Table 2). The first and second classes of soil heavy metal samples are considered to in-volve no soil heavy metal pollution. The soil heavy metal concentrations of the third class are defined as back-ground values. The fourth and fifth soil heavy metal classes require intensive monitoring and consideration of reme-dial action. Among these 194 samples, 142, 53, 35, 152 and 23 samples fell into the fourth class for Cd, the fourth class for Cr, the fifth class for Cr, the fourth class for Cu and the fifth class of Cu respectively.
The land use map, as a vector and polygon database, was digitized in scale 1:5,000 and developed by the Food and Agriculture Department of the Council of Agriculture of the Republic of China and the Department of Land Ad-ministration of the Ministry of the Interior of Republic of China in 1984. The irrigation channel map, as a vector and polyline database, was digitized by our research group
Fig. 1 Sampling points
based on 1:1,200 scale cadastral maps provided by the Irrigation Association of Changhua County in 1998. Figure 3a depicts the industrial and built-up areas and the irrigation system in the study area from GIS software Arcview 3.0. The suspected pollution sources map is the point database at scale 1: 5,000 and was digitized by Taiwan EPA (Fig. 3b). Most of the industrial plants in this study area served metalwork, electroplating and metal surface treatment industries (Fig. 3b). These industrial plants have been suspected of discharging wastewater into irrigation channels in this study area.
In this study, the variogram and cross-variogram models of normalized Cd, Cr and Cu and the indicator variograms of these three soil heavy metals with various thresholds were fitted within VARIOWIN 2.2 (Pannatier 1996). The values of the fourth and fifth classes (Table 2) were selected as the thresholds for these soil heavy metals in indicator kriging. Factorial kriging and indicator kriging were performed by modified kriging programs in GSLIB (Geostatistical Software Library; Deutsch and Journel 1992). Each spatial component and probability map based on thresholds of soil heavy metals determined in kriging, was estimated with over 1,064 (50·50 m) cells in 38 col-umns and 28 rows. The estimates of the spatial compo-nents and the exceeding threshold probabilities of these
Fig. 2
Linkage between the GIS database, kriging methods and cartography
Table 1
Statistics of soil heavy metals. Units in mg/kg
Cd Cr Cu Mean 0.743 10.433 54.455 SD 0.575 13.070 67.099 Minimum 0.120 0.650 11.370 Maximum 3.600 76.690 556.000 25th 0.380 2.708 23.000 50th 0.570 5.120 32.845 75th 0.878 12.000 61.043 Table 2
Soil class. Units in mg/kg Soil heavy metals 1 2 3 4 5
As – <4 4–9 10–60 >60 Cd – <0.05 0.05–0.39 0.40–10 >10 Cr – <0.10 0.10–10 11–16 >16 Cu <1 1–11 12–20 21–100 >100 Hg – <0.10 0.10–0.39 0.40–20 >20 Ni – <2 2–10 11–100 >100 Pb – <1 1–15 16–120 >120 Zn <1.5 1.5–10 11–25 26–80 >80 Fig. 3
heavy metals were input into Arcview 3.0 (ESRI 1998) to map and identify the patterns of spatial variation and sources of soil heavy metals, Cd, Cr and Cu, at different scales over the area of interest. Moreover, for each soil heavy metal, the different combinations of soil variations estimated by factorial kriging and the probability of ex-ceeding thresholds (fourth and fifth classes as showed in Table 2) estimated by indicator kriging were mapped and extracted using Arcview 3.0 to delineate and verify pollu-tion sources and the area of soil pollupollu-tion.
Geostatistics
Geostatistical methods are based on the regionalized variables theory that states that variables in an area have both random and spatial properties (Journel and
Huijbregts 1978). A geostatistical variogram of data should first be determined. A variogram quantifies the commonly observed relationship between the values of data pertain-ing to the samples and the samples’ proximity. The vari-ogram, c(h), is defined as,
cð Þ ¼ 1=2h ð ÞVar Z x½ ð Þ Z x þ hð Þ ð1Þ where h denotes the lag distance between pairs of points, Var represents the variance of the argument, Z(x) is the value of the regionalized variable of interest at location x, and Z(x + h) represents the value at location x + h. An experimental variogram, c(h), is given by, cð Þ ¼ 1= 2n hh ½ ð ÞX
n hð Þ
i¼1
Z xið þ hÞ Z xið Þ
½ 2 ð2Þ
where c(h) denotes the variogram for an interval lag distance class h, and n(h) represents the number of pairs separated by lag distance h.
Similarly, the spatial correlations, called cross-variograms, between two variables (a and b) can be written as, cabð Þ ¼ 1= 2n hh ½ ð ÞX n hð Þ i¼1 Zaðxiþ hÞ Zað Þxi ½ Zb ðxiþ hÞ Zbð Þxi ð3Þ An indicative goodness of fit (IGF), a quantitative visual fit technique, is calculated every time the two-dimensional nested model is modified. IGF is calculated with the fol-lowing equation (Pannatier 1996):
IGF ¼ 1 N XN k¼1 Xn kð Þ i¼0 P ið Þ P nðkÞ j¼0 P jð Þ D kð Þ d ið Þ cð Þ ^i ccð Þi r2 2 ð4Þ
where N denotes the number of directional variograms, n(k) represents the number of lags relative to variogram k, D(k) is the maximum distance relative to variogram k, P(i) denotes the number of pairs for lag i of variogram k, d(i) represents the mean pair distance for lag i of variogram k, c(i) is the experimental measure of spatial continuity for lag i, ^ccð Þ denotes the modelled measure of spatial conti-i
nuity for d(i), r2represents the covariance of data for the cross-variogram and the cross-variance, the maximum absolute experimental value of all measures for the cross-madogram, the correlation coefficient for the cross-cor-relogram and the cross-standardized variogram. The IGF is a number of standardization without units; a value close to zero indicates a good fit (Pannatier 1996). Therefore, the IGF value is comparable from one to another, allowing the user of VARIOWIN 2.2 to numerically check how well models fit the experimental measures (Pannatier 1996).
Factorial kriging
The multivariate regionalization of a set of random func-tions can be represented in a spatial, multivariate linear model, which allows easy manipulation of multivariate data (Wackernagel 1995). The nested cross-variogram can thus be modelled as linear combinations:
cabð Þ ¼h X Ns u¼1 cuabð Þ ¼h X Ns u¼1 buabguð Þh ð5Þ
where Nsis the number of the spatial scale, buab are coef-ficients, and gabuð Þ are variogram functions.h
A set of random functions, Zi(x), can be decomposed into sets of spatially uncorrelated factors (Rouhani and Wa-ckernagel 1990; Goovaerts 1992; WaWa-ckernagel 1995). The co-kriging estimator of the spatial component Zku at point x0 is, Zuk ð Þ ¼x0 X m i¼1 Xn j kjiZi xj ð6Þ
The kriging system can be solved using, Xm g¼1 Xn s¼1 ksgcig xj xs li ¼ buikg u xj xo andX n s¼1 ksi¼ 0 ð7Þ
where li is the Largrange multiplier; gu(xj,x0) is the value taken by the uth basic variogram function, gu(h), between the ath sampling point and x0; i = 1,..., p, and j = 1, ..., n.
Indicator kriging
Indicator kriging estimates the probability of exceeding specific threshold values, zk, at a given location. In indi-cator kriging, the data, z(x), are transformed into an in-dicator variable as follows.
i x; zkð Þ ¼ 1; if z xð Þ zk 0; otherwise
ð8Þ
At an unsampled location, x0, the probability that z(x) £ zk can be estimated using a linear combination of neigh-bouring indicator variables. This ordinary indicator kri-ging estimator is,
Pr ob z x0½ ð Þ zk=ð Þn ¼X n
a¼1
where i(xa;zk) represents indicator values at xa, a=1,ÆÆÆ,n, and ka, determined by solving the following kriging sys-tem, is the kriging weight of i(xa;zk) used in estimating Prob[z(x0) £ zk/(n)].
An ordinary indicator kriging system can be solved using, Xn b¼1 kbci xa xb;zk þ l ¼ ciðxa x0;zkÞ andX n b¼1 kb¼ 1 ð10Þ
where l is the Largrange multiplier; ci(xa–xb;zk) is the indicator variogram between indicator variables at the ath and bth sampling points ; ci(xa–x0;zk) is the variogram between indicator variables the ath sampling point and x0, and a = 1,..., n.
Results and discussion
Spatial correlation and modelFollowing Lin’s work (2002) on variography, this study determines experimental variograms, cross-variograms for normalized (zero mean and unit variance) Cd, Cr and Cu soil heavy metals, and indicator variograms, according to the thresholds (fourth and fifth soil classes in Table 2) for these heavy metals with a single active lag (1,200 m) and lag interval (120 m). Moreover, a relatively consistent set of good fit models was established to fit these experimental variograms, using VARIOWIN 2.2 (Pannatier 1996). The variogram model with the lowest value of IGF is selected as a good fit model from all sought models. The good fit variogram and cross-variogram models of these three heavy metals (Lin 2002), listed in Table 3, were sums of three structures, according to the following equation [Eq. (11)]. The equation includes a nugget effect term, and two spherical models with ranges of 288 and 1,200 m for short and long spatial scales respectively.
where a, b = 1,..., 7. b0
abis the nugget effect; b1abis the sill of the short-range spherical structure and b2
abis the sill of the long-range spherical structure.
These results demonstrate that experimental variograms and cross-variograms of soil heavy metals Cd, Cr and Cu can be modelled by three spatial structures that include the nugget effect and two spherical models – short-range and long-range (Lin 2002). These nest variogram (short-range and long-range) structures for soil heavy metals may be dominated by the output of waste from industrial plants into nearby fields and, through irrigation channels, into far fields.
These co-regionalization analyses also reveal that the variogram models of normalized Cd, Cr and Cu with higher nugget effects (0.370–0.480) and higher nugget effect ratios (34.30–46.20%) display higher small-scale variations, possibly because of waste from the industrial plants. The cross-variogram models also reveal that spatial correlations among these three heavy metals showed large nugget effects (0.364–0.384) and high nug-get effect ratios (46.50–61.20%). Variograms and cross-variogram models of Cd, Cu and Cd–Cu are shown as examples in Fig. 4.
Similarly, the indicator variogram models of Cd, Cr and Cu, based on the forth and fifth classes, were fitted in VARIOWIN 2.2, as listed in Table 4. These variogram models revealed three spatial structures, including the nugget effect and two spherical models – short-range and long-range. These indicator variogram models included high nugget effect ratios, ranging from 41.52 to 87.04%. Moreover, the nugget effects of the indicator variogram models of Cd, Cr and Cu, based on the fourth class thresholds, were less than those based on the fifth class thresholds, indicating that higher pollution values cor-responded to higher small-scale variations of Cd, Cr and Cu. These variations might be primarily determined by the locations of the irrigation channels and industrial plants in the study area. Figure 5 shows indicator vari-ogram models of Cd and Cr as examples of indicator variogram. cabðhÞ ¼ bo abþ b1ab 32 h 288 1 2 h 288 3 h i þ b2 ab 32 h 1200 1 2 h 1200 3 h i for 0 < h < 1200 cabðhÞ ¼ b0 abþ b1abþ b2ab for h > 1200 ( ð11Þ Table 3
Variogram and cross-variogram models of standardized Cd, Cr and Cu
Standardized Short range Long range Nugget effect Short-range sill Long-range sill Ratio of nugget IGF
variable (m) (m) (c0) (c1) (c2) c0/(c0+c1+c2) Cd 288 1,200 0.370 0.240 0.470 0.343 1.86E-03 Cd–Cr 288 1,200 0.364 0.210 0.021 0.612 4.58E-03 Cd–Cu 288 1,200 0.368 0.152 0.272 0.465 1.56E-03 Cr 288 1,200 0.480 0.360 0.200 0.462 1.63E-03 Cr–Cu 288 1,200 0.384 0.272 0.160 0.471 1.55E-03 Cu 288 1,200 0.480 0.200 0.370 0.457 9.00E-04
Fig. 4
Variogram and cross-variogram models of normalized a Cd, b Cu, c Cd–Cu
Table 4
Indicator variogram models of Cd, Cr and Cu
Indicator Soil class Short range Long range Nugget effect Short range sill Long range sill Ratio of nugget IGF
variable (m) (m) (c0) (c1) (c2) c0/(c0+c1+c2) Cd 4 288 1,200 0.088 0.070 0.054 0.415 3.08E-03 Cr 4 288 1,200 0.148 0.054 0.008 0.705 1.89E-03 Cr 5 288 1,200 0.131 0.022 0.010 0.803 1.14E-03 Cu 4 288 1,200 0.092 0.054 0.012 0.581 2.87E-03 Cu 5 288 1,200 0.094 0.010 0.004 0.870 1.48E-03 Fig. 5
Indicator variograms of indicator variables based on a Cd>0.4 mg/kg, bCr>11 mg/kg, c Cr>16 mg/kg
Spatial variation and sources of pollution Figures 6 and 7 show the short-range and long-range spatial components of Cd, Cr and Cu. The areas with high Cd short-range variation were in the east and south parts of the study area. The areas with large short-range varia-tions of Cr and Cu were northeast and northwest of the study area. These spatial component maps of Cd, Cr and Cu, exhibiting spatial patterns similar to those of short-range variations, reveal that high spatial variations occur around the industrial plants, as shown in Fig. 6, revealing that sources of short-range variations of soil Cd, Cr and Cu closely correspond to the positions of local industrial plants.
Figure 7 illustrates that the long-range variations of Cd, Cr and Cu all showed similar spatial patterns, except for Cr south of the study area. These spatial patterns of compo-nents Cd, Cr and Cu were related not only to the location of industrial plants but also to that of the irrigation channels shown in Fig. 7. The overlaid maps (Figs. 6 and 7) of spatial variations of Cd, Cr and Cu, industrial plants and irrigation channels also reveal that these soil heavy metals were contaminated by waste water from industrial (metalwork, electroplating and metal surface treatment) plants, transported to the fields by irrigation systems. Figures 6 and 7 also confirm locations of suspected pol-lution sources in northeast and south of this study area.
Figures 8 and 9 show indicator kriged probability maps based on the thresholds of Cd, Cr and Cu. Figure 8a in-dicates that the areas with a probability greater than 0.6 for exceeding the threshold of the fourth class (0.4 mg/kg) of Cd are strongly related to the location of irrigation chan-nels and industrial plants. Most areas south of the study area showed a probability greater than 0.6 for exceeding the thresholds of the fourth and fifth classes (11 and 16 mg/kg) of Cr (Fig. 8b, c). Figure 9 shows that the areas with a probability above 0.6 for exceeding the threshold of the fourth class (21 mg/kg) of Cu almost covered this entire study area. The probability map for the threshold of the fifth class (100 mg/kg) of Cu shows that the area with probability above 0. 8 is located in the northeast as shown in Fig. 9.
Areas of high variation and pollution
GIS was used to find the areas of high variation (>0.5) and high probability (>0.8) for exceeding the thresholds of Cd, Cr and Cu. Most areas with high short-range variation and high probability for exceeding 0.4 mg/kg of Cd in the south and northeast are highly correlated with the location of industrial plants of the study area as shown in Fig. 10a. Moreover, the high long-range variation and high proba-bility for exceeding 0.4 mg/kg of Cd areas are located northeast and northwest of the study area as shown in Fig. 10a. The areas showing high short-range variation, long-range variation and probability are located in the
Fig. 6
Short-range spatial components of a Cd, b Cr, c Cu
Fig. 7
northeast as shown in Fig. 10a. Thus, the areas, irrigation channels and industrial plants in the northeast, south and northwest may be considered as first and second priorities for further environmental monitoring for soil Cd. Fig-ure 10a also indicates that these three areas confirm most locations of suspected pollution sources, except the northwest.
Figure 10b, c shows a few areas in which the concentration of Cr exceeded the fourth and fifth classes of Cr, with high variations. The few areas with short-range variation >0.5, long-range variation >0.5 and high probabilities (>0.8) for concentrations to exceed those of the fourth (11 mg/kg) and fifth (16 mg/kg) classes of Cr are randomly distributed across the northeast and south as shown in Fig. 10b, c. Figure 11a shows that the areas with a high probability of Cu exceeding the threshold for the fourth class (21 mg/kg) of Cu and a high short-range variation are in the northeast and northwest. Figure 11a reveals that most areas in the northeast showed high probability of Cu exceeding the thresholds for the fourth class (21 mg/kg) of Cu and high long-range variation. The areas with a short-range varia-tion, long-range variation and high probabilities of Cu exceeding the fourth (21 mg/kg) class of Cu are primarily distributed in the northeast as shown in Fig. 11a.
Fig. 8
Kriged probabilities of a Cd>0.4 mg/kg, b Cr>11 mg/kg,
cCr>16 mg/kg Areas with high variation and high probabilities of a Cd>0.4 mg/kg,Fig. 10 bCr>11 mg/kg, c Cr>16 mg/kg (SLV short- and long-range variation; LV long-range variation; SV short-range variation)
Fig. 9
Figure 11a also reveals that these two areas verify most locations of suspected pollution sources, except the northwest.
Few areas in the northeast and northwest showed both a high probability of Cu exceeding the thresholds of the fifth class (100 mg/kg) of Cu and a high short-range variation as shown in Fig. 11b. In the northeast, few areas showed both a high probability of Cu exceeding the thresholds for the fifth class (100 mg/kg) of Cu and high variation as depicted in Fig. 11b. These results reveal that the northeast must be intensively monitored to remedy soil Cu. Fur-thermore, the soil, irrigation channels and industrial plants in the northeast of the study area should be prior-itized for further Cd, Cr and Cu monitoring and remedi-ation. A new monitoring programme should be established in the northwest of this study area for identifying pollution sources.
Conclusions
This study has determined the spatial variations in mea-sured Cd, Cr and Cu concentrations at various spatial scales, and the probability of soil pollution by these heavy metals over the area of interest. Factorial kriging estab-lished spatial component maps of soil heavy metals at short (288 m) and long (1,200 m) range. The maps clearly demonstrated that, at these two scales, variations of Cd, Cr and Cu were strongly related to the location of industrial plants and irrigation systems in the study area. The maps also confirm most locations of suspected pollution
sources. Indicator kriging mapped the areas in which soil had a high probability of exceeding the thresholds of Cd, Cr and Cu. Moreover, a GIS can be applied to combine factorial kriging and indicator kriging to determine the areas of a high variation at different scales and a high probability of soil heavy metal concentrations exceeding thresholds, to provide various scenarios to support further monitoring and remediation of pollution by soil heavy metal Cd, Cr and Cu.
References
Batista AC, Sousa AJ, Batsta MJ, Viegas L (2001) Factorial kriging with external drift: a cause study on the Penedono Region, Portugal. Appl Geochem 16:921–929
Bocchi S, Castrignano` A, Fornaro F, Maggiore T (2000) Appli-cation of factorial kriging for mapping soil variation at field scale. Eur J Agron 13:295–308
Carlon C, Critto A, Marcomini A, Nathanail P (2001) Risk based characterisation of contaminated industrial site using multivariate and geostatistical tools. Environ Pollut 111(3): 417–427
Castrignano` A, Giugliarini L, Risaliti R, Martinelli N (2000a) Study of spatial relationships among some soil physico-chemi-cal properties of a field in central Italy using multivariate geo-statistics. Geoderma 97:39–60
Castrignano` A, Goovaerts P, Lulli L, Bragato G (2000b) A geo-statistical approach to estimate probability of occurrence of Tuber melanosporum in relation to some soil properties. Geo-derma 98:95–113
Chang TK, Shyu GS, Lin YP, Chang NC (1999) Geostatistical analysis of soil arsenic content in Taiwan. J Environ Sci Health Part A 34(7):1485–1501
Couto EG, Stein A, Klamt E (1997) Large area spatial variability of soil chemical properties in central Brazil agriculture. Ecosyst Environ 66(2):139–151
Cressie N (1990) The origins of kriging. Math Geol 22(2):239–252 Cressie N (1993) Statistics for spatial data. Wiley, New York Deutsch CV, Journel AG (1992) Geostatistical software library and
user’s guide. Oxford University Press, New York
Dobermann A, Goovaerts P, George T (1995) Sources of soil variation in an acid Ultisol of the Philippines. Geoderma 68:173–191
Einax JW, Soldt U (1998) Multivariate geostatistical analysis of soil contaminations. Fresenius J Anal Chem 361:10–14 ESRI (Environmental systems Research Institute) (1998) Arcview
3.0. ESRI, Redlands, California
Facchinelli A, Sacchi E, Mallen L (2001) Multivariate statistical and GIS-based approach to identify heavy metal sources in soils. Environ Pollut 114 313–324
Goovaerts P (1992) Factorial kriging analysis: a useful tool for exploring the structure of multivariate spatial soil information. J Soil Sci 43:597–619
Goovaerts P (1994) Study of spatial relationships between two sets of variables using multivariate geostatistics. Geoderma 62:93–107
Goovaerts P (1998) Geostatistical tools for characterizing the spatial variability of microbiological and physico-chemical soil properties. Bio Fertil Soils 27:315–334
Jime´nez-Espinosa R, Chica-Olmo M (1999) Application of geo-statistics to identify gold-rich areas in the Finisterre-Fervenza region, NW Spain. Appl Geochem 14(1):133–145
Journel AG, Huijbregts CJ (1978) Mining geostatistics. Academic Press, San Diego
Fig. 11
Areas with high variation and high probabilities of a Cu>21 mg/kg, bCu>100 mg/kg (SLV short- and long-range variation; LV long-range variation; SV short-range variation)
Juang KW, Lee DY (1998a) A comparison of 3 kriging methods using auxiliary variables in heavy-metal contaminated soils. J Environ Qual 27:355–363
Juang KW, Lee DY (1998b) Simple indicator kriging for esti-mating the probability of incorrectly delineating hazardous ar-eas in a contaminated site. Environ Sci Technol 32:2487–2493 Lin YP, Chang TK (2000) Geostatistical simulation and estimation
of the spatial variability of soil zinc. J Environ Sci Health Part A 35(3):327–347
Lin YP, Chang TK, Teng TP (2001) Characterization soil lead by comparing of sequential Gaussian simulation, simulated an-nealing simulation and kriging methods. Environ Geol 41(1/ 2):189–199
Lin YP (2002) Multivariate geostatistical methods to identify and map spatial variations of soil heavy metals. Environ Geol (in press)
Meuli R, Sculin R, Webster R (1998) Experience with the repli-cation of regional survey of soil pollution. Environ Pollut 101(3):311–320
Pannatier Y (1996) VARIOWIN Software for spatial data analysis in 2D. Springer, Berlin Heidelberg New York
Rouhani S, Wackernagel H (1990) Multivariate geostatistical approach to space-time data analysis. Water Resour Res 26:585–591
Smith M, Williams R (1996) Examination of methods for evalu-ating remining a mine waste site. Part II. Indicator kriging for selective remediation. Eng Geol 43(1):23–30
Smith JL, Halvorson JJ, Papendick RI (1993) Using multiple variable indicator kriging for evaluating soil quality. Soil Sci Soc Am J 57:743–749
Van Meirvenne M, Goovaerts P (2001) Evaluating the probability of exceeding a site-specific soil cadmium contamination threshold. Geoderma 102:75–100
Wackernagel H (1995) Multivariate geostatistics: an introduction with applications. Springer, Berlin Heidelberg New York Wang XJ (1998) Kriging and heavy-metal pollution assessment in
waste water irrigated agricultural soil of Beijing eastern farming regions. J Environ Sci Health Part A 33(6):1057–1073
Wang XJ, Zhang ZP (1999) A comparison of conditional simu-lation, kriging and trend surface analysis for soil heavy metal pollution pattern analysis. J Environ Sci Health Part A 34(1):73– 89
Webster R, Atteia O, Dubois JP (1994) Coregionalization of trace metals in the soil in the Swiss Jura. Eur J Soil Sci 45:205–218 White JG, Welch RM, Norvell WA (1997) Soil zinc map of the USA using geostatistics and geographic information systems. Soil Sci Soc Am J 61:185–194
Zhang R, Shouse P, Yates S (1999) Estimates of soil nitrate dis-tributions using cokriging with pseudo-crossvariogram. J En-viron Qual 28:424–428
Zhu HC, Charlet JM, Poffijn A (2001) Randon risk mapping in southern Belgium: an application of geostatistical and GIS techniques. Sci Total Environ 272:203–210
