**A Novel Approach to Harmonize Vulnerability Assessment in Carbonate and Detrital Aquifers at Basin Scale**

### **Leticia Baena-Ruiz \* and David Pulido-Velazquez**

Spanish Geological Survey (IGME), Urb. Alcázar del Genil, 4. Edificio Zulema, bajo, 18006 Granada, Spain; d.pulido@igme.es

**\*** Correspondence: l.baena@igme.es

Received: 23 September 2020; Accepted: 21 October 2020; Published: 23 October 2020

**Abstract:** The DRASTIC (D: Depth to water; R: Net recharge; A: Aquifer media; S: Soil media; T: Topography; I: Impact of vadose zone; C: Hydraulic conductivity) index is usually applied to assess intrinsic vulnerability in detrital and carbonate aquifers, although it does not take into account the particularities of karst systems as the COP (C: Concentration of flow; O: Overlying layers above water table; P: precipitation) method does. In this paper we aim to find a reasonable correspondence between the vulnerability maps obtained using these two methods. We adapt the DRASTIC index in order to obtain reliable assessments in carbonate aquifers while maintaining its original conceptual formulation. This approach is analogous to the hypothesis of "equivalent porous medium", which applies to karstic aquifers the numerical solution developed for detrital aquifers. We applied our novel method to the Upper Guadiana Basin, which contains both carbonate and detrital aquifers. Validation analysis demonstrated a higher confidence in the vulnerability assessment provided by the COP method in the carbonate aquifers. The proposed method solves an optimization problem to minimize the differences between the assessments provided by the modified DRASTIC and COP methods. Decision trees and spatial statistics analyses were combined to identify the ranges and weights of DRASTIC parameters to produce an optimal solution that matches the COP vulnerability classification for carbonate aquifers in 75% of the area, while maintaining a reliable assessment of the detrital aquifers in the Basin.

**Keywords:** groundwater vulnerability; carbonate aquifers; optimized DRASTIC; COP; decision trees; nitrate validation

### **1. Introduction**

Groundwater pollution is a widespread problem affecting most aquifers all over the world due to the increasing agricultural and industrial human activity [1,2]. It is the result of the interaction between the anthroposphere and the hydrosphere, where substances from the different land uses penetrate the groundwater, leading to impacts on environmental water quality and human health. The protection of the groundwater resources has become a priority due to the importance of groundwater for human supply, irrigation and dependent ecosystems, especially in semi-arid regions [3–5]. The degree of protection depends on intrinsic groundwater vulnerability, which is defined as the susceptibility of aquifers to pollution arising from anthropogenic activity [6]. Several methods have been proposed to assess intrinsic vulnerability by using conceptual groundwater flow models and taking transport processes into account [7,8]. The most frequently employed methods are the index-based approaches [9,10] of which DRASTIC (D: Depth to water; R: Net recharge; A: Aquifer media; S: Soil media; T: Topography; I: Impact of vadose zone; C: Hydraulic conductivity) [9] is the most common. DRASTIC has been

applied to various types of aquifer, though some authors have revealed that it is not suitable for assessing vulnerability in some karstic aquifers [10–13].

Karstic aquifers present physical particularities that, in general, make them more sensitive to contamination [14,15]. Specific approaches, such as the COP (C: Concentration of flow; O: Overlying layers above water table; P: precipitation) method [10], have been developed specifically to assess groundwater vulnerability in carbonate (karstic) aquifers. The COP method has been extensively applied in different research studies to this type of aquifer [16–21].

However, despite significant physical differences between detrital and karstic aquifers, the same approaches are often applied to simulate groundwater processes in both types of aquifers. In numerical groundwater flow models, the "equivalent porous medium" assumes the validity of the Darcy's law in karstic aquifers, which is the most frequent numerical approach deduced for detrital aquifers. We find many examples of application of this approach both for detrital aquifers [22–24] and karstic aquifers [25,26]. In the same way, index methods such as DRASTIC could be suited to assessing vulnerability in both detrital and karstic aquifers.

In this context, recent research studies have attempted to adapt the DRASTIC method to karst aquifers in recent years [12,27–30]. Some of these [12,28,29] modified DRASTIC by including and/or removing parameters to account for karst characteristics. Other studies [27,30] adapted DRASTIC by modifying the classic assessment according to the distribution of nitrate concentration, although the most contaminated areas do not always imply higher vulnerability [31,32]. A detailed review of different approaches proposed to modify the DRASTIC method is given in [13].

These modifications of the DRASTIC method were performed using different methodologies/approaches: single-parameter sensitivity analysis [3,5,33], calibration by correlation analysis with pollutants [4,34,35], analytic hierarchy process [36–38] amongst others. Different data mining algorithms have also been applied to improve DRASTIC performance [39–42], predict areas vulnerable to groundwater contamination or to identify hydrogeological factors influencing groundwater contamination [39,43,44]. Data mining algorithms aim to extract knowledge from previously unknown and indistinguishable data, and are used as operational tools to find optimal solutions in high-dimension problems. Although many authors have employed these techniques for different purposes, more research is needed to explore the potential of the statistical techniques—including data mining—to deal with the uncertainties in the weights and ratings of index-based methods in the groundwater vulnerability assessments [45].

In general, the aforementioned index-based methods have been proven to provide satisfactory vulnerability assessment in a variety of regions. However, different index methods often provide dissimilar results in a single study area, making it difficult to compare and validate results [46–48]. In complex groundwater systems containing various aquifers types, an integrated vulnerability assessment is needed to homogenize criteria and compare results at basin scale and between different case studies. In general, a harmonized method that is applicable to all aquifer types would be more likely to be implemented worldwide [12,49]. Previous work [12] has aimed at developing a new method (DRISTPI; Depth to water, Recharge, Impact of vadose zone, Soil, Topography, Preferential Infiltration) by adding and removing certain parameters from the classic DRASTIC method that can be applied to any type of aquifer. It was tested in karstic aquifers, but not in detrital aquifers.

In this paper we propose a novel approach to standardize the DRASTIC method to assess vulnerability in a basin composed by both, detrital and carbonate aquifers. We proposed an adaptation of the most commonly applied method for detrital aquifers—the DRASTIC method—to obtain a reliable vulnerability assessment for carbonate aquifers. The DRASTIC method can be applied in both detrital and carbonate aquifers, though it does not always provide appropriate vulnerability assessment in carbonate aquifers. In contrast, the COP method is specific to carbonate aquifers and it yielded better results in the validation analysis in our case study. Therefore, we proposed to adapt DRASTIC to harmonize the vulnerability assessment for different aquifer typologies, in order to make the results comparable in a basin that contains a wide variety of geological formations. Harmonization of criteria in the vulnerability assessment would allow the dissimilarities provided by different vulnerability methods to be dealt with. Although in the literature we can find several attempts to adapt the DRASTIC method to the particularities of each case study, none of them have proved the applicability of the modified DRASTIC in different aquifer typologies. *Water* **2020**, *12*, x FOR PEER REVIEW 3 of 26 literature we can find several attempts to adapt the DRASTIC method to the particularities of each case study, none of them have proved the applicability of the modified DRASTIC in different aquifer

In this study, an optimization problem is solved to establish a correspondence between DRASTIC and COP by maintaining the original formulation of the DRASTIC method. It aims to find the weights and ranges of the DRASTIC parameters that maximize the correspondence between the qualitative vulnerability classes obtained using DRASTIC and COP methods. The optimization problem is solved through a new approach that combines spatial statistics analysis and data mining (decision trees). Although decision trees have been applied to predict the sensitivity to contaminants based on groundwater vulnerability [50,51], they have not been used to adapt/improve the DRASTIC method. typologies. In this study, an optimization problem is solved to establish a correspondence between DRASTIC and COP by maintaining the original formulation of the DRASTIC method. It aims to find the weights and ranges of the DRASTIC parameters that maximize the correspondence between the qualitative vulnerability classes obtained using DRASTIC and COP methods. The optimization problem is solved through a new approach that combines spatial statistics analysis and data mining (decision trees). Although decision trees have been applied to predict the sensitivity to contaminants based on groundwater vulnerability [50,51], they have not been used to adapt/improve the DRASTIC

The optimization methodology was applied to five carbonate aquifers in the Upper Guadiana Basin, where the validation analysis demonstrated that COP method produces better vulnerability assessment than the original DRASTIC one. The optimal solution obtained for carbonate aquifers (O-DRASTIC) was also tested in the three detrital aquifers within the basin and the validation analysis also shows significant improvement in the results comparing with the original DRASTIC. A sensitivity analysis of the changes introduced to define the optimum DRASTIC reveals the influence of the various intrinsic characteristics on the severity of vulnerability for different aquifer typologies. method. The optimization methodology was applied to five carbonate aquifers in the Upper Guadiana Basin, where the validation analysis demonstrated that COP method produces better vulnerability assessment than the original DRASTIC one. The optimal solution obtained for carbonate aquifers (O-DRASTIC) was also tested in the three detrital aquifers within the basin and the validation analysis also shows significant improvement in the results comparing with the original DRASTIC. A sensitivity analysis of the changes introduced to define the optimum DRASTIC reveals the influence of the various intrinsic characteristics on the severity of vulnerability for different aquifer typologies.

#### **2. Materials and Methods 2. Materials and Methods**

#### *2.1. Study Area and Data 2.1. Study Area and Data*

The Upper Guadiana Basin (Figure 1) is located in the central part of Spain. It is composed of eight groundwater bodies including five unconfined mixed (carbonate and detrital) and three unconfined detrital aquifers extending over approximately 14,000 km<sup>2</sup> . The Upper Guadiana Basin (Figure 1) is located in the central part of Spain. It is composed of eight groundwater bodies including five unconfined mixed (carbonate and detrital) and three unconfined detrital aquifers extending over approximately 14,000 km2.

**Figure 1.** Location of study area within the Upper Guadiana Basin (Spain). **Figure 1.** Location of study area within the Upper Guadiana Basin (Spain).

The climate in this area is typically continental and semiarid. Mean annual precipitation over the period 1974–2015 was 445 mm and the mean annual temperature is 14.7 °C. Mean potential evapotranspiration is 700 mm/year. The climate in this area is typically continental and semiarid. Mean annual precipitation over the period 1974–2015 was 445 mm and the mean annual temperature is 14.7 ◦C. Mean potential evapotranspiration is 700 mm/year.

The area is predominantly flat, bounded by mountain landscapes to the north (Sierra de Altomira) and south (Campo de Montiel). The area is predominantly flat, bounded by mountain landscapes to the north (Sierra de Altomira) and south (Campo de Montiel).

Connectivity between groundwater bodies is structurally complex, with strong natural interaction between groundwater and surface water. Under natural conditions groundwater flow discharges into the central aquifer of the system (Mancha Occidental Aquifer) [52] which flows eastwards. Connectivity between groundwater bodies is structurally complex, with strong natural interaction between groundwater and surface water. Under natural conditions groundwater flow discharges into the central aquifer of the system (Mancha Occidental Aquifer) [52] which flows eastwards.

The predominantly dry climate and the prevalence of irrigated agriculture means that the Upper Guadiana Basin has been intensively pumped, and this has led to some of its aquifers being declared overexploited [53]. The water table is highly variable, lying at less than 1 m to more than 250 m, though over most of the basin it lies below 15 m.

Rainfall is the main source of aquifers' recharge. Mean annual recharge varies between 45 and 70 mm/year, although there is disagreement about these values [52–54].

The soils in the basin mainly belong to the cambisol group according to the FAO (Food and Agriculture Organization of the United Nations) classification [54,55]. Regosol and others, such as luvisol and podzol, can be found in the southeast [54]. Soil texture in the northern part of the basin is predominantly silty-loam, whereas the soil in the southern part it is peaty.

The geology is complex, including mixed carbonate–detrital aquifers. In the southern half of the Upper Guadiana Basin, the aquifers are predominantly composed of limestones, with many karstified zones [53]. In general, the karst is not very developed and there are no swallow holes in these aquifers. Many parts of the central aquifer are formed by Tertiary detrital deposits [53]. The northern aquifers are more heterogeneous. There are no large karstified areas and other formations of metamorphic materials can be found. The detrital aquifers are mainly composed of Tertiary and Quaternary alluvial materials.

The unsaturated zone is formed by poorly permeable lithologies in the northern part of the basin, with higher permeabilities in the southern part [54].

Conductivity in the Upper Guadiana Basin varies widely. To the north, conductivity is low (below 1.5 <sup>×</sup> <sup>10</sup>−<sup>4</sup> <sup>m</sup>/s); in the central part there are zones with values higher than 5 <sup>×</sup> <sup>10</sup>−<sup>4</sup> <sup>m</sup>/s, while conductivity to the southern is mainly in the range 3.5 <sup>×</sup> <sup>10</sup>−4–5 <sup>×</sup> <sup>10</sup>−<sup>4</sup> <sup>m</sup>/s [56].

### *2.2. DRASTIC and COP Vulnerability Maps*

DRASTIC and COP vulnerability maps were calculated in the five mixed aquifers following the proposal made in [9,10], respectively (Figure 2). A detailed explanation of these methods can be found in the Appendix A. Tables A4 and A5 (Appendix B) summarizes the data sources and the methodology applied to calculate the parameters of DRASTIC (hereinafter, we will call it "original DRASTIC") and COP, respectively.

The "original DRASTIC" values vary between 41 and 171 within the study area. The final index was reclassified into five qualitative classes (Very low: <100; Low: 100–120; Moderate: 120–160; High: 160–180; Very high: ≥180) by grouping the original categories proposed in [9] to obtain the same number of classes as the COP vulnerability map. The values of parameters Aquifer media, Soil media, Impact of vadose zone and Conductivity (Figure 2(a3,a4,a6,a7)) show significant heterogeneity, with clear differences between the northern and central–southern areas. These differences are finally reflected in the vulnerability map (Figure 2(a8)).

Figure 2b shows the results of applying the COP method as proposed by [10]. The protection conferred by Overlaying the layers (Figure 2(b2)) generally decreases from north to south. It is very low in the southern part due to the presence of limestones and dolomites outcrops. The surface features related to the Concentration of flow (Figure 3(b1)) produce only a slight reduction of protection over 80% of the area, except in the southern Campo de Montiel aquifer, where the protection from Overlying layers is greatly reduced.

While the highest values of the DRASTIC index are located in the center and southern part of the basin (classified as "Moderate vulnerability"), the highest values of the COP index appear in the southern zone (Campo de Montiel aquifer). Both indices give their lowest values in the north area (Sierra de Altomira and Lillo Quintanar aquifers).

The percentage overlap between the vulnerability classes obtained using COP and DRASTIC is 55.75%. There is a significant coincidence in the "Very low" and "Low" vulnerability classes, which cover nearly 80% of the area in the COP map. However, the coincidence in "High" and "Very high" classes is almost null. A misclassification of the highest vulnerability areas of groundwater dependent ecosystems such as the Upper Guadiana Basin would lead to erroneous planning and management decisions, possibly leading to significant environmental impacts.

**Figure 2.** (**a**) DRASTIC (D: Depth to water; R: Net recharge; A: Aquifer media; S: Soil media; T: Topography; I: Impact of vadose zone; C: Hydraulic conductivity) maps: (**a1**) Depth to water; (**a2**) Net recharge; (**a3**) Aquifer media; (**a4**) Soil media; (**a5**) Topography; (**a6**) Impact of vadose zone; (**a7**) Hydraulic conductivity; (**a8**) Vulnerability map; (**b**) COP(C: Concentration of flow; O: Overlying layers above water table; P: precipitation) maps: (**b1**) Reduction of protection due to Concentration of flow; (**b2**) Degree of protection from Overlying layers; (**b3**) Reduction of protection due to the precipitation factor; (**b4**) Vulnerability map. **Figure 2.** (**a**) DRASTIC (D: Depth to water; R: Net recharge; A: Aquifer media; S: Soil media; T: Topography; I: Impact of vadose zone; C: Hydraulic conductivity) maps: (**a1**) Depth to water; (**a2**) Net recharge; (**a3**) Aquifer media; (**a4**) Soil media; (**a5**) Topography; (**a6**) Impact of vadose zone; (**a7**) Hydraulic conductivity; (**a8**) Vulnerability map; (**b**) COP(C: Concentration of flow; O: Overlying layers above water table; P: precipitation) maps: (**b1**) Reduction of protection due to Concentration of flow; (**b2**) Degree of protection from Overlying layers; (**b3**) Reduction of protection due to the precipitation factor; (**b4**) Vulnerability map. *Water* **2020**, *12*, x FOR PEER REVIEW 6 of 26

**Figure 3.** Validation of vulnerability maps in mixed aquifers: (**a**) DRASTIC + 5 × land use (LU); (**b**) COP + 5/LU; (**c**) mean nitrate concentration map (1975–2015). **Figure 3.** Validation of vulnerability maps in mixed aquifers: (**a**) DRASTIC + 5 × land use (LU); (**b**) COP + 5/LU; (**c**) mean nitrate concentration map (1975–2015).

*2.3. Validation of Vulnerability Maps* 

DRASTIC pollution index by other authors [27,60]:

The rates assigned to LU [3,27,33,34] are shown in Table 1.

The DRASTIC vulnerability map is validated by adding an LU factor to the DRASTIC index (Equation (1)), as proposed in other studies [27,58,59]. Making an analogy, the index defining protection against pollution could be based on the COP vulnerability (Equation (2)). In this case the inverse of the LU factor is considered since a higher COP index indicates lower vulnerability levels. Since the principal land use is agriculture, the pollution index and the protection-against-pollution index will maintain the factor of 5 used to weight the LU term, as was originally included in the

Protection െ against െ Pollution index (COP) = COP index +

**Table 1.** Rates of LU.

Pollution index (DRASTIC) = DRASTIC index + 5 × LU (1)

**LU Rate** 

Agro-forestry areas 7 Airports 2

Broad-leaved forest 2 Burnt areas 2 Complex cultivation patterns 8

Annual crops associated with permanent crops 8

5

LU (2)

Contaminant loads are linked to human activities and land use (LU). Therefore, in addition to

### *2.3. Validation of Vulnerability Maps*

Contaminant loads are linked to human activities and land use (LU). Therefore, in addition to aquifer vulnerability, these are a key factor in determining groundwater contamination [3,57]. In the study area, the intensive agriculture and the use of nitrogen fertilizers has provoked high levels of nitrates in many groundwater areas. Since nitrate is not naturally found in groundwater, it is considered to be a good indicator of contamination by human impact, especially in agricultural zones. The DRASTIC vulnerability map is validated by adding an LU factor to the DRASTIC index (Equation (1)), as proposed in other studies [27,58,59]. Making an analogy, the index defining protection against pollution could be based on the COP vulnerability (Equation (2)). In this case the inverse of the LU factor is considered since a higher COP index indicates lower vulnerability levels. Since the principal land use is agriculture, the pollution index and the protection-against-pollution index will maintain the factor of 5 used to weight the LU term, as was originally included in the DRASTIC pollution index by other authors [27,60]:

$$\text{Pollution index} \left( \text{DRASTIC} \right) = \text{DRASTIC index} + 5 \times \text{LU} \tag{1}$$

$$\text{Projection} - \text{againt} - \text{Position index} \left( \text{COP} \right) = \text{COP} \text{ index} + \frac{5}{\text{LU}} \tag{2}$$

The rates assigned to LU [3,27,33,34] are shown in Table 1.


**Table 1.** Rates of LU.

The correlations (R-squared coefficient) with the pollution index derived from DRASTIC and COP methods (Figure 3) were determined from the 214 observation points where nitrate concentration data was available (1974 to 2015 provided by the Spanish Geological Survey and the River Basin Authority). The DRASTIC pollution index (Figure 3a) is weakly correlated with the mean nitrate concentration (R<sup>2</sup> = 0.04) in the carbonate aquifers. However, the protection-against-pollution index linked to the COP map (Figure 3b) shows strong correlation (R<sup>2</sup> = 0.78). Accordingly, in this validation analysis, the COP vulnerability assessment for the study area provides a more reliable assessment than the DRASTIC index one.

As can be observed in Figures 2 and 3, both DRASTIC and COP indices suggest low vulnerability in some zones with high nitrate concentration (in the western part of the basin) that results from the intensive agricultural exploitation since early 1970's. Thus, groundwater contamination not only depends on intrinsic vulnerability, but also the land use, type of crops, type of irrigation, etc. For this reason, the most contaminated areas do not always correspond to the most vulnerable ones [31,32].

### *2.4. Methodology: Optimization of DRASTIC Method*

In this paper we adapt the DRASTIC method to minimize differences with the COP vulnerability map for the five carbonate aquifers in the Upper Guadiana Basin, which has been proven to provide better results in those aquifers according to the validation analysis. The main objective is to obtain a harmonized method that allows vulnerability in carbonate and detrital aquifers to be assessed in a homogenous way, in a system comprising varying geological formations. The harmonization of criteria to assess groundwater vulnerability will make the results comparable at basin scale. To this end, the DRASTIC index is recalculated by varying the ranges and weights of its parameters (fulfilling certain constraints to maintain a rational definition), in order to minimize the differences with the COP vulnerability map for the five mixed (carbonate and detrital) aquifers in the Upper Guadiana Basin. Data mining techniques are then applied to identify the ranges and weights that provide an optimal solution.

Figure 4 shows the flowchart of the proposed optimization problem and the methodology applied. Two objective functions (O.F.) were tested: (A) to maximize the percentage of spatial coincidence (area, Si) between DRASTIC and COP vulnerability classes; and (B) to minimize the distance (di; see Equation (3)) between the vulnerability classes in DRASTIC and COP. The decision variables in this optimization problem are (1) the ranges (of the DRASTIC index and its parameters) and (2) the weights of the DRASTIC parameters.

coincidence (area, Si) between DRASTIC and COP vulnerability classes; and (B) to minimize the distance (di; see Equation (3)) between the vulnerability classes in DRASTIC and COP. The decision

**Figure 4.** Flowchart of methodology. **Figure 4.** Flowchart of methodology.

The ranges of the DRASTIC parameters and index are proposed based on the available data in the study area and covering a wide range of hypothetical cases. The following constraints regarding the ranges proposal are imposed: The ranges of the DRASTIC parameters and index are proposed based on the available data in the study area and covering a wide range of hypothetical cases. The following constraints regarding the ranges proposal are imposed:




\* Original classification of DRASTIC parameters.


*Water* **2020**, *12*, 2971

For each combination of parameter ranges and weights used to generate a DRASTIC map, the distance (*d<sup>i</sup>* ) or difference between vulnerability classes reported by DRASTIC and COP is calculated by Equation (3):

$$d\_{\bar{l}} = \sum \alpha\_{\bar{j}k} \times x\_{\bar{j}k} \tag{3}$$

where


In order to establish a dimensionless threshold to quantify the distance (*d<sup>i</sup>* ), it is expressed as a percentage of the maximum calculated distance (*dmax*) over all calculated DRASTIC indices:

$$d\_i(\%) = \frac{d\_i}{d\_{\text{max}}} \times 100\tag{4}$$

In order to reduce the number of calculations, we employ data mining techniques (decision trees) to select the values of the variables domain to be tested. The optimal solution is sought following two steps:

	- i. First, DRASTIC vulnerability maps are calculated modifying the ranges of parameters and the classification of the DRASTIC index. The weights of parameters are the same as proposed in [9].
	- ii. All the DRASTIC indices are evaluated through the objective functions and the results are classified intro three categories (Table 3).
	- iii. Decision trees are applied in order to find out the ranges for each parameter that gives the highest coincidence (*Max*(*S<sup>i</sup>* )) and a lowest distance (*Min*(*d<sup>i</sup>* )) between vulnerability classes assigned using DRASTIC and COP.
	- i. In this second step, the weights of parameters are introduced as new variables to compute all the feasible combinations of weights and parameter ranges selected in the previous step. The DRASTIC index is calculated for all the combinations of weights and selected classifications in step 1.
	- ii. The new set of DRASTIC maps are evaluated by means of the objective functions.
	- iii. For each parameter, decision trees are applied again to determine the weight that yields greatest similarity between the DRASTIC and COP maps.


**Table 3.** Classification criteria of objective functions for the decision trees algorithm.

The main objective of decision trees in this study is to identify the ranges and weights for each DRASTIC parameter involved in any combination that yields the maximum spatial coincidence (*S<sup>i</sup>* ) and the minimum distance between vulnerability classes (*d<sup>i</sup>* ) (*S<sup>i</sup>* = 3 or *d<sup>i</sup>* = 1 according to Table 3). We establish these threshold values according to the distribution of results in the first set of combinations. Decision trees reduce the computational cost of the optimization problem and allow the most relevant variables to be identified in the vulnerability assessment in carbonate aquifers.

The CHAID (Chi-square Automatic Interaction Detection) algorithm [62,63] is applied in decision trees, considering the objective functions of the optimization problem (*S<sup>i</sup>* and *d<sup>i</sup>* ) as dependent variables. The proposed ranges and weights of parameters and classifications of DRASTIC are the independent variables, and the chi-squared test of significance is used as the splitting criterion in the CHAID algorithm. Each dataset generated in steps 1 and 2 is partitioned into a training set (70%) and a testing set (the remaining 30%) in order to assess the performance of the models. The goodness-of-fit of the classification is evaluated using the precision index [51,64]:

$$Precision = \frac{\sum\_{n=1}^{n=3} TP\_n}{\sum\_{n=1}^{n=3} (TP\_n + FP\_n)} \tag{5}$$

where:


Decision trees can represent the relationship between variables and output class using specific rules following the pathway from the root node to the terminal node [63,65,66]. A priori, the number of terminal nodes on the tree can determine the number of rules but it may generate a large number of irrelevant pieces of information. Only the most relevant variables (rules with the highest population for each decision tree with precision above 50%) whose terminal nodes in the tree are classified as *S<sup>i</sup>* = 3 or *d<sup>i</sup>* = 1 are selected to generate all the feasible combinations to compute the DRASTIC indices.

### **3. Results**

### *3.1. Optimization of the DRASTIC Method*

### 3.1.1. Ranges Optimization

In the first optimization step, we obtained a total of 11,520 different DRASTIC maps. The spatial coincidence (*S<sup>i</sup>* ) of the vulnerability classes between the DRASTIC and COP maps rose to 61.31% (from 55.75% in "original DRASTIC") and the minimum distance between vulnerability classes (*d<sup>i</sup>* ) fell to 20.85% (doriginal DRASTIC = 24.72%).

The selected ranges (from Table 2) for each parameter extracted from decision trees were as follows:


The ranges proposed in [9] were selected for parameters D and T in the decision tress. For conductivity, the selected classifications (C7 and C8) assign lower rates to conductivity values.

### 3.1.2. Weights Optimization

The selected ranges for parameters and the DRASTIC index were combined, varying the weights of parameters between one and five. Due to a large number of generated combinations, we first considered weights one, three and five. This resulted in 35,700 new DRASTIC maps. The spatial coincidence (*S<sup>i</sup>* ) between vulnerability classes increased to 70.34% and the minimum distance between vulnerability classes (*d<sup>i</sup>* ) fell to 8.45%.

• T\* and T1 for Topography; • C7 and C8 for Conductivity.

3.1.2. Weights Optimization

Decision trees were applied again to determine the optimum weight for each parameter. We aimed to determine if the value of the weight for each parameter provides relevant information in the optimization problem. We found that a weight equal to five did not appear in relevant rules in parameters D and T, whereas a weight equal to one did not appear for parameters R and S. All weights (one, three, and five) were found in rules for parameters A, I and C. Decision trees were applied again to determine the optimum weight for each parameter. We aimed to determine if the value of the weight for each parameter provides relevant information in the optimization problem. We found that a weight equal to five did not appear in relevant rules in parameters D and T, whereas a weight equal to one did not appear for parameters R and S. All weights (one, three, and five) were found in rules for parameters A, I and C.

*Water* **2020**, *12*, x FOR PEER REVIEW 12 of 26

The ranges proposed in [9] were selected for parameters D and T in the decision tress. For

The selected ranges for parameters and the DRASTIC index were combined, varying the weights of parameters between one and five. Due to a large number of generated combinations, we first considered weights one, three and five. This resulted in 35,700 new DRASTIC maps. The spatial

conductivity, the selected classifications (C7 and C8) assign lower rates to conductivity values.

A new set of combinations of ranges and weights were computed to find an optimum between the gaps left due to constraints. We introduced the mid-weights for each parameter, discarding those not found in the relevant rules. We included the following weights for each parameter: W<sup>D</sup> = 2; W<sup>R</sup> = 4; W<sup>A</sup> = 2 and 4; W<sup>S</sup> = 4; W<sup>T</sup> = 2; W<sup>I</sup> = 2 and 4; W<sup>C</sup> = 2 and 4. The total of combinations provide two optimal solutions: A new set of combinations of ranges and weights were computed to find an optimum between the gaps left due to constraints. We introduced the mid-weights for each parameter, discarding those not found in the relevant rules. We included the following weights for each parameter: WD = 2; WR = 4; WA = 2 and 4; WS = 4; WT = 2; WI = 2 and 4; WC = 2 and 4. The total of combinations provide two optimal solutions:

• Optimum of O.F. *Min*(*d<sup>i</sup>* ): *d<sup>i</sup>* = 8.45; *S<sup>i</sup>* = 42.91; • Optimum of O.F. *Min*(*di*): *di* = 8.45; *Si* = 42.91; • Optimum of O.F. *Max*(*Si*): *di* = 13.05; *Si* = 70.34;

between vulnerability classes (*di*) fell to 8.45%.

• Optimum of O.F. *Max*(*S<sup>i</sup>* ): *d<sup>i</sup>* = 13.05; *S<sup>i</sup>* = 70.34;

The second solution was considered the best solution because the gain in *S<sup>i</sup>* was higher than the loss in *d<sup>i</sup>* . Finally, the best solution was refined by adjusting the classification of the DRASTIC index to increase the spatial coincidence with the COP map, obtaining the optimum DRASTIC. The classification of the optimum DRASTIC does not match with the ranges proposed originally in [9]. The objective functions take the following values for the optimum DRASTIC: *S<sup>i</sup>* = 76.75%; *d<sup>i</sup>* = 10.92%. The second solution was considered the best solution because the gain in *Si* was higher than the loss in *di*. Finally, the best solution was refined by adjusting the classification of the DRASTIC index to increase the spatial coincidence with the COP map, obtaining the optimum DRASTIC. The classification of the optimum DRASTIC does not match with the ranges proposed originally in [9]. The objective functions take the following values for the optimum DRASTIC: *Si* = 76.75%; *di* = 10.92%. Figure 5 shows the dot-plot of the all the DRASTIC indices calculated. It reveals the efficacy of

Figure 5 shows the dot-plot of the all the DRASTIC indices calculated. It reveals the efficacy of using decision trees in the methodology to reduce the number of combinations to be tested when seeking the optimal solution. Each set of combinations improves the objective functions. Black dots show those including the mid-weights for each parameter that produce improvement in the objective function "*di*", but not in "*Si*". Red dots indicate the DRASTIC indices that provide *Max*(*S<sup>i</sup>* ) and *Min*(*d<sup>i</sup>* ). using decision trees in the methodology to reduce the number of combinations to be tested when seeking the optimal solution. Each set of combinations improves the objective functions. Black dots show those including the mid-weights for each parameter that produce improvement in the objective function "*di*", but not in "*Si*". Red dots indicate the DRASTIC indices that provide *Max*(*Si*) and *Min*(*di*).

**Figure 5.** Results of objective functions for all DRASTIC combinations.

### **Figure 5.** Results of objective functions for all DRASTIC combinations. *3.2. Analysis of Optimum DRASTIC (O-DRASTIC)*

O-DRASTIC keeps the ranges proposed in [9] for parameters D and T. Only ranges of parameter C change in O-DRASTIC (C8 from Table 2). In general, the value of the conductivity rates is reduced in O-DRASTIC. The spatial distribution of conductivity ranges shows slight changes. The weights (W) of parameters of O-DRASTIC are shown in Table 4, compared with the original DRASTIC weights.



Results of O-DRASTIC reveal that Depth to water and Conductivity have no a significant impact on vulnerability for our case study area. The ranges of parameter C in O-DRASTIC also support this conclusion given that the new classification of this parameter reduces the rate assigned to conductivity values. This conclusion is confirmed by the single-parameter sensitivity analysis of O-DRASTIC, in which the parameter C takes a mean effective weight of 1.7% (compared to its empirical weight 4.3%). In general, in the study area, the carbonate aquifers have a scarcely developed karst. In fact, approximately half the study area shows lower rates in the conductivity map. Although conductivity takes high values in some zones, the ranges of this parameter change in O-DRASTIC with respect to the original DRASTIC. The new classification in O-DRASTIC assigns lower rates to the conductivity values. Therefore, this parameter losses influence, which is also reflected in the weight. Moreover, conductivity can be considered as implicit in the aquifer media parameter (A), which increased in weight in O-DRASTIC in our case study. The reduced weight of parameter D may be due to the large depth to water table over most of the study area. The mean effective weight obtained in the sensitivity analysis was 3.2% (empirical weight = 4.3%) Aquifer media and Soil media gained great significance in the vulnerability assessment, as well as Recharge, albeit by a reduced amount. Impact of vadose zone continues to be an important factor in O-DRASTIC.

The weights in O-DRASTIC are consistent with the concept of COP methodology, where the vulnerability assessment is based mainly on the degree of protection afforded by overlying layers and the way in which the water percolates (recharges) into the aquifer.

The main changes between the original and optimum DRASTIC vulnerability maps mostly occur in areas of "Very low" and "Moderate" vulnerability of the original DRASTIC, since these are the predominant classes in this vulnerability map. Remarkably, those changes do not always occur in one direction. "Moderate" vulnerability shifts towards "Low" or "Very high" vulnerability depending on the zone, whereas other "Moderate" vulnerability areas remain with the same class. In the same way, some "Very low" vulnerability zones jump up a vulnerability class to "Low", while other maintain the same class. We analyzed the differences in the distribution of the parameter rates in these areas.

Figure 6a shows that the main factor causing the vulnerability to change from "Moderate" to "Low" is Depth to water. Since the weight of this parameter changed from five to one in O-DRASTIC, zones with higher rates of D report a greater fall in the vulnerability value, leading the vulnerability class to drop by one level. This graph also shows how some areas with rates of three, five and seven change from "Moderate" to "High" vulnerability, which demonstrates that other parameters influence the "Very high" vulnerability class.

On the other hand, we can observe in Figure 6a that the rate of parameter S (soil media) is equal to eight over more than 80% of the area where vulnerability increased from "Moderate" to "Very high", corresponding with soils with a high organic content, and outcrops of limestone and dolomites. Furthermore, the rate of S is equal to four over nearly 90% of the area where vulnerability fell from "Moderate" to "Low". Zones where no change in vulnerability class was observed have soils with medium values.

A similar conclusion regarding Soil media is drawn in the areas where vulnerability rose from "Very low" to "Low" (Figure 6b). In these zones we can also observe higher rates of parameters A (aquifer media) and I (impact of vadose zone). These zones correspond to karstified areas and highly permeable layers.

The results highlight the influence of Aquifer media, Soil media and Impact of vadose zone in the vulnerability of carbonate aquifers.

These conclusions are supported by the single-parameter sensitivity analysis carried out for O-DRASTIC, which yields higher effective weights in parameters A, S and I and lower effective weights for parameters D and C.

*Water* **2020**, *12*, x FOR PEER REVIEW 14 of 26

**Figure 6.** Distribution of rates of the different parameters within the area where significant changes are observed from (**a**) "Moderate" or from (**b**) "Very low" classes in the original DRASTIC to other vulnerability classes in O-DRASTIC. **Figure 6.** Distribution of rates of the different parameters within the area where significant changes are observed from (**a**) "Moderate" or from (**b**) "Very low" classes in the original DRASTIC to other vulnerability classes in O-DRASTIC.

On the other hand, we can observe in Figure 6a that the rate of parameter S (soil media) is equal to eight over more than 80% of the area where vulnerability increased from "Moderate" to "Very high", corresponding with soils with a high organic content, and outcrops of limestone and dolomites. Furthermore, the rate of S is equal to four over nearly 90% of the area where vulnerability fell from "Moderate" to "Low". Zones where no change in vulnerability class was observed have In general, the new ranges of Conductivity in O-DRASTIC contribute to an overall drop in vulnerability class, but the weights of parameters cause a sharp jump in vulnerability class, especially in the southern part of the basin, where Aquifer media and Soil media have the greatest influence (higher rates in this area).

soils with medium values. A similar conclusion regarding Soil media is drawn in the areas where vulnerability rose from The values of the optimum DRASTIC (O-DRASTIC) vary between 52 and 178 and the optimal vulnerability classes are the following:


Figure 7a shows a better distribution of vulnerability values in O-DRASTIC within the COP vulnerability classes (compared to the original DRASTIC). Close similarities can also be appreciated between the O-DRASTIC and COP vulnerability maps (Figure 7b). The vulnerability classes in for parameters D and C.

(higher rates in this area).

• "Very low": 52–107; • "Low": 107–130; • "Moderate": 130–138; • "High": 138–146; • "Very high": ≥146;

vulnerability classes are the following:

O-DRASTIC overlap with COP over 76.75% of the basin, representing an improvement of 20% relative to the original DRASTIC. This spatial coincidence is particularly improved in the class of "Very high" vulnerability, with a 90% coincidence between COP and O-DRASTIC. DRASTIC overlap with COP over 76.75% of the basin, representing an improvement of 20% relative to the original DRASTIC. This spatial coincidence is particularly improved in the class of "Very high" vulnerability, with a 90% coincidence between COP and O-DRASTIC.

between the O-DRASTIC and COP vulnerability maps (Figure 7b). The vulnerability classes in O-

*Water* **2020**, *12*, x FOR PEER REVIEW 15 of 26

These conclusions are supported by the single-parameter sensitivity analysis carried out for O-DRASTIC, which yields higher effective weights in parameters A, S and I and lower effective weights

In general, the new ranges of Conductivity in O-DRASTIC contribute to an overall drop in vulnerability class, but the weights of parameters cause a sharp jump in vulnerability class, especially in the southern part of the basin, where Aquifer media and Soil media have the greatest influence

The values of the optimum DRASTIC (O-DRASTIC) vary between 52 and 178 and the optimal

Figure 7a shows a better distribution of vulnerability values in O-DRASTIC within the COP

**Figure 7.** Distribution of vulnerability values of original DRASTIC and O-DRASTIC within vulnerability classes in COP (**a**); vulnerability maps for original DRASTIC, O-DRASTIC and COP (**b**). **Figure 7.** Distribution of vulnerability values of original DRASTIC and O-DRASTIC within vulnerability classes in COP (**a**); vulnerability maps for original DRASTIC, O-DRASTIC and COP (**b**).

O-DRASTIC was validated using the pollution index (Equation (1)) and the correlation with nitrate concentration (R<sup>2</sup> = 0.653) in carbonate aquifers improved with respect to the original DRASTIC).

Lastly, we assessed the vulnerability of the three detrital aquifers in the Upper Guadiana Basin by applying original DRASTIC and O-DRASTIC, and we validated the maps following Equation (1). Original DRASTIC showed a good correlation (R<sup>2</sup> = 0.765) in detrital aquifers but O-DRASTIC gave a significantly better correlation (R<sup>2</sup> = 0.862). Both methods (original DRASTIC and O-DRASTIC) perform a better vulnerability assessment in detrital aquifers than in carbonate aquifers.

Figure 8 shows the change that O-DRASTIC produces (compared to original DRASTIC) in terms of distance between vulnerability classes. Negative values mean the vulnerability class drops in O-DRASTIC compared to original DRASTIC, while positive values mean the vulnerability class in O-DRASTIC is higher than in original DRASTIC. Three zones are thus distinguished in the Upper Guadiana Basin: the southern area (Campo de Montiel) where vulnerability increases by one or two classes. This area is characterized by a higher recharge rate and a karstified aquifer. In the mid-basin (including Mancha Occidental I, Mancha Occidental II, Consuegra Villacañas and Rus-Valdelobos), vulnerability decreases generally by one class, although most of the area remains unchanged; in the northern zone (Lillo-Quintanar, Sierra de Altomira and La Obispalía) only small areas jump up a vulnerability class, and by only one level.

DRASTIC).

O-DRASTIC was validated using the pollution index (Equation (1)) and the correlation with nitrate concentration (R2 = 0.653) in carbonate aquifers improved with respect to the original

Lastly, we assessed the vulnerability of the three detrital aquifers in the Upper Guadiana Basin by applying original DRASTIC and O-DRASTIC, and we validated the maps following Equation (1). Original DRASTIC showed a good correlation (R2 = 0.765) in detrital aquifers but O-DRASTIC gave a significantly better correlation (R2 = 0.862). Both methods (original DRASTIC and O-DRASTIC)

Figure 8 shows the change that O-DRASTIC produces (compared to original DRASTIC) in terms of distance between vulnerability classes. Negative values mean the vulnerability class drops in O-DRASTIC compared to original DRASTIC, while positive values mean the vulnerability class in O-DRASTIC is higher than in original DRASTIC. Three zones are thus distinguished in the Upper Guadiana Basin: the southern area (Campo de Montiel) where vulnerability increases by one or two classes. This area is characterized by a higher recharge rate and a karstified aquifer. In the mid-basin (including Mancha Occidental I, Mancha Occidental II, Consuegra Villacañas and Rus-Valdelobos), vulnerability decreases generally by one class, although most of the area remains unchanged; in the northern zone (Lillo-Quintanar, Sierra de Altomira and La Obispalía) only small areas jump up a

perform a better vulnerability assessment in detrital aquifers than in carbonate aquifers.

**Figure 8.** Spatial distribution of change in O-DRASTIC vulnerability classes in comparison to the **Figure 8.** Spatial distribution of change in O-DRASTIC vulnerability classes in comparison to the original DRASTIC.

The largest difference in vulnerability class assigned is for carbonate aquifers. More than 17% of the carbonate aquifers area rises by two classes or more vulnerability classes. Only 7% of the detrital aquifers area jumps by this much, while 62% of the detrital aquifers show no change in the The largest difference in vulnerability class assigned is for carbonate aquifers. More than 17% of the carbonate aquifers area rises by two classes or more vulnerability classes. Only 7% of the detrital aquifers area jumps by this much, while 62% of the detrital aquifers show no change in the vulnerability class.

### **4. Discussion and Conclusions**

vulnerability class.

original DRASTIC.

**4. Discussion and Conclusions**  This paper demonstrates that the DRASTIC method can be adapted to assess vulnerability in carbonate aquifers by undertaking some simple modifications of the weights and ranges of the parameters. Most recent studies to adapt the DRASTIC method to carbonate aquifers have aimed to transform DRASTIC by including and/or removing parameters that take the karst characteristics into account [12,28,29]. Other studies modified the classic assessment according to nitrate concentrations [27,30], though the most highly contaminated groundwater does not always imply higher vulnerability [31,32]. Instead of using either of these previous approaches, our study establishes a This paper demonstrates that the DRASTIC method can be adapted to assess vulnerability in carbonate aquifers by undertaking some simple modifications of the weights and ranges of the parameters.Most recent studies to adapt the DRASTIC method to carbonate aquifers have aimed to transform DRASTIC by including and/or removing parameters that take the karst characteristics into account [12,28,29]. Other studies modified the classic assessment according to nitrate concentrations [27,30], though the most highly contaminated groundwater does not always imply higher vulnerability [31,32]. Instead of using either of these previous approaches, our study establishes a correspondence between DRASTIC and a vulnerability method that was specifically developed for karstic aquifers, the COP method, and therefore we avoid making conceptual changes in the original definition of DRASTIC method.

Our methodology is based on an optimization approach that identifies the ranges and weights of DRASTIC parameters that minimize the differences in the vulnerability assessment compared with the COP method, which was developed specifically for karstic aquifers. It allows us to identify the significance of the different DRASTIC parameters in the vulnerability assessment in karstic aquifers. Decision trees and spatial statistics analyses are combined to identify the ranges and weights of parameters that provide the optimal DRASTIC classes in terms of coincidence and distance with the COP vulnerability map. This optimization approach could be applied by minimizing the differences with respect to any other reference method for assessing vulnerability, whose results has been previously validated and considered as reasonable. Although many data mining techniques have been applied in groundwater vulnerability assessments [39,40,42], only a few studies have used decision trees in vulnerability studies [50,51]. They have been mainly used to assess other groundwater quality problems [65,67,68].

Our proposed method was applied in the Upper Guadiana Basin, an agricultural area that overlies carbonate and detrital aquifers. The socio-economic and hydrogeological particularities of the basin highlight the need to establish unified management measures at basin scale, not only regarding groundwater exploitation but also in terms of protecting the good quality of the groundwater resource. Therefore, harmonization of criteria to assess groundwater vulnerability would allow comparison at basin scale and overcome the issue of dissonant results provided by contrasting vulnerability methods. Our approach assumes that the COP method provides a better approximation of the vulnerability in the case study. This assumption was tested by means of a validation analysis in which we show that the "protection-against-pollution" index (derived from COP and LU data) more closely correlates (R<sup>2</sup> = 0.78) with nitrate concentration than the DRASTIC "pollution index". The results of the validation show that the pollution index derived from the original DRASTIC for carbonate aquifers is not correlated with nitrates (R<sup>2</sup> = 0.04), whereas in detrital aquifers there is a close correlation ((R<sup>2</sup> = 0.76). Other authors have previously pointed out that the original DRASTIC method does not significantly correlate with nitrate concentration in agricultural areas [45,59,69]. In contrast, the pollution index derived from O-DRASTIC shows significant correlation with nitrates for both, carbonate and detrital aquifers (R<sup>2</sup> = 0.65 and R<sup>2</sup> = 0.86, respectively).

The optimal solution (O-DRASTIC) obtained in this optimization problem shows that changing the range of the Conductivity parameter produces a small drop in the vulnerability class compared with the original DRASTIC but does not lead to a significant improvement in the objective functions (coincidence and distance between vulnerability classes). The reduced significance of Conductivity is also confirmed by the much-reduced weight assigned to this parameter in O-DRASTIC. Other DRASTIC adaptation for carbonate aquifers also pointed the reduced significance of Conductivity to assess the potential "protection-against-pollution" in karstic systems [12]. We also observed a reduced significance of the Depth to water table in our case study (WD = 1). Other research studies that modified the weights of DRASTIC parameters [3,69,70] also found the Depth to water parameter to be insignificant. Specifically in karstic aquifers Depth to water is not so relevant in protecting an aquifer from contamination because of the high transit velocity through the vadose zone [12]. Moreover, the low significance of in our case study may be due to fact that the water table lies below 30 m over most of the basin. The same argument was also stated in another case study [3]. The reduced significance of Depth to water and Conductivity is confirmed by a sensitivity analysis. Topography and Impact of vadose zone maintain their weights (W<sup>T</sup> and W<sup>I</sup> ) as defined in the original DRASTIC, while the remaining parameters (Recharge, Aquifer media and Soil media) are given maximum weights (W<sup>R</sup> = W<sup>A</sup> = W<sup>S</sup> = 5). The large increments in the weights given to Aquifer media and Soil media in O-DRASTIC show that they are the most significant factors controlling the vulnerability in karstic aquifers. Other authors concur that these parameters are the most significant [3,70,71]. The principal change in O-DRASTIC with respect to the original DRASTIC is in the weights of parameters, which highlights that these parameters embrace most of the uncertainty in the DRASTIC vulnerability assessment [72].

Our optimal solution provides an improvement of 20% in terms of coincidence between the vulnerability classes assigned by DRASTIC and COP. This improvement was achieved by applying spatial statistical analyses and decision trees, which allowed potential solutions to be obtained by exploring only a 0.1% of the total dimensionality of the defined optimization problem. The proposed method also helps to achieve a better understanding of the parameters and variables of the "equivalent detrital approach" that really influence vulnerability in this karstic system. This optimal solution was tested for the carbonate and detrital aquifers in our case study, but it should also be tested in other different aquifers with similar hydrogeological characteristics in order to prove its applicability in a broader context under different management framework.

In summary, results show that COP and O-DRASTIC report higher vulnerability classes than the original DRASTIC method over 36% of the total area overlying carbonate aquifers. The greatest differences between the original DRASTIC and O-DRASTIC are produced for the carbonate aquifers rather than the detrital aquifers. This confirms that the reliability of DRASTIC vulnerability assessment is significantly better for detrital aquifers than for karst aquifers. These underestimations of vulnerability in karstic aquifers when applying the classic DRASTIC is due to the physical particularities of these aquifers and their greater sensitivity to pollution [1,12,31].

### *Hypotheses, Limitations and Future Works*

We have demonstrated the applicability of the method in the case of carbonate aquifers where the karst is not highly developed. Its applicability to well-developed karst aquifers also needs to be tested. The main adopted assumptions and limitations of the general methodology adopted are:


Future work should be developed to verify the applicability of O-DRASTIC in other case studies, including aquifers with different physical (climatic and hydrogeological settings) and management particularities in order to withdrawal more generalized conclusions.

**Author Contributions:** Conceptualization, L.B.-R. and D.P.-V.; Data curation, L.B.-R.; Investigation, L.B.-R.; Methodology, L.B.-R. and D.P.-V.; Software, L.B.-R.; Supervision, D.P.-V.; Validation, L.B.-R.; Visualization, L.B.-R.; Writing—original draft, L.B.-R. and D.P.-V.; Writing—review and editing, L.B.-R. and D.P.-V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research did not receive any specific grant from funding agencies in the public, commercial, or not-for-profit sectors.

**Acknowledgments:** This paper was partially supported by the SIGLO-AN (RTI2018-101397-B-I00) project from the Spanish Ministry of Science, Innovation and Universities (Programa Estatal de I+D+I orientada a los Retos de la Sociedad), the GeoE.171.008-HOVER and the GeoE.171.008-TACTICproject from GeoERA organization funded by European Union's Horizon 2020 research and innovation program.

**Conflicts of Interest:** The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

### **Appendix A. DRASTIC and COP Methods**

### *Appendix A.1. DRASTIC Method*

DRASTIC method was developed by [9] to assess intrinsic groundwater vulnerability in any type of aquifer.

This method considers that there are seven parameters/variables influencing the vulnerability to contamination: Depth to water table (D), Net recharge (R), Aquifer media (A), Soil media (S), Topography (T), Impact of vadose zone (I) and Hydraulic conductivity (C). A rate of importance is assigned to the parameters according to the value or characteristics of each parameter (Table A1).

The values of the parameters are weighted to obtain the DRASTIC index, which is calculated following Equation (A1):

$$\text{DRAST index} = 5 \times D + 4 \times R + 3 \times A + 2 \times S + 1 \times T + 5 \times I + 3 \times C \tag{A1}$$

The DRASTIC index was originally classified into eight vulnerability levels according to some color codes [9] (Table A2) although it is usually grouped into five vulnerability classes that do not match with the original proposal in [9].




**Table A1.** *Cont.*

**Table A2.** Color codes for the DRASTIC index.


### *Appendix A.2. COP Method*

The COP method was developed by [10] to assess intrinsic groundwater vulnerability in carbonate aquifers.

This method considers the properties of layers overlying the water table (O factor), the concentration of flow (C factor) and precipitation (P factor) as the main parameters influencing groundwater vulnerability in carbonate aquifers. The concept of this method is to assess the natural protection of groundwater determined by the overlying soils and the unsaturated zone, which may be modified by the infiltration process and climatic conditions.

Each factor is divided into subfactors, whose formulation is detailed explained in [10]. The factors are classified in ranges and the COP index calculated as the product of the three factors following the Equation (A2):

$$\mathbf{COP} \text{ index} = \mathbf{C} \times \mathbf{O} \times \mathbf{P} \tag{A2}$$

The ranges of factors and the COP index are shown in Table A3.


**Table A3.** Values for COP parameters and vulnerability classes for the COP index.

### **Appendix B. Data Source and Methodology to Calculate DRASTIC and COP**



**Table A5.** Data source and methodology to calculate COP parameters.


### **References**


**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article*
