A meta-analysis of published semivariograms to determine sample size requirements for assessment of heavy metal concentrations at contaminated sites

Additional keywords: confidence interval, guideline exceedance, mapping, simulation, soil contamination.

Soil contamination is a leading global issue and is becoming more prominent with past and present industrial activity, increasing population density, and resulting urban expansion (Andronikov et al. 2000; Lee et al. 2006; Lacarce et al. 2012). An area is classed as ‘contaminated’ when substances, such as heavy metals, exceed background concentrations and this is often due to anthropogenic influences (National Environment Protection (Assessment of Site Contamination) Measure 1999, ASC NEPM)). Heavy metals are of particular concern as they are influenced by human activity and remain in the soil for a considerable amount of time (Markus and McBratney 1996). Soil contamination, especially heavy metal contamination, can have substantial impacts upon human and environmental health (Helios Rybicka 1996; Mielke et al. 1999; Khan et al. 2008) and it is therefore crucial that contaminated areas are identified and managed.

If substances occur above established guideline values at a test site, there is a need for further investigation and potential remediation (ASC NEPM; VROM 2000; BC MoE 2014). Site assessment and remediation is costly, with a report by the NSW Environment Protection Authority (NSW EPA 2013) estimating the cost to be AU$100–$200 million per year for testing and clean-up of contaminated sites. Testing methods vary depending on the site and contaminant, and are often determined using an ad hoc approach based on a priori knowledge of past activities in the area, the site-specific risk assessment, and the associated conceptual model developed in the preliminary investigation of the site (ASC NEPM; US EPA 2002a; Glavin and Hooda 2005).

Government organisations may regulate the sampling conducted by establishing recommended sampling procedures (VROM 2000; Theocharopoulos et al. 2001). In Europe, where many different guidelines are established, sampling methods may instead vary from organisation to organisation (Theocharopoulos et al. 2001; de Zorzi et al. 2008). Such a variety of methods can have significant impacts upon the accuracy of analyses and in turn affect the precision of the remediation process, resulting in false positive or false negative results (Cattle et al. 2002). Contaminated site analysis needs to be timely and cost-effective, and one way to enhance this is by improving the sampling scheme, which relates to collection of sufficient samples and locating them so they identify contamination levels appropriately.

In Australia, sampling is conducted in the preliminary investigation and the results used to inform subsequent detailed investigation (ASC NEPM). Preliminary sampling is conducted with the purpose of estimating mean contamination at the site and comparing this, along with the 95% upper confidence interval, to the developed health investigation levels (HILs) and if these levels are exceeded, further investigation is required (ASC NEPM). Further investigation involves conducting detailed sampling at the site to provide a better picture or distribution of contaminants, with the number of samples to collect guided by the conceptual site model and results of the initial sampling (ASC NEPM).

Delineating contamination at a site can use either a design- or model-based approach, depending on the purpose of the investigation (i.e. preliminary or detailed). Design-based, or probability-based, sampling approaches select points based on probability and randomisation to predict a mean or another statistic, whereas model-based approaches, also known as purposive or judgemental approaches, rely on a model for predicting means (De Gruijter et al. 2006). Many studies adopt a model-based, or systematic, sampling scheme where a grid is selected and samples are taken at set intervals (US EPA 2002b). As a result, the samples are not independent of each other and will provide a biased variance. In situations in which we wish to estimate the mean contamination or the 95% confidence interval, a probabilistic sampling scheme is best (ASC NEPM; US EPA 2002b).

In contrast, grid-based sampling is highly useful and has been recommended for mapping and kriging variables (Pettitt and McBratney 1993). Studies have also complemented grid designs with shorter scale samples taken so as to model the short-range spatial variation that may not otherwise be detected (Lark 2002; Karunaratne et al. 2014). Detecting shorter-range variation and reducing bias of samples assist in increasing precision in mapping, thus preventing site misclassification.

Aside from using a suitable sampling scheme, the number of samples (herein referred to as the ‘sample size’) can affect the precision and reliability of contaminated site assessment. Error may arise from insufficient sampling, or failing to sample in areas where contaminant concentrations are higher (i.e. ‘hotspots’) and in turn increase risk of misclassification (Tiller 1992; Cattle et al. 2002). Studies have sought to quantify sample size; however, they have been applied to assessment of soil properties such as loss on ignition, soil texture, and pH, rather than specifically for contaminants, which can be highly variable (McBratney and Webster 1983; Kerry and Oliver 2004). Past studies have also sought to determine suitable grid spacing for soil spatial analyses (Chang et al. 1998; McBratney and Pringle 1999); McBratney and Pringle (1999) used published semivariograms in the analyses. There are several studies that have conducted sampling at contaminated sites and so the question arises as to whether it would be possible to use these published studies in as part of a meta-analysis to estimate sample size requirements and provide broad guidelines for investigation.

Many published studies assessing heavy metal distribution use model-based sample approaches, which, as mentioned earlier, will have biased variances. There are three general types of studies available to predict variation at a site: those that adopt a probabilistic design, those that adopt a systematic sampling scheme but do not provide a semivariogram, and those that adopt a systematic sample design and do provide semivariograms. The latter of these were used throughout the current study. If summary statistics and variogram parameters are provided by the study, it is possible to extract an unbiased variance using a method proposed by Domburg et al. (1994).

By conducting a meta-analysis utilising variogram parameters from several compiled contamination studies, this study aimed to: (1) determine the optimal sample size for estimating mean heavy metal concentration at a site and calculation of the proportion of the site exceeding the set guideline value, and (2) determine suitable sample sizes for mapping heavy metal distribution by simulating the spatial variation described by the semivariograms from each of the studies. The outcomes of this study will provide broad indications of the range of samples required for contaminated site assessment, which may improve accuracy of decision making and improve efficiency of assessment.

A literature search was performed to find peer-reviewed research on the spatial variation of heavy metal contaminants. Studies were chosen if they provided experimental semivariograms, topsoil samples (up to 10 cm) and if the study area was less than 5 km². The studies chosen covered a variety of land uses including industrial, agricultural, mine sites, and residential areas (Table 1). Data were obtained for a range of heavy metals, with the most commonly provided being lead (Pb), zinc (Zn), and cadmium (Cd), and so these were used for determining sample sizes for mean estimates of concentrations and for mapping contaminant distribution.

Table 1.  Summary of location, land use, and sampling designs for each chosen study

The majority of the studies added to the database used logarithmic transformations on skewed heavy metal concentrations and the transformed values were used in favour of raw values for later analyses. Some studies used transformed values in their analyses, but only provided untransformed means; to overcome this, the raw means were transformed using an equation derived from arithmetic moments and the lognormal distribution:

where μ_ln is the normalised mean, μ_raw is the skewed (raw) mean, and σ_raw is the skewed standard deviation.

Throughout this study, heavy metal concentrations were compared with the most conservative HILs provided by the National Environment Protection (Assessment of Site Contamination) Act 1994 (Cth), which are used in Australia for contaminated site assessment. Guideline values for Pb, Zn, and Cd are 300 mg kg^–1, 7400 mg kg^–1 and 20 mg kg^–1 respectively.

After compiling the variogram database we considered two situations: one where we wish to identify how many samples are required to estimate whether the mean contaminant concentration is below or above a guideline value, and another for mapping the concentration of the contaminant to identify whether the concentration at each spatial location was above or below the guideline value.

Equation 2 was used to calculate the 95% confidence interval (CI) around the mean. Given an estimate of the mean () or the variance (s²), the sample size (n) impacts on the width of CI due to the t_crit and the standard error of the mean, , decreasing with increasing sample size:

It was hypothesised that the mean would not be equal to the guideline value, whereas if the null hypothesis were proven true, the mean would be equal to the guideline value. Therefore, if the 95% confidence interval did not overlap with the contaminant guideline value, the null hypothesis would be rejected.

However, to be able to calculate the required sample size, it needs to be ensured that the variances were unbiased. Many of the selected studies used systematic sample designs, for example a grid or transect design (Table 1). Since the designs had no probabilistic component the estimated variances may be biased. To extract an unbiased estimate of the variance for each study included in the meta-analysis, the method established by Domburg et al. (1994) was implemented (Eqn 3). This method estimates the variance of the sample mean of an area using the semivariogram parameters:

where is the mean semivariance calculated using variogram parameters provided by the study and n is the sample size used in the study.

Calculation of the mean semivariance () was accomplished by simulating a matrix of 1000 sets of random pairs of coordinates within a simulated study area in the shape of a square and equal in size to each published study. All analyses were performed using R (R Core Team 2016).

Once the mean semivariance was calculated for each study, unbiased variance across each study area was obtained using Eqn 3. Subsequently, sample size requirements were quantified using Eqn 2, with a t_crit of 0.025, the degrees of freedom (n − 1) obtained in each study, and substituting the newly calculated unbiased variance of the mean into the equation. The sample size was increased until the 95% CI was above or below the guideline value.

Although the mean and confidence interval may be under the guideline value, some part of the site could still be contaminated. We estimated this for each study based on the cumulative upper probability above the guideline value using the t-distribution in Genstat 16 (VSN International Ltd 2013), based upon the optimal sample size.

As a case study, data from a subset of studies quantifying Pb were selected for simulation. Unconditional Gaussian simulation incorporating variogram parameters from each study was used to simulate contaminant distribution via the ‘gstat’ package in R (Pebesma 2004). A total of 1000 simulations were created for each study, with the set grid area closely resembling the original study area to maintain proportions.

From these simulated fields, the sample design consisted of 90% of locations being on a grid and the remaining 10% taken a short distance from randomly selected grid points. The latter portion of samples would provide short range samples for variogram modelling, similar to the approach suggested by Lark (2002). The range of sample sizes tested included 40, 60, 80, 100, 120, 140, 160, 180, 200, 250, 300, 350, 400, 450, and 500 samples. Each tested sample size was used to map soil contamination using ordinary kriging for mapping onto the grid.

The simulated grid values were assigned a class based on whether the Pb concentration exceeded (thus classed as ‘contaminated’) or was less than the guideline value (and so classed as ‘uncontaminated’). These classifications of the simulated values formed the ‘observed’ classes, or ‘truth’. The predicted values underwent a similar process; however, the guideline values were compared with the prediction intervals of each point, and used an additional classification of ‘high uncertainty’ where the guide value lay within the prediction interval (Fig. 1). The allocated classes of corresponding point locations within the observed and predicted datasets were then compared, and if they matched, they were assigned the class of ‘correct’, whereas if they mismatched, they were classed as ‘error’ (Including Type I or II); if the prediction interval indicated high uncertainty, the point was classed as ‘uncertain’.

Fig. 1.  Classification guide for predictions (adapted from Johnson et al. (2017)). PI refers to prediction intervals.

The compiled studies incorporated a variety of sample sizes, mean concentrations, and semivariogram parameters (Table 2). Sample sizes ranged from 48 to 665 samples and these were taken over areas from 0.0004 to 3500 km². Some studies presented mean values as logged, others presented the raw values and logged them before variogram analyses, and some used the raw values without transformation. Variogram models included linear, spherical, double spherical, Gaussian, exponential, and stable. The majority of nugget values were less than 1.0; however, some studies obtained a nugget semivariance of over 151000 as the data were not transformed before analysis. Distance parameters were often related to area size (i.e. the larger the area, the larger the distance parameter), yet this was not always the case. Mean values for many of the studies exceeded the most conservative ASC NEPM guideline values (20 mg kg^–1 for Cd, 6000 mg kg^–1 for Cu, and 300 mg kg^–1 for Pb).

Table 2.  Summary of parameters from each compiled study for each metal. Mean is as provided in original study, transform describes the type of (if any) transformation used by original study
n, sample size; c0, nugget semivariance; c1, c2, structural semivariance; d1, d2, distance parameters

The predicted sample sizes required for each study were very different to the sample sizes actually used (Table 3); this is justified as the purpose of these studies was to map contamination rather than estimate mean content. More samples were required for Pb compared with Cd and Cu. Overall, the majority of sample sizes required to determine whether the mean concentration was significantly different from the guideline value was four samples or less, which is quite a small number of samples in relation to the size of the study area. Some studies required a much greater sample size, such as that by Shi et al. (2008) for Cd, which required 108 samples; and the study by Bourennane et al. (2006), which required 122 samples to determine whether Pb exceeded the guideline value. Key drivers of the sample size requirement were the variability and how close the mean was to the guideline value.

Table 3.  Sample size (n), mean and transformation described in the original studies, the associated ASC NEPM health investigation levels (HILs), calculated unbiased variance, predicted sample size requirement for estimating the mean and proportion of each site exceeding the established HILs.

Although we could calculate the mean and show that it was significantly greater or less than the guideline value, some of the site may still contain concentrations that exceed the guideline. To detect whether or not areas of elevated concentrations were missed, the proportion of each area exceeding the guideline was calculated based on the estimated sample size (Table 3). Overall, the proportion of each site exceeding the guideline values for Cd and Cu were very low, except for the study by Shi et al. (2008) observing Cd. There were two studies in which >90% of the site exceeded the guideline value for Pb; however, it was studies in which the mean did not exceed the guideline but showed a proportion of values exceeding the guideline that would be at greatest risk of misclassification. Contaminant distribution is most often heterogeneous and obtaining a mean value does not provide an accurate representation, increasing risk of misclassification. More accurate representation could be obtained by mapping contamination at a site.

A subset of study parameters was taken from the Pb dataset and used in a simulation exercise to determine optimal sample size for mapping. Study area size (and therefore size of simulated fields) varied among studies, and as a result, grid intervals also varied to ensure the total number of grid points in the simulated fields remained consistent (Table 4).

Table 4.  Grid specifications for simulations of Pb.

Comparison between observed (simulated) and predicted classifications indicated that increasing sample size increased the number of points classified as correct while decreasing the amount of uncertainty (Fig. 2). Increasing sample size also resulted in an increase in the amount of error; however, the rate of increase was relatively small (maximum 5% increase overall). Regardless of sample size, uncertainty of predictions remained much higher than the number of samples classified correctly, where the majority of the studies had >50% of points classed as uncertain. The rate of change varied between studies, with larger study areas having lower rates, but there was generally no more than 10% increase in correct classification. There was an inflection in the results from simulation of Ersoy et al. (2008) (Fig. 2) and substantial variation throughout the categories for that specific study.

Fig. 2.  Comparison between classification and sample size for each subset Pb study.

Rather than classifying based upon a prediction that is common among studies, the current study classified based upon calculated prediction intervals. Based on the calculations, the prediction intervals presented in the study are quite large and sample size does not greatly affect them. In terms of selecting a suitable sample size, this varied for each study and the number of points correctly classified stabilised at around 200 samples (except the first set of results from Ersoy et al. (2008)), but was slightly better around 500 samples.

To understand the underlying mechanisms for the trends, Pearson correlations were used to compare the % values for each category with variogram parameters of each study, mean, area and sampling density (Table 5). Nugget semivariance (c₀) influenced all categories negatively, but the greatest indicated the higher the nugget semivariance, the lower the number of points classified as ‘correct’. As the structural semivariance (c₁) increased, error also increased, but this did not seem to influence uncertainty (r = 0.183). The distance parameter for the variogram was also related to error (r = 0.533) as well as the mean, and error decreased with increasing mean (r = −0.462). Sample area also shared a relationship with the magnitude of error and uncertainty (r = 0.681 and 0.445 respectively), both parameters increasing as area increased.

Table 5.  Pearson correlations between original features of subset studies (variogram parameters, mean, area, and density) and the simulation classifications. Areas highlighted in bold type indicate stronger correlations.
n, sample size; c0, nugget semivariance; c1, c2, structural semivariance; d1, d2, distance parameters

This study used a novel method to provide an indication of broad sample size requirements for assessment of heavy metal contamination by obtaining variogram parameters from a variety of studies. Using previous studies provides more realistic values, rather than simulating an ideal (and often less likely) scenario, as contamination can be highly variable both spatially and in magnitude (Glavin and Hooda 2005). The studies obtained encompassed a variety of land uses and sources of contamination, helping make this study more widely applicable to several situations.

The number of samples required to show whether the 95% CI of the mean concentration exceeded the guideline values was relatively small. The majority required four samples or less; however, some exceeded this, most likely due to the variance in relation to the sample size, inducing a wider confidence interval compared with other studies. This could be attributed to various influencing factors including the original contamination source, its location, historical land use, soil morphology, climate, and anthropogenic factors affecting environmental dynamics.

Calculation of the proportion of each site exceeding the guideline was valuable as it demonstrated that although the mean may not exceed the guideline, there were still points on the site that did exceed it. The inverse was also true, with some means exceeding the guideline, yet some proportion of the site did not exceed it. Detailed investigation may negate this issue if sufficient samples are taken. Furthermore, it was useful to observe the proportion of each site exceeding the guideline values and compare these with the sample size requirements as studies with higher sample requirements possessed greater variability at the site, which may have been missed if only a small number of samples were taken. The outcomes of this research demonstrate the importance of calculating proportion of the site exceeding the guideline values, rather than relying upon the mean value alone.

Although it is useful to use the mean to obtain an overall indication of the concentration at a site, especially if samples are limited; there is much variation that may be missed. This is especially true in terms of soil contamination as there may be a mix of diffuse and point sources; these point sources may be missed, potentially undermining scientific findings and resulting in failure to remediate where it is required (Marchant et al. 2011). It is therefore more reliable if a site were mapped using interpolation to obtain more precise predictions.

A suitable sample size was suggested to be around 200 samples as the rate of change in classification error and uncertainty stabilised around this value. Based upon prediction intervals and determined classes, the lowest amount of uncertainty was detected around n = 500, indicating the more samples, the less the uncertainty, which may not be realistic in many cases. Obtaining 500 samples over any study area would be very costly both economically and temporally, and regardless of sample size, uncertainty was still quite high with little decrease with increasing sample size. With increasing sample size, the amount of classification error increased, rather than decreased, so this would also need to be weighed up when deciding upon the number of samples to collect. Contaminated sites, especially following environmental catastrophes, require urgent and timely assessment and so unless the resources can be afforded, collection of so many samples would not be economically sound. Complementary methods in both physical assessment and mapping are worth exploring; for example, the use of proximal sensing would allow the collection of more samples (Horta et al. 2015), and the use of covariates as supplementary data may provide a more accurate map of contaminant distribution (Johnson et al. 2017).

It is essential to use unbiased variances when calculating CIs for comparison with guideline values. To ensure unbiased variances it is recommended to use design-based sampling schemes, or, if not practical due to site access and cost, use correction methods such as the Domburg equation (Domburg et al. 1994). It could also be useful to explore sampling methods that are able to take available information, such as that obtained in initial site investigation, into account. Building of a conceptual model for contaminated site assessment depends upon collection of information such as site history, layout, and topography (ASC NEPM). This site information may also be useful in development of a less biased sampling scheme. Conditioned Latin Hypercube sampling derives a sampling scheme that takes ancillary variables into account, thus reducing bias (Minasny and McBratney 2006).

It is good to estimate the mean contaminant concentration of a site; however, to delineate areas that exceed contamination guideline values and inform remediation, greater detail needs to be used to reconcile cost and efficiency. The research findings in this paper and those presented in the compiled studies provide further evidence of the importance of mapping and interpolation.

The number of samples required for estimating whether the mean exceeds the guideline value were very low. However, evident through the calculation of proportion of site contamination, estimating the mean may miss a large portion of the variation at the site, especially as heavy metal contamination is highly variable. Therefore, it is better to use interpolation methods such as kriging to detect this variation. Estimates of plausible sample sizes for mapping a site was estimated at 200, with 500 samples resulting in the lowest amount of uncertainty. Collecting so many samples may be unrealistic in many cases and therefore shows that sample sizes and schemes are site-specific. To improve accuracy, it would be worth exploring improving efficiency in other facets of contamination assessment, such as in detection and in the reporting stage.

The authors declare no conflicts of interest.