Journal of Environmental Quality 31:1875-1884 (2002)
© 2002 American Society of Agronomy, Crop Science Society of America, and Soil Science Society of America
TECHNICAL REPORTS
Heavy Metals in the Environment
Mapping the Probability of Exceeding Critical Thresholds for Cadmium Concentrations in Soils in the Netherlands
D. J. Brus*,a,
J. J. de Gruijtera,
D. J. J. Walvoortb,
F. de Vriesa,
J. J. B. Bronswijkc,
P. F. A. M. Römkensa and
W. de Vriesa
a W. de Vries, Alterra, Green World Research, P.O. Box 47, 6700 Wageningen, the Netherlands
b Agricultural University Wageningen, Duivendaal 10, P.O. Box 37, 6700 AA Wageningen, the Netherlands
c National Institute of Public Health and the Environment, P.O. Box 1, 3720 BA Bilthoven, the Netherlands
* Corresponding author (d.j.brus{at}alterra.wag-ur.nl)
Received for publication April 24, 2001.
 |
ABSTRACT
|
|---|
The probability of exceeding critical thresholds of Cd concentrations in the soil was mapped at a national scale. The critical thresholds in soil were based on food quality criteria for Cd in crops or in organs of cattle (Bos taurus), and were calculated by inverting a regression model for the Cd concentration in the crop, with the Cd concentration in soil, soil organic matter (SOM) content, clay content, and pH as predictors. The probability of exceeding the critical threshold for Cd in soil per node of a 500- x 500-m grid was approximated by Monte Carlo simulation, using the estimated cumulative distribution functions (cdf) of SOM, clay, pH, and Cd as input. The cdfs were estimated by simple indicator kriging with local prior means. For SOM, clay, and pH, detailed maps of soil type and land use were used to define subregions with assumed constant local means of the indicators (a priori distributions). The cdfs were sampled by Latin hypercube sampling. We accounted for correlation between the actual and critical Cd concentrations in soil by drawing Cd values from cdfs conditional on SOM and clay. The estimated probability for grassland is negligible, even in areas with high Cd concentrations in soil, and for maize (Zea mays L.) land the probability is almost everywhere smaller than 5%. For arable soils, however, these probabilities commonly are larger than 5% when sugar beet (Beta vulgaris L.) or wheat (Triticum aestivum L.) is taken as a reference crop, and locally exceed 50%.
Abbreviations: cdf, cumulative distribution function ccdf, conditional cumulative distribution function LGN3+, landuse database of the Netherlands SOM, soil organic matter content
 |
INTRODUCTION
|
|---|
THERE IS WORLDWIDE concern about the dispersion of heavy metals such as Cd, Pb, Cu, and Zn into the rural environment through atmospheric deposition and the use of agrochemicals and manure. For the Netherlands, Van Driel and Smilde (1981) reported high heavy metal concentrations in recent river sediments that grossly exceeded the concentrations considered acceptable for cultivated soils. Somewhat later, Wiersma et al. (1986) reported high concentrations of Cd and Pb in cereals, often exceeding the proposed maximum allowable concentrations for human consumption in the Netherlands. In 1991 critical thresholds for heavy metal concentrations in agricultural soils were published (Landbouwadviescommissie milieukritische stoffen, 1991). These thresholds, referred to as LAC values, serve as a warning system: if the heavy metal concentrations in the soil exceed these thresholds, then there is a substantial probability that the quality standards for heavy metal concentrations in the crop or organs of cattle are exceeded. Four land use types are defined, and for each land use type two LAC values are presented, one for sand and one for clay or peat. Within these soil parent material categories the LAC is a constant.
Recently, a new method has been proposed to calculate critical thresholds of heavy metal concentrations in soils derived from quality standards for Cd in crops used for human consumption or as fodder, or for Cd in organs of cattle. The core of the method is the use of regression models to predict the heavy metal concentration in crops from heavy metal concentrations in soils, soil organic matter content, clay content, and pH as predictors. In many soils, these basic soil properties govern the transfer of heavy metals from soils to crops. For Cd, these regression models describe the variation of the metal concentration in various crops reasonably well. Therefore, the new method is a promising tool to derive soil quality standards related to health risks from Cd.
Contrary to the LAC values, the new critical thresholds also vary within soil parent material categories, because the basic soil properties SOM, clay, and pH vary within these categories. The aim of this study was to make maps of the new critical thresholds for Cd concentration in soils, and of the probability that the actual Cd concentration in soils exceeded these new critical thresholds. This was done for the dominant land use types in the Netherlands: grassland, maize land, and arable land. The critical thresholds for Cd in soils and the probability of exceeding the critical threshold for Cd in soil at a given location was valuated for the actual land use at that location.
 |
MATERIALS AND METHODS
|
|---|
Basic Soil Properties and Land Use
The uptake of heavy metals in soils by a crop is controlled by soil organic matter (SOM) content, clay content, and pH. Basically, these soil properties largely determine the concentrations of heavy metals in pore water, and therefore the amount of heavy metals available to crops. Measurements of these basic soil properties have been derived from the Soil Information System of the Netherlands. From this database we selected all points with a laboratory measurement of these properties on samples taken from the topsoil, beneath the sod layer, if present, from the surface, to about 25 cm below the surface. Soil organic matter was determined by loss on ignition, clay was determined gravimetrically, and pH was determined in a 1 M KCl solution. The total number of points (measurements) was 4105 for clay, 4051 for SOM, and 4187 for pH. These measurements have been used to estimate frequency distributions for grouped units of the Soil Map of the Netherlands at a scale of 1:50 000.
We also required information about the actual land use because uptake of heavy metals varies for different crops, and because land use may affect the concentrations of heavy metals in soils. Previous research has shown, for instance, that Cd concentrations in uncultivated soils are less than in soils under cultivation (Van Drecht et al., 1996). In this study we used the land use database LGN3+, which is derived from satellite imagery from 1995 and 1997, topographic maps, and statistics on agriculture of the Centraal Bureau voor Statistiek (Thunnissen and de Wit, 2000). The 25- x 25-m pixels were aggregated to pixels of 500- x 500-m by calculating the land use class with the largest area.
Cadmium Data
Various datasets on Cd concentrations were used (Table 1). The first five datasets of Table 1 were also used by Tiktak et al. (1998) to estimate means of 500- x 500-m blocks. These datasets were supplemented by the data of the National Soil Quality Monitoring Network (sampling rounds 1993, 1994, and 1995), data of the Soil Quality Monitoring Networks of the provinces of Groningen, Friesland, Drenthe, Gelderland (sampling round 1998), Flevoland, Utrecht, Noord-Holland, Zuid-Holland, and Noord-Brabant, and a dataset collected in an occasional project in the province of Zeeland. The total number of measurements was 4094. The size and shape of the sampled plots (geostatistical support) differs between the datasets. Common supports are squares of 20 x 20 m2 (IB, IB/RIKILT, CCRX), 100 x 100 m2 (PSM-DRa), 150 x 150 m2 (ZLD), fields or forest stands (BLGG, PSM-GR, PSM-DR, PSM-GLD, SC-DLO), or all fields of a farm (NSM, PSM-NBR). Also, the depth and thickness of the sampled layer differs between the datasets. Commonly used layers are 0 to 20 cm for arable soils or 0 to 5 cm for grassland (BLGG, IB, IB/RIKILT), and 0 to 10 cm (regardless of land use). Soil samples were destructed with either HNO3, HCl, or aqua regia, and analyzed with flame atomic absorption spectrometry or the inductive coupled plasma technique. For the IB/RIKILIT dataset the soil samples were ashed at 450°C (Wiersma et al., 1986). Similar to Tiktak et al. (1998), we did not correct for differences in sampling depth or extraction method.
For the BLGG dataset consisting of 1086 points (Table 1), the geographical coordinates were unknown and therefore these data could not be used in a conventional interpolation method. However, as will be shown below, these data could be used to estimate the local prior means (a priori distributions) of the heavy metals.
All datasets were collected to assess baseline concentration levels, that is, concentration levels commonly found in soils at that time. Places that are known or suspected to be heavily polluted by point sources were either not sampled or omitted afterward. We also omitted samples taken from the recent flood plains along the rivers Maas and Rhine, and two outliers having a Cd concentration larger than the mean plus four times the standard deviation.
Critical Thresholds for Cadmium Concentrations in Soils Derived from Quality Standards for Crops
Quality standards in crops are the starting points for calculating critical thresholds of heavy metals, in this case cadmium, in soils. Once a model is derived that relates Cd concentration in crops to soils, then the soil Cd concentrations that result in exceeding the quality standards in crops can be estimated.
For arable land, wheat and sugar beet were taken as reference crops. For potatoes (Solanum tuberosum L.), a third commonly grown crop on arable soils in the Netherlands, there was no relation between Cd concentrations in soil and in potatoes. The Cd concentration in grass, maize, wheat, and sugar beet was modeled by the Freundlich equation:
 | [1] |
where Qplant and Qsoil are the Cd concentrations (mg kg-1 dry matter) in the crop and the soil, respectively, A is the Freundlich constant, and B is a nonlinearity coefficient. We assumed that log(A) is a linear function of pH, log(SOM), and log(clay):
 | [2] |
where pH is the soil pH (pHKCl), SOM is the soil organic matter content (in %), and clay is the soil clay content (in %). Combining Eq. [1] and [2] gives:
 | [3] |
Moreover, we took account of different values for B at different Cd concentrations in soil, by fitting a continuous piecewise linear model with one knot (Montgomery and Peck, 1992):
 | [4] |
where Qknot is the Cd concentration in the soil at the knot, and Iknot is an indicator having a value of 1 when Qsoil > Qknot, and 0 in other cases. Eq. [4] can also be written as:
 | [5] |
For each crop, all possible models were fitted, and the best model was selected on the basis of Mallow's Cp, under the following constraints (Montgomery and Peck, 1992): Qsoil must be part of the model, B and C must be positive, and A1 and A3 must be negative.
Note that when a piecewise linear model is not significantly better than a single linear model without a knot, a model is selected without the last term of Eq. [4]. The best location of the knot was estimated by fitting models for a range of values for the knot, then selecting the model with the smallest residual variance. Table 2 shows the estimated regression coefficients, the estimated knot, the percentage of variance accounted for, and the residual variance of the selected models.
View this table:
[in this window]
[in a new window]
|
Table 2. Estimated regression coefficients, percentage of variance accounted for (R2adj), and residual variance ( 2) of (piecewise) linear regression model for log(Cd) in four crops. n, Number of measurements; A0A3, B, C, see Eq. [4]; , not fitted.
|
|
The Cd concentration in the soil at which the quality standard for a given crop is exceeded, that is, the critical threshold for Cd concentrations in the soil, was calculated by rearranging the two equations of Eq. [5], and predicting the Cd concentration in the crop when the Cd concentration in the soil equals Qknot, Qplant,knot:
 | [6] |
Note that the critical threshold concentration in the soil depends on pH, SOM, and clay, and as a result varies in space, because the coefficient A depends on pH, SOM, and clay. For Qplant,crit we inserted the quality standard for Cd in fodder, which equals 1.0 mg kg-1 for grass and 0.5 mg kg-1 for maize and sugar beet, or the quality standard for Cd in food for human consumption, which equals 0.15 mg kg-1 for wheat (Projectgroep veterinaire milieuhygiene, 1997).
Critical Threshold for Cadmium Concentrations in Soils Derived from Quality Standards for Organs of Cattle
For grasssland, we also calculated critical thresholds for Cd in soil from the food quality standard for Cd in the kidney of cattle. Kidney appeared to be the most critical organ for Cd. In this case not only the transfer of Cd from soil to crop must be considered, but also the transfer from crop to cattle and Cd uptake through direct ingestion of soil. The internal Cd concentration of the kidney (mg kg-1) was thus calculated as (Projectgroep FBS, 1999):
 | [7] |
where Mgrass and Msoil are the amounts of grass and soil consumed, respectively (kg d-1), fpa is the transfer coefficient of Cd in grass to animal via roots, and fsa is the transfer coefficient of Cd in soil to animal. Uptake of Cd from drinking water was ignored because its contribution is negligibly small. The Cd concentration in grass was modeled in the same way as in maize and wheat (see previous section). Table 2 shows the model selected. Note that for grass the piecewise linear model was not significantly better than the model without the knot, Eq. [3]. Substitution of Qgrass by Qsoil and rewriting gives the following equation for calculating the critical threshold in soil from the critical threshold in kidney:
 | [8] |
Note that Qsoil,crit is on both sides of the equation, and therefore must be calculated by iteration. The values used for Mgrass and Msoil were 15 and 0.5 kg d-1, respectively (Projectgroep veterinaire milieuhygiene, 1997). Following Projectgroep FBS (1999) we assumed that the transfer of Cd from soil or grass into the kidney is equal, that is, fsa = fpa, and equals 2.99. Again, Qsoil,crit varies in space because the Cd concentration in grass is determined by the pH and the SOM of the soil (Table 2). The food quality standard for Cd in kidney is 2.5 mg kg-1 (Projectgroep veterinaire milieuhygiene, 1997).
Estimating the Actual and Critical Cadmium Concentrations in Soils by Simple Indicator Kriging with Local Prior Means
General Approach
If we know the values of the three basic soil properties in the regression model at a given point, then we can calculate the critical threshold concentration of Cd in soil, and we can establish if the actual Cd concentration in the soil at this point is above or below the critical threshold. The Cd concentrations and the soil properties were interpolated to map the probability that the critical threshold for Cd concentrations in the soil is exceeded. In establishing the probability of exceeding the critical threshold concentration, the uncertainty in both the interpolated Cd concentration and in the critical threshold were taken into account. Conventional kriging techniques focus on deriving optimal estimates of the contaminant concentration and the associated variance of the estimation error. To estimate the probability that the Cd concentration exceeds the critical threshold, we need the entire probability distributions of Cd and the soil properties at points. We may assume that the estimation errors follow normal distributions, but this is unrealistic in many situations. Other solutions are to use disjunctive kriging (Von Steiger et al., 1996) or indicator kriging (Goovaerts et al., 1997; van Meirvenne and Goovaerts, 2000). We chose indicator kriging because it can accommodate measurements less than the detection limit.
The cumulative distribution functions (cdfs) for SOM, clay, pH, and Cd were estimated for all 500- x 500-m grid nodes with agriculture according to LGN3+. All cdfs were estimated by simple indicator kriging with local prior means (Goovaerts, 1997, p. 284331). The method is described in four steps: (i) choosing a number of cutoff values, (ii) indicator coding of measurements, (iii) estimating the local prior means, (iv) updating the local prior means.
Earlier investigations (e.g., Bak et al., 1997; Tack et al., 1997; McGrath and McCormack, 1999) have shown that the Cd concentration in the soil is related to the clay and SOM content. The previous section showed that for a given crop the critical threshold of the Cd concentration in soil is a function of SOM, clay, and pH. Therefore, to calculate the probability that the actual Cd concentration exceeds the critical threshold, the correlation between the actual and critical Cd concentrations was taken into account. We did this by drawing Cd values from conditional cumulative distribution functions (ccdfs), that is, cumulative distribution functions conditional on soil parent material defined in terms of SOM and clay. Basically, all updated cumulative distribution functions are conditional on the measurements in the neighborhood of the grid node. Nevertheless, we only use the term conditional cumulative distribution function for Cd, which is conditional on the SOM and clay value at the grid node, whereas we use the term cumulative distribution function for SOM, clay, and pH. The exact values for SOM and clay at a grid node are unknown, and therefore we have to estimate all possible ccdfs of Cd at each node.
Choosing Cutoff Values
The indicator approach starts with choosing cutoff values that provide a reasonable discretization of the frequency distribution. We chose the nine deciles of the sample histogram as cutoff values. As we have special interest in large concentrations, we included two extra cutoff values in the upper tail of the distribution, the 95th percentile (P95) and 99th percentile (P99) (Table 3).
View this table:
[in this window]
[in a new window]
|
Table 3. Estimated percentiles of the frequency distribution for three soil properties and the soil Cd concentration used as cutoff values in indicator kriging.
|
|
Indicator Coding of Measurements
Each measurement was transformed into a K vector of indicators where K is the number of cutoffs such that i(x;zk) = 1 if z(x)
zk, or 0 otherwise. Here, i(x;zk) is the indicator value for cutoff value zk at the location with coordinates x. Note that for each point there are 11 indicators because we have 11 cutoff values. For example, a Cd concentration of 0.25 mg kg-1 would result in a vector of indicators equal to [0 0 0 1 1 1 1 1 1 1 1]' (Table 3). For cutoff values less than the detection limit, we do not know if the Cd concentration is less than these cutoff values or not, and consequently the indicator values are missing: i(x;zk) = 1 if zd(x)
zk, and missing otherwise. Here, zd(x) is the detection limit at location x. For example, for a Cd measurement less than a detection limit of 0.5, the vector of indicators becomes [* * * * * * * 1 1 1 1]'. Note that in this way these so-called left-censored measurements are not omitted completely in the interpolation but are used as far as they contain information.
Estimating Local Prior Means
The indicators are used to estimate cumulative frequencies, that is, areal fractions with values less than the cutoffs. The estimated cumulative frequencies are used as a priori estimates of the cumulative probabilities at points. For example, the estimated cumulative frequency of Cd concentrations
0.25 mg kg-1 is 0.40, and consequently at any point in the Netherlands the a priori cumulative lower probability of 0.25 Cd mg kg-1 can be estimated by 0.40. For SOM, clay, and pH, it is unrealistic to assume constant cumulative frequencies throughout the Netherlands. In geostatistical terms, it is unrealistic to assume that the means of the indicators are constant throughout the Netherlands. The soil map of the Netherlands at a scale of 1:50 000 was thus used to define subregions with assumed constant local means of the indicators. Soil organic matter and pH are influenced by land use, so for these two properties the soil map was combined with the land use database LGN3+. The resulting map units were grouped into 14, 13, and 13 categories for SOM, clay, and pH, respectively. The cumulative frequencies (means of indicators) for these categories were estimated by the measurements derived from the soil information system.
For Cd we considered four categories of parent material (Table 4) as well as two land use practices (cultivated and uncultivated). For uncultivated clay, uncultivated peat, and uncultivated peaty clay there were too few measurements (<100) to obtain reliable estimates of the cumulative frequencies. For these soil parent material categories no distinction was made in land use, resulting in five categories: uncultivated sand, cultivated sand, clay (undifferentiated), peat (undifferentiated), and peaty clay (undifferentiated). The local means of the indicators for Cd were estimated with the parent material (SOM and clay) and land use as recorded at the sample locations, and not as depicted on maps because of impurities and spatial variation within map units. Note that we could use the BLGG data (Table 1) with unknown geographical coordinates to estimate the local means of the indicators.
To estimate the local means of the indicators, special attention was paid to the 707 measurements less than detection limits. We used the following procedure to estimate the local means of the indicators:- Estimate the parameters of the lognormal distribution by maximum likelihood estimation, using all data including the censored data (Akritas et al., 1994);
- For any censored measurement, independently draw one random value from the part of the distribution that is less than the detection limit;
- Backtransform all values, including the randomly drawn values;
- Calculate the percentage of measurements less than the 11 cutoffs (means of indicators);
- Repeat Steps 2 to 4 one-hundred times, leading to 100 estimates of the indicator means;
- Calculate the averages of the 100 estimates of the indicator means.
Figure 1
shows the local prior means for the five categories.

View larger version (22K):
[in this window]
[in a new window]
|
Fig. 1. Cumulative frequency distributions of Cd (mg kg-1) in the topsoil for categories defined in terms of soil parent material and land use.
|
|
Updating the Local Prior Means
A first estimate of the discretized cumulative distribution function at a given node was provided by the local means of the indicators. Note that this a priori distribution is not location specific, that is, within a subregion the estimated cdf is constant. This is counter-intuitive because if a node is surrounded by sample locations with small values, then we expect a relatively small probability of exceeding a given cutoff value, and with high neighboring values we expect a relatively large probability. Therefore the local means were updated by interpolating the residuals of the indicators, that is, the indicators minus the local prior means. For example, for a node within the category "cultivated sand" and a Cd concentration of 0.25 mg kg-1, the vector of residuals equals:
For all sample locations except the BLGG data (geographical coordinates unknown), the vector with residuals was calculated. Next, the residuals were interpolated one by one to the grid nodes, that is, first the residuals of the first indicator were interpolated, then those for the second indicator, etc. For interpolation we used the program kt3d with the option simple kriging and expected value 0 of GSLIB (Deutsch and Journel, 1998). Experimental variograms were also estimated with GSLIB. Parameters of variogram models were estimated by weighted least squares, using the number of pairs as weights (Table 5). Fitted models were spherical with nugget (22 times), double spherical with nugget (16 times), and pure nugget (6 times). A double spherical model can be described by:
 | [9] |
where
(h) is the semivariance between points with mutual distance h, c0 is the nugget variance, c1 is the first structured variance component, a1 is the first range, c2 is the second structured variance component, and a2 is the second range.
The interpolated residuals were added to the a priori cumulative probabilities. The lower cumulative probability for a cutoff zk at the grid node with coordinates x0, the left-hand term in Eq. [10], was estimated by:
 | [10] |
where F(zk|s0) and F(zk|si) are the a priori lower cumulative probabilities for zk at the grid node x0 and at the measurement location xi, respectively,
(xi) is the simple kriging weight of the measurement location xi, and I(xi;zk) is the random indicator variable for the cutoff value zk at location xi.
The exact values of SOM and clay at the nodes are unknown, and this implies that the a priori distributions of Cd are unknown. We only know the cdfs of SOM and clay at the nodes, thus we can sample these cdfs a large number of times and determine the soil parent material category for each pair of SOM and clay so that we obtain a probability for each parent material category and the associated a priori distribution. Ignoring the uncertainty about land use, we have in principle four possible a priori distributions at each node. All four have been updated by adding the interpolated residuals leading to four series of K conditional cdf (ccdf) values per node.
Each (conditional) cdf value must lie in the interval [0,1] and each series of K (conditional) cdf values must be nondecreasing. We checked for these constraints, and order-relation deviations have been corrected by averaging the upward and downward corrections (Deutsch and Journel, 1998, p. 8186).
To simulate values by sampling the (conditional) cdf we need a complete (conditional) cdf, not only the cumulative probabilities at cutoff values. For soil property values between the smallest and largest cutoff, the (conditional) cumulative probabilities were interpolated linearly. For Cd and pH, probabilities less than the smallest cutoff were extrapolated by a power function (Goovaerts, 1997, p. 280). The exponent of the power function was estimated from the data by linear regression. Table 6 shows the estimated powers and arbitrarily chosen minima. There were very few measurements less than the smallest cutoff for SOM and clay. This implies that the a priori probabilities for the smallest cutoff were nearly 0. Consequently, we extrapolated linearly to a minimum of 0 for these properties. For all soil properties and Cd, the (conditional) cumulative probabilities beyond the largest cutoff were extrapolated linearly to a chosen maximum (Table 6).
View this table:
[in this window]
[in a new window]
|
Table 6. Minima, maxima, and exponent ( ) of power function used to extrapolate beyond the lower and upper tail of conditional cumulative distribution function of Cd. When = 1, the power function is equivalent to linear extrapolation.
|
|
Estimating the Probability of Exceeding the Critical Threshold Concentrations for Cadmium in Soils by Monte Carlo Simulation
Estimates of the probability of exceeding the critical threshold for Cd concentrations in soils at each node were made using Monte Carlo methods. The following algorithm describes the process:- Draw 100 values of SOM, clay, and pH values by Latin hypercube sampling (LHS) (McKay et al., 1979; Stein, 1987) from the respective cdfs of SOM, clay, and pH; this results in 100 triples of SOM, clay, and pH;
- Calculate the critical threshold, Qsoil,crit, from the first triple of SOM, clay, and pH;
- Determine the ccdf of Cd from SOM and clay of the first triple, and from the land use according to LGN3+;
- Draw 100 Cd values from this ccdf by LHS sampling;
- Compare the 100 Cd values of Step 3 with the Qsoil,crit value of Step 2, and count the number of times the Cd value is larger than Qsoil,crit value;
- Repeat Steps 2 to 5 for Triples 2100, and sum the number of times Qsoil,crit is exceeded;
- Estimate the probability by the total number of times Qsoil,crit is exceeded, divided by 10000.
At each node we thus obtained 100 drawn values for the critical threshold of Cd in soil and 10000 drawn values for the actual Cd concentration in soil. We used the median of these values as estimates of the critical and actual Cd concentrations in soil to make maps.
In Step 1, SOM, clay, and pH values were drawn independently, that is, we assumed that these properties were independent. To check the validity of this assumption we calculated correlation coefficients between these three properties in intersections of the SOM, clay, and pH subregions. It appeared that the correlation was generally weak.
 |
RESULTS AND DISCUSSION
|
|---|
Actual Cadmium Concentrations in Soil
The geographical pattern with estimated Cd concentrations resembles, more or less, the pattern with clustered units of the Soil Map of the Netherlands 1:50 000 (Fig. 2)
. In broad outlines, peat soils have the largest concentrations of Cd, fluvial and marine clay soils have intermediate concentrations, and sandy soils have the smallest concentrations. This illustrates the effect of clay and SOM on the actual Cd concentration. Sandy, cultivated areas have somewhat larger concentrations than their uncultivated counterparts, such as the ice-pushed ridges covered with forest and heath in the center of the Netherlands (Veluwe and Utrechtse Heuvelrug), and the coastal dunes. However, we cannot conclude from this that the Cd concentrations in agricultural soils are increased by the input of agrochemicals, because the natural background values of the cultivated sandy soils can also be relatively large (high weathering potential), and because less Cd may have been leached due to relatively high pH values in cultivated areas.

View larger version (61K):
[in this window]
[in a new window]
|
Fig. 2. Estimated Cd concentration (mg kg-1) in the topsoil. Median of 10 000 values simulated with the updated conditional distribution functions of Cd are taken as estimates. The class boundaries are chosen such that the number of cells in each class is approximately equal.
|
|
To find locations where the Cd concentrations are greater or less than the concentrations that can be expected on the basis of the local SOM and clay content, we repeated the Monte Carlo simulation using the a priori distributions of Cd as input instead of the updated distributions. For SOM and clay we used the updated distributions. In general, the differences between the "a priori" Cd concentrations and the "updated" Cd concentrations are small (Fig. 3)
. Three regions where the updated Cd concentrations are markedly different (absolute difference > 0.2 mg kg-1) from the a priori Cd concentrations are (i) the southeastern part of Brabant and the adjacent part of Limburg west of the river Maas; (ii) South Limburg, and (iii) the so-called "Groene Hart" in the west, that is, the area between Rotterdam, Amsterdam, and Utrecht. The first anomaly is caused by airborne pollution originating from smelters in Belgium and Brabant that emit Cd and Zn. (Van Wensem and Vegter, 1998). The large Cd concentrations in the "Groene Hart" (third anomaly) are possibly related to the use of compost from the cities such as Amsterdam. Areas with concentrations somewhat smaller than expected occur mainly in the northern part of the country.
Critical Threshold for Cadmium Concentrations in Soils
The critical threshold for Cd in soils differs considerably between the land use types (Fig. 4)
, and increases from arable soils, to maize land, to grassland. When sugar beet is taken as a reference crop this leads to critical thresholds somewhat smaller than for wheat as a reference crop, despite the fact that for sugar beet the quality standard (0.5 mg kg-1) is larger than for wheat (0.15 mg kg-1). For grassland, the critical thresholds for Cd concentrations in soils are on average very large, regardless of whether the quality standard for the kidney of cattle or the quality standard for grass was taken as a starting point. The quality standard based on kidneys of cattle leads to the smallest critical thresholds for Cd in soils. For this quality standard, the critical threshold for Cd in soils varies from 3.2 to 18 mg kg-1. The large spread is related to differences in soil parent material (clay and SOM). The critical threshold for maize land with Cd in soils varies from 1.6 to 4.2 mg kg-1, for arable land it varies from 0.06 to 2.1 mg kg-1 when sugar beet is taken as a reference crop, and from 0.02 to 5.0 mg kg-1 for wheat as a reference crop.

View larger version (13K):
[in this window]
[in a new window]
|
Fig. 4. Cumulative frequency distributions of estimated critical threshold for Cd concentrations in soils (mg kg-1) (A) and of estimated probability of exceeding this critical threshold (B). For grassland the probability of exceeding the critical threshold for Cd concentrations in soil is 0.0 nearly everywhere, when the quality standard for Cd in kidney (2.5 mg kg-1) or in grass (1.0 mg kg-1) is used to calculate the critical threshold for Cd in soils.
|
|
For arable soils, the critical threshold for Cd in soils for tillage of sugar beet is depicted (Fig. 5)
. For grassland the critical thresholds for Cd in soils were derived from the quality standard for kidneys of cattle. The geographical pattern is largely determined by the actual land use. For instance, the spottiness of the areas with cover sand in the eastern and the southern parts of the Netherlands is related to the occurrence of maize land, arable land, and grassland. Red pixels are pixels with predominantly maize or arable crops, both having relatively small critical thresholds in soil for Cd (<0.38 mg kg-1, see also Fig. 5), whereas the yellow and green pixels have grassland with much larger critical thresholds (6.6 to 8.1 mg kg-1). Uncultivated areas such as the ice-pushed ridges of the Veluwe and the Utrechtse Heuvelrug and the coastal dunes, are not mapped because the critical threshold for Cd in soils is not defined for these uncultivated areas.

View larger version (61K):
[in this window]
[in a new window]
|
Fig. 5. Estimated critical threshold for Cd (mg kg-1) in soils for the actual land use. For arable soils the reference crop is sugar beet. The class boundaries are chosen such that the number of cells in each class is approximately equal.
|
|
The effect of SOM, clay, and pH on the critical threshold in soils can also be seen on the map. For instance, the fluviatile clay soils in the Betuwe (in the middle of the Netherlands, along the river Rhine) with grassland have critical thresholds that are considerably larger than the adjacent sandy soils with grassland. The peat and peaty clay soils with grassland such as in the "Groene Hart" have very large critical thresholds (>10 mg kg-1). A similar effect can be seen for arable soils. The critical thresholds on marine clays with arable land (e.g., in Zeeland in the southwest, and in Friesland and Groningen in the north) are larger (1.7 to 2.0 mg kg-1) than on the sandy soils with arable land (e.g., in Drenthe and Groningen in the northeast) (0.06 to 0.38 mg kg-1).
Two remarks must be made. First, for all crops the regression models must be extrapolated for many nodes to derive the critical threshold for Cd in soil from the quality standards for the crops. For grass, maize, wheat, and sugar beet this was the case for 42, 33, 33, and 41% of the nodes, respectively (see Appendix). These large percentages imply that we are rather uncertain about the critical threshold for Cd in soils in considerable parts of the agricultural area. The second remark concerns the model for Cd in the kidney of cattle. We are rather confident of the model for the uptake of Cd in soil by grass, but we need more data on the transfer of Cd in soil and in grass to the kidney of cattle. Therefore, the critical thresholds for Cd in grassland soils derived from the quality standard for Cd in kidney are tentative only. On the other hand, we believe that the estimated critical thresholds are safe because we assumed that the transfer coefficient from soil to animal equals that of grass to animal, whereas it is probably smaller (Eq. [7]).
Probability of Exceeding the Critical Threshold for Cadmium Concentrations in Soil
The probability that the actual concentrations of Cd exceed the critical threshold in soil are negligibly small for the entire grassland area. For 98% of the grassland area the probability that the critical threshold for Cd in soil is exceeded is less then 0.1% using the critical threshold in grassland soils derived from the quality standard for kidneys of cattle. The probabilities of exceeding the thresholds are negligible for the peat and clay soils in the "Groene Hart" with estimated Cd concentrations in soil up to 2.7 mg kg-1 (Fig. 6)
. The large Cd concentrations in these soils are compensated by the large critical thresholds here. For maize land the probabilities of exceeding the critical threshold for Cd in soils are also small: only 3% of the present maize land area has a probability larger than 5%. However, for arable soils, these probabilities are much larger. For example, with sugar beet as a reference crop, 75% of these soils have probabilities larger than 5%, and 38% of the soils have a probability larger than 50%. For wheat as the reference crop, 82% of these soils have probabilities larger than 5%, and 12% have a probability larger than 50%. For the arable soils where the regression model was not extrapolated to derive critical thresholds, 57% have probabilities larger than 5%, and 7% has probabilities larger than 50% when sugar beet is the reference crop. Thus, based on the area where the regression model need not be extrapolated, the estimated areal percentages with exceedence probabilities larger then 5 or 50% are considerably reduced. For wheat as the reference crop, these areal percentages are reduced to 80% with a probability of exceeding the threshold of >5% and 5% with a probability of exceeding the threshold of >50%.

View larger version (65K):
[in this window]
[in a new window]
|
Fig. 6. Estimated probability of exceeding critical threshold for Cd in soils. For arable soils the reference crop is sugar beet.
|
|
 |
SUMMARY AND CONCLUSIONS
|
|---|
The critical threshold for Cd concentrations in soils differs considerably between the three land use types: grassland, maize land, and arable land. This is partly due to the differences in the quality standards for Cd in cattle kidney, maize and wheat, or sugar beet that are used as starting points for calculating the critical thresholds in grassland soils, maize land soils, and arable soils, respectively. For arable soils, the critical thresholds for Cd, derived from the quality standard for sugar beets, are somewhat smaller than for wheat, even though the quality standard for sugar beet is larger. Due to spatial variation in pH, SOM, and clay, the critical threshold for Cd in soils also varies in space within a given land use type. This spatial variation is considerable especially for grassland: the critical threshold for Cd in grassland soils varies from 3.2 to 18 mg kg-1.
For grassland the probability of exceeding the critical threshold for Cd in soils is negligible (<0.1%) everywhere, even in areas with high baseline concentrations of Cd in soil. For maize land the probabilities are small (<5%) almost everywhere. For arable soils, however, these probabilities are usually greater than 5% when sugar beet or wheat is taken as a reference crop, and locally exceed 50%.
To derive the critical thresholds for Cd in soils from the quality standards for Cd in grass, maize, wheat, and sugar beet, or the kidney of cattle, the regression models must be extrapolated in large parts of the Netherlands. For arable soils where the regression models need not be extrapolated, the areal percentages with high exceedence probabilities are much smaller than for the total area with arable soils.
Extrapolation with Regression Model
To determine whether the regression model is suitable for predicting the Cd concentration in the crop at a given prediction point (node), we might evaluate if all predictors individually fall within the range of the values at the original data points used to fit the model. However, for multiple regression models this method is inappropriate because there can be prediction points that satisfy the individual condition, but still fall outside the multidimensional cloud of original data. Therefore it is proposed to calculate, for all prediction points, the quantity h0 defined as (Montgomery and Peck, 1992):
 | [A1] |
where x0' is the p vector with values of the predictors at the prediction point (p is the number of regression coefficients), and X is the (n x p) matrix with values of the predictors at the n data points used to fit the model. When h0 > 2p/n, the model must be extrapolated, and the predictions and prediction variances are not reliable.
If we insert the estimated critical threshold for Cd in the soil at a given node in vector x0 of Eq. [A1], and h0
2p/n, then we can safely derive the critical threshold for Cd in soil from the quality standard for Cd in the crop at this node.
 |
ACKNOWLEDGMENTS
|
|---|
We are very grateful to the authorities of the provinces of Groningen, Friesland, Drenthe, Flevoland, Gelderland, Noord-Holland, Zuid-Holland, Utrecht, Zeeland, and Noord-Brabant for giving permission to use the data of the Provincial Soil Quality Monitoring Networks.
 |
REFERENCES
|
|---|
- Akritas, M.G., T.F. Ruscitti, and G.P. Patil. 1994. Statistical analysis of censored environmental data. p. 221240. In G.P. Patil and C.R. Rao (ed.) Handbook of statistics. Vol. 12. Environmental statistics. North-Holland, Amsterdam.
- Bak, J., J. Jensen, M.M. Larsen, G. Pritzl, and J. Scott-Fordsmand. 1997. A heavy metal monitoring programme in Denmark. Sci. Total Environ. 207:179186.
- Deutsch, C.V., and A.G. Journel. 1998. GSLIB geostatistical software library and user's guide. 2nd edition. Oxford Univ. Press, New York.
- Goovaerts, P. 1997. Geostatistics for natural resources evaluation. Oxford Univ. Press, New York.
- Goovaerts, P., R. Webster, and J.P. Dubois. 1997. Assessing the risk of soil contamination in the Swiss Jura using indicator geostatsitics. Environ. Ecol. Stat. 4:3148.
- Landbouwadviescommissie milieukritische stoffen. 1991. LAC-signaalwaarden. Ministry of Agric., Nature Conservation and Fishery, the Hague, the Netherlands.
- McGrath, D., and R.J. McCormack. 1999. The significance of heavy metal and organic micropollutants in soils. Rural Environ. Ser. 23. Johnstown Castle Res. Centre, Wexford, Ireland.
- McKay, M.D., W.J. Conover, and R.J. Beckman. 1979. A comparison of three methods for selecting values of input variables in the analysis of output from a computer code. Technometrics 21:239245.
- Montgomery, D.C., and E.A. Peck. 1992. Introduction to linear regression analysis. Wiley, New York.
- Projectgroep FBS. 1999. Function-specific soil quality evaluation system 2: Function-specific soil quality target values. Ministry of Agric., Nature Conservation and Fisheries, IKC-Agriculture, Ede, the Netherlands.
- Projectgroep veterinaire milieuhygiene. 1997. Veterinaire milieuhygienewijzer 1997. Veterinaire Hoofdinspectie van de Volksgezondheid, Rijswijk, the Netherlands.
- Stein, M. 1987. Large sample properties of simulations using latin hypercube sampling. Technometrics 29:143151.
- Tack, F.M.G., M.G. Verloo, L. Vanmechelen, and E. Van Ranst. 1997. Baseline concentration levels of trace elements as a function of clay and organic carbon contents in soils in Flanders (Belgium). Sci. Total Environ. 201:113123.
- Thunnissen, H., and A. de Wit. 2000. The national land cover database of the Netherlands. Proc. 14th Congr. of the Int. Soc. Photogrammetry and Remote Sensing (ISPRS), Amsterdam. 1623 July 2000. GITC, Lemmer, the Netherlands.
- Tiktak, A., J.R.M. Alkemade, J.J.M. van Grinsven, and G.B. Makaske. 1998. Modelling cadmium accumulation at a regional scale in the Netherlands. Nutr. Cycling Agroecosyst. 50:209222.
- Van Drecht, G., L.J.M. Boumans, D. Fraters, H.F.R. Reijnders, and W. van Duijvenbooden. 1996. National images of metal load from non-point sources to the soil and of metal concentrations in the topsoil, as well as the relation between concentrations and load. (In Dutch, with English abstract.) Rep. 714801006. RIVM, the Netherlands.
- Van Driel, W., and K.W. Smilde. 1981. Heavy-metal contents of Dutch arable soils. p. 305313. In Sonderdruck aus Landwirtschaftliche Forschung. Zeitschrift des Verbandes Deutscher Landwirstschaftlicher Untersuchungs-und Forschungsanstalten. Sonderheft 38. Kongressband 1981. J.D. Sauerlander's Verlag, Frankfurt am Main, Germany.
- Van Meirvenne, M., and P. Goovaerts. 2000. Evaluating the probability of exceeding a site-specific soil cadmium contamination threshold. Geoderma 102:75100.
- Van Wensem, J., and J.J. Vegter. 1998. Large scale heavy metal pollution: The case for an integrated risk and land-use management approach. Volume 1. Proc. 6th Int. FZK-TNO Conf. Cont. Soil, Edinburgh, UK. 1721 May 1998. Thomas Ltd., Telford, UK.
- Von Steiger, B., R. Webster, R. Schulin, and R. Lehmann. 1996. Mapping heavy metals in polluted soil by disjunctive kriging. Environ. Pollut. 94:205215.[Medline]
- Wiersma, D., B.J. Van Goor, and N.G. Van der Veen. 1986. Cadmium, lead, mercury, and arsenic concentrations in crops and corresponding soils in the Netherlands. J. Agric. Food Chem. 34:10671074.
This article has been cited by other articles:

|
 |

|
 |
 
W. Li, C. Zhang, J. E. Burt, A.-X. Zhu, and J. Feyen
Two-dimensional Markov Chain Simulation of Soil Type Spatial Distribution
Soil Sci. Soc. Am. J.,
September 1, 2004;
68(5):
1479 - 1490.
[Abstract]
[Full Text]
[PDF]
|
 |
|

|
 |

|
 |
 
D. J. Brus and M. J. W. Jansen
Uncertainty and Sensitivity Analysis of Spatial Predictions of Heavy Metals in Wheat
J. Environ. Qual.,
May 1, 2004;
33(3):
882 - 890.
[Abstract]
[Full Text]
[PDF]
|
 |
|