Correlation between Geochemical and Multispectral Patterns in an Area Severely Contaminated by Former Hg-As Mining

Abstract

In the context of soil pollution, plants suffer stress when exposed to extreme concentrations of potentially toxic elements (PTEs). The alterations to the plants caused by such stressors can be monitored by multispectral imagery in the form of vegetation indices, which can inform pollution management strategies. Here we combined geochemistry and remote sensing techniques to offer a preliminary soil pollution assessment of a vast abandoned spoil heap in the surroundings of La Soterraña mining site (Asturias, Spain). To study the soil distribution of the PTEs over time, twenty-seven soil samples were randomly collected downstream of and around the main spoil heap. Furthermore, the area was covered by an unmanned aerial vehicle (UAV) carrying a high-resolution multispectral camera with four bands (red, green, red-edge and near infrared). Multielement analysis revealed mercury and arsenic as principal pollutants. Two indices (from a database containing up to 55 indices) offered a proper correlation with the concentration of PTEs. These were: CARI2, presenting a Pearson Coefficient (PC) of 0.89 for concentrations >200 mg/kg of As; and NDVIg, PC of −0.67 for >40 mg/kg of Hg. The combined approach helps prediction of those areas susceptible to greatest pollution, thus reducing the costs of geochemical campaigns.

1. Introduction

Mining and industrial activities commonly release contaminants, which in turn affect the quality of environmental compartments [1,2]. In this regard, pollution has been reported not only in urban zones [3,4] but also in more rural or natural areas, where the quality of sensitive ecosystems is being seriously compromised [5]. Although industries are not usually located in rural areas, this is not true for mining operations as the ores of interest are commonly found in these areas. In the case of mining activities in hilly regions and subsequent waste disposal, residues are exposed to hard weathering and a challenging geomorphology. Therefore, pollutants can be widely dispersed, and the costs and uncertainties of geochemical campaigns may be very high [6].

Given these considerations and in the context of soil pollution studies, there is a need to find appropriate non-invasive approaches that reduce the cost and time required for these campaigns. Remote sensing techniques are widely used in natural environment research, such as studies on forest studies [7,8,9], relative climate change [10,11,12], topography [13,14,15,16,17], and precision agriculture [18,19], among others. Of note, in recent years, some studies have combined techniques of remote sensing and geochemistry, especially in areas affected by heavy metal pollution, either with satellites or high-altitude manned aircrafts [20,21,22,23] or unmanned aerial vehicles (UAVs) [24,25,26,27,28]. Remote sensing uses electromagnetic radiation to obtain information about an object or phenomenon without involving physical contact [29]. Data are gathered through measurements of electromagnetic reflectance and different frequencies [30]. Within this discipline, UAVs are one of the supports used for the development of remote sensing tools for soil applications, with most of the required information being collected by means of optical sensors, which have also been recently used in geological applications [31,32].

To date, studies addressing the electromagnetic spectrum have been based on the use of proximal or airborne sensors. However, several authors have recently used proximal aerial sensors, also known as UAVs, for this purpose, due to the advantages they present over satellite, and high-altitude manned aircrafts, such as short revist time and high spatial resolution [33]. This approach has also found applications in soil contamination research [33,34]. The direct study of contaminants in soil through remote sensing is complex; therefore, indirect techniques are used, among which vegetation as a bioindicator of soil condition stands out [30]. One of the most common ways to analyze vegetation with multispectral remote sensors is through vegetation indices, which are simple algebraic operations between the different bands of the electromagnetic spectrum [35,36,37]. In this regard, Shabou et al. [38] used the NDVI and MIR indices over test fields to achieve bare soil measurements from satellite images.

With respect to target pollutants, metal mining is a frequent source of potentially toxic elements (PTEs). These elements are characterized not only by their toxicity, but also by their capacity to bioaccumulate, biomagnify and persist in the food chain [39], giving rise to multiple health and environmental problems [40]. Among PTEs, mercury (Hg) is highly toxic to organisms, even at very low concentrations, especially when present in methyl-Hg form [41]. In this context, in 2005, the EU promoted an initiative to reduce the emissions of this element. Later, in 2017, the Minamata Convention, a global treaty to regulate the emissions and distribution of this toxic element, was approved. This treaty encourages countries to reduce and monitor the different sources of Hg pollution [42]. Since then, various studies have addressed the environmental impact of Hg mining [43,44,45].

Historically, Spain has been one of the world’s largest producers of Hg, especially the Almadén district [46], the largest Hg deposit in the world. In the case of Asturias (NW Spain), numerous cinnabar deposits have been mined, the remains of which are now scattered throughout the region mainly as spoil heaps enriched in Hg-As waste. The impact of some of these mining operations has been examined in detail. As shown in the studies carried out at the mines of El Terronal [47] and Caunedo [43], and at the Muñón-Cimero complex (Soterraña), which is the focus of the present study [48,49]. Here we describe an innovative approach that combines a classical geochemical campaign and a UAV-based multispectral airborne camera to monitor the impact of Hg-As spoil heaps. The main objective of this study is a first approach to the possibility of evaluating As or Hg levels in soil polluted areas through low-cost multispectral sensors.

2. Materials and Methods

In order to summarize and synthetize the procedures followed, the different steps carried out along the study (data collection, geochemistry and remote sensing techniques, statistical analysis, etc.) are shown in Figure 1.

2.1. Context of the Study Area: Activity in La Soterraña, Geology and Vegetation

The main activity of the Hg mountain mine of La Soterraña (Asturias, Spain. Figure 2) took place between the 1940s and the mid-1970s, the latter decade marking the crisis of the Hg sector. During its time of operation, ore mining and processing were carried out onsite. The room and pillar mining method was used on a low temperature hydrothermal epigenetic type deposit, in which the main mineral paragenesis comprised cinnabar, realgar and oropiment, with important deposits of arsenopyrite, marcasite and pyrite hosted in fractured limestones and shales [48]. Once the materials had been extracted, the ore was treated at La Soterraña mine itself, together with the material of other smaller mines located in the province. For this purpose, rotary kilns, condensers and coolers were used, as well as other equipment typical of Hg metallurgy.

The tailings generated by ore mining and processing were indiscriminately dumped in a spoil heap, which is the principal focus of pollution. At present, the mechanical dispersion of this material is the main cause of pollutant distribution in the surrounding soils, in addition to the processes of oxidation and leaching. This main spoil heap is located over layers of limestone and shale, along an average slope of 30°. There is virtually no vegetation on the heap due to the high concentrations of metals and its steep slope. Downstream of the heap, herbaceous colonies of Centaurea nigra and Equisetum telmateia predominate. These species are of special interest due to their capacity to bioaccumulate Hg and As [49,50]. In fact, E. telmateia and C. nigra showed notable Bioaccumulation Indices (BAC); BAC_As = 14.39 and BAC_Hg = 76.95, and BAC_As = 2.30 and BAC_Hg = 4.84 respectively [51]. In the surroundings of the spoil heap, there are also notable tree species, including Betula celtibérica, Quercus robur and Castanea sativa.

2.2. UAV Flight, Sampling, Preparation and Chemical Determinations

The first step consisted of designing the coverage of UAV flights. For the purpose of this study, we selected a surface area of 25 ha downstream of the spoil heap (Figure 1), which is the area most affected by erosion and washing by surface runoff.

The UAV flight was carried out with a BlueGrass Multiespectal model UAV (Parrot SA, Paris, France) which incorporates a Parrot Sequoia plus camera. This camera includes a multispectral sensor, with four spectral bands: green, red, red edge and near infrared, and a sunlight sensor, which provides absolute reflectance measurements without the need of calibration. The use of this camera with the software PIX4D Mapper (Pix4D S.A. Prilly, Switzerland) allowed to obtain the optimal results, as reflected by other authors [52,53,54]. With this procedure, a reduction in the time of the photogrammetric survey and data processing is achieved. The photogrammetric survey is carried out with a transversal and longitudinal reception of 70%, due to the strong orography, and with an average height of 90 m. The date of the flight was chosen based on the vegetation conditions during the early spring of 2019 in order to obtain radiometric information of the herbaceous vegetation when younger; in addition, in this period of the year soil was not affected by extreme weather conditions such as droughts, or alterations due to agro-livestock activities such as mowing or cattle grazing.

In turn, a geochemical campaign composed of a total of 27 soil samples was designed. The samples were gathered randomly. Ten of these samples were collected within the limits of the area covered by the UAV flight (Figure 2), while the others were gathered outside this area in order to assess the degree of pollution in the surroundings for future research purposes. Each sample comprised five increments of the first 20 cm of the soil by means of a Dutch Edelman probe, forming a 1-m edge square and its central point to improve representativeness.

Sample preparation started with in-situ sieving of the soil through 20-mm mesh to remove large particles. The finer fraction was then placed in sterilized plastic bags and stored at 4 °C until being taken to the laboratory. Here, the fractions were air-dried at room temperature to prevent the evaporation of Hg and quartered to maintain sample representativeness. Finally, for each sample, 1 g was ground in an agate mortar to a particle size of < 100 μm and sent to an ISO 9002 and ISO-17025 accredited laboratory in Vancouver (BC, Canada). The concentrations of the following PTEs were determined after aqua regia extraction and ICP-MS (in parenthesis, their respective detection limits): As (0.5 ppm); Cd (0.1 ppm); Co (0.1 ppm); Cr (1 ppm); Cu (0.2 ppm); Hg (0.01 ppm); Mn (1 ppm); Mo (0.1 ppm); Ni (0.1 ppm); Pb (0.1 ppm); Sb (0.1 ppm); Tl (0.1 ppm); V (2 ppm); and Zn (1 ppm).

2.3. Statistical Analysis and Graphical Representations

The geochemical data were analyzed by means of descriptive statistics and compared with the official Risk-Based Soil Screening Levels (RBSSLs) [55]. Moreover, Pearson’s Coefficient between elements was also calculated. These operations allowed us to identify the main pollutants in the study area. Additionally, the concentration of As and Hg for the whole area was mapped by means of ordinary kriging, applying log-transformation to improve normality [56]. Apart from providing an estimation of the distribution of elements in the soil, the interpolations are also useful to establish graphical comparisons between distribution and the maps showing the vegetation indices. All map representations and graphical information were treated through Geographical Information Systems (GIS), using a licensed copy of ArcMap version 10.2.1. (ESRI, Redlands, CA, USA).

2.4. Multi-Spectral Analysis Tools: Vegetation Index

The multispectral images used were obtained from a Sequoia + multispectral camera mounted on a UAV Parrot Bluegrass Fields. The camera had a spectral resolution of four bands—green, (550 nm ± 40 nm), red (660 nm ± 40 nm), red edge (735 nm ± 10 nm) and near infrared (790 ± 40 nm)—all with a high spatial resolution of 16 cm.

The first step of multispectral analysis was the achievement of a raster dataset from each band, which was carried out with software PIX4D Mapper using all photograms per each spectral resolution. The second step consisted in obtaining a new raster dataset with information derived from an algebraic operation between the bands obtained in first step. These new raster datasets are known as vegetation indices, due to their widespread use for determining vegetation status. Up to a total of 55 different vegetation indices were tested with the aim to obtain the best correlations.

The list of the 55 vegetation indices used was extracted from a remote sensing index database [57]. Initially all of them were studied, adapted to the available data and tested. Those found with the greatest correlation with the available analytical data were the Chlorophyll Absorption Ratio Index 2 (CARI-2), and the Normalized Difference Vegetation Index-green (NDVIg).

The CARI2 index was defined by Kim [58] to study leaf area index (LAI), absorbed photosynthetically active radiation, and primary production. It is described by Equation (1):

$$CARI-2=\frac{\left|\left(a·\left|Red\right|+\left|Red\right|+b\right)\right|}{\sqrt{a^2+1}}·\frac{Red Edge}{Red}$$

(1)

where:

$$\begin{gathered} a=\left(\left[Red Edge\right]-\left[Green\right]\right)/150 \\ b=\left[Green\right]-\left(a ·\left[Green\right]\right) \end{gathered}$$

In contrast, NDVIg is a modification of the classical NDVI index that uses absorbance of green light (550 nm), which is more sensitive to a wide range of chlorophyll-a concentrations than red light [56]. It is defined using Equation (2):

$$NDV{I}_g=\frac{Near Infrared-Green}{Near Infrared+Green}$$

(2)

3. Results and Discussion

3.1. Geochemical Characterization

A previous study by Fernández et al. [59] reported maximum total concentrations of 23,000 mg/kg of As and 6000 mg/kg of Hg, corresponding to the mining waste of the spoil heap. As a consequence of spoil heap weathering, the soils of the area showed extremely high concentrations of both elements, as reflected in the descriptive statistics (

Table 1. Descriptive statistics for the 27 soil samples analyzed. Concentrations expressed as mg/kg.

PTE	RBSSL	Range	Mean	Median	Typical Deviation	CV
As	40	44–9920	632.1	150.0	1876.4	296.8
Cd	2	0.2–0.9	0.4	0.5	0.2	40.3
Co	25	2.7–17	13.4	13.4	3.0	22.1
Cr	10,000	18–57	26.9	25.0	8.0	29.6
Cu	55	11–98	33.2	28.8	18.8	56.6
Hg	1	1.95–860	68.0	29.2	158.0	232.4
Mn	2135	251–1270	707.2	704.0	248.3	35.1
Mo	6	0.4–3.6	1.1	1.0	0.6	61.6
Nor	65	10–42	28.9	28.0	6.5	22.5
Pb	70	24–74	37.5	34.0	11.2	29.8
Sb	5	0.2–44	3.4	2.0	8.1	238.0
Tl	1	0.1–2	1.2	2.0	0.9	74.1
V	50	20–83	38.9	35.0	12.6	32.4
Zn	455	65–214	121.3	113.0	29.7	24.5

) corresponding to the chemical analyses of the 27 soil samples taken downstream of the heap. Table 1 also includes the official Risk-Based Soil Screening Levels (RBSSLs) or threshold values in force for the use of soil [55].

The differences between the mean and the median, or the high coefficients of variation (CV), strongly suggest Hg and As as the main pollutants in the area as they present irregular distributions. All the samples exceeded the RBSSL for Hg (1 mg/kg) and As (40 mg/kg) [46], values that correspond to the category “Other Uses”, to which the soil of study currently belongs because it is located in agricultural areas dedicated to pasture or crop cultivation.

Furthermore, most of the samples were also above the soil-chemical background levels of the area determined in pristine areas by Loredo et al. (4.18 mg/kg for Hg and 39 mg/kg for As) [48]. Although the enrichment of these elements was consistent with the mineral paragenesis of the zone described in Section 2.1, the extremely high concentrations detected revealed a gradual degradation of the spoil heap over time due to weathering processes. These are attributable mainly to water erosion by surface and subsurface runoff, which circulates throughout the heap during the months of flooding (Figure 2 and Figure 3). In this regard, higher concentrations of Hg and As were observed in areas close to the runoff, thereby contributing to the dispersion of these pollutants downstream.

With respect to the rest of the PTEs, Sb also appeared in an anomalous form, although to a lesser extent than Hg and As, as revealed by its range or high CV, although the similarity between average and median reveals more normality of the Sb distribution. In this regard, Sb appeared to be linked to As (metalloids usually geochemically associated) [60]. The rest of the PTEs presented stable values and were not hazardous for the environment in this case. The highest concentrations of As and Hg were recorded within the UAV flight zone and in the two samples to the north, near the mining area (Figure 3). Outside the area covered by the UAV, although the concentrations of both contaminants were substantially lower for most samples, they continued to be high. Given this observation, it would be advisable for future studies to collect more samples outside this area in order to better delimit the extent of contamination. Finally, bivariate correlations revealed high correlations between Hg and Cr (0.777), V (0.736) and Pb (0.697). In contrast, As showed an almost full correlation only with Sb (0.991). Of note, there was a low correlation between As-Hg (0.270), reflecting independence between variables. This observation can be explained by the highly irregular distribution of As and Hg and their very different behaviors.

3.2. Mutispectral Analysis

The results presented in Section 3.1 revealed the severe contamination present within the UAV flight area, the main contaminants being Hg and As. Therefore, it is feasible that such high concentrations affect the physiology of plants in the zone. Of the 55 indices tested, those with the greatest correlation with the available analytical data were CARI-2 and NDVIg for As and Hg, respectively. These indices are shown in Figure 4.

The correlations between the concentrations of Hg and As obtained by chemical analysis and the vegetation indices were established by means of Pearson’s coefficient and are shown in Figure 5. The results obtained showed that the higher the concentration, the higher the correlation, in absolute value. In the case of As, it was revealed a significant correlation starting with Pearson value’s of 0.732 (corresponding to concentrations above 117 ppm) up to 0.970 (close to 400 ppm). Hg also showed a significant correlation above 40 ppm (−0.670 Pearson coefficient), up to −0.730 (70 ppm).

The CARI-2 index correlated better with geochemical values for As, maintaining a logarithmic (increasing) correlation that reached stabilization at 200 ppm (Pearson’s Coef. = 0.89). The positive correlation reflects a direct relationship, that is, the higher the concentration, the greater the stress present in plants. However, when the concentration reached 200 ppm, the plants reach the limit of tolerance for the vegetation described in Section 2.1.

Similarly, the NDVIg index presented an optimal correlation with Hg by means of an exponential (decreasing) function whose stabilization was achieved at 40 ppm, a concentration under the industrial limit and above the limits for residential and recreational uses, but unlike As, with a negative correlation (Pearson’s Coef. = −0.67). This implies that the higher the concentration, the lower the value of the NDVIg index. This finding allowed us to identify a direct alteration of the vegetation with increasing concentrations of Hg, since the definition of the index itself shows that the lower its value, the less photosynthetic activity occurs [61].

In any case, it can be said that the results of the vegetation indices presented high reliability for As and acceptable reliability for Hg at the points where the concentrations of these two pollutants were very high. All things considered, the mapping of the CARI2 and NDVIg indices (Figure 4) allows prediction of where the concentrations of As and Hg, respectively, will be highest, thus reducing the costs of sampling campaigns.

3.3. Geostatistical Analysis

Figure 6 shows the geostatistical interpolations obtained by ordinary kriging. As anticipated in the statistical analysis, the distribution of As and Hg was highly irregular along the UAV flight area. Extremely polluted zones near the spoil heap alternated with others downstream presenting lower concentrations of these heavy metals, although always exceeding RBSSLs. However, for academic purposes, such alternation is not necessarily negative as the differences in concentrations allowed us to test the reliability of the vegetation indices at a range of concentrations.

Our results also revealed that Hg showed less mobility than As in soils, a property that was addressed in previous studies [62,63,64]. This was appreciable since changes in concentration were far greater for Hg, or in other words, the bullseye was more easily identifiable. The lower mobility of Hg also contributes to explaining why Pearson’s Coefficient was lower for this heavy metal. In this case, low normality hampered the estimation. Nevertheless, even under this circumstance, the correlation obtained was acceptable.

3.4. Synergy between Multispectral Data and Geostatistical Analysis

Figure 7 aims to synergize the geostatistical results (Figure 6) and the multispectral analysis (Figure 4). Areas in red correspond to concentrations of As and Hg exceeding the thresholds described in the previous section (>200 mg/kg for As and >40 mg/kg for Hg). Our findings imply that, in these locations, it is possible to use multispectral techniques to predict the extent to which the pollution affects plants and subsequently predicting the degree of pollution. More specifically, around 89% of reliability was achieved for As and 67% for Hg on the basis of Pearson’s Coefficient (Figure 5).

We also observed that most of the area lent itself to study; although results are more reliable for As than for Hg. This finding could be attributable to the previously mentioned higher mobility of As, which reached higher concentrations than Hg far from the spoil heap. In contrast, the areas with low concentrations of Hg and As are shown in green, where there was no correlation between the pollution levels and the CARI2/NDVIg results.

4. Conclusions

Hg-As mining spoil heaps pose an environmental risk due to their highly toxic load. Therefore, it is important to monitor their evolution in order to evaluate the dispersion of the contaminants. Here we examined the impact of weathering on the As- and Hg-rich spoil heap of the former La Soterraña mine. To this end, the synergy between classical geochemical tools and remote sensing methodologies using a UAV-based airborne multispectral camera was explored.

First, the geochemical campaign confirmed severe As and Hg pollution of soils downstream of the spoil heap, and to a lesser extent the presence of Sb, in concentrations several orders of magnitude above RBSSLs. On the basis of the high level of pollution found, we reasoned that the vegetation of the area would be affected. This notion was confirmed by remote sensing techniques involving a multispectral camera mounted on a UAV. Algebraic operations with remote sensing data, also called vegetation indices, revealed reliable results between As and the CARI-2 index, and between Hg and the NDVIg index, the former giving better results.

In this regard, our results demonstrated a very high correlation between plant stress and high concentrations of As (0.89 correlation when As concentration > 200 ppm) and Hg (0.67 correlation when Hg concentration > 40 ppm). These findings pave the way for new studies that combine classical geochemical tools and remote sensing.

The dispersion of the contaminants is attributable to the physical-chemical weathering of the spoil heap, mainly through solid waste leaching processes, which results in the distribution of As and Hg towards lower levels of the heap. However, thanks to the high correlations between the concentrations of the two contaminants and the CARI-2 and NDVIg indices, the mapping of these indices facilitates the identification of pollution hotspots in areas particularly affected by contamination, thereby reducing the costs of geochemical characterization. On the whole, this study emphasizes the potential of combining geochemical tools and remote sensing technology for soil pollution research. However, future research efforts are expected to bring about the fine tuning of both the geochemical aspects and the multispectral information.