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

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.

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,4,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,11,12,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,17,18,19,20,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,23,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,28,29,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,37,38] amongst others. Different data mining algorithms have also been applied to improve DRASTIC performance [39,40,41,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,47,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.

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.

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.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².

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).

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,53,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 × 10⁻⁴ m/s); in the central part there are zones with values higher than 5 × 10⁻⁴ m/s, while conductivity to the southern is mainly in the range 3.5 × 10⁻⁴–5 × 10⁻⁴ m/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.

Table A4. Data source and methodology to calculate DRASTIC parameters.

Factor	Data Source	Methodology
D	Data from simulation flow model (River Basin Authority) (mean of the simulated data for each grid point from 1974 to 2015).	Spatial interpolation using IDW (Inverse Distance Weighted) of groundwater level data and reclassification into D index values.
R	Recharge time series calculated from SACRAMENTO model. Mean recharge value in the period 1974–2015.	Estimation of the mean net recharge taking into account the different hydrology cycle variables in the period 1974–2015.
A	Hydrogeological map of Spain 1:200,000.	Direct assignment of A index values for each hydrogeological unit delimitated.
S	Soil Map of Spain 1:1,000,000.	Direct assignment of S index values for each soil type.
T	Digital Terrain Model at 100 × 100 m cell size.	Calculation of the slope raster file and reclassification of values into T index values.
I	Lithostratigraphic map of Spain 1:200,000.	Direct assignment of I index values for each lithostratigraphic unit.
C	Flow model at 1000 × 1000 m cell size.	Spatial interpolation using IDW of conductivity data and reclassification into C index values.

and

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

Factor	Subfactor	Data Source	Methodology
C	-	-	There are no catchment areas to swallow holes in these aquifers.
	Scenario 2	Vegetation from CORINE LAND COVER and slope from Digital Terrain Model (100 × 100 m cell size) Karstic features from previous research works & fieldwork and litthostratigraphic map 1:200,000.	Carbonate lithologies with low karstification are considered as fissured formations. Limestones and dolomites with high or very high permeability are considered as scarcely developed karst. For sv factor, the vegetation cover is considered high when more than 30% of the surface is covered. Assignment of the values for the karstic features, vegetation and slope according to COP methodology.
O	OS—Soil	Soil Map of Spain 1:1,000,000.	Assignment of the Os values after classify the different types of soil according to the COP methodology.
	OL—Lithology	Lithostratigraphic map of Spain 1:200,000.	Classification of each lithology according to COP methodology and determination of thickness of vadose zone from 3D flow model.
P	PQ—Precipitation quantity	Precipitation data from SPAIN02 [73] in the grid within the Upper Guadiana Basin (mean rainfall taking into account data above 0.5 mm/day).	Reclassification of the precipitation values into the PQ subfactor values, taking into account the average rainfall in the wet years. Precipitation series from SPAIN02 between 1974 and 2015 were used to extract the mean annual precipitation for wet years.
	PI—Temporal distribution	Precipitation data from SPAIN02 (number of rainy days in the grid within the Upper Guadiana Basin).	Counting of the number of rainy days above 0.5 mm for each cell in the SPAIN02 grid. For the estimate of the rainy days per year, meteorological historical series between 1974 and 2015 from SPAIN02 were analysed.

(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.

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]:

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

(1)

$$\mathrm{Protection}-\mathrm{against}-\mathrm{Pollution} \mathrm{index} \left(\mathrm{COP}\right)=\mathrm{COP} \mathrm{index}+\frac{5}{\mathrm{LU}}$$

(2)

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

Table 1. Rates of LU.

LU	Rate
Agro-forestry areas	7
Airports	2
Annual crops associated with permanent crops	8
Broad-leaved forest	2
Burnt areas	2
Complex cultivation patterns	8
Coniferous forest	2
Construction sites	2
Continuous urban fabric	10
Discontinuous urban fabric	10
Dump sites	9
Fruit trees and berry plantations	7
Industrial or commercial units	8
Inland marshes	1
Land principally occupied by agriculture, with significant areas of natural vegetation	5
Mineral extraction sites	3
Mixed forest	3
Natural grasslands	3
Non-irrigated arable land	5
Olive groves	6
Pastures	5
Permanently irrigated land	8
Road and rail networks and associated land	2
Sclerophyllous vegetation	3
Sparsely vegetated areas	3
Sport and leisure facilities	2
Transitional woodland-shrub	2
Vineyards	5
Water bodies	1

.

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² = 0.04) in the carbonate aquifers. However, the protection-against-pollution index linked to the COP map (Figure 3b) shows strong correlation (R² = 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 Figure 2 and Figure 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.

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 non-continuous parameters (A, S and I) were assumed to be invariant and we classify these parameters as proposed in [9]. We modify the DRASTIC index classification (DR) and the ranges of only three numerical parameters in DRASTIC: Depth to water (D), Topography (T) and Conductivity (C). The proposed ranges are shown in

  Table 2. Proposed classifications for the parameters and DRASTIC index.

  Vulnerability	DR1	DR2	DR3	DR4	DR5	DR6	DR7	DR8	DR9	DR10	DR11	DR12	DR13	DR14	DR15	DR16
  Very low	<70	<70	<70	<80	<80	<90	<90	<90	<100	<100	<100	<100	<110	<120	<120	<130
  Low	80	90	100	90	100	100	100	110	110	115	120	130	120	130	130	140
  Moderate	90	110	130	100	120	110	115	130	120	125	140	140	130	140	145	150
  High	100	130	160	110	140	120	130	150	130	140	160	150	140	150	160	160
  Very high	>100	>130	>160	>110	>140	>120	>130	>150	>130	>140	>160	>150	>140	>150	>160	>160
  D (m) rate	D *	D1	D2	D3	D4	D5	D6	D7	D8	D9	D10	D11	D12	D13	D14	D15
  10	1.5	5	5	10	10	13	15	15	20	20	25	25	30	30	30	35
  9	5	10	10	15	20	21	20	30	25	40	30	50	35	55	60	70
  7	10	15	20	20	30	30	25	45	30	60	35	75	40	80	90	105
  5	15	20	40	25	40	37	30	60	35	80	40	100	45	100	120	140
  3	23	25	80	30	50	54	35	75	40	100	45	125	50	127	150	175
  2	30	30	160	35	60	88	40	90	45	120	50	150	55	160	180	210
  1	>30	>30	>160	>35	>60	>88	>40	>90	>45	>120	>50	>150	>55	>160	>180	>210
  T (%) rate	T *	T1	T2	T3	T4
  10	2	1	<3	4	9
  9	6	2	7	8	18
  5	12	4	14	12	27
  3	18	8	22	20	36
  1	>18	>8	>22	>20	>36
  C (m/day)	C *	C1	C2	C3	C4	C5	C6	C7	C8
  1	0.04–4	≤0.1	≤2	≤2	≤3	≤3	≤5	≤6	≤15
  2	4–12	3	5	6	6	12	12	15	30
  4	12–28	6	12	15	15	25	25	30	45
  6	28–40	15	25	30	30	40	40	75	60
  8	40–80	30	40	40	75	80	80	80	75
  10	>80	>30	>40	>40	>75	>80	>80	>80	>75

  * Original classification of DRASTIC parameters.

  . Due to the narrow variability of the data in the study area the recharge ranges were not modified.
• The number of classes and rates adopted for all parameters are as proposed in [9]. We change only the distribution of numerical data within the ranges in order to adapt them to the characteristics of the case study.
• The modification of ranges is based on statistical criteria according to the distribution of data in the study area (quantile method, equal intervals, natural clusters in ArcGis and increasing the limit of the ranges by a constant value) [32,61].
• We established a minimum range amplitude of 5 for the DRASTIC vulnerability classes.
  We also established certain constraints in the weights of parameters:
• Only weights between 1 to 5 are considered;
• The sum of the weights had to be 23 for each combination (as the original proposal in [9]).

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

$$d_i={\displaystyle \sum }{\alpha }_{jk}\times x_{jk}$$

(3)

where

• α_jk is the total area of “j” vulnerability class of COP overlapping with the “k” vulnerability class of DRASTIC;
• x_jk is a weight from 0 to 4 depending on the number of jumps from one vulnerability class to another.

In order to establish a dimensionless threshold to quantify the distance (d_i), it is expressed as a percentage of the maximum calculated distance (d_max) over all calculated DRASTIC indices:

$$d_i(\%)=\frac{d_i}{d_{max}}\times 100$$

(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:

• Ranges optimization:
  • 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].
  • All the DRASTIC indices are evaluated through the objective functions and the results are classified intro three categories (

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

    Objective Function	Value	Class
    Si (spatial coincidence)	<30%	1
    	30–50%	2
    	>50%	3
    di (distance)	<30%	1
    	30–50%	2
    	>50%	3

    ).
  • Decision trees are applied in order to find out the ranges for each parameter that gives the highest coincidence (Max(S_i)) and a lowest distance (Min(d_i)) between vulnerability classes assigned using DRASTIC and COP.
• Weights optimization:
  • 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.
  • The new set of DRASTIC maps are evaluated by means of the objective functions.
  • For each parameter, decision trees are applied again to determine the weight that yields greatest similarity between the DRASTIC and COP maps.

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_i) and the minimum distance between vulnerability classes (d_i) (S_i = 3 or d_i = 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_i and d_i) 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{{{\displaystyle \sum }}_{n=1}^{n=3}T{P}_n}{{{\displaystyle \sum }}_{n=1}^{n=3}\left(T{P}_n+F{P}_n\right)}$$

(5)

where:

• n = number of classes;
• TP_n = number of correctly recognized class examples in the class n;
• FP_n = number of examples incorrectly assigned to the class n;

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_i = 3 or d_i = 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_i) 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_i) fell to 20.85% (d_{original DRASTIC} = 24.72%).

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

• DR11, DR12, DR14, DR15 and DR16 for DRASTIC classification;
• D*, D1, D2, D3 and D4 for Depth to water;
• T* and T1 for Topography;
• C7 and C8 for Conductivity.

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_i) between vulnerability classes increased to 70.34% and the minimum distance between vulnerability classes (d_i) fell to 8.45%.

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.

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_D = 2; W_R = 4; W_A = 2 and 4; W_S = 4; W_T = 2; W_I = 2 and 4; W_C = 2 and 4. The total of combinations provide two optimal solutions:

• Optimum of O.F. Min(d_i): d_i = 8.45; S_i = 42.91;
• Optimum of O.F. Max(S_i): d_i = 13.05; S_i = 70.34;

The second solution was considered the best solution because the gain in S_i was higher than the loss in d_i. 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_i = 76.75%; d_i = 10.92%.

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 “d_i”, but not in “S_i”. Red dots indicate the DRASTIC indices that provide Max(S_i) and Min(d_i).

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. Weights of parameters of original DRASTIC and optimum (O)-DRASTIC.

DRASTIC Parameters	D	R	A	S	T	I	C
W (original DRASTIC)	5	4	3	2	1	5	3
W (O-DRASTIC)	1	5	5	5	1	5	1

, 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.

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).

The values of the optimum DRASTIC (O-DRASTIC) vary between 52 and 178 and the optimal vulnerability classes are the following:

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

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 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.

O-DRASTIC was validated using the pollution index (Equation (1)) and the correlation with nitrate concentration (R² = 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² = 0.765) in detrital aquifers but O-DRASTIC gave a significantly better correlation (R² = 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.

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

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² = 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² = 0.04), whereas in detrital aquifers there is a close correlation ((R² = 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² = 0.65 and R² = 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_T and W_I) as defined in the original DRASTIC, while the remaining parameters (Recharge, Aquifer media and Soil media) are given maximum weights (W_R = W_A = W_S = 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:

• The ranges of categorical non-continuous parameters (Aquifer media, Soil media and Impact of vadose zone) are not modified in this optimization procedure. We consider that the Delphi criteria proposed in [9] can be applied to establish the relative significance of each range with respect to potential pollution.
• Other algorithms and/or techniques (for example, a Random Forest algorithm) could be employed to achieve the goal in a more efficient way.
• We have not studied the whole domain of potential solutions, and a wider spectrum of parameter ranges could be tested to find other optimal solutions. Moreover, the optimization procedure provides local optimum solutions.
• Although decision trees help to reduce the dimensionality of the optimization problem, the methodology involves a large number of calculations, which might handicap extending the method to other case studies.
• The proposed methodology requires a previous validated vulnerability assessment in the study area in order to optimize the DRASTIC method.

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.