JEQ Grow Your Career With ASA
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 3 January 2006
Published in J Environ Qual 35:21-36 (2006)
DOI: 10.2134/jeq2004.0389
© 2006 American Society of Agronomy, Crop Science Society of America, and Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in ISI Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via ISI Web of Science (2)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Vanderlinden, K.
Right arrow Articles by Giráldez, J. V.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Vanderlinden, K.
Right arrow Articles by Giráldez, J. V.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by Vanderlinden, K.
Right arrow Articles by Giráldez, J. V.
Related Collections
Right arrow Mining Wastes
Right arrow Spatial Variability
Right arrow Ecosystem Restoration
Right arrow Heavy Metals
Right arrow Soil Pollution

TECHNICAL REPORTS

Heavy Metals in the Environment

Mapping Residual Pyrite after a Mine Spill Using Non Co-Located Spatiotemporal Observations

Karl Vanderlindena,*, Rafaela Ordóñezb, Maria J. Poloc and Juan V. Giráldezc

a CIFA "Las Torres-Tomejil", I.F.A.P.A.-C.I.C.E. (Junta de Andalucía), Ctra. Sevilla-Cazalla, km 12,2, 41200 Alcalá del Río, Sevilla, Spain
b CIFA "Alameda del Obispo", I.F.A.P.A.-C.I.C.E. (Junta de Andalucía), 14080 Córdoba, Spain
c Dep. of Agronomy, University of Córdoba, Avda. Menéndez Pidal, s/n, 14080 Córdoba, Spain

* Correfsponding author (karl.vanderlinden.ext{at}juntadeandalucia.es)

Received for publication October 18, 2004.

   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Monitoring of soil chemical properties for pollution assessment generally requires destructive soil sampling and results in spatiotemporal datasets where data from different sampling dates are non co-located. The objective of this study was to assess the spatial distribution of residual pyrite sludge at a reclaimed site, using temporally non co-located data on pH; soil oxidizable fraction (SOF); and EDTA-extractable Fe, Zn, and Cu from six different sampling dates over a period of 2 yr. During this period spatially averaged pH and Zn concentrations ranged, respectively, from 4.4 to 6.6 and from 60 to 140 mg kg–1, with minimum pH values of below 2.7. The data were merged into a single dataset for each chemical property after applying a normal score (ns) transform. Normal score pH was significantly negatively correlated with the ns metal concentrations. A principal component analysis (PCA) showed that normal score pH, Zn, and Fe were associated with the residual contamination, while ns Cu, SOF, and elevation were related with historic contamination. The spatial dependence between the properties was found to be scale-dependent. The best ns estimates were produced by ordinary kriging with an anisotopic variogram model, for the properties related with Principal Component (PC) I, while those associated with PC II were best estimated using simple kriging with varying local means. A classified ns pH map showed that 33% of the study area reached at least once values of below 4 during the 2-yr period. This part of the area should be excluded to ensure successful revegetation.

Abbreviations: m, stationary mean • MAE, mean absolute error • ME, mean error • ml, varying local mean • ns, normal score • OKa, ordinary kriging with anisotropic variogram model • OKi, ordinary kriging with isotropic variogram model • PC, principal component • PCA, principal component analysis • RF, random function • RMSE, root mean squared error • SBF, sugar beet foam • SKlm, simple kriging with varying local means • SOF, soil oxidizable fraction


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
ON 28 Apr. 1998, the tailings reservoir at the Aznalcóllar pyrite mine (southwestern Spain) collapsed and released approximately 6 Mm3 of acid water and sludge into the Agrio River and, subsequently, into the Guadiamar River, causing a huge flood and the posterior deposition of approximately 2 Mm3 of pyritic sludge on the river bed and banks. A 64-km-long fringe of the river bank was affected by a sludge layer with a depth ranging from a few millimeters to more than one meter near the river, and which occupied a mean width of 600 to 700 m. The total area affected covered approximately 4634 ha of agricultural land on both sides of the river.

Shortly after the spill, a cleanup program was set up to remove the toxic residue from the river bed and banks. Between May and December 1998, a total volume of 7 Mm3 of sludge and contaminated soil was removed. In an effort to further reduce the residual pollution nearly 1 Mm3 of contaminated material was removed during 1999. Subsequently, after applying limestone, the cleaned soil was tilled to mix the remaining sludge with the topsoil layer.

The tailings were characterized by a very fine granulometry, with a particle diameter below 100 µm and a high pyrite (FeS2) content of between 680 and 780 g kg–1 (Cabrera et al., 1999; Simón et al., 1999; Vidal et al., 1999). Important quantities of other sulfide minerals were present too, yielding a broad spectrum of metals such as Zn, Fe, As, Cu, Sb, Co, Tl, Bi, Cd, Ag, Hg, and Se (listed in order of their importance) with concentrations ranging from 8 g kg–1 (Zn) to 0.01 g kg–1 (Se).

A project was set up to monitor the soil chemical properties of the reclaimed land on which remained a small—but highly spatially variable—quantity of pyritic sludge (Antón-Pacheco et al., 2001), mixed with the topsoil layer. This mixture still represents a potential risk of contamination, depending on local soil characteristics and climatological conditions. In this paper, we focus on one, intensively sampled monitoring site, situated about 8 km downstream from the tailings reservoir, where soil samples were periodically taken at different locations and analyzed, yielding a non co-located spatiotemporal dataset.

To assess the spatiotemporal distribution of the pollutants, it is necessary to understand the dynamics of heavy metals in soils, which are complex and influenced by numerous factors such as pH, soil organic matter, soil texture, redox potential, and temperature, to name only a few (Alloway, 1990). Generally, heavy metals become increasingly mobile and available as the pH decreases. The acidification starts when the pyrite sludge (FeS2) or other sulfide minerals are exposed to atmospheric oxygen and water, whereby the sulfidic component is oxidized to sulfate and H+ ions are liberated. Vanderlinden et al. (2006) give an overview of the chemical reactions involved in the pyrite oxidation. The formation of sulfuric acid decreased the pH of the residual tailings on the land and increased metal mobility and availability.

Wunderly et al. (1996) and Xu et al. (2000) have modeled the oxidation and leaching process under both saturated and unsaturated conditions. They showed that laboratory and field estimates of pyrite oxidation rates cover a wide range of values and that oxygen availability is a very important factor. Förstner and Calmano (1998) pointed out that acidity production in pyrite holding dredged material can develop many years after storage, when its buffering capacity diminishes. They also warned about the reliability of long-term prognoses on the behavior of these materials, since redox and pH controlled systems are associated with nonlinear and delayed processes. Meyer et al. (1999), from a column experiment, found that the application of sewage sludge and compost increased the oxidation rate of pyrite in lignite mine spoils, compared to the control or mineral fertilizer.

Simón et al. (2002) also summarized the chemical processes related to the oxidation of the tailings and analyzed their effects on a calcareous soil in the affected zone. They observed that after 15 mo only the top 20 mm of the soil were acidified and that most of the pollutants did not penetrate deeper than 100 mm, which could be attributed to the scant rainfall during that period. Domènech et al. (2002) determined the rate of sludge oxidation and the existence of metal retention processes from a column experiment. They found that the pH of the leaching water dropped to values of near 2 after 260 d and that Zn, Cd, and Co leached completely since no retention processes existed. Further details on the accident and related topics can be found in Grimalt and MacPherson (1999) and Ayora et al. (2001).

The analysis of spatiotemporal environmental datasets has received considerable attention in geostatistical literature. Kyriakidis and Journel (1999) reviewed existing geostatistical space–time models and distinguished three conceptual viewpoints: (i) approaches involving a single spatiotemporal random function, (ii) approaches involving vectors of space random functions, and (iii) approaches involving vectors of time series. De Iaco et al. (2003) extended the use of the Linear Coregionalization Model from spatial to spatiotemporal processes. The Bayesian Maximum Entropy (BME) framework (Christakos, 1990, 2000; Christakos et al., 2002) is another spatiotemporal estimation method that can incorporate a wide variety of hard or soft secondary information (knowledge bases). These spatiotemporal models are normally used to model environmental processes with large data sets obtained from automated monitoring networks of air quality or meteorological or hydrological variables (Christakos and Serre, 2000; De Iaco et al., 2003). These networks provide generally high frequency observations (hourly, daily) at a relatively small number of locations. When dealing with destructive sampling techniques (e.g., soil auger sampling) no further sampling can be done at the same location, giving rise to temporally non co-located spatiotemporal data sets (e.g., Papritz et al., 1993; Papritz and Flühler, 1994; Papritz and Webster, 1995a, 1995b). In addition, this kind of sampling is costly and time consuming, so that the data are usually only available at a small number of locations at a small number of instances. For these reasons it becomes under these circumstances practically impossible to infer reliable (marginal) space–time covariance models and to use one of the above mentioned spatiotemporal modeling approaches. However, Douaik et al. (2004), used the BME approach to map the spatiotemporal distribution of soil salinity using data from 419 field observations and 13 to 20 laboratory observations at 19 instances.

For agricultural studies Stein et al. (1998) distinguished between independent and dependent spatial random fields at different instances. In the first case, the spatial variability on each instance could be modeled with a unique variogram, whereas in the second case, the spatial random fields were considered replicates, with variograms proportional to a single basic variogram. In the latter case, it was postulated or decided that the range of spatial dependence was constant in time. Generally, the basic variogram can be estimated using observations collected at each of the instances, making this methodology especially useful when dealing with non co-located spatiotemporal data, where there are too few observations to estimate the variogram for each instance. For the interpolation on each instance, traditional well-known geostatistical tools can be used (Goovaerts, 1997; Deutsch and Journel, 1998). This approach was successfully applied to assess the spatiotemporal development of downy mildew epidemics in cabbage (Stein et al., 1994) and to map wind-blown mass transport (Sterk and Stein, 1997).

The main objectives of this study are to (i) analyze the relations between soil chemical properties and evaluate which ones best represent the residual pollution, (ii) adapt and apply the Stein et al. (1998) approach to infer the spatial distribution of the residual pollution at the site using non co-located observations of soil chemical properties on different instances, and (iii) compare the performance of three kriging techniques and evaluate their capacity to map the residual pyrite sludge to assess further remediation and to develop successful revegetation schemes. In a companion paper (Vanderlinden et al., 2006) the temporal dynamics of the soil chemical properties at the site are addressed.


   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
Site Description and Sampling Scheme
The Aznalcóllar mining area is situated within the Iberian pyrite belt, approximately 25 km west of Seville in southwestern Spain (Fig. 1 ). Early evidence of mining activity dates from ancient Roman times. Due to this activity and to natural processes, high heavy metal concentrations had already been reported in the Guadiamar River basin before spillage occurred in April 1998 (Cabrera et al., 1987).



View larger version (64K):
[in this window]
[in a new window]
  		Fig. 1. Location of the Iberian pyrite belt within Andalusia, Spain, and of the study zone in the province of Seville. The ortho-photograph gives a detailed view of the meander and the monitoring site along the Guadiamar River during the cleanup activities in July 1998.

 
In this paper, we focus on an intensively sampled monitoring site (37°26' N, 6°13' W), situated at a meander in the Guadiamar River, left bank, about 8 km downstream from the tailings reservoir (see Fig. 1). The terrain is rather flat (see Fig. 2 ) and was flooded by the sludge from north to south. The study area was cleaned at the end of 1998 (see also Fig. 1 for an ortho-photograph taken during the cleanup actions; Junta de Andalucía, 2003), and was then tilled to incorporate previously applied pH correctors such as sugar beet foam (SBF) or lime. At the site, recovered tailings were temporarily stored during the cleaning operations, to be removed later by bulldozers filling larger trucks that carried the load back to the open mine pit. The recharging process of the recovered tailings augmented the remaining contamination at the site in the bands left by the bulldozer tracks.



View larger version (78K):
[in this window]
[in a new window]
  		Fig. 2. Location of the 503 sampling points within the monitoring site depicted in Fig. 1, according to the sampling date. The background gray-scale map represents local relative topography.

 
Classified multispectral images of the area (Antón-Pacheco et al., 2001) showed a pattern of parallel remaining sludge bands, from northwest to southeast. The soil is a Typic Xerofluvent (Soil Survey Staff, 1999) with a prevailing loam to sandy loam texture, although coarser materials such as gravel and sand may occasionally appear, in accordance with the alluvial character of the site. These coarse-textured spots form a potential gateway to the underlying aquifer for heavy metals.

Periodical sampling took place on six different dates between June 1999 and December 2000 (June 1999, September 1999, December 1999, March 2000, May 2000, and September 2000). The locations of the sample points are shown in Fig. 2. There were a total of 503 sample points, partly clustered and partly on a pseudo-regular grid, covering an area of approximately 9 ha. The June 1999 sampling was intended as a reconnaissance survey and covered only a small part of the study area. Only the May 2000 survey covered the whole area, while the others covered adjacent or partly overlapping areas.

All soil samples were taken with a 7-cm-diameter auger from the 0- to 15-cm horizon and stored in marked plastic bags for their transport to the laboratory. Sampling point coordinates were measured on each date with a Topcon (Livermore, CA) topographic station and georeferenced to a common coordinate system with the x axis parallel to the direction of the remaining sludge bands. Soil moisture conditions ranged from very dry (June 1999) to wet (December 1999). During February 2000, SBF was applied to the site and consecutively mixed with the topsoil layer by tillage which temporally affected soil pH values, as the results show.

Analytical Methods
All samples were air-dried and sieved (< 2 mm). The pH was measured in both water and CaCl2 solutions (Lab pH-meter Basic 20; Crison SA, Alella, Barcelona, Spain). The soil oxidizable fraction (SOF) of the contaminated soil was determined through oxidation by Cr2O7K2 following a procedure similar to that used for the determination of organic matter (Nelson and Sommers, 1996). Iron, Zn, and Cu were extracted with a solution of CH3COONH4 0.5 M + CH3COOH 0.5 M + EDTANa2 0.02 M and subsequently measured by atomic absorption spectrophotometry (Atomic Absorption Spectrometer 3110; PerkinElmer, Wellesley, MA).

Coefficient of Correlation and Principal Component Analysis
To control for linear relationships between variables, the coefficient of correlation was calculated and used to group the variables that were most related to the residual sludge. Principal component analysis (PCA) (Davis, 1973; Webster and Oliver, 1990; Wackernagel, 1998) is an ordination method that can elucidate relations among variables by identifying common underlying processes. The resulting principal components (PCs) or factors are orthogonal linear combinations of the variables. The first PC explains a major part of the total variance of the data set, and each successive PC explains a smaller part of the remaining variance. The different PCs are then related with common processes that affect the variables. Varimax rotation can than be used to make the interpretation of the PCs easier (Webster and Oliver, 1990). Principal component analysis was used here to confirm the results obtained from the comparison of the coefficients of correlation and to find one or more PCs that could be associated with the residual contamination.

Geostatistical Framework
Outset
We use the term dataset to denote the entire spatiotemporal dataset, and subdataset to denote each spatial dataset that corresponds to a single date. On each sampling date a different number of samples was collected in different subareas of the monitoring site, as shown in Fig. 2. For example, the area sampled in March 2000 lies completely outside previously sampled areas. However, classified multispectral images (Antón-Pacheco et al., 2001) showed that the overall remaining sludge pattern was similar in the entire monitoring site.

The spatial distribution of the residual pyrite sludge did not change in time, but those of the soil chemical properties that reflected this contamination did, due to, for example, pyrite oxidation or the application of limestone or SBF, processes that occur homogeneously throughout the monitoring site. Taking this into account we consider that their spatial patterns remained constant in time.

The dataset, consisting of data on different soil chemical parameters, from the six sampling dates, can be considered as, (i) non co-located spatiotemporal datasets, one for each analyzed property or (ii) as six co-located spatial subdatasets, one for each date, containing each data for all properties analyzed. In the first case, data from other sampling dates should be taken into account to assess the spatial distribution of each soil chemical property on a single sampling date. This is a difficult task since there are no co-located observations on different dates. In the second case, the spatial distribution for a single soil chemical property on a certain date would be assessed taking into account data on other chemical properties, measured on the same date and at the same location. Since the additional information (secondary variables) is co-located with the primary variable, it is of little value for the spatial interpolation (Goovaerts, 1997). In addition, on each sampling date, only the surveyed part of the area can be mapped in this way and the information from other sampling dates is not taken into account.

The use of pseudo-crossvariograms (Rouhani and Myers, 1990; Myers, 1991) with cokriging could offer another alternative, but was not adopted here after the criticisms of Papritz and Flühler (1994), who evaluated the usefulness of pseudo-crossvariograms with a spatiotemporal non co-located dataset, and of Wackernagel (1998). Indicator kriging (Goovaerts, 1997) is commonly used to deal with skewed data but could also be valid in this context. Using the deciles of the distributions of each variable for each date as threshold values leads to standardization and makes it possible to produce maps of the probability of exceeding a threshold value on a certain date. We did not consider this possibility in this study due to a large number of indicator variograms that would necessarily have to be inferred and modeled.

Sterk and Stein (1997) standardized their wind-blown mass transport measurements for four wind storms and calculated a stratified variogram to overcome the inconveniences of a small data set. Although for each storm the wind velocity and average mass transport were significantly different, the same basic spatial correlation structure was considered in the spatial interpolation of the transported mass of each storm. From here on we will consider the first option and treat our data as non co-located spatiotemporal data sets, one for each considered soil chemical property.

Normal Score Transform
Inspired by Sterk and Stein (1997) we standardized the subdatasets using the normal score transform y(x) = {phi}[z(x)], with z(x) as the original data, &{phi}(·) the transformation function, and y(x) the transformed data, or normal scores (ns), with a standard normal histogram (Deutsch and Journel, 1998; Goovaerts, 1997). Tied data values were randomly ordered during the transformation. For each soil chemical property, the transform was applied to the six subdatasets. The conditional cumulative distribution function for each soil chemical property can be written, after Goovaerts (1997), as:

[1]


where N = {sum}ni with i = 1, ... , 6, and n is the available data in each subdataset. The terms Zi(x and Y(x) are the random functions (RF) of the original and transformed data, respectively, and G(y) is the standard normal cumulative distribution function.

Through this transform the subdatasets are standardized (zero mean and unit variance). Note that assuming temporal dependence (consistency) of the spatial pattern of the residual pollution during the monitoring period, the RFs of the ns-transformed subdatasets, Yi(x), can be considered identical, which means that the pooled ns-transformed dataset can be represented by a single RF, Y(x). Therefore, it is very important that the data that constitute the different subdatasets, zi(x), represent well the corresponding RF, Zi(x), so that accurate transformation functions, &{phi}i(·), are defined.

This approach makes it possible to treat and analyze the spatiotemporal datasets (one for each soil chemical property) as single ns spatial datasets and to interpret the results in terms of the ns-transformed properties. If conclusions have to be drawn in terms of the original data scales on a certain date, then the results will have to be back-transformed, according to the corresponding back-transformation function, &{phi}–1(·).

Variography
The residual sludge is thus quantified through ns-transformed soil chemical properties, which are assumed constant in time. This offers the possibility of inferring a single time-independent sample variogram for each ns-transformed chemical property. At the same time, it may overcome problems related with the small number of observations of some subdatasets and the different sampling areas covered during each sampling campaign. The pooled sample auto-variogram was computed from the merged ns-transformed dataset using the traditional equation (Goovaerts, 1997):

[2]
where Np(h) is the number of data pairs separated by a distance h and y(xi) is the ns-transformed sample value at location xi. The stratified sample auto-variogram was computed from the ns-transformed subdatasets using the following expression:

[3]

where Ns(h) = nj(h) is the total number of data pairs separated by a distance h, and nj(h) is the number of data pairs belonging to the jth subdataset, separated by a distance h. Using Eq. [3] only data values from the same subdataset are paired and possible inhomogeneities of the merged ns-transformed dataset do not affect the results. We assumed that the spatial auto-correlation structure of the ns-transformed data did not change between the different sampling dates and that it was valid for the total area covered by the different data sets. This decision was motivated by the consideration that the total area was subject to the same reclaiming actions and that it was sufficiently small to ignore meteorological variations within it. The homogeneity of the merged ns-transformed dataset was controlled by comparing the values of the pooled and stratified sample variograms at the first lags (not shown). If the merged ns dataset were not homogeneous, then the pooled variogram (Eq. [2]) would be expected to produce larger values at the first lags as a consequence of the loss of spatial continuity at short distances. Note that Eq. [2], in contrast to Eq. [3], allows data pairs between data from different sampling dates.

To analyze the joint spatial variation between the variables we also calculated pooled sample cross-variograms from the merged ns-transformed dataset:

[4]
where the subscripts u and v denote two regionalized variables.

A shortcoming of the product-moment correlation coefficient is that it does not take not into account the location of the sample points, so it is unable to detect a changing correlation with changing scale. The sample codispersion function (Matheron, 1965; Goovaerts and Webster, 1994; Webster and Oliver, 2001):

[5]
a lag-dependent correlation, overcomes this disadvantage and can be used to characterize the scale-dependence of the relationships between the variables. It can be interpreted as the correlation coefficient between the spatial increments of Yu(x) and Yv(x) and can be calculated directly from the sample auto- and cross-variograms. If the relation between both variables does not change with spatial scale, then {rho}uv(h) will be constant for all h and the variables are said to be intrinsically correlated. Under the condition of second-order stationarity, {rho}uv(0) equals the ordinary product-moment correlation coefficient.

Spatial Interpolation
By kriging the ns-transformed data followed by a back-transformation, it is possible to obtain maps of the entire area for each date. In this way data from other sampling dates are taken into account without having to recur to cokriging. Cokriging would require the estimation and modeling of cross-variograms with observations from different sampling dates, which is not possible in this case since the data on different sampling dates are not co-located.

The random function Y(x) can be broken down as:

[6]
where m(x) is the deterministic component, {xi}'(x) is the stochastic spatially dependent term, and {xi}'' represents a normally distributed spatially uncorrelated noise term which includes all the variation that our geostatistical model cannot explain. This component contains the nugget variance, but accounts also for errors that are due to the violation of the underlying hypothesis of kriging. Simple kriging with varying local means (SKlm) (Goovaerts, 1997, 1999) offers a possibility of estimating the expected value of Y at an unsampled point, while considering the deterministic component of Eq. [6]:

[7]
where m(xo) and m(xi) are the assumed known local varying means, at the unsampled point xo and at the n neighboring points xi, respectively, y(xi) are the n neighboring observations of Y, and {lambda}iSKlm are the SKlm weights. These weights were obtained by solving the simple kriging system using the covariance function of the residual random function:

[8]

Ordinary kriging (OK):

[9]
assumes that the local mean is unknown, but stationary within on xi centered search neighborhoods, and filters it out, imposing the condition that the weights should add up to unity to ensure unbiasedness. The corresponding system of equations can be solved either in terms of isotropic variogram models that do not take into account possible anisotropy or with anisotropic variogram models that capture the anisotropy of the spatial correlation structure.

Cross-validation (CV) is used here to compare the performance of OK with an anisotropic variogram model and SKlm, to interpolate the ns-transformed soil chemical properties. The results obtained for OK with an isotropic variogram model are used as a reference, since it does not take into account extra information on the spatial distribution of the residual pollution (deterministic component of the RF). Cross-validation is done by eliminating consecutively each data value from the data set and then estimating it. Different statistical parameters are then used to compare estimated and observed values. Here we will use:

(i) the mean error (ME), which quantifies the global bias of the estimators, y*(xi):

[10]

(ii) the mean absolute error (MAE):

[11]
which gives the precision or dispersion of the estimates against the observations, y(xi), and

(iii) the root mean square error (RMSE):

[12]
which is similar to the MAE, but gives more weight to the largest errors.

To assess the temporal evolution of the soil chemical properties at the site (Vanderlinden et al., 2006), the interpolated ns values have to be backtransformed to the original data scales of each sampling date. This is not at all straightforward since the optimal character of the kriging estimates is lost and unbiasedness can no longer be guaranteed. The nonlinear back-transform, {phi}–1(·), magnifies the interpolation errors of the ns estimates, yielding biased back-transformed estimates. In this case, especially, underestimation (negative bias) of the high metal concentrations and, although less serious, overestimation (positive bias) of the low pH values are important. A thorough analysis of this issue is given in Vanderlinden et al. (2006).

Using the interpolation method that performed best in the cross-validation, ns maps were produced for each transformed soil chemical property. These maps give a time-independent representation of the spatial distribution of the residual contamination at the site in terms of ns-transformed soil chemical properties. To avoid the back-transform and biased back-transformed estimates each ns map can be classified according to a ns threshold value that corresponds to a certain value of the soil chemical property. There are two types of misclassification errors; false positives and false negatives. The probability of a Type I error (false positive) is represented by {alpha} and ß is the probability of a Type II error (false negative).


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
CONCLUSIONS
 REFERENCES
 
Exploratory Data Analysis
Descriptive Statistics
To control for the effect of the clustering of sample points and possible preferential sampling in remaining sludge spots, declustering weights (Deutsch and Journel, 1998) were calculated, but no significant influence on the descriptive statistics could be detected. Table 1 shows the descriptive statistics of pH, SOF, and the EDTA-extractable concentrations of Fe, Zn, and Cu for the six sampling dates from June 1999 to September 2000, with 42, 83, 74, 96, 129, and 79 sampling points, respectively.


View this table:
[in this window]
[in a new window]
  		Table 1. Descriptive statistics for pH; soil oxidizable fraction (SOF); and EDTA-extractable Fe, Zn, and Cu for the six sampling dates.

 
Temporal variations could clearly be observed for pH and metal concentrations, while SOF seemed to remain constant. These fluctuations in time could be best appreciated in the case of pH, with mean values ranging from 4.4 in June 1999 to 6.6 in March 2000, due to the application of SBF during February 2000. Figure 3 shows the histograms of pH and Cu, before and after the SBF application. Note how the bimodal histogram of pH in December 1999 shifted to the right in March 2000, while the proportion of Cu observations of <20 mg kg–1 increased from 46 to 68%, after adding SBF. The other metal concentrations also followed this pattern because they become less mobile and available as the pH increases. The statistical distributions of the metal concentrations were highly positively skewed, especially for Fe and Zn, which was, together with the wide range of values, a typical characteristic of metal-contaminated soils. Note also the large differences between the mean and the median for these properties, hence the need for a transformation of the metal concentrations. It can be concluded from Table 1 that the situation in September 2000 was quite similar to that observed in June 1999, despite the application of SBF. This observation demonstrates the fact that, more than 2 yr after the cleanup, pyrite oxidation was still occurring in the areas with residual sludge. A detailed analysis of the temporal evolution of the soil chemical properties at this site is reported in Vanderlinden et al. (2006).



View larger version (45K):
[in this window]
[in a new window]
  		Fig. 3. Histograms of pH and EDTA-extractable Cu, before and after the sugar beet foam application in February 2000.

 
Correlation Coefficient
The correlation matrix of the ns-transformed variables is shown in Table 2. It was calculated using the entire ns-transformed merged data set. The elevation of the sample points, elev, was also considered, to evaluate possible relationships with the local topography, bearing in mind that lower areas were covered by a thicker sludge layer and that cleaning at these sites could have been less effective. Negative correlation coefficients between pH and metal concentrations were expected since increasing acidity enhances metal mobility and availability. It can be seen from Table 2 that Fe and Zn were especially correlated with pH. The SOF was less correlated with these variables, although this variable was first thought of as being a good indicator of the presence of pyrite, given the large amount of oxidizable material in the sludge. Only the Cu concentration was related to SOF. Local topography had no influence on the variables studied, except for the Cu concentration. According to this information, the variables considered could be split up into two groups: pH and Fe and Zn concentration on one hand and elev, SOF, and Cu concentration on the other. The scale-dependence of the relationships between the variables was analyzed with the sample codispersion function in the Variography section.


View this table:
[in this window]
[in a new window]
  		Table 2. Correlation matrix for elevation (elev); normal score (ns)–transformed properties of pH; soil oxidizable fraction (SOF); and EDTA-extractable concentrations of Fe, Zn, and Cu.

 
Principal Component Analysis
Principal component analysis of the six variables showed that the first PC explained 51% (= 3.07/6) and the first two PCs together explain 74% [= (3.07 + 1.38)/6]. The eigenvalue that corresponded to the third PC equaled 0.69 and explained a smaller part of the total variance (0.11%) than each variable on its own. For this reason, only the first two PCs were retained. Figure 4a shows the position of the six variables in the plane defined by the first two PCs, called the circle of correlations, where the coordinates represent the correlation coefficients between the variables and the PCs. This representation reveals the proximity of the different variables within a circle with a unit radius.



View larger version (12K):
[in this window]
[in a new window]
  		Fig. 4. Circles of correlation, showing the configuration of elevation (elev); normal score (ns) pH; ns soil oxidizable fraction (SOF); and ns EDTA-extractable Fe, Zn, and Cu, in the plane defined by the first two principal components (PC): (a) without rotation and (b) after Varimax rotation.

 
To improve the contribution of the different variables to the PCs, a Varimax rotation of the PC axis was performed, as shown in Fig. 4b. The six variables could be separated into two groups. A first group, with a major contribution to PC I, represented the residual contamination within the study area. The antagonism between the pH and the metal concentration (Fe and Zn) could be clearly observed. A second group, with a major contribution to PC II, represented the effect of elevation on the spatial distribution of the SOF and the Cu concentration. Aerial photographs from the 1950s revealed a former course of the Guadiamar River at the lowest parts of the study area. This could explain the accumulation of organic matter (SOF) and Cu in the former river bed and banks, since Cu is generally found to be bound on organic matter (Sposito, 1989). This source of contamination dates from before the spillage of 1998 (Cabrera et al., 1987). It is probable that the concentrations of Fe and Zn also showed this spatial pattern before the spillage, but the offsetting effect of the greater contribution by the sludge and the subsequent cleaning masked it. This was confirmed by the multivariate spatial structural analysis in the next section.

Finally, each data value was represented against the first two PCs to detect data grouping and differences between data that stem from different sampling dates (not shown). No data grouping could be observed, which confirms the hypothesis of the uniformity of the ns-transformed observations originating from different sampling campaigns.

Geostatistical Analysis
Variography
Sample auto-variogram surfaces of the ns-transformed soil chemical properties and elevation (see Fig. 5 ) were calculated, using the Variowin software (Pannatier, 1996), to assess the anisotropy of the underlying spatial correlation structure of the remaining pyritic sludge. In the case of pH, Fe, and Zn, the pattern consisted of alternating bands of high and low sample variogram values in the x direction. This direction corresponded to the orientation of the remaining sludge bands that could be observed from the multispectral images. A different spatial pattern was observed for SOF and Cu, which could be associated with previous depositions of Cu, and with elev.



View larger version (130K):
[in this window]
[in a new window]
  		Fig. 5. Sample variogram surfaces for normal score (ns) pH; ns soil oxidizable fraction (SOF); ns EDTA-extractable Fe, Zn, and Cu; and elevation, elev.

 
For all attributes, the direction of the major continuity was the x direction. In the y direction, harmonic directional sample auto-variograms were obtained for pH, Fe, and Zn. The directional and omnidirectional sample auto-variograms were calculated using Eq. [2]. These variograms, with fitted models, are shown in Fig. 6 . In the case of SOF and Cu, anisotropy in the spatial correlation structure was mainly observed at the largest separation distances. The cyclic pattern, that characterized the pH, Fe, and Zn variograms in the y direction, did not appear in those attributes. This harmonic behavior could be modeled using the hole effect model (Cressie, 1991; Deutsch and Journel, 1998) with a range, or wavelength, of approximately 20 m. This kind of pattern is common in geology and has also been observed in the context of the spatial variability of agronomical attributes, where they are induced by fertilizer applicators (Solie et al., 1999) or irrigation machines.



View larger version (49K):
[in this window]
[in a new window]
  		Fig. 6. Directional and omnidirectional sample variograms, with their respective fitted models for normal score (ns) pH; ns soil oxidizable fraction (SOF); and ns EDTA-extractable Fe, Zn, and Cu. The fitted parameters are represented in Table 3.

 

View this table:
[in this window]
[in a new window]
  		Table 3. Fitted parameters of the isotropic and anisotropic variogram models for normal score (ns)–transformed properties of pH; soil oxidizable fraction (SOF); and EDTA-extractable Fe, Zn, and Cu concentrations.

 
All sample variograms were calculated and those without harmonic behavior were modeled using Variowin (Pannatier, 1996). A mixed automatic–visual approach using current graphing software (Golden Software, 2002) was employed to fit models that contained a hole effect model. The sample variograms in Fig. 6 were modeled using three basic permissible models (Goovaerts, 1997; Deutsch and Journel, 1998). The exponential model is defined by an effective range, a, and a positive contribution to the sill variance, c:

[13]

The spherical model is given by:

[14]
with an actual sill, a. The hole effect model is given by the following expression:

[15]
where a the distance to the first peak. All linear combinations of these basic models also constitute permissible models. In this study, we used a maximum of four structures, including the nugget effect, co. Table 3 shows the parameters of the fitted models. The hole effect model is only positive semi-definite in one dimension (Journel and Huijbregts, 1978) and could, therefore, only be used in the y direction. The model is annulated in the x direction, choosing a large range of 105 m in that direction.

The spatial cross dependence between the ns-transformed chemical properties was quantified by means of sample cross-variograms (Eq. [4]). Using these cross-variograms and the corresponding sample auto-variograms, the sample codispersion function was calculated using Eq. [5], to assess the scale dependence of these spatial relationships. The omnidirectional cross-variograms and codispersion functions are shown in Fig. 7 . It can be observed how the codispersion functio napproached the ordinary product-moment correlation coefficient (see Table 2) at the largest separation distances. As in the correlation and PC analysis, again a different behavior could be observed between the two groups of soil properties. The strongest spatial dependence was observed between pH, Fe, and Zn. The sample codispersion functions for these properties showed that the correlation between the increments changed with distance; for increments between nearby points, up to approximately 10 m, the spatial relationships were weaker, while their strength remained more or less constant for increments of data points separated by more than 10 m. This could be due to relatively large analytical errors for Fe and Zn in comparison with the nugget effect, especially for the largest values, associated to locations within the residual sludge spots. This could also explain the larger nugget variance of these properties (approximately 50% of the total sill, see Table 3). Following Goovaerts and Webster (1994), conclusions on the origin of the nugget variance of the auto-variograms could be derived from the nugget variance of the corresponding cross-variograms. The nugget variance consists of unexplained microscale variation that is not captured by the sampling design and measurement or procedural errors. The measurement or procedural error of both properties could be assumed independent (i.e., the covariance of these variations was zero), so that the nugget variance of the cross-variograms derived only from unexplained microscale variance. For these properties the sample cross-variograms had a small nugget, indicating that there existed little correlation between increments at the microscale. This suggested that the nugget variances in the auto-variograms of Fe and Zn arose from analytical and procedural errors. The spatial dependence between properties of the second group (SOF, Cu, and elev) was weaker, in accordance with the correlation coefficients in Table 2. The inverse relation of SOF and Cu with elev reached a maximum strength for increments between data points separated by more than 30 m and increased gently toward the largest separation distances. Although the SOF–Cu cross-variogram showed a clear spatial cross dependence, the correlation between the increments (codispersion function) was not scale independent. No spatial dependence was observed between the properties of Group 1 and elev. The remaining combinations of the properties of Group 1 and 2 showed a weak spatial cross dependence. Note that the strength of the spatial relationships between the metals of both groups increased with distance, suggesting a common large scale pattern that can be associated with the historical pollution at the site.



View larger version (31K):
[in this window]
[in a new window]
  		Fig. 7. Omnidirectional sample cross-variograms and codispersion functions for the pooled normal score (ns)–transformed chemical properties: pH; EDTA-extractable Fe, Zn, and Cu; soil oxidizable fraction (SOF); and elevation (elev).

 
Varying Local Means
The multispectral images and the sample variogram surfaces (see Fig. 5) showed a clear cyclic spatial pattern in the residual pollution. This knowledge was incorporated into the estimation process using SKlm, considering this discontinuity or trend as varying local means. These local mean surfaces were inferred for each variable, using OK with a pure nugget effect and an extremely stretched elliptical search neighborhood, with the major axis along the x direction. This is the same as local averaging and is similar to moving-windows averaging. Local mean surfaces with a strong discontinuous pattern along the y direction and almost continuous bands along the x direction were obtained (not shown), which represented the deterministic component of the spatial pattern of the residual pollution, consisting of parallel residual sludge strips. These estimated varying local mean surfaces are used in SKlm as if they were known exactly. For ns-transformed SOF and Cu a circular search neighborhood was used, since these variables did not exhibit this anisotropic spatial pattern. Residuals were obtained according to Eq. [8] and used to calculate residual sample variograms.

Performance of the Different Interpolation Algorithms
Cross-validation was used to compare OK with isotropic variogram models (OKi), OK with anisotropic variogram models (OKa), and SKlm using the inferred varying local means and the covariance functions of the residual RFs. The error was evaluated statistically in terms of the ME, MAE, and RMSE, and represented in Fig. 8 . The ME was only negative for the ns-transformed pH. Compared to the MAE and RMSE, the ME values are rather small and could be ignored. In terms of the MAE, ns-transformed pH and Zn behave similarly, and so did ns-transformed SOF and Cu. The ns-transformed Fe had the largest MAE values, in accordance with the large nugget effect observed for this variable. The minimum MAEs for ns-transformed pH, Zn, and Fe were reached for OKa, while SKlm yields the lowest MAEs for ns-transformed SOF and Cu. Approximately the same pattern was observed for the RMSE, except for Zn, which showed a slightly smaller RMSE for OKi. Based on this analysis we decided to use OKa for the interpolation of the sludge-related ns-transformed variables and SKlm for ns-transformed SOF and Cu. All calculations were performed using the KT3D program of the GSLIB software (Deutsch and Journel, 1998).



View larger version (14K):
[in this window]
[in a new window]
  		Fig. 8. Performance of three interpolation methods based on the kriging algorithm for the interpolation of normal score (ns)–transformed soil chemical properties. OKi, ordinary kriging with an isotropic variogram model; OKa, OK with an anisotropic variogram model; and SKlm, simple kriging with local varying means.

 
Normal Score Maps
Figure 9 shows the block-kriged ns maps for the five variables, obtained with OKa and SKlm using the GSLIB software (Deutsch and Journel, 1998). Block dimensions were 2.5 x 2.5 m, and four discretization points were used. The spiky effect of the ns maps for pH, Fe, and Zn was caused by the large differences between the nugget effect of the variograms in the x and y direction (see Fig. 6). The model of the latter was forced through the nugget effect of the variogram model in the x direction. Near the origin, a steep-sloping variogram in the y direction was obtained, reflecting a large variability at short distances in that direction. This is not the case for ns-transformed SOF and Cu, which show a smoother spatial variability (see Fig. 9).



View larger version (95K):
[in this window]
[in a new window]
  		Fig. 9. Block kriged normal score maps for pH; soil oxidizable fraction (SOF); and EDTA-extractable Fe, Zn, and Cu. Block dimensions are 2.5 x 2.5 m.

 
Both variables are interpolated using SKlm with previously calculated varying local mean surfaces and the covariance function of the residual RFs. These maps constituted a time-independent representation of the spatial distribution of the analyzed soil chemical properties and of the residual contamination at the site.

Normal score maps are, however, of little use for decision makers or managers. What really matters is whether a block is estimated as being contaminated or safe, to design the cleanup or reforestation strategy, and to calculate its total cost (Goovaerts and Saito, 2000). Considering a threshold value above which remediation is mandatory or revegetation becomes problematic is an alternative to back-transformation of the ns maps. Here, as an example, we considered the pH, because of its important contribution to the first PC, which represented residual contamination. We selected a threshold value of pH = 4 at which the growth of plants and trees is limited and below which the oxidation rate of Fe2+ to Fe3+ becomes very low, constituting a pyrite oxidation rate-limiting step (Xu et al., 2000). For each sampling date, the ns map was then classified according to the equivalent ns threshold values for pH = 4. Here, we considered only the worst case, yielding an equivalent ns value of –0.28 for September 1999.

In the classified pH map, shown in Fig. 10 , about 38% of the pH values observed for September 1999 have a pH < 4, while, respectively, 36 and 33% of the cross-validation and grid block estimates were below the threshold value. Twelve percent of the CV estimates were estimated erroneously as pH < 4 (Type I error, false positive), while 14% of the CV estimates were classified erroneously as pH > 4 (Type II error, false positive).



View larger version (39K):
[in this window]
[in a new window]
  		Fig. 10. Classified normal score pH map for September 1999, according to a threshold value of pH = 4.

 
The map shown in Fig. 10 constituted a valuable tool for designing a successful reforestation scheme and thus reducing the number of trees that die after planting. About 33% of the study area should be excluded for reforestation. The Type II error can be taken into account by considering a buffer zone around the contaminated areas of Fig. 10. Further research will be undertaken to estimate the local uncertainty, using indicator kriging and stochastic simulation, to improve the assessment of the delineation of appropriate areas for successful revegetation.


   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
REFERENCES
 
Pyrite oxidation resulted in acidification and increasing metal mobility and availability. During the 2-yr monitoring period spatially averaged pH fluctuated between 4.4 to 6.6, with minimum values of near 2.5 for each sampling date. The application of SBF after the first year contributed significantly to this fluctuation, although similar average pH values were observed during the first and last sampling campaign. The median of the EDTA-extractable Fe and Zn concentrations ranged, respectively, from 109 to 509 and from 17 to 94 mg kg–1, inversely proportional to the pH.

The underlying spatial distribution of the residual pyrite sludge contamination at a reclaimed site could be inferred from the ns-transformed soil chemical properties observed at different locations during consecutive sampling campaigns. For each property and sampling date, the ns-transform standardized the data and transformed them into a Gaussian distribution, which was especially useful when dealing with skewed distributions (e.g., EDTA-extractable Fe, Zn, and Cu). These ns-transformed datasets were then merged into a single dataset for each property, assuming that the underlying spatial correlation structure did not change in time. This hypothesis was supported by the clear spatial pattern of the residual sludge at the study site, as observed from multispectral images.

The correlation matrix and a PC analysis of the ns-transformed properties showed that the considered variables grouped together according to the underlying features that they represented. On one hand, ns pH, ns Fe, and ns Zn yielded a major contribution to PC I, and were related to the residual pyrite sludge contamination at the site. On the other hand, ns SOF, ns Cu, and elevation had a major contribution to PC II and represented the historic contamination at the lower parts of the site due to flooding. Properties of both groups showed scale-dependent relationships. For the metal concentrations the lack of correlation between spatial increments at short distances was associated with the analytical and procedural errors, as concluded from the cross-variance analysis.

Ordinary kriging with anisotropical variograms performed best for the interpolation of ns pH and ns Fe and Zn concentrations, while SKlm was selected for the interpolation of ns SOF and ns Cu. These maps give a time-independent representation of the residual contamination at the site in terms of ns-transformed variables. Applying the corresponding back-transforms to the ns interpolated maps, the spatial distribution of the phenomenon studied can be reconstructed for each sampling date, in terms of the representative soil chemical properties. However, the back-transform implies the loss of the unbiasedness and minimum estimation variance properties of the estimates.

Since we were basically interested in assessing the reforestation scheme, it was sufficient to classify the pH ns map using as a threshold value the ns equivalent of pH = 4 in the worst case (i.e., the sampling date for which the largest ns value for pH = 4 was observed). It was found that approximately 33% of the study area should be excluded to guarantee successful revegetation.

These residual sludge-contaminated zones are not only a problem for revegetation but also constitute a serious treat for the environmental recovery of the area since prolonged episodes of intense rainfall could trigger the surface (runoff and erosion) and subsurface movement (leaching) of residual contaminants in the future.


   ACKNOWLEDGMENTS
 
The authors are especially grateful to M. Armenteros for his indispensable help with the field work, and to I. Ordóñez, C. del Moral, C. Lara, and A. García for the laboratory analysis. We also acknowledge the constructive comments of O. Wendroth and two anonymous reviewers. Financial support was supplied by the co-financed FEDER-CICYT Projects 1FD97-0765 and 1FD97-0765, and the DGIFA CAP Junta de Andalucía PIA-03-050 project.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES