Incorporating Remote Sensing Techniques To The DRASTIC Index To Assess Groundwater Vulnerability In The Mining Area of the Rio Sonora Aquifer in Northwestern Mexico

Groundwater metal pollution is a major concern for societies, especially in areas where the mining industry is important. Index-based techniques, as the DRASTIC index, are often used to assess the intrinsic groundwater vulnerability and could be modied to evaluate the aquifer vulnerability to specic contaminants. Mines, mining wastes and related features are detectable with remote sensing techniques. In this work we evaluate the vulnerability of the Rio Sonora Aquifer to metallic pollution by the traditional DRASTIC method and by the addition of a land use (Lu) parameter in which possible sources of metals (detected by remote sensing) were considered (DRASTIC+Lu). The methodology allowed us to locate possible sources of metallic contamination. The Sonora River channel showed the higher vulnerability in both calculated vulnerability indices (DRASTIC and DRASTIC+Lu). Generally, the addition of the land use parameter caused a decrease in vulnerability but also a local increase where possible sources of metals were found. Thus, the modied method facilitated the identication of highly vulnerable areas which is relevant to better protect the studied aquifer.


Introduction
Groundwater is a valuable resource for human life and economic development. Its quantity and quality are of vital importance in arid and semi-arid areas, where the climatic conditions are characterized by low rainfall and high evapotranspiration, impacting sur cial water resources and the aquifer recharge. The concept of groundwater vulnerability was rst introduced by J. Margat in 1968; today it is of importance for the protection of groundwater resources. Assessing the vulnerability of an aquifer permits the identi cation of areas that are more susceptible to being contaminated, allowing to carry out effective protection measures and management plans for pollutants or wastes. Vulnerability assessment is really relevant as remediation of aquifers would be di cult and expensive (Aydi, 2018;Yin et al., 2013).
The intrinsic vulnerability de nes the vulnerability of an aquifer to a variety of pollutants, independently of their nature, and it is related to the aquifer's features (hydrological, geological, and hydrogeological) (Oke, 2020). In that sense the vulnerability depends on the resistance of the aquifer itself when receiving pollutants from outside, the lower the resistance, the greater the vulnerability. In uencing factors are: depth to groundwater; net recharge rate; aquifer media; topography; vadose zone; hydraulic conductivity; aquifer thickness; and, pumping density rate in case of over-pumping (Abu-Bakr, 2020).
Aquifers will have different reactions to different pollutants due to their physicochemical characteristics.
In those cases, it is more appropriate to talk about the speci c vulnerability which de nes the vulnerability to a speci c contaminant or group of contaminants considering the contaminants' properties and its interaction with the aquifer (Gogu and Dassargues, 2000;Voutchkova et al., 2021).
Many methods have been developed to assess the groundwater vulnerability; they can be classi ed into three types: simulation methods, statistical methods, and index methods. The index-based techniques have the advantage that they do not depend on data availability or similarities (Barbulescu, 2020), being widely used. One of the most widely used index-based methods is the DRASTIC index, developed by the United States Environmental Protection Agency (EPA) to assess the potential for groundwater contamination (Aller, 1987). DRASTIC considers seven parameters: Depth to water table (D), net Recharge (R), Aquifer media (A), Soil media (S), Topography (T), Impact of the vadose zone (I), and hydraulic Conductivity (C) which, together, form the acronym.
Frequently, new parameters are added to the seven main hydrogeological parameters of the DRASTIC index. Additional parameters used by authors include: land-use (Kozłowski and Sojka, 2019), lineament (Abdullah et al., 2015), proximity to rivers, residential areas and roads (Aydi, 2018), hydraulic parameters (Lappas, I and Matiatos, I, 2014), redox state of the aquifer (Voutchkova et al., 2021), adsorption capacity of soils (Jr and Viero, 2006), contamination index (Cd) and heavy metal pollution index (HPI) ( Groundwater contamination related to the mining industry is an important global issue. Sulphide oxidation and the associated acid mine drainage (AMD) or acid rock drainage (ARD) is considered as one of the main water pollutants in many countries that have historic or current mining activities. AMD is prominent in both active and abandoned mining sites (Simate and Ndlovu, 2014 (Wang et al., 2021), it is calculated as the level of greenness using imagery. NDVI is also a useful tool for distinguishing the boundaries of vegetated terrain from tailings impoundments, which the NDVI primarily assigns negative pixel values (Firozjaei et al., 2021;Schimmer, 2008;Zeng et al., 2017).
To our knowledge, this is the rst attempt to evaluate groundwater vulnerability to metallic pollution by the addition of a land use parameter in which possible sources of metals are considered to provide greater certainty to the vulnerability assessment at mining areas. We applied a supervised classi cation method to detect possible sources of metals on the area based on the NDVI values of known mining sites, mining wastes and mineralized areas of the study site.
The aims of the present work were: (i) propose a new method combining remote sensing and the DRASTIC procedure (modi ed DRASTIC method), (ii) identify possible sources of heavy metals (active and inactive mines, mining wastes and mineralized areas) by a remote sensing work using a supervised classi cation procedure based on NDVI, (iii) assessing the ground water vulnerability to metal pollution at the mining area of The Rio Sonora basin and, (iv) comparing and validating the results obtained by the DRASTIC method and the proposed modi ed DRASTIC method.

Study site and context
The Rio Sonora aquifer is located in northwestern Mexico at the state of Sonora (Fig. 1), covering an area of about 12,615 km2. The area accounts for a population of 23,261 inhabitants (CIAD, 2013) from Aconchi, Arizpe, Banámichi, Baviácora, Huépac, San Felipe de Jesús and Ures municipalities. Groundwater major use is for agriculture, followed by industrial and domestic usage.
It is an uncon ned aquifer; its lateral limits correspond to intrusive igneous rocks of the granitic type and extrusive rocks of the rhyolitic and andesitic type. In some areas, such as the Ures and San Felipe de Jesús valleys, the Báucarit formation emerges, which is a conglomerate complex with medium to low permeability. The Rio Sonora Aquifer is a porous-medium consisting of unconsolidated cobbles, gravels, and sands with good granular porosity; thus, good permeability is restricted to the riverbed and tributary occurs in the region. The main mining operations correspond, from north to south, to the El Gachi mine (Pb-Zn) located east of Arizpe, the Santa Elena mine (Au) localized east of Banámichi, the San Felipe mine (Cu, Pb, Zn and Au), as well as, El Jaralito (W) and the Washington mine (Cu, W and Mo) situated west and east of Baviácora, respectively (Calmus et al., 2018). The El Gachi mine is currently inactive and corresponds to a distal replacement deposit related to a porphyry copper-like environment, lodged in sediments that correspond to the Early Cretaceous (Zuñiga-Hernández, 2010). The Santa Elena mine is the largest active mining site within the study area; it is considered of hydrothermal origin, the elements with the highest abundance are Au, Ag, Zn and Pb (Calmus et al., 2018). The San Felipe mining area is currently abandoned; it has several mineral deposits in the form of hydrothermal veins. The main extraction minerals were pyrite, chalcopyrite, galena and sphalerite. This area is also characterized by the presence of abandoned tailings (Del Rio-Salas et al., 2019). The Jaralito mine is of the skarn-type deposit and has records of pegmatite bodies with W-Be mineralization (Roldán-Quintana, 1991). The Washington mine is currently abandoned and corresponds to a breccia sedimentary rock located in a sequence of volcanic rocks of the Tarahumara Formation (Simmons and Sawkings, 1983; Zuñiga-Hernández, 2010).
There are records of metal contamination in the region since the 1980s (Gomez-Alvarez et al., 1990). Del Rio-Salas et al., (2019) observed that e orescence minerals and mine tailings from the area have a high toxicity and potential to affect the quality of water (groundwater and sur cial) in the region. The mobility and accessibility of some potentially toxic elements (Zn, Pb, and As) were investigated from the mining tailings of San Felipe de Jesús and adjacent agricultural soils located within the aquifer. Zinc was mainly recovered from labile fractions in oxide-rich tailings (~ 60%) and in a lower amount from sul de-rich tailings (~ 30%). The percentage of mobile fractions (sum of water-soluble, exchangeable, and bound to carbonate fractions) in agricultural soils was as follows: Zn ~ 60%, Pb ~ 15%, and As ~ 70% The Quaternary is de ned by the presence of alluvial and uvial deposits located in the channels of the rivers, streams and ood plains (CONAGUA, 2015).

Evaluation of the intrinsic aquifer vulnerability
The DRASTIC method was used to determine the intrinsic aquifer vulnerability at the Rio Sonora aquifer.
The model was developed by the Environmental Protection Agency (EPA) in 1987 to evaluate the potential for groundwater contamination (Aller, 1987). Today it is the most frequently used method to It is based on four major assumptions: i) the contaminant is introduced at the ground surface; ii) the contaminant is ushed into the groundwater by precipitation; iii) the contaminant has the mobility of water; and iv) the area evaluated is 40 ha or larger ( Aller, 1987). The method produces index numbers derived from the rating (r) and weights (w) assigned to each parameter, the higher the DRASTIC index, the greater is the groundwater pollution potential. The DRASTIC index (D) was calculated as follows: The DRASTIC index varies from 23 to 230. The categories to interpret the DRASTIC index are: very low vulnerability (23-64), low vulnerability (64-100), medium vulnerability (106-147), high vulnerability (147-188), very high vulnerability (188-230).
The maps corresponding to the seven parameters were constructed with available hydrogeological data in the GIS-ArcView software. The depth to groundwater was estimated from 2014 piezometric data from 168 wells. Piezometric data was obtained from the Water National Commission (CONAGUA) from Wells located along the Sonora River. For the areas where there was no data, it was interpolated using the Inverse distance weighted interpolation based on existing values. The recharge was estimated from daily precipitation data (1925 to 2012) obtained from CONAGUA weather stations at Arizpe, Sinoquipe, Banámichi, Huépac, Aconchi, Mazocahui, Ures and Topahue (Fig. 1) using the formula: Where RN is the net recharge (mm), P is the mean annual precipitation (mm) and Er is the real annual evapotranspiration (mm). The Er was estimated using the Turc method (Belmonte-Jiménez et al., 2005).
Recharge data was interpolated to the whole study area using the Inverse distance weighting (IDW) method. The aquifer media and the Impact of the vadose zone thematic maps were constructed based on existing geological maps from the Sonora Water Commission; data was reclassi ed according to the types of lithology recommended by Aller, (1987). The soil media thematic map was constructed based on the World Soil Information (2017) data base; data was reclassi ed according to the soil types In this work, remote sensing was used in order to locate potential sources of metals related to mining activity (mining tailings, active/inactive mines) or areas with exposed natural geochemical anomalies in an e cient and economical way covering large areas affected by mining activities. Freely-available satellite images were obtained from SENTINEL 2 at a resolution of 10 m; they were processed with the free access Quantum GIS software (QGIS Development Team, 2019). We used two May 22, 2019 images covering the study area presenting 0.1% cloudiness, in a UTM / WGS84 projection (Universal Transversal Mercator) of 100km x 100km. We selected the bands 2, 3, 4 and 8 with a spatial resolution of 10 meters. They were processed using QGIS and the Semi-Automatic Classi cation Plugin (SCP) to carry out the atmospheric correction using the dark pixel subtraction method (DOS1). Selected bands were combined to produce a mosaic dataset using the Mosaic Raster Layer de SAGA (System for Automated Geoscienti c Analyses) by the Nearest Neighbor method.
Spectral characteristics of tailings and other background features (water bodies, vegetation and bare soil) were studied by the photointerpretation of 5 points (selected manually) of each class located inside the study area.
NDVI was employed to distinguishing the boundaries of vegetated terrain from tailings impoundments, mining sites and areas with exposed natural geochemical anomalies, which are characterized by negative pixel values. The NDVI assigns values in a range from -1 to 1, negative values and close to -1 correspond to bodies of water and negative values close to zero correspond to areas with bare or rocky soil, sand or snow (Saravanan et al., 2019); thus, we focus on negative values close to 0. NDVI was calculated as follows:

NDVI = NIR − Red NIR + Red
To increase the precision of the method and effectively distinguish from bare soils and built-up areas, the optimal NDVI threshold was de ned based on the minimum and maximum NDVI values exhibited by 5 training elds (areas of known identity) including mining wastes and active/inactive mines of the study site (Table 1). The polygons generated as described above were subsequently veri ed and validated manually to ensure the selection of polygons corresponding to the class "possible sources of metallic contamination" counting for mining tailings, active/inactive mines and areas with exposed natural geochemical anomalies. We considered areas of at least 400 m 2 or four pixels. Detection was validated by comparing visually the identi ed polygons to existing information (Zeng et al., 2017).

DRASTIC + Lu
The DRASTIC method was modi ed by introducing a Land Use (Lu) parameter (100m 2 resolution) in which the potential sources of metals (identi ed using remote sensing) were integrated as one class type. The original Lu data was obtained from the Mexican National Institute of Statistic and Geography (INEGI, 2016). Ratings were de ned based on the potential to lead to metallic contamination of each class. Assigned ratings and the Lu parameter weight are shown in Table 2. Most of the study area is covered by secondary vegetation, grasslands, forest and scrub, with assigned rating of 1, considered with little in uence on metal pollution. Ratings of 5, 6 and 7 were assigned to areas with human activity. Ratings to agriculture activity depends on irrigation type as the water in ltration promotes the mobility of contaminants to the aquifer. Finally, a rating of 10 was assigned to the mining and mineralized areas due to their polluting potential. The DRASTIC + Lu index was calculated as follows: D Lu = D r D w + R r R w + A r A w + S r S w + T r T w + I r I w + C r C w + Lu r Lu w In order to compare both methods, the results obtained in the original and modi ed DRASTIC method were normalized to scale from 0 to 100. The index was divided into ve categories: very low vulnerability (<20), low vulnerability (20-40), medium vulnerability (40-60), high vulnerability (60-80), very high vulnerability (80-100).

Map removal sensitivity analysis of the DRASTIC + Lu
To investigate the effects of adding the Lu parameter on the vulnerability map, we performed the map removal sensitivity analysis which show the sensitivity of the vulnerability map by removing one or more parameters from the analysis. It is computed in the following way (Babiker et al., 2005): where S is the sensitivity measure expressed in terms of variation index, V and V´ are the unperturbed (DRASTIC) and the perturbed (DRASTIC + Lu) vulnerability indices respectively, and N and n are the number of data layers used to compute V and V´.

Correlation between groundwater metal concentrations and calculated indexes
The groundwater metal concentrations from governmental surveys (public data) carried out by the National Water Commission (CONAGUA) were used to associate and correlate the metal pollution in the groundwater to the vulnerability indexes obtained with the original DRASTIC and modi ed methods. Used data is freely available at http://www. deicomisoriosonora.gob.mx/ deicomiso.html and corresponds to Al, Cu, Zn, Mn y Fe average concentrations (measurements between august and september 2014 in mg L −1 ) from 28 wells located following the Sonora River (Figure 1). We calculated the spearman correlation due to the non-normality of data.

Assessment of the intrinsic aquifer vulnerability with DRASTIC
The seven thematic maps representing the parameters of the DRASTIC method were prepared using existing input data based on GIS.
Groundwater Depth (D) According to the input data obtained from CONAGUA, the parameter varies between 0 and 23 m. The ratings assigned according to the DRASTIC method (Aller, 1987) ranged between 3 and 10 ( Fig. 2). Input data was obtained from points located following the Sonora River from north to south of the studied aquifer and not well distributed on the study area; thus, some imprecision could be attributed. Nevertheless it is largely accepted that in the area groundwater depth ranges from 1 to 30 m, increasing in the direction of the Sonora River, which is shaped by the basement topography underlying alluvial and uvial deposits, as well as by the narrowing of the Rio Sonora that lifter the water level, causing it to emerge (Archundia et al., 2021; CONAGUA, 2015).
Recharge (R) This parameter was calculated based on precipitation input data from 8 climatological stations dating between 1925 -2012 and IDW interpolation. Some imprecision could be expected as climatological stations were located following the Sonora River, where the municipalities of the area are located. Recharge values of 0 to 51 mm were obtained in the central part of the southern portion, values between 51 and 102 mm were obtained on the margins of the southern portion and values from 102 to 178 mm were found in the upper and more extensive part. The ratings assigned according to the DRASTIC method (Aller, 1987) were between 1 and 9 (Fig. 3).
Aquifer media (A) The lithology of the study area is mainly of igneous and metamorphic rocks, which are poorly permeable and non-porous unless fractured (rating 3). Sequences of sedimentary rocks, such as sandstones and shales, appear sporadically within the study area (rating 6). In some areas located near of San Felipe de Jesús and Arizpe conglomerates and basalt occur, and a high rating was assigned (7 and 8 respectively) as they present interconnected pores and fractures. Conglomerate are also common in the surrounding of the Sonora River channel. The alluvium (rating 10) is mainly located in the Sonora streambed and its tributaries (Fig. 4).
Soil (S) According to the World Soil Information (2017), four soil texture types occur in the study area (Fig.  5). Soils with a large fraction of clay occur surrounding the Sonora River (rating 3). Sandy-clay soils are present in the lower proportions and especially in the middle and southern portion of the study area (rating 4). Clay-loam soils are located on high altitudinal areas (rating 5). Sandy clay loam soils are found in the highest altitudinal areas where igneous and metamorphic rocks are located (rating 6) (Fig. 5). It is important to note that the World Soil Information combines existing regional and national updates of soil information worldwide so it may lack precision at a local scale.
Topography (T) Most of the study area has slope values greater than 18%, near the Sonora River and at the southeastern part very low and null slope values occur. Rating varies between 1 and 10 ( Figure 6).
Impact of the vadose zone (I) Igneous and metamorphic rocks with reduced permeability are located at the eastern and western parts of the study area (rating 4). Sedimentary and permeable rocks occur from north to south following the Sonora River channel (rating 6). Basaltic rocks, which could show signi cant porosity and permeability, are in spots at the north and the center (rating 7). Finally, the alluvium, highly permeable, is found from north to south associated to the Sonora River channel (rating 8) (Fig. 7).
Hydraulic Conductivity (C) Hydraulic conductivity values were obtained from high resolution global lithology existing maps; thus some imprecision could be expected. Most of the study area present values between 0 and 28.55 cm day −1 . Some spots which coincide with the presence of carbonate rocks show values greater than 81 cm day −1 , corresponding to an average conductivity in free aquifers (Fig. 8). The area with the lowest hydraulic conductivity (rating 1) corresponds to the igneous body known as the Aconchi batholith.
DRASTIC index. The DRASTIC vulnerability index was computed overlaying the seven hydro-geological parameter maps presented above: values range from 64 to 171. The groundwater vulnerability of the Rio Sonora aquifer uctuates from low, at some places in the center and the south, to high (Fig. 9). The results of this study show that 23.61 % of the study area has low pollution potential, 72.95 % has moderate pollution potential and 3.43 % has high pollution potential. The distribution of the Rio Sonora aquifer intrinsic vulnerability is mainly conditioned by the lithology and the impact of the vadose zone. The yellow and green colors correspond to areas with igneous rocks, which were classi ed as low permeable rocks. In the center of the study area, the orange color corresponds to the conglomerate rocks and the red color represents the channel of the Sonora River where the alluvium is found. In the southeast zone there are also colorations that indicate high vulnerability, most likely due to the lithology (alluvial deposits) and a low slope. This is of great importance since most of the drinking water supply wells are located near the channel of the Sonora River of high vulnerability.

Evaluation of the Rio Sonora aquifer vulnerability to metal pollution (DRASTIC + Lu index)
The DRASTIC method was modi ed by introducing an extra parameter considering the Land Use (Lu) in which the potential sources of metals, identi ed by remote sensing, were incorporated.
3.2.1 Detection of possible sources of metallic contamination by remote sensing Figure 10 show the averaged spectral characteristics of tailings, water bodies, vegetation and bare soil in the study site. When comparing spectral signatures, we notice the different behaviors of background features, which validate the use of remote sensing for the identi cation of possible sources of metallic contamination as mine tailings. The group of "Water body" and "Vegetation" show typical behaviors. The re ectance value of a water body decreases as the wavelength increases, due to low re ection in the near infrared. Vegetation presents a peak in the infrared spectrum. Both the Mine Tailings and the Bare Soil groups present high re ectance values in the studied bands. However, the rst tend to have a re ectance peak in the red band and higher re ectance values.
The NDVI values threshold de ned based on considered training elds was between -0.2504 and 0.1946 it allowed detecting non-vegetated areas including considered training elds (Fig. 11).
Base on de ned NDVI threshold and after the manual ltration, 60 areas representing the possible sources of metallic pollution were located ( Fig. 12

DRASTIC + Lu index
Identi ed possible sources of metallic pollution were integrated to the Land use information of the Mexican National Institute of Statistic and Geography (INEGI, 2016). The DRASTIC +Lu index was computed as described in the methodology section. Figure 13 shows the DRASTIC+Lu index which varied from 69 to 201. The study area presents moderate and high vulnerability values, according to the vulnerability scale proposed by Aller (1987

Comparison between calculated indexes.
Normalized vulnerability maps (on a scale of 1 to 100) obtained by the DRASTIC and DRASTIC+Lu methods are shown in Figure 14. We notice that in both maps the areas of lower vulnerability are in upper parts, corresponding to igneous rocks and areas with pronounced slope. Greatest vulnerability values (dark red color) correspond to conglomerate rocks, the presence of agricultural areas and low slope. The mining and mineralized areas are almost imperceptible on the scale of the maps, but they increased the vulnerability at a local level. Comparing both maps, a general decrease in vulnerability can be observed in the DRASTIC +Lu map. The area of very high vulnerability (80-100) decreased by 15%, the extent of high vulnerability (60-80) decreased by 81.3%, the zone of medium vulnerability increased slightly (8.4%), the area of low vulnerability increased in 94% and the area of very low vulnerability increased by 154.9%.
This analysis allowed the identi cation of high vulnerability areas on which protection efforts need to be focused and to localize small sites with strong negative local effects.

Validation of the proposed DRASTIC + Lu method
The map removal sensitivity analysis was preferred because it tests the sensitivity of operations between map layers (Thapa, 2018) helping to evaluate its in uence on the vulnerability assessment. Table 3   To better understand the metal (and others) pollutants fate and dynamics in the study area it is important to perform studies on the prediction of groundwater ow, super cial soil and subsurface (vadose zone) properties and reactivity against metals. In that sense and in order to improve the vulnerability assessment it is suggested to characterize soil types and properties as well as ground water depth at a local scale.

Conclusion
Page 15/31 The proposed remote sensing methodology allowed to locate 60 possible sources of metallic contamination, corresponding to mining tailings, open-pit wastes, areas with exposed natural mineralization and areas with current mining activity. Within the identi ed sources, only three correspond to exposed natural mineralization areas, possibly indicating that the method might not be sensitive enough. The use of higher resolution satellite images could improve the precision to locate them but signi cantly increasing the costs of the method.
The Sonora River channel has the higher vulnerability in both calculated vulnerability indices (DRASTIC and DRASTIC+Lu) indicating that the aquatic resources of the populations located within the study area are threatened. To improve the precision of the modi ed method results, an effort in the characterization of the water depth and soil type information at the local scale is needed. With the application of the modi ed DRASTIC index and, in comparison with the results obtained with the traditional DRASTIC index, a decrease in vulnerability was observed, mainly in conglomerate rocks and higher elevation areas located parallel to the Sonora River channel. This general decrease in vulnerability coupled with a local increase of vulnerability where possible sources of metals are located, facilitated to highlight highly vulnerable areas which is relevant to better protect the studied aquifer. Figure 1 Study area corresponding to the Sonora River aquifer. The gure contains the municipalities, known mineralizes zones, mine tailings, mining sites and hydrogeological sampling points used to validate the results.

Figure 11
NDVI at a portion of the study site and the San Felipe de Jesús training eld.

Figure 12
Possible sources of metallic pollution detected by remote sensing and training elds.

Figure 13
The Sonora River DRASTIC+Lu index. Figure 14