Next Article in Journal
A Call for More Snow Sampling
Next Article in Special Issue
Reservoir Quality of Upper Jurassic Corallian Sandstones, Weald Basin, UK
Previous Article in Journal
Fluid-Rock Interactions in a Paleo-Geothermal Reservoir (Noble Hills Granite, California, USA). Part 2: The Influence of Fracturing on Granite Alteration Processes and Fluid Circulation at Low to Moderate Regional Strain
Previous Article in Special Issue
Two-Scale Investigation of the Retention Behavior of a Well-Graded Mixed Soil
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Evaluating Spatial Regression-Informed Cokriging of Metals in Soils near Abandoned Mines in Bumpus Cove, Tennessee, USA

1
Colliers Engineering & Design, Charlotte, NC 28217, USA
2
Department of Geosciences, East Tennessee State University, Johnson City, TN 37614, USA
*
Author to whom correspondence should be addressed.
Submission received: 20 September 2021 / Revised: 12 October 2021 / Accepted: 14 October 2021 / Published: 20 October 2021
(This article belongs to the Collection Early Career Scientists’ (ECS) Contributions to Geosciences)

Abstract

:
Inorganic contaminants, including potentially toxic metals (PTMs), originating from un-reclaimed abandoned mine areas may accumulate in soils and present significant distress to environmental and public health. The ability to generate realistic spatial distribution models of such contamination is important for risk assessment and remedial planning of sites where this has occurred. This study evaluated the prediction accuracy of optimized ordinary kriging compared to spatial regression-informed cokriging for PTMs (Zn, Mn, Cu, Pb, and Cd) in soils near abandoned mines in Bumpus Cove, Tennessee, USA. Cokriging variables and neighborhood sizes were systematically selected from prior statistical analyses based on the association with PTM transport and soil physico-chemical properties (soil texture, moisture content, bulk density, pH, cation exchange capacity (CEC), and total organic carbon (TOC)). A log transform was applied to fit the frequency histograms to a normal distribution. Superior models were chosen based on six diagnostics (ME, RMS, MES, RMSS, ASE, and ASE-RMS), which produced mixed results. Cokriging models were preferred for Mn, Zn, Cu, and Cd, whereas ordinary kriging yielded better model results for Pb. This study determined that the preliminary process of developing spatial regression models, thus enabling the selection of contributing soil properties, can improve the interpolation accuracy of PTMs in abandoned mine sites.

1. Introduction

Mining industry contamination is a major environmental problem and may arise in soils as a result of associated blasting, transportation, ore extraction, and materials processing, and the adverse effects may be exacerbated by associated deforestation [1]. Perhaps most notably, mining activities discharge potentially toxic metals (PTMs) into the environment. PTMs that commonly accumulate in soils include Cadmium (Cd), Chromium (Cr), Copper (Cu), Mercury (Hg), Nickel (Ni), Lead (Pb), Manganese (Mn), and Zinc (Zn) [2]. Therefore, locating the source and distribution of PTMs in soil can provide a scientific basis for environmental risk assessment and soil reclamation.
Natural manifestations of toxic metals in soils are significantly related to the intricate transport of metalliferous parent bedrock and accumulate in soil by means of geochemical weathering and/or pollutants from atmospheric deposition [3,4]. Unlike the majority of pollutants from organic origin, metals do not biodegrade and may disrupt the environment by accumulating in biological organisms and/or are biologically transformed into organic complexes [5,6]. Such activity increases toxicity and potential for remobilization and may directly endanger human health, agriculture, and ecotoxicology [6,7].
Mining of metalliferous ores enriches concentrations of PTMs in soils from the presence of unmanaged waste rock, mine tailings, and slag. Soils in the vicinity of mines and tailings are often low in organic matter and macro- or micro-nutrients due to the presence of high PTMs [8,9]. Therefore, accretion of PTMs in soils at un-restored mine sites may endanger the environmental ecosystem and public health. Except in highly acidic soils, they occur primarily as insoluble compounds [10]. Insolubility affects PTM transport in soils, with above 90% found in surficial soils (to 15 cm depth) [11]. Groundwater contamination may arise due to slow metal transport through soil and saprolite [10]. A range of factors may interfere with PTM pollution near abandoned mine sites, including current and historic mineral resource inventory, extraction, and processing methods and chronology, as well as the surrounding geomorphology and ecology [12]. PTM accumulation is further complicated because some metals demonstrate environmental idiosyncrasies [12]. Bioavailability, movement, and spatial distribution of PTMs in soils is also affected by total concentration [13,14], pH [15], organic matter [16], texture [17,18], cation exchange capacity (CEC) [19], and redox reactions [5,20]. In this regard, researchers evaluated various soil physico-chemical properties in relation to PTM spatial accumulation, determined their influence on and correlation with the PTMs, and estimated the concentration of PTMs based on the soil properties [21,22,23]. In particular, clay content, soil pH, and CEC were known to be highly correlated with Co, Cu, Cr, Ni, Pb, and Zn concentrations in Spain’s Northern Plateau [21], while physico-chemical soil properties contributed to the spatial distribution of heavy metals in and surrounding an industrial plant in China [23]. In prior research by authors of the current study, spatially weighted multivariate regression models created for Zn, Pb, Cu, Mn, and Cd using soil physico-chemical properties yielded better results over ordinary least squares regression models in an abandoned mine site in Bumpus Cove, TN [22].
Another commonly used approach to model environmental data is spatial interpolation, which may include such techniques as inverse distance weighting (IDW), artificial neural networks (ANN), and kriging. These methods have been used to interpolate groundwater geochemistry [24], depth to a geologic marker [25,26], and geologic reservoir parameters [27]. Similarly, interpolation can help predict the spatial distribution of PTM pollution in unsampled areas that can aid to refine environmental site assessment. Therefore, many researchers have interpolated the spatial distribution of PTMs and other elements of concern [28,29,30,31,32,33]. Kriging offers statistical advantages and is favored to IDW when variogram parameters are known [34], and ordinary kriging was shown to outperform ANN in mapping depth to a geologic marker [25].
Kriging interpolation is a broadly used spatial interpolation technique that can provide a linear optimized unbiased estimate at unsampled locations within a study area [35,36]. In particular, the ordinary kriging method, based on measured data and a semivariogram model to predict unknown points, is frequently used to study the spatial distribution characteristics of PTMs in soils [37,38,39]. It is important to obtain the most realistic prediction model on the distribution of PTMs. For this purpose, researchers have found that compared to univariate kriging method, multivariate cokriging methods that incorporate pertinent soil physico-chemical variables such as soil organic matter, cation exchange capacity, and presence of iron oxide as covariates can provide better interpolation accuracy [40,41,42,43].
There is a need to evaluate efficiency of different geospatial statistical methods such as spatially weighted multivariate regression, kriging, and cokriging to predict the spatial distribution of PTMs in soil. In this regard, the objective of this research is to compare the results of ordinary kriging and spatial modeling informed cokriging, and to assess their prediction accuracy.

2. Materials and Methods

2.1. Study Area

Bumpus Cove is located in the Unaka Range of the Blue Ridge physiographic province in northeast Tennessee and spans both Washington and Unicoi Counties. The cove forms a synclinal valley between northeast-southwest trending Embreeville Mountain (~885 m) and Rich Mountain (~1035 m) [44]. The valley consists of well-jointed Cambrian aged Shady Dolomite and the surrounding ridges are formed by an older Cambrian aged more resistant formation, known as Chilhowee Group sandstones [44]. Due to its regional geology, Bumpus Cove is one of the richest mineralized areas in East Tennessee and contains extensive iron, zinc, lead, and manganese deposits. Zinc, lead, and iron sulfides are present within the Shady Dolomite geologic formation, and oxidized deposits of zinc, lead, iron, and manganese develop as a result of in situ weathering within the residual clay [44]. Bumpus Cove contains 47 abandoned mines which were operational from the late 1700s to the early 1950s. The Peach Orchard Mine was an important supplier of zinc in the United States during its operation from 1916 to 1926 and 1931–1943 [44].
The study area is 0.67 km2 and is located in Unicoi County within the southwest corner of the ~19.5 km2 Bumpus Cove Creek Watershed, which drains northeastward into the Nolichucky watershed (HUC 06010108) (Figure 1). The study area is forested and surrounded by United States Forest Service land to the north, south, and west, and bordered by residential properties of the Embreeville community to the east. The region experiences a humid temperate climate (Köppen climate classification Cfa), typically with wet summers and dry winters. The average annual temperature of Unicoi County is 14 °C (57 °F) and the average annual precipitation is 104 cm (41 in) [45].
The abandoned Peach Orchard Mine is situated in the middle of the study area at the southwest end of the cove in the vicinity of Bumpus Cove Creek. Additionally, five other abandoned Pb, Zn, and Mn mines are centralized in the study area and operated simultaneously during the same time (Figure 1). No additional mapped mines are located upstream of the study area.

2.2. Sampling Procedure

The methods are outlined in Figure 2 and described in the following paragraphs. Soil samples were collected from 10 November to 5 December 2016, using a semi-random sampling strategy, from a 51-cell equal area grid (each cell measured ~130 × ~115 m). After removal of the humus layer, at minimum one soil sample was collected from a random location within each grid cell, from a depth of 0–15 cm (n = 52 samples). Samples were transported to the East Tennessee State University (ETSU) Geosciences Soils Laboratory in a cooler with ice. There, physical parameters (grain size distribution, moisture content, bulk density, and porosity) were measured. Metals (Mn, Zn, Pb, Cu, and Cd) were measured using flame atomic absorption spectrometry (FAAS) at the ETSU Environmental Health Sciences Laboratory. A portion of each soil sample was sent to an external laboratory for additional chemical analyses (total organic content (TOC), cation exchange capacity (CEC), pH). Sampling and lab methods are described in detail in Magno et al., 2019 [22].

2.3. Data Analysis and Kriging

Magno et al. (2019) [22] describe statistical analyses and development of spatial regression models for the metals (Mn, Zn, Pb, Cu, and Cd) using physico-chemical characteristics of the soil samples. This research builds upon that analysis using the findings to select variables for cokriging interpolation models and to select an appropriate neighborhood size. The findings of the prior research are summarized in Table 1 and list the explanatory variables retained in spatial regression models and the neighborhood size of the spatial weights matrix.
For each metal, normality was assessed by examining histograms, and where appropriate, a log transform was applied. Two kriging interpolation models were produced in ArcGIS Pro 2.7.0 [46] and compared: (1) ordinary kriging using the optimize utility; (2) ordinary cokriging using the covariates retained in the spatial regression models in the prior research [22]. The number of neighbors for each cokriging model was calculated by overlaying a circle with radius equal to the optimal neighborhood size identified in the prior study and counting the number of neighbors (nearby sampling points) falling within the circle. The geometric interval classification scheme was selected for each kriging surface and the same classification was used for the sampled data for easy visual comparison.
Five metrics were used to compare the two models for each metal: mean error (ME), root mean square (RMS), mean error standardized (MES), root mean square standardized (RMSS), and average standard error (ASE). MES is the mean of the prediction error divided by the standard error, which should be near zero. RMSS is the standard deviation of the residuals, divided by the standard error, which should be close to one. If RMSS > 1, variability is underestimated and if RMSS < 1, variability is overestimated. ASE is the average of the prediction standard errors. The difference between ASE and RMS was calculated to assess whether variability was overestimated (ASE − RMS > 0) or underestimated (ASE − RMS < 0) in each model. To select a best model, we relied more heavily on the standardized MES and RMSS, and ASE-RMS. All metrics were calculated using cross validation in ArcGIS Pro 2.7.0. Distribution of the predicted and measured concentrations were compared for each model, as was a scatterplot of predicted versus measured concentrations. Maps were compared to assess how well each was able to capture the overall pattern in metal concentration as well as outlier values. The best model selected for each metal was superior in a majority of the metrics and plots, and appropriately captured the spatial pattern of measured data, including outliers.

3. Results

3.1. Descriptive Statistics

Descriptive statistics for all variables were published previously and are included here for context (Table 2) [22]. Soil moisture content fluctuated between 8 and 53%, and generally increased near the east-central part of the study area. Soil bulk density ranged from 1 to 2% and increased eastward. Soil pH ranged from very acidic (3.6) around the northwest and southeast edges, to neutral (7.6) especially towards the middle valley area. The neutral pH values run parallel to the fold axis of the Shady Dolomite formation. The dolomite possibly increases the soil pH. CEC varied in the study area, but the majority ranged from 2 to 8 meq/100 g and generally increased near the central valley. Most soil had TOC ranging from only 3 to 18%, and gradually decreased from northwest to southeast. Overall low CEC and TOC is expected at and around a disturbed, un-reclaimed mine site.
The overall soil texture was well-graded sands (94.5% sand, 0.04% silt, and 0.01% clay) as described by the Unified Soil Classification System (USCS), where higher sand concentration was found along the floodplains of the stream and lowest in the southwestern ridge. Silt content showed the opposite spatial distribution to the sand texture, with silt less predominant along stream floodplain areas and at a maximum in the southwest. Clay showed a decreasing trend from the northeast to southwest.
The Zn, Mn, and Pb concentrations show two orders of magnitude differences between samples (12–1354 mg/kg for Zn, 6–2574 mg/kg for Mn, and 33–2271 mg/kg for Pb). Cu and Cd showed more steady distribution, with ranges of 1–65 and 7–40 mg/kg, respectively. Among the PTMs, only Pb exceeded EPA acceptable limits of 420 mg/kg in soils [47].

3.2. Kriging and Cokriging Model Comparison

For all metals, model output diagnostics (Table 3) were compared to identify the superior model. Optimally, diagnostics for the superior model have ME and MES closest to zero, RMSS closest to one, and the smallest values for RMS and ASE. Results were mixed such that for all metals, some diagnostics suggested kriging was better, while other suggested that the cokriging model was superior. By evaluating diagnostics, cross validation plots of predicted versus measured values, distribution of predicted and measured values, and map products (Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7), we determined that ordinary cokriging produced the better model for Zn, Mn, Cu, and Cd and ordinary kriging produced the better model for Pb.

3.2.1. Manganese

The cokriging model was determined to be superior for Mn, with ME and MES closer to zero. RMS and ASE were slightly smaller for the kriging model, but RMSS was closer to one for the cokriging model. For Mn, the kriging model over estimated variability (ASE > RMS) while the cokriging model underestimated variability (ASE < RMS), but less so (Table 3). The nugget for both variograms was near zero, indicating that the sampling scale was appropriate for the scale of local spatial variability. After additionally reviewing the cross-validation graph and interpolated maps, for Mn, the cokriging model was selected as the better model (Figure 3).

3.2.2. Zinc

The cokriging model was determined to be the better model for zinc because, even though ME and MES were larger in the cokriging model indicating that the model overestimated Zinc concentration on average, the difference was slight (Table 3). Furthermore, RMS and ASE were closer to zero and RMSS was closer to one in the cokriging model (Figure 4). The difference between ASE and RMS was positive in both models, indicating that both overestimated variability. The nugget for the variogram on the cokriging model only was near zero. Examination of the two interpolated maps shows minimal differences between the two interpolation models. The cokriging model was determined to be superior to the kriging model due to the predicted versus measured cross validation plot, which demonstrated a poorer model fit for the kriging model.

3.2.3. Lead

The kriging model was selected as the superior model for lead because the two standardized diagnostics, MES and RMSS were more favorable. While both models overestimated variability (ASE > RMS), the kriging model did so by a smaller amount. The cross-validation plots show an overall underestimation of Pb concentrations, concurring with the positive value for mean error. Only the variogram for the kriging model had a nugget near zero although the nugget of the cokriging model semivariance appears to have been overestimated by the optimization process. Given the mixed results, the final determination of the superior model was decided by weighing the diagnostics with qualitative evaluation of the interpolated maps.

3.2.4. Copper

The cokriging model with log transformation was selected as the better model for copper. The log transformation was used because Cu concentrations were log normally distributed. Both models showed good results from the five cross validation metrics, with the cokriging model having better values for RMS and ASE (Table 3). Further, ASE − RMS was smaller for the cokriging model, indicating it better captured variability in the data. Only the nugget on the variogram for the cokriging model was zero and the variogram for the kriging model was flat, likely due to a high variability/short distance outlier visible in both variograms. We further compared cross-validation scatterplots and distributions of measured and predicted values, which confirmed that the cokriging model was superior (Figure 6).

3.2.5. Cadmium

The cokriging model was selected as the superior model for cadmium because all model diagnostics agreed, with smaller ME, MES, RMS, and ASE, and RMSS closest to one for the cokriging model. Notably, the variograms were very similar, but the nugget for the cokriging model was smallest due to an outlier at short distance evident on the variogram for the kriging model (Figure 7a). Further, the interpolated map produced by cokriging appeared to better capture the influence of drainages in the valley bottom (Figure 7).

4. Discussion

For the purpose of abandoned mine reclamation, it is essential to obtain the most realistic spatial distributions of PTMs. In this regard, the present research used results of spatial regression models published previously [22] to inform cokriging interpolation models. The covariates retained in spatial regression models and the associated optimal spatial neighborhood were used to select cokriging variables and neighborhood size. This was necessary because spatial regression models cannot be used as predictive models due to the auto-regressive spatial lag and spatial error variables in the model. Therefore, we developed interpolation models to estimate metal concentrations at locations where samples were not collected. Indeed, for this research, the value of the spatial regression models was to identify important explanatory variables under conditions of spatial autocorrelation. Once identified, in the present research we sought to improve upon optimized ordinary kriging models by systematically selecting cokriging variables and settings from the prior spatial regression models.
As presented in results, findings were mixed. For four metals (Mn, Zn, Cu, and Cd), the cokriging models were superior to the kriging models, and it appears that the development of spatial regression models was a worthwhile preliminary step. For Pb, this was not the case. Further discussion of model diagnostics presented in Table 3 and the maps and charts in Figure 3, Figure 4, Figure 5, Figure 6 and Figure 7 follow.
We elected to consider six diagnostics (ME, RMS, MES, RMSS, ASE, and ASE − RMS), recognizing that for comparison purposes, the standardized diagnostics (MES and RMSS) and the composite diagnostic ASE-RMS were preferred. MES is preferred over ME [48] with values near zero indicating small errors. MES assesses whether the model overpredicts (MES < 0) or underpredicts (MES > 0) values, in this case metal concentrations. RMSS assesses overall error, with larger errors contributing more to the diagnostic than smaller errors. A small number near one indicates consistently small errors with few or no large outlier errors. ASE-RMS provides information on how well the model captures variability in the data, indicating whether it overestimates (ASE − RMS > 0) or underestimates (ASE − RMS < 0) variability in values [49]. We had concerns regarding the nugget effect in the variograms, because this suggests variability in metal concentrations at a geographic scale larger than the sampling scale, or sampling error. After selection of the best model, only the variogram for Cd had a nonzero nugget, although it was small. Unfortunately, there is no way to predict the variogram a priori, and development of the field plan including sampling grid size was made by emulating sample spacing in other soil studies reviewed from the literature. In future studies, we recommend that soil sampling for heavy metals, and Cd in particular, use a finer sampling resolution to capture local scale variability suggested by the nugget effect.

4.1. Manganese

For Mn, MES and RMSS were very similar, with the cokriging model having a slight edge. The cokriging model was selected because, while the ASE − RMS value of −27.887 indicated that the model underestimated variability, the underestimation was less severe than the overestimation of variability present in the kriging model (ASE − RMS = 39.669). Further, the cokriging model reasonably captured outliers without the obvious bulls-eye effect present in the kriging map (Figure 3a).
The best Mn prediction model used the highest distance threshold of 300 m, i.e., 16–24 neighbors, to estimate meaningful spatial correlation. The larger neighborhood size suggests that Mn was transported from the known mine locations, and primarily concentrated downstream in lower elevation valley areas. The metal transport seems to be affected by drainage and showed better model accuracy when % silt and moisture content were added as covariates. Higher moisture content was found in downstream soil, with low silt content. The previous study [22] in the same study area indicated surface soils exhibited a well-graded sand textural class where most PTMs beside Zn showed a significant positive correlation with % sand, and a significant negative correlation with % silt. Moreover, as anticipated, the turbulent stream in a hilly landscape deposited sands as bed load material from the surrounding slopes of Chilhowee Group sandstones. Mn concentration in lower elevation valleys could also be related to localized ponding, which creates wetland conditions. This finding was supported by the field survey and higher measured moisture content in the valley. In an oxygen-depleted environment, facultative bacteria utilize Mn as terminal electron acceptors to make energy. As a result, Mn is reduced and becomes soluble in soil water. As water mobilizes the soluble Mn, it eventually may encounter oxygenated air through root channels or other macropores, where it oxidizes and redeposits as concentrated Mn [50].

4.2. Zinc

For Zn, the cokriging model was selected as the best model, but there was little difference in the diagnostics of the two models. Likewise, there was little difference between the maps, although the cokriging map better captured the outliers on the southeastern periphery. The variogram for the kriging model (Figure 4a) is flat, with a maximal nugget, while that of the cokriging model (Figure 4b) has a very short range. These suggest that kriging interpolation is not optimal for the Zn data. The best Zn prediction model from the spatial regression modeling used a distance threshold of 150 m, i.e., 6–8 neighbors, to estimate meaningful spatial correlation, which means Zn concentration at the study area is spatially dependent on fewer neighboring data points at a shorter distance.
In the cove, Zn concentration was overall higher than the common 10–300 mg/kg range. Unlike Mn, Zn was concentrated upgradient of the mines, and showed lesser degree of transportation through drainage, which supports the use of the smaller distance threshold (fewer neighbors) in the interpolation. The resulting distribution could be because of upgradient mine tailings and waste rocks or due to localized bioaccumulation. Unmapped mines in the study area may have also impacted Zn concentrations, as the mine locations acquired from the USGS are approximate. The cokriging result indicated that Zn concentration was sensitive to pH and CEC, where Zn was found in higher concentration in soil with near neutral pH and low CEC. The solubility of Zn decreases with an increase in pH. When soil pH exceeds 5.5, Zn becomes less mobile as it is adsorbed in the clay minerals [51] and our findings agreed, whereby a higher concentration of Zn was found in the clay rich soil along the ridges of the western part of the study area. CEC is also known to affect PTM concertation in soil. High CEC has been reported to impact retention of PTMs [52,53,54]. For soils of the Bumpus Cove, CEC was a significant covariate in the cokriging model, however, the result indicated a negative relation between Zn and CEC [22], in contrast to the significant positive correlation between Zn and CEC found at a different site by Navas and Machin [52].

4.3. Lead

For Pb, the diagnostics suggest that the kriging model was superior, and it better captured outlier values albeit at the expense of increased spatial variability (ASE − RMS = 56.086).
Pb was concentrated around the abandoned mine areas in carbonate bedrock and did not show much evidence of transportation. At pH > 6.0, in soils with high calcite, iron and phosphorus content, Pb is not transported and is captured in soil, which might explain limited transportation of Pb. The model results did not improve with addition of covariates like pH and CEC, though other studies have shown that concentration of Pb could be sensitive to soil physico-chemical properties [21,22,52].

4.4. Copper

For Cu, the cokriging model with log transform was selected as the better model with a smaller RMS and ASE, and a low variability (ASE − RMS = −1.290). The cokriging map also better captured high outliers.
Cu was not a minable ore at Bumpass Cove, and likely came from the northwest Chilhowee Group Sandstone bedrock. The general trend of the Cu prediction map showed a decrease from northwest to southeast direction, validating its source. It was interesting to see a northwest–southeast trending linear patch with Cu concentration, which correlated with a mining related disturbance at the specific location. Like Pb, concentration of Cu is also known to correlate with soil physcio-chemical properties like clay content, soil pH, and cation exchange capacity, TOC, and BD [21,22,55]. However, the addition of these covariates did not produce improvement in the model results for distribution of Cu at the cove.

4.5. Cadmium

Lastly, Cd had small errors and RMSS close to one for both models. Cokriging had better diagnostics (smaller MES, RMSS closer to one, and a smaller ASE-RMS) and better captured outliers; it was therefore selected as the best model.
Like Cu, Cd also had no ore sources and was possibly found as trace metal in carbonate bedrock that concentrated primarily downstream, in the lower elevation valley. The Cd prediction map improved when pH, % sand, and % silt were added as covariates. Studies found that Cd concentration in soil is pH dependent; with lower pH, Cd becomes soluble and can be transported, whereas as pH increases, Cd becomes trapped in the soil [56]. Based on the distribution of pH in the soil of Bumpass Cove, we propose that soluble Cd was transported from high elevation low pH areas, and was deposited downstream at the valley region as the soil pH increased. The reason for the presence of a high concentration of Cd in sandy soil is unclear, but at the cove, sandy soil was prevalent in the valley.

5. Conclusions

To obtain the most realistic spatial distributions of PTMs we used optimized ordinary kriging and cokriging models by systematically selecting cokriging variables and neighborhood settings from prior statistical analyses. The PTMs analyzed in this study are Zn, Mn, Pb, Cu, and Cd. While each model successfully produced PTM concentration maps, the findings were mixed in terms of whether kriging or cokriging models were deemed superior for each metal. For Mn, Zn, Cu, and Cd, the cokriging models improved interpolation accuracies, where for Pb the kriging model yielded better results.
The interpolation maps indicated that Mn, Pb and Cd were concentrated around the vicinity of the mine area. Mn and Cd were further transported downstream possibly due to seasonal wetland condition and variable pH. Zn and Pb soil concentration were elevated, with Pb exceeding the USEPA soil regulatory limit.
The study established that soil physico-chemical properties influence the spatial distribution of PTMs in mineralized areas. Moreover, prior to applying spatial interpolation method, development of spatial regression models is an important preliminary step, and should be performed to select the contributing soil properties as this process can improve the final outcome.

Author Contributions

Conceptualization, M.M., I.L. and A.N.; methodology, M.M., I.L. and A.N.; formal analysis, M.M., I.L. and A.N.; investigation, M.M., I.L. and A.N.; resources, A.N.; data curation, M.M.; writing—original draft preparation, M.M.; writing—review and editing, M.M., I.L. and A.N.; visualization, M.M.; supervision, I.L. and A.N.; project administration, M.M.; funding acquisition, M.M., I.L. and A.N. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Association of Environmental and Engineering Geologists (AEG) Stout Scholarship and East Tennessee State University (ETSU) Honors College Office of Undergraduate Research and Creative Activities.

Acknowledgments

The authors would like to thank Brian Evanshen and Phil Scheuerman for use of the East Tennessee State University Department of Environmental Health Sciences laboratory and guidance with laboratory safety and procedures. The authors also acknowledge Rick Turner and the Tennessee Wildlife Resources Agency (TWRA) for permitting access to the study area and for showing interest in this work. Additional thanks are due to Jamie Kincheloe, Mick Whitelaw, Alex McClain, Macon Magno, Samuel Hauser, and Isaac Shockley for their assistance in data collection and processing.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Weiss, F.T.; Leuzinger, M.; Zurbrügg, C.; Eggen, R.I.L. Chemical Pollution in Low and Middle-Income Countries; Eawag: Dubendorf, Switzerland, 2016. [Google Scholar]
  2. Dowdy, R.H.; Volk, V.V.; Nelson, D.W.; Elrick, D.E.; Tanji, K.K. Movement of heavy metals in soils. In Chemical Mobility and Reactivity in Soil Systems; American Society of Agronomy: Madison, MI, USA, 1983; pp. 227–240. [Google Scholar]
  3. Reeder, R.J.; Schoonen, M.A.A.; Lanzirotti, A. Metal Speciation and Its Role in Bioaccessibility and Bioavailability. Rev. Miner. Geochem. 2006, 64, 59–113. [Google Scholar] [CrossRef]
  4. Pierzynski, G.M.; Vance, G.F.; Sims, J.T. Soils and Environmental Quality; CRC Press: Boca Raton, FL, USA, 2005. [Google Scholar]
  5. Jain, C. Metal fractionation study on bed sediments of River Yamuna, India. Water Res. 2004, 38, 569–578. [Google Scholar] [CrossRef]
  6. Förstner, U.; Wittmann, G.T. Metal Pollution in the Aquatic Environment, 2nd ed.; Springer Science and Business Media: Berlin/Heidelberg, Germany, 2012. [Google Scholar]
  7. Alloway, B.J. Heavy Metals in Soils; Springer Science and Business Media: Berlin/Heidelberg, Germany, 1995. [Google Scholar]
  8. He, Z.L.; Yang, X.E.; Stoffella, P.J. Trace elements in agroecosystems and impacts on the environment. J. Trace Elem. Med. Biol. 2005, 19, 125–140. [Google Scholar] [CrossRef]
  9. Vega, F.; Covelo, E.; Andrade, M. Competitive sorption and desorption of heavy metals in mine soils: Influence of mine soil characteristics. J. Colloid Interface Sci. 2006, 298, 582–592. [Google Scholar] [CrossRef]
  10. Li, Z.; Schuman, L.M. Heavy metal movement in metal-contaminated soil profiles. Soil Sci. 1996, 161, 656–666. [Google Scholar] [CrossRef]
  11. Chang, A.C.; Warneke, J.E.; Page, A.L.; Lund, L.J. Accumulation of Heavy Metals in Sewage Sludge-Treated Soils. J. Environ. Qual. 1984, 13, 87–91. [Google Scholar] [CrossRef]
  12. Ersoy, A.; Yunsel, T.Y.; Cetin, M. Characterization of Land Contaminated by Past Heavy Metal Mining Using Geostatistical Methods. Arch. Environ. Contam. Toxicol. 2004, 46, 162–175. [Google Scholar] [CrossRef]
  13. Alloway, B.; Davies, B. Trace element content of soils affected by base metal mining in Wales. Geoderma 1971, 5, 197–208. [Google Scholar] [CrossRef]
  14. Kabala, C.; Singh, B.R. Fractionation and Mobility of Copper, Lead, and Zinc in Soil Profiles in the Vicinity of a Copper Smelter. J. Environ. Qual. 2001, 30, 485–492. [Google Scholar] [CrossRef] [Green Version]
  15. Qin, F.; Ji, H.; Li, Q.; Guo, X.; Tang, L.; Feng, J. Evaluation of trace elements and identification of pollution sources in particle size fractions of soil from iron ore areas along the Chao River. J. Geochem. Explor. 2014, 138, 33–49. [Google Scholar] [CrossRef]
  16. Tessier, A.; Campbell, P.G.C.; Bisson, M. Sequential extraction procedure for the speciation of particulate trace metals. Anal. Chem. 1979, 51, 844–851. [Google Scholar] [CrossRef]
  17. Hernandez, L.; Probst, A.; Probst, J.L.; Ulrich, E. Heavy metal distribution in some French forest soils: Evidence for atmospheric contamination. Sci. Total Environ. 2003, 312, 195–219. [Google Scholar] [CrossRef] [Green Version]
  18. Yao, Q.; Wang, X.; Jian, H.; Chen, H.; Yu, Z. Characterization of the Particle Size Fraction associated with Heavy Metals in Suspended Sediments of the Yellow River. Int. J. Environ. Res. Public Health 2015, 12, 6725–6744. [Google Scholar] [CrossRef]
  19. Cabral, A.R.; Lefebvre, G. Use of Sequential Extraction in the Study of Heavy Metal Retention by Silty Soils. Water Air Soil Pollut. 1998, 102, 329–344. [Google Scholar] [CrossRef]
  20. Rodríguez, L.; Ruiz, E.; Alonso-Azcárate, J.; Rincón, J. Heavy metal distribution and chemical speciation in tailings and soils around a Pb–Zn mine in Spain. J. Environ. Manag. 2009, 90, 1106–1116. [Google Scholar] [CrossRef]
  21. Santos-Francés, F.; Martínez-Graña, A.; Zarza, C.; Sánchez, A.G.; Rojo, P.A. Spatial Distribution of Heavy Metals and the Environmental Quality of Soil in the Northern Plateau of Spain by Geostatistical Methods. Int. J. Environ. Res. Public Health 2017, 14, 568. [Google Scholar] [CrossRef]
  22. Magno, M.A.; Nandi, A.; Luffman, I.E. Using Spatial Regression to Model Potentially Toxic Metal (PTM) Mobility Based on Physicochemical Soil Properties. Appl. Environ. Soil Sci. 2019, 2019, 6432571. [Google Scholar] [CrossRef] [Green Version]
  23. Gu, G.-Q.; Wan, X.-M.; Zeng, W.-B.; Lei, M. Analysis of the Spatial Distribution of Heavy Metals in Soil from a Coking Plant and Its Driving Factors. Huanjing Kexue 2021, 42, 1081–1092. [Google Scholar]
  24. Mirčovski, V.; Gičevski, B.; Dimov, G. Hydrochemical characteristics of the groundwaters in Prilep’s part of Pelagonia valley–Republic of Macedonia. Rud. Geol. Naft. Zb. 2018, 33, 111–119. [Google Scholar] [CrossRef] [Green Version]
  25. Šapina, M. A comparison of artificial neural networks and ordinary kriging depth maps of the lower and upper Pannonian stage border in the Bjelovar Subdepression, Northern Croatia. Rud. Geol. Naft. Zb. 2016, 31, 75–85. [Google Scholar] [CrossRef] [Green Version]
  26. Kiš, I.M. Comparison of ordinary and universal kriging interpolation techniques on a depth variable (a case of linear spatial trend), case study of the Šandrovac Field. Rud. Geol. Naft. Zb. 2016, 31, 41–58. [Google Scholar] [CrossRef]
  27. Ivšinović, J.; Malvić, T. Application of the Radial Basis Function interpolation method in selected reservoirs of the Croatian part of the Pannonian Basin System. Min. Miner. Depos. 2020, 14, 37–42. [Google Scholar] [CrossRef]
  28. Huo, X.-N.; Li, H.; Sun, D.-F.; Zhou, L.-D.; Li, B.-G. Combining Geostatistics with Moran’s I Analysis for Mapping Soil Heavy Metals in Beijing, China. Int. J. Environ. Res. Public Health 2012, 9, 995–1017. [Google Scholar] [CrossRef] [Green Version]
  29. Amini, M.; Afyuni, M.; Fathianpour, N.; Khademi, H.; Flühler, H. Continuous soil pollution mapping using fuzzy logic and spatial interpolation. Geoderma 2005, 124, 223–233. [Google Scholar] [CrossRef]
  30. Davis, H.T.; Aelion, C.M.; McDermott, S.; Lawson, A.B. Identifying natural and anthropogenic sources of metals in urban and rural soils using GIS-based data, PCA, and spatial interpolation. Environ. Pollut. 2009, 157, 2378–2385. [Google Scholar] [CrossRef] [Green Version]
  31. Ding, Q.; Cheng, G.; Wang, Y.; Zhuang, D. Effects of natural factors on the spatial distribution of heavy metals in soils surrounding mining regions. Sci. Total Environ. 2017, 578, 577–585. [Google Scholar] [CrossRef]
  32. Qiao, P.; Lei, M.; Yang, S.; Yang, J.; Guo, G.; Zhou, X. Comparing ordinary kriging and inverse distance weighting for soil as pollution in Beijing. Environ. Sci. Pollut. Res. 2018, 25, 15597–15608. [Google Scholar] [CrossRef]
  33. Mahmoudabadi, E.; Sarmadian, F.; Savaghebi, G.H.; Alijani, Z. Accuracy assessment of geostatistical methods for zoning of heavy metals in soils of urban-industrial areas. Int. Res. J. Appl. Basic Sci. 2012, 3, 991–999. [Google Scholar]
  34. Kravchenko, A.N. Influence of Spatial Structure on Accuracy of Interpolation Methods. Soil Sci. Soc. Am. J. 2003, 67, 1564–1571. [Google Scholar] [CrossRef]
  35. Li, J.; Heap, A.D. A review of comparative studies of spatial interpolation methods in environmental sciences: Performance and impact factors. Ecol. Inform. 2011, 6, 228–241. [Google Scholar] [CrossRef]
  36. Oliver, M.A.; Webster, R. Kriging: A method of interpolation for geographical information systems. Int. J. Geogr. Inf. Syst. 1990, 4, 313–332. [Google Scholar] [CrossRef]
  37. Wang, Y.; Duan, X.; Wang, L. Spatial distribution and source analysis of heavy metals in soils influenced by industrial enterprise distribution: Case study in Jiangsu Province. Sci. Total Environ. 2019, 710, 134953. [Google Scholar] [CrossRef]
  38. Reza, S.K.; Baruah, U.; Singh, S.K.; Das, T.H. Geostatistical and multivariate analysis of soil heavy metal contamination near coal mining area, Northeastern India. Environ. Earth Sci. 2014, 73, 5425–5433. [Google Scholar] [CrossRef]
  39. Li, X.; Yang, H.; Zhang, C.; Zeng, G.; Liu, Y.; Xu, W.; Wu, Y.; Lan, S. Spatial distribution and transport characteristics of heavy metals around an antimony mine area in central China. Chemosphere 2016, 170, 17–24. [Google Scholar] [CrossRef]
  40. Chen, T.; Liu, X.; Li, X.; Zhao, K.; Zhang, J.; Xu, J.; Shi, J.; Dahlgren, R. Heavy metal sources identification and sampling uncertainty analysis in a field-scale vegetable soil of Hangzhou, China. Environ. Pollut. 2008, 157, 1003–1010. [Google Scholar] [CrossRef]
  41. Chen, T.; Chang, Q.; Liu, J.; Clevers, J.; Kooistra, L. Identification of soil heavy metal sources and improvement in spatial mapping based on soil spectral information: A case study in northwest China. Sci. Total Environ. 2016, 565, 155–164. [Google Scholar] [CrossRef]
  42. Smichowski, P.; Gómez, D.; Frazzoli, C.; Caroli, S. Traffic-Related Elements in Airborne Particulate Matter. Appl. Spectrosc. Rev. 2007, 43, 23–49. [Google Scholar] [CrossRef]
  43. Liu, D.; Wang, Z.; Zhang, B.; Song, K.; Li, X.; Li, J.; Li, F.; Duan, H. Spatial distribution of soil organic carbon and analysis of related factors in croplands of the black soil region, Northeast China. Agric. Ecosyst. Environ. 2006, 113, 73–81. [Google Scholar] [CrossRef]
  44. Rodgers, J. Geological and Mineral Deposits of Bumpass Cove, Unicoi and Washington Counties; Bulletin 54; Division of Geology: Nashville, TN, USA, 1948. [Google Scholar]
  45. National Weather Service National Weather Service Climate. Available online: http://w2.weather/gov/climate/xmacis/php?wfo=mrx (accessed on 10 June 2016).
  46. ArcGIS Pro 2.7.0 2020; Environmental Systems Research Institute Inc.: West Redlands, CA, USA, 2020.
  47. United States Department of Agriculture. Natural Resources Conservation Service Heavy Metal Soil Contamination. In Soil Quality—Urban Technical Note No. 3; United States Department of Agriculture: Auburn, AL, USA, 2000. [Google Scholar]
  48. Oliver, M.; Webster, R. A tutorial guide to geostatistics: Computing and modelling variograms and kriging. Catena 2014, 113, 56–69. [Google Scholar] [CrossRef]
  49. Interstate Technology Regulatory Council. Geospatial analysis for optimization at Environmental Sites (GRO-1). In Evaluate Geospatial Method Accuracy; ITRC: Washington, DC, USA, 2016. [Google Scholar]
  50. Vepraskas, M.J.; Sprecher, S.W. Overview of aquic conditions in hydric soils. In Aquic Conditions and Hydric Soils: The Problem Soils; Vepraskas, M.J., Sprecher, S.W., Eds.; Soil Science Society of America: Madison, WI, USA, 1997; pp. 1–22. [Google Scholar]
  51. Fageria, N.K.; Nascente, A.S. Management of Soil Acidity of South American Soils for Sustainable Crop Production. Adv. Agron. 2014, 128, 221–275. [Google Scholar] [CrossRef]
  52. Navas, A.; Machín, J. Spatial distribution of heavy metals and arsenic in soils of Aragón (northeast Spain): Controlling factors and environmental implications. Appl. Geochem. 2002, 17, 961–973. [Google Scholar] [CrossRef] [Green Version]
  53. Orhue, E.R.; Frank, U.O. Fate of some heavy metals in soils: A review. J. Appl. Nat. Sci. 2011, 3, 131–138. [Google Scholar] [CrossRef] [Green Version]
  54. Ghosh, M.; Singh, S.P. A review of phytoremediation of heavy metals and utilization of its by products. Asian J. Energy Environ. 2005, 6, 214–231. [Google Scholar]
  55. Imperato, M.; Adamo, P.; Naimo, D.; Arienzo, M.; Stanzione, D.; Violante, P. Spatial distribution of heavy metals in urban soils of Naples city (Italy). Environ. Pollut. 2003, 124, 247–256. [Google Scholar] [CrossRef]
  56. Yun, S.-W.; Yu, C. Immobilization of Cd, Zn, and Pb from Soil Treated by Limestone with Variation of pH Using a Column Test. J. Chem. 2015, 2015, 641415. [Google Scholar] [CrossRef] [Green Version]
Figure 1. Bumpus Cove study area (0.67 km2) in northeast TN, USA, showing topographic contours and mine and sample locations.
Figure 1. Bumpus Cove study area (0.67 km2) in northeast TN, USA, showing topographic contours and mine and sample locations.
Geosciences 11 00434 g001
Figure 2. Description of methods.
Figure 2. Description of methods.
Geosciences 11 00434 g002
Figure 3. Variograms and interpolated Mn concentration using ordinary kriging (a) and cokriging (b). The cokriging model, marked with an asterisk, was selected as the superior model.
Figure 3. Variograms and interpolated Mn concentration using ordinary kriging (a) and cokriging (b). The cokriging model, marked with an asterisk, was selected as the superior model.
Geosciences 11 00434 g003
Figure 4. Variograms and interpolated Zn concentration using ordinary kriging (a) and cokriging (b). Cokriging model, marked with an asterisk, was selected as the superior model.
Figure 4. Variograms and interpolated Zn concentration using ordinary kriging (a) and cokriging (b). Cokriging model, marked with an asterisk, was selected as the superior model.
Geosciences 11 00434 g004
Figure 5. Variograms and interpolated Pb concentration using ordinary kriging (a) and cokriging (b). Kriging model, marked with an asterisk, was selected as the superior model.
Figure 5. Variograms and interpolated Pb concentration using ordinary kriging (a) and cokriging (b). Kriging model, marked with an asterisk, was selected as the superior model.
Geosciences 11 00434 g005
Figure 6. Variograms and interpolated Cu concentration using ordinary kriging (a) and cokriging (b). Cokriging model, marked with an asterisk, was selected as the superior model.
Figure 6. Variograms and interpolated Cu concentration using ordinary kriging (a) and cokriging (b). Cokriging model, marked with an asterisk, was selected as the superior model.
Geosciences 11 00434 g006
Figure 7. Variograms and interpolated Cd concentration using ordinary kriging (a) and cokriging (b). Cokriging model, marked with an asterisk, was selected as the superior model.
Figure 7. Variograms and interpolated Cd concentration using ordinary kriging (a) and cokriging (b). Cokriging model, marked with an asterisk, was selected as the superior model.
Geosciences 11 00434 g007
Table 1. Summary of spatial regression model output used to inform choice of cokriging variables and settings [22].
Table 1. Summary of spatial regression model output used to inform choice of cokriging variables and settings [22].
MetalCovariatesNeighborhood SizeNumber of Neighbors
Manganese (Mn)Moisture content, silt300 m16–24
Zinc (Zn)pH, cation exchange capacity150 m6–8
Lead (Pb)pH, cation exchange capacity150 m6–8
Copper (Cu)Sand, silt, bulk density200 m8–12
Cadmium (Cd)pH, sand, silt150 m6–8
Table 2. Descriptive statistics of soil properties [22].
Table 2. Descriptive statistics of soil properties [22].
VariablesMinMaxMeanStandard DeviationSkewnessKurtosis
Soil Properties
moisture content (%)8.2975.2124.5913.21.43.1
bulk density (g/cm3)0.791.791.250.20.2−0.07
cation exchange capacity (meq/100 g)1.5216.584.833.51.82.9
total organic carbon (%)3.231.58.64.72.610.6
sand (%)8899950.03−0.5−0.8
silt (%)11140.030.6−0.8
clay (%)0310.010.3−0.2
Metal Concentration
Mn (mg/kg)6.292574.93344.49652.22.44.9
Zn (mg/kg)11.801354.16302.52402.71.71.6
Pb (mg/kg)33.432271.43326.69529.62.76.4
Cu (mg/kg)1.1464.6713.9613.22.35.9
Cd (mg/kg)7.1440.0011865.13.617.8
Table 3. Ordinary kriging (OK) and ordinary cokriging (OCK) diagnostics from interpolation models for five metals. The better model is in bold font.
Table 3. Ordinary kriging (OK) and ordinary cokriging (OCK) diagnostics from interpolation models for five metals. The better model is in bold font.
MnZnPbCuCd
ModelOKOCKOKOCKOKOCKOKOCKOKOCK
Cokriging covariates-MC, silt-pH, CEC-pH, CEC-sand, silt, BD-pH, sand, silt
Neighbors 16–24 6–8 6–8 8–12 6–8
ME−29.315−19.71412.81414.31623.28020.7340.086−0.1930.2370.0350
RMS611.194611.253409.981378.532435.264405.22412.79812.1665.0053.954
MES−0.0287−0.02120.02920.03580.03660.0412−0.0017−0.1090.03980.00857
RMSS1.0221.0210.9670.9660.8760.8240.8341.2130.9120.947
ASE571.525583.366424.396392.328491.351494.57414.42010.8755.5524.187
ASE−RMS39.669−27.88714.41713.79756.08689.3501.621−1.2900.54650.2333
Covariates: MC = moisture content; CEC = cation exchange capacity, BD = bulk density. Diagnostics: ME = mean error; RMS = root mean square; MES = mean error standardized; RMSS = root mean square standardized; ASE = average standard error.
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Magno, M.; Luffman, I.; Nandi, A. Evaluating Spatial Regression-Informed Cokriging of Metals in Soils near Abandoned Mines in Bumpus Cove, Tennessee, USA. Geosciences 2021, 11, 434. https://0-doi-org.brum.beds.ac.uk/10.3390/geosciences11110434

AMA Style

Magno M, Luffman I, Nandi A. Evaluating Spatial Regression-Informed Cokriging of Metals in Soils near Abandoned Mines in Bumpus Cove, Tennessee, USA. Geosciences. 2021; 11(11):434. https://0-doi-org.brum.beds.ac.uk/10.3390/geosciences11110434

Chicago/Turabian Style

Magno, Melissa, Ingrid Luffman, and Arpita Nandi. 2021. "Evaluating Spatial Regression-Informed Cokriging of Metals in Soils near Abandoned Mines in Bumpus Cove, Tennessee, USA" Geosciences 11, no. 11: 434. https://0-doi-org.brum.beds.ac.uk/10.3390/geosciences11110434

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop