First Steps into an Integrated Karst Aquifer Vulnerability Approach (IKAV). Intrinsic Groundwater Vulnerability Analysis of the Yucatan Karst, Mexico

Abstract

Karst groundwater vulnerability maps are important tools for the development of groundwater management and protection strategies. However, current methodologies do not always match regional characteristics and parameter adaptations are necessary. In addition, other important processes such as dilution and aquifer residence time are not included in vulnerability analysis for the complications of evaluating two or more criteria simultaneously. The integrated karst aquifer vulnerability approach (IKAV) project aims to develop an integrated approach to include these parameters and estimate global change implications in current and future scenarios. As a first step, intrinsic vulnerability methodologies are studied in order to highlight important parameters and the congruence with regional characteristics of the Yucatan karst. Results demonstrate agreement between methods for the evaluation of high and very high vulnerabilities and their relation with fissures and dolines. Moderate vulnerabilities are assigned to more than 50% of the area. However, moderate vulnerabilities, assigned to the coastal area and the Southern hill, are highly questionable. Intrinsic features affecting moderate classes vary according to the method. Parameter sensitivity analysis and overlap analysis demonstrate the influence of depth to the unsaturated zone, soils, precipitation, and slope on moderate values. Therefore, such parameters must be re-evaluated and discretized according to the characteristics of the study area to match Yucatan regional characteristics.

1. Introduction

Karst groundwater resources are at constant pollution risk due to their own natural features like dolines, fissures, and conduits that can allow a fast infiltration of pollutants into the aquifer. These fast-flowing pathways increase groundwater exposure to pollution since a contaminant would reach groundwater faster, experiencing little to no degradation process in comparison with unconsolidated aquifers. As karst aquifers are the source of water supply for almost 25% of the world population, as estimated by Ford and Williams [1], their protection is urgent, principally in areas where groundwater is the only freshwater source. From the necessity to evaluate and characterize how endangered a given aquifer is, either consolidated or not, the vulnerability concept emerged as an essential part of groundwater management strategies. Vulnerability refers to the aquifer’s susceptibility to contamination by human activities [2].

To date, several methods to estimate either source (spring or well) or resource (groundwater) vulnerability in karst areas have been proposed and tested on different areas of the world. An accurate overview of karst vulnerability methods developed during the last three decades can be found in the work of Iván and Mádl-Szònyi [3]; nevertheless, more literature introducing regional adaptations or new methods have also been published, but in languages other than English [4,5].

The majority of current methodologies use a multi-parameter approach to represent the variables affecting a theoretical pollutant’s travel time, being discretized using score intervals according to a relative degree of protection [6]. In general, development of groundwater vulnerability maps can be seen as the evaluation of multiple layers representing geology, pedology, meteorology, hydrology, and others for their influence on the vertical travel time from surface to the water table. Nevertheless, among the proposed methodologies, parameters are categorized and evaluated with different scores and weights, displaying contrasting outcomes when several methods are applied over the same area [7,8,9,10].

To work out vulnerability maps, it is necessary to deal first with some factors that complicate the process, such as data availability and the determination of an adequate method to be used for a given karst area. Some vulnerability methods need specific data, which in some cases is sparse or even null; for example, Malík and Svasta propose point classifications based on spring discharge data to estimate the degree of karstification as part of the REKS method [11]. This data is not easy to obtain, such as in the case of coastal karst areas with underwater springs, like the Yucatan karst aquifer, thus not allowing the usage of such approach. Current methodologies have been validated via field tests, mostly tracer tests, on their respective areas of development; nevertheless, high subjectivity and personal interpretation regarding intrinsic features behavior also influence results. For example, the EPIK method [12] considers epikarst as a feature which increase vulnerability due to probable infiltration though vertical shafts; this interpretation is opposite to that presented in the PaPRIKa method, where epikarst behaves as a protective feature for its hydrogeological function as perched aquifer; therefore, delaying infiltration [13]. Despite the fact that vulnerability maps can be interpreted by stakeholders without a broad hydrogeological knowledge, outcomes must be taken as theoretical approximations and further analysis is always encouraged.

Global change elements (anthropogenic practices, extreme events) are important factors to be included into a vulnerability analysis. If we consider precipitation as an external stressor, tropical areas benefit from excessive rainfall events since substantial water volumes, infiltrating the subsurface, would promote dilution [14]. Methodologies like the Slovene approach take into account precipitation analyzed in terms of volume, intensity and extreme events [15]. Other methodologies like KARSTIC and DRISTPi consider recharge as a parameter derived from precipitation regimes [16,17]. Nevertheless, vulnerability assessment in tropical areas, with pronounced seasonal precipitation patterns, can be over/under estimated since precipitation ratings from the before mentioned methods are based on year averages. Additionally, dilution potential will depend on pollutant volumes which can be estimated in current scenarios (untreated wastewater, landfills leaching, etc.), while its residence time in the aquifer would be shorter depending on the degree of karstification. Despite effects of global change and its influence on real vulnerability scenarios, there is no method which evaluates such influence together in a vulnerability analysis. Inclusion of dilution and pollutant residence time must be considered carefully in a vulnerability analysis, as the evaluation of two or more independent criteria can lead to ambiguous situations [18].

Combination of travel time, aquifer residence time, and pollutant concentration would resolve into a new, integrated groundwater vulnerability model, increasing the confidence when considering influence of global change factors on regional vulnerability. The IKAV project, affiliated to the INOWAS group (research unit under the Department of Hydrosciences of the Faculty of Environmental Sciences at Technische Universität Dresden, Germany), aims to develop such an integrated methodology to deal with and evaluate current and future vulnerability scenarios.

As a first step, analysis of selected intrinsic vulnerability methodologies, their congruence with regional features, and parameter sensitivity analysis are contemplated to reach an adaptive karst groundwater vulnerability methodology (Figure 1). According to data availability, a total of eight methodologies were chosen and applied on the Yucatan State, Mexico. This paper presents the analysis of four of them (DRISTPi, KARSTIC, RISKE and the Slovene approach), the other being presented in a previous work [10].

Outcomes from this work will serve as a basis to highlight critical parameters, weights adaption, and further considerations to be implemented in an adaptive karst vulnerability approach. Due to its unique characteristics and the increasing aquifer contamination, the chosen study area is Yucatan, Mexico. This work aims to highlight advantages, disadvantages, and the relationship between regional features and vulnerability rates from the four methods here applied. Additionally, this work extends the study and analyzes the outcomes and spatial match of vulnerability classes between the four methods, including a map layer-removal sensitivity analysis with “final vulnerability class” sensitivity, based on vulnerability class changes.

2. Materials and Methods

2.1. Study Area

The Yucatan Peninsula is a trans-boundary limestone platform which embraces areas of Mexico, Belize, and Guatemala (Figure 2). In general, the aquifer is considered as a well-developed karst with large conduit systems at variable scale ranges [19]. The study area, the Yucatan State (39,524 km²), is settled north of the Peninsula; some of the main features in the area are the nearly flat topography in most of the region, lack of surface flow, and “Cenotes” alignment (dolines exposing water).

Dolines are dispersed in the area and show a characteristic semi-circular alignment of about 180 km of diameter at northwest [21]; this alignment, named the “Cenote ring”, is considered as the surface expression of an asteroid impact at the end of the Cretaceous [21,22,23]. The northeastern Yucatan also displays a high doline density [22,24,25]; nonetheless current doline maps, based on contours maps at 1:50,000 scale, underestimate a high percentage of such features (Figure 3a). This underestimation was demonstrated by the application of a new methodology to estimate dolines from digital elevation models [26].

In the southern Yucatan, terrain elevates up to 204 masl in a hilly area named “Sierrita de Ticul”; nevertheless, most of the area is considered flat with elevation gradually increasing from 0 at the coastline to approximately 30 m along the base of the hill [27]. Combination of flatness and karstic development do not allow surface streams to generate, infiltration is; therefore, diffuse at regional scale [28]. The flat topography is reflected on slope, which demonstrates a very low change in elevation for most of the area with exception of the Sierrita de Ticul hill (Figure 3b). Yucatan experiences wet (May–October) and dry seasons (November–April) with a yearly average precipitation fluctuating from 550 mm/y northwest to 1500 mm/y southeast [29]. Recharge has been estimated between 14% and 20% of the mean annual precipitation, meaning high evapotranspiration rates [30,31]. For this study, precipitation values were obtained from 62 stations and records covering a 25-year period were used (from 1990 to 2015) to estimate mean annual precipitation volumes (Figure 3c) [32]. Yucatan also experiences natural phenomena like tropical storms and hurricanes, events that generate extreme precipitation and increase the recharge [33].

Hydraulic gradients are low ranging from 7 to 10 mm/km, suggesting high hydraulic conductivity [34]; this is also reflected on water table depth, where the hill area reflects the deeper aquifer in comparison with the flat plain (Figure 3d). Groundwater flows in a northwest direction [30,35], discharging via submarine springs located in coastal zones [23], with an estimated higher discharge during the wet season [22]. Hydraulic conductivities have been estimated for the Sierrita de Ticul upland and the rest of the area, with values of 5.5 × 10⁻³ and 1.11 m/s, respectively [36,37]. However, since this study focuses on intrinsic resource vulnerability, the hydraulic conductivity parameter was not utilized.

Soils distribution is displayed in terms of Edaphology classification. Rendzinas, Lithosols, and chromic Luvisols predominate in approximately 46%, 29%, and 10% of the area according to public datasets [38]. Point sampling data from INEGI (National Institute of Statistics and Geography) shows agreement regarding Lithosols and Luvisols, with less than 7 cm of depth [39,40]. Medium soil textures seem to predominate in central areas, while coarser soils are found in the coastal rim (Figure 3e).

Due to regional anthropogenic practices, mostly related with untreated wastewater disposal, and groundwater as the only source of water supply on the region, aquifer protection becomes critical with an urgent establishment of protection strategies to maintain water quality, making Yucatan an important area for aquifer vulnerability studies. Currently there are no protection strategies in Yucatan to deal with groundwater contamination; however, some areas of Yucatan have been declared as protected under the title of “hydrogeological reserve” [41,42].

2.2. Methods and Data

Four rating methods to estimate water resource vulnerability were applied in Yucatan. Outcomes from DRISTPi, KARSTIC, RISKE, and the Slovene approach are analyzed to highlight intrinsic features leading to different vulnerability classes and to study the congruence of results with regional characteristics of the Yucatan karst. High subjectivity; however, is expected from the assignation of some features not clearly defined by some methods.

Datasets to develop the necessary map layers were gathered from different public sources and literature review (

Table 1. Data utilized to create map layers for the four methods presented here and those presented previously by Moreno-Gómez [10].

Map Layers	DRISTPi	KARSTIC	RISKE	SA	PaPRIKa	COP	Pi	EPIK	Data
Soil thickness			X	X	X	X	X	X	Borehole data a
Soil texture	X	X	X	X	X	X	X		Borehole data a
Lithology	X	X	X	X	X	X	X		Lithology map a
Fracturing	X	X	X	X	X	X	X		Fracturing maps a
Unsaturated zone (Depth)	X	X		X		X	X		Digital elevation model b/groundwater levels a
Epikarst			X		X		X	X	Sinkholes density map
Confined conditions				X		X			Literature review
Slope	X	X	X	X	X	X	X	X	DEM b
Vegetation				X		X	X	X	Vegetation maps a
Karstification		X	X	X	X	X	X	X	Sinkhole/fissure density maps
Rainfall volume				X		X			25 years historic data c
Rainfall intensity				X		X			25 years historic data c
Recharge	X	X					X		Precipitation c, APLIS d method
Surface features							X		Sinkholes and fracturing maps a
Effective field capacity							X		Settled with minimum values
Hydraulic conductivity (soil)							X		Soil maps a, Saxton equations
Hydraulic conductivity (aquifer)		X							Na
Rock reservoir		X	X		X				Na

Note: ^a Data publicly available at http://www.beta.inegi.org.mx (Vector maps at 1:50,000 scale); ^b digital elevation model (ASTER GDEM, 30 m resolution) from https://www.earthexplorer.usgs.gov; ^c data publicly available at http://clicom-mex.cicese.mx/; ^d multi-parameter method to estimate recharge based on altitude (A), precipitation (P), lithology (L), infiltration (I) and soil (S).

). Some parameters were neglected (see KARSTIC and RISKE methods section) since they are directly related to aquifer lateral flow, hence source vulnerability.

2.2.1. DRISTPi Method

This method is based on the well-known DRASTIC method [43] with addition of one scenario to evaluate effects of non-karstic areas. This method also includes a new parameter to evaluate preferential infiltration (Pi), which can be specific for each scenario. This multi-parameter model defines vulnerability rates from given values and weights according to:

V_I = D(2) + R(4) + I(5) + S(2) + T + PI(5),

(1)

where V_I—vulnerability index; D—depth to water table; R—recharge; I—lithology of the vadose zone; S—soil; T—topography (slope); Pi—preferential infiltration. The numbers in brackets represent the assigned weights. DRISTPi allows the re-assignment and modification of weights. More detailed information about this method can be found in the work of Jiménez-Madrid et al. [17].

2.2.2. KARSTIC Method (KRSTI)

KARSTIC is also an adaptation for karst of the DRASTIC method; it is an extensive vulnerability index since it evaluates nine parameters, four of them combined to create a complex variable (K and I).

V_I = K × F(10) + A(3) + R(4) + S(2) + T + I × D(20) + C(3),

(2)

where V_I—vulnerability index; K—karst development; F—fracturing degree; A—aquifer medium; R—recharge; S—soil; T—topography (slope); I—lithology of the unsaturated zone; D—depth to water table; C—hydraulic conductivity of the aquifer. For the purposes of this research, A and C parameters (aquifer media and aquifer hydraulic conductivity) where neglected since they aim to estimate source vulnerability; this modification was followed by the percentage redistribution of the original weights from the original formula to obtain:

V_I = K × F(11) + R(5) + S(3) + T(2) + I × D(22).

(3)

Applications and study cases related with this method can be consulted in the work of Davis et al. [16].

2.2.3. RISKE Method (ISKE)

RISKE is also an index-based vulnerability mapping method inspired by the Swiss EPIK method. It aims at characterizing the intrinsic vulnerability of the karstic system to vertical infiltration on a catchment scale. The method considers the most prominent intrinsic characteristics of karstic aquifers that condition the behavior of any pollutant entering the system by considering five criteria that describe both the structure of the aquifer and its functionality. Vulnerability rate is assessed by:

V_I = R(0.1) + I(0.4) + S(0.1) + K(0.2) + E(0.2),

(4)

where V_I—vulnerability index; R—aquifer rock; I—infiltration conditions (slope); S—soil; K—degree of karstification; E—epikarst. In (4), the R parameter is related to water circulation in the phreatic zone; therefore, it was excluded and weights were distributed in percentage to finally get:

V_I = I(0.45) + S(0.11) + K(0.22) + E(0.22).

(5)

A precise summary of this method and its development is found in the work of Petelet-Giraud et al. [44].

2.2.4. The Slovene Approach

The Slovene Approach (SA) is based on the Spanish COP method, in accordance with the COST Action 620 conceptual framework for assessment of vulnerability. This methodology comprises three elements; in addition to intrinsic vulnerability, hazards and contamination risk maps are also evaluated. In this study, just the intrinsic vulnerability section is applied. The SA encompasses three basic factors (and sub-factors) that condition water infiltration and the possible transport of contaminants from the surface to the saturated zone. Another particularity of the method is the consideration of the temporal hydrologic variability which is prominent in contaminant transport processes. Calculation of vulnerability rates follows:

V_I = (os + ol × cn) × (sv × sf) × (rd × se),

(6)

where V_I—vulnerability index; os—soil texture and thickness; ol—lithology and fracturing; cn—confinement; sv—slope and vegetation; sf—surface features; rd—rainy days; se—extreme events. This is the equation for areas which do not generate surface flow towards a swallow hole. Parameters inside brackets correspond to overlying layers, concentration of flow, and precipitation factors, respectively, (O, C, P). A detailed description about the parameters and assigned rates are presented by Ravbar and Goldscheider [15].

The methods presented here were applied and validated in karst areas with distinctive characteristics. Final vulnerability ranges vary according to the method, while vulnerability classes are related to specific ranges. Notice the inverse range of the SA method since it evaluates the protection degree, an antonym for vulnerability (

Table 2. Vulnerability classes and index variation according different methods.

Vulnerability Description	Vulnerability Class	Original Test Site
		Sierra de Cañete, Spain and Neblón River Basin, Belgium	Rapid Creek Basin, South Dakota, USA	Fontanilles and Cent-Fonts Karstic Aquifers, France	Podstenjšek Springs Catchment, Slovenia
		DRISTPi Index	KARSTIC Index	RISKE Index	SA Index
Very low	1	17–50	0–200	0–0.8	4–15
Low	2	50–80	200–350	>0.8–1.6	2–4
Moderate	3	80–110	350–600	>1.6–2.4	1–2
High	4	110–140	600–800	>2.4–3.2	0.5–1
Very high	5	>140	>800	>3.2–4	0–0.5

).

For dataset analysis, management, and development of the necessary layers to obtain vulnerability maps, Arc Catalogue and Arc Map version 10.5 were used. Vector data was transformed into raster with a resolution of 900 m² in agreement with the digital elevation model (DEM) used in this study; maps were projected with the World Geodetic System 1984. Aiming to evaluate the influence of single features (e.g., soil texture, karst features, etc.) on final vulnerability classes, vulnerability maps were overlapped against each feature (I, S, K, and E, as in the RISKE case) using the tabulate function of ArcGIS. Hereinafter, vulnerability will be described as VLV for very low; LV for low; MV for moderate; HV for high; and VHV for very high.

2.3. Sensitivity Analysis

The sensitivity analysis was carried out by computing Yucatan unique condition sub-areas to avoid processing large numbers of map grids. A unique condition sub-area (UCS) represents a unique combination of parameters having an extension from one to several map grids. For comparison, theoretical UCSs were also developed for each one of the methods using all rates from each parameter obtaining all possible combinations. Map or parameter sensitivity is obtained from the removal of one layer at time [45], which determines how each parameter affects the vulnerability index without the influence of the weighting factor [46]; It is expressed as:

S_xi = |V_I/N − V_xi/N − 1|,

(7)

where S_xi is the sensitivity associated with the removal of parameter X; V_I the original vulnerability index corresponding to the ith UCS; V_xi the vulnerability index of the ith UCS computed after removing parameter X; N correspond to the number of maps used to estimate vulnerability.

The variation index (effective weight) is also used to estimate index change in terms of percentage of the original index; this is calculated for each grid map (or UCS) by:

V_xi = ((V_I − V_xi)/V_I) × 100

(8)

where V_xi is the variation index of parameter X for a given UCS; V_I and V_xi are the original and perturbed vulnerability index, respectively, of the ith UCS. The variation index will explain the effective relative importance of each parameter in the study area. The equation is presented according to Marín [9] and is identical to the effective weight presented by Napolitano [47]. The previous equations assist in analyzing the index change after removal of one layer; however, a high sensitivity does not always reflect a change of vulnerability class which is dependent on index ranges. In order to estimate the effect of each parameter on final vulnerability classes two equations based on areas are introduced as:

Δ_C = Cl − Cl_xi,

(9)

and

C_S = Ʃ ((A_ΔC)/A_T) × 100   if Δ_C > 0,

(10)

where Δ_C is the vulnerability class change for a given UCS after removal of one map layer; Cl is the original vulnerability class and Cl_xi is the class assigned after removal of one map; Cs represents the area, in percentage, experiencing a vulnerability class change; A_ΔC is the area of grids which experiences a change in vulnerability class and A_T is the summed area of all UCSs.

3. Results

Vulnerability and Parameter Sensitivity

The DRISTPi method (Figure 4a) displays MV covering 59% of the study area followed by HV and VHV at 33% and 8%, respectively. VHV is associated with areas of high sinkhole density where infiltration is expected to be high, while HV covers the area of considerable fracturing; the rest of the area is denoted with MV without further discretization. From visualization, it seems that no other features but those representing preferential infiltration conditions are leading to high vulnerability outcomes. This is not reflected from sensitivity analysis since, for the Yucatan case, Pi display lower values in comparison with the all possible UCS scenarios. According to

Table 3. Average sensitivity values for the Yucatan case and all possible UCS scenarios for the DRISTPi method.

Sensitivity	D	R	I	S	T	Pi
Yucatan Sxi	1.1	2.7	3.8	2.6	2.4	1.6
All UCS Sxi	1.6	2.1	2.8	1.4	2.4	2.6
Yucatan Vxi	12.1	29.7	34.0	5.0	5.6	13.6
All UCS Vxi	10.3	22.1	28.0	11.6	5.5	22.5
Yucatan Cs	45.2	72.5	100.0	9.3	19.2	24.9
All UCS Cs	34.7	70.7	83.2	38.8	18.4	66.2

Note: Yucatan UCSs = 920; all UCSs = 123,200.

, the most influential layer in Yucatan is related lithology (I). Removal of this layer affects the most final outcomes according to S_xi, V_xi and C_s. Lithology in the area, being characterized solely as limestone, explains the influence on MV for most of the area without further discretization.

The KRSTI method displays vulnerabilities ranging from VHV to LV. As shown in Figure 4b, MV is assigned to 57% of the area while VHV, HV, and LV cover 18%, 12%, and 14%, respectively. LV is displayed on the hill area were depth of the unsaturated zone is above 40 m. This reflects the impact of the unsaturated zone (lithology and thickness) due to a high weight assigned to this parameter. Being the area characterized solely as limestone, depth to groundwater is significant for this method. Sensitivity analysis in

Table 4. Parameters sensitivity of the KARSTIC method without the A and C layers.

Sensitivity	K	R	S	T	I
Yucatan Sxi	4.1	2.1	6.2	5.5	14.0
UCS Sxi	28.7	13.5	15.4	16.9	25.7
Yucatan Vxi	12.3	22.1	5.7	7.3	52.6
UCS Vxi	37.1	9.4	6.3	4.0	40.1
Yucatan Cs	23.0	24.6	0.0	0.1	25.0
UCS Cs	55.9	12.3	7.7	4.5	57.7

Note: Yucatan UCSs = 444; all UCSs = 118,800.

supports the outcomes, since removal of I map demonstrates the highest influence of this layer for both the Yucatan case and all possible combination scenarios, with exception of the UCS Sxi.

The ISKE method displays four vulnerability classes. VHV are not visible in Figure 4c due to a low 4% of coverage. Areas with high sinkhole and fracturing density are denoted as HV, covering 31%. This method presents the highest MV area covering 63% of Yucatan. Just 1% of LV is displayed in the hill area, matching grids of high slope where vertical infiltration could be diminished. The I parameter, or infiltration (for this method it refers to slope percentage), is the most influential according sensitivity (

Table 5. Average sensitivity values for the Yucatan case and all possible combination scenarios for the ISKE method without the R parameter.

Sensitivity	I	S	K	E
Yucatan Sxi	0.17	0.10	0.08	0.08
UCS Sxi	0.12	0.08	0.05	0.05
Yucatan Vxi	39.6	37.6	23.3	24.1
UCS Vxi	37.1	35.1	26.2	23.2
Yucatan Cs	68.0	25.3	45.0	44.8
UCS Cs	98.8	14.6	15.7	14.4

Note: Yucatan UCSs = 62; all UCSs = 625.

). The fact that Yucatan is mostly a flat plain explains the area percentage being characterized with MV. Table 5 shows also a high variation index derived from the removal of the soil layer (S); however, the final map does not reflect such influence.

The SA is the only method displaying VLV areas; however, in a low percentage. Since Yucatan is not being analyzed as swallow hole catchment area, locations with a considerable density of karst surface features (dolines, fissures) are categorized as HV, covering 16% of the area. This method also displays a high percentage of MV with 62%, including areas with both shallow water tables at the coast and deeper saturated zones at the Sierrita de Ticul uplands. The 19% classified as LV correspond to fine texture soils with thickness above 30 cm. The influence of the Os parameter is clear from visualization as a consequence of the homogeneity displayed by the Rd and Se parameters (P factor). Average sensitivity for the Yucatan case supports the analysis previously mentioned. Yucatan S_xi, V_xi, and C_s demonstrate the high influence of the Os layer. Nevertheless, comparison of sensitivity from the Yucatan case against all the possible UCSs does not show correspondence as presented by previous methods (

Table 6. Average sensitivity values for the SA.

Sensitivity	Os	Ol	Sv	Sf	Rd	Se
Yucatan Sxi	0.29	0.08	0.12	0.23	0.15	0.09
UCS Sxi	0.13	0.18	0.17	0.19	0.16	0.17
Yucatan Vxi	65.98	34.02	5.56	38.10	11.11	0.00
UCS Vxi	31.86	68.14	30.47	52.78	25.69	30.56
Yucatan Cs	100	97.7	0	48.7	3.01	0
UCS Cs	99.9	99.9	99.6	96.5	91.3	66.7

Note: Yucatan UCSs = 30; all UCSs = 207,346.

).

A resume of vulnerability percentages is displayed in Figure 5. Despite MV is the most remarkable class, covering most of the area, it does not match the Yucatan regional characteristics. Display of MV on coastal areas, with shallow water tables (less than 1 m of depth) and sands as soil cover, can be considered inaccurate from a theoretical basis. Results for MV are triggered by the approaches themselves on how to define HV and VHV classes, making swallow holes catchment areas the most vulnerable to pollution, automatically decreasing the vulnerability for the rest of the area. This leads to a clear misclassification, making necessary the adjustment of parameters, index ranges, and weights to deal with areas where solely diffuse infiltration conditions exist. In order to establish a regionalization, it is necessary to highlight the parameters affecting MV by each method and to define further classifications.

To go further into the analysis, overlap between final vulnerability maps and single features was carried out in order to estimate the spatial match with vulnerability classes (

Table 7. Important parameters leading specific vulnerability classes derived from spatial match with final maps.

Method	Vulnerability Class
	VHV	HV	MV	LV	VLV
DRISTPi	Fissuring	Fissuring	Lithology	Not displayed	Not displayed
KRSTI	Sinkhole density and fissures	Sinkhole density and fissures	Shallow water table	Deep water table	Not displayed
ISKE	Negligible	Sinkhole density and fissures	Low slope percentage	High slope percentage	Not displayed
SA	Not displayed	Sinkhole density and fissures	Soil thickness and texture	Soil thickness and texture	Soil thickness and texture

). Results from the overlap analyses reiterate outcomes from sensitivity analysis, giving a base to develop a new or adapted intrinsic vulnerability method to match the characteristics of Yucatan.

The influence of sinkholes, in terms of karst feature density, and fissuring leading VHV and HV is clear. However, MV can be further discretized to match the regional characteristics. Having multiple factors affecting MV, considerations to match regional features are further evaluated.

4. Discussion

Each methodology evaluated here display, to some extent, agreement with Yucatan characteristics (low slopes, shallow water table, and a regional diffuse infiltration). Nevertheless, they also show areas where the assigned vulnerability can be arguable. This questionable regional classification depends on the method and the values given to specific parameters, which are also dependent on the characteristic of the area where the original methods were developed and tested. The most significant vulnerability class in Yucatan according to this study is MV; this shows agreement with results obtained in previous studies for the EPIK, PI, and COP methods [10]. Moderate values in the area could be explained by different factors, like a regional diffuse infiltration, lithology mostly composed by limestone, low slopes, and a shallow water table in most of the region. These intrinsic features are the most important according to sensitivity analysis and the overlap of final maps with individual parameters. The four methods agree that high sinkhole density areas and extensive fissuring indicate either VHV or HV, of course varying depending on the vulnerability rates given for each one of them. Having this as a base to adapt or generate a new intrinsic vulnerability map, analysis of MV and their relation with regional characteristics is necessary to further classify the area. To adapt or create a new intrinsic groundwater vulnerability approach, considerations presented further could help to develop a new adaptive methodology.

The first process for karst groundwater vulnerability evaluation must be the separation between areas promoting point infiltration at the surface (swallow holes and the catchment feeding it) and the areas with diffuse surface infiltration. This is highly recommended since categorizing diffuse areas as “rest of the area” will inherently decrease vulnerability classes for those areas without surface point infiltration. This will lead to a different evaluation of soils, slope, and other features affecting runoff generation (vegetation, precipitation intensity). For example, a clay-rich soil and high slopes are more likely to promote runoff, hence will increase vulnerability towards a swallow hole if this condition exists (Figure 6a). However, having the same intrinsic characteristics in a diffuse infiltration area leads to an inverse analysis, since the scenario will theoretically minimize the possibility of vertical infiltration; therefore, vulnerability would decrease for the given grid map (Figure 6b). This evaluation is not considered by multiple methodologies since they mostly focus on slope as the sole feature promoting surface runoff.

Estimation of MV in coastal areas, as shown by the four final maps, could be considered inaccurate. Despite considering a regional diffuse infiltration, which is clearly less vulnerable than surface-focused infiltration, the shallow water table in the coastal rim could be reached by pollution faster in comparison with southern areas, such as the Sierrita de Ticul hill. Therefore, a method must contemplate differences in unsaturated zone thickness based on these ranges. To achieve this classification, it is necessary to increase the importance (or weight) of the unsaturated zone thickness and re-classify rates according to regional groundwater table intervals. For example, a Jenks classification could be used to reduce the variance between values contained in a given interval and to increase the variance between different intervals, hence defining regional ranges statistically. The KARSTIC method presents a high weight for this parameter, classifying the hill area as LV; however, shallow water table regions remain with MV.

The area being mostly a flat terrain, categorized with a regional diffuse infiltration, an approach to classify topography (slope) must consider solely the influence on vertical infiltration. For example, the SA is the only method presented here where a high slope means an increase in vulnerability contemplating its influence on runoff generation. Not having significant surface streams and without surface point infiltration in Yucatan, high slopes must be considered as promoting low vulnerability, since vertical infiltration, at those given map grids, would be decreased. For this case, the RISKE method represents such behavior since areas with a strong change in elevation display LV, mainly in the hill area at South of Yucatan.

The Yucatan lithology is characterized as limestone without further discretization. Despite the Yucatan limestone corresponding to different ages, as displayed in Figure 2, the majority of vulnerability methodologies consider solely the lithological material and its karstic development for a vulnerability evaluation. Most of the methods analyzed in this work consider several lithological materials (marls, breccia, shale, sandstone, etc.) for their influence on infiltration capacities (see DRISTPi, KARSTIC, and SA methods). However, when applying these methodologies in Yucatan, lithological materials analyzed by the methods but inexistent in the study area are also evaluated, thus affecting the final outcomes. To better adapt lithology influence on regional vulnerability, this could be coupled with fissuring degree to discretize areas from low to high possible infiltration. This could be similar to the PI method evaluation, where lithology and fracturing are directly related.

When methodologies focus on point infiltration, the parameters affecting runoff towards a swallow hole, generally the slope, are the most influential. For diffuse infiltration areas, soils are the first layer affecting pollutant transport and must be evaluated accordingly. The influence of soils on vulnerability is generally analyzed according to its texture and thickness; however, low weights are assigned to this parameter in most of the methods. Despite soils in Yucatan being considered as absent or for general purposes, due to their shallowness and karst outcrops, they must be given a more important role when diffuse infiltration areas are being studied. For example, the SA method displays the importance of soils since they are the most influential parameter to define MV and LV areas. This is supported by the sensitivity analysis as a consequence of the regional homogeneity of the P and C layers in the Yucatan case. If we analyze the hill area in Figure 4d, we notice that LV and VLV match fine texture soils with a considerable thickness.

Regarding precipitation as an external stressor, or medium for pollutant transport, it is important to consider the fact of Yucatan being a tropical area. Either as recharge or precipitation volumes, this parameter must be also re-classified according to regional precipitation quantities. In this case, effect of precipitation on vulnerability must be contemplated regarding its seasonal influence due to the well-defined precipitation regimes and extreme events period. Additionally, this parameter must focus solely on terms of pollutant travel time and not for dilution potential (see COP method).

5. Conclusions

In this work four methodologies for intrinsic groundwater resources vulnerability mapping where applied in the Yucatan state, northern part of the Yucatan Peninsula, as a first step in the development of an integrated karst groundwater vulnerability method. Results show a predominance of MV for more than 50% of the study area according to the four methods. After overlap of final vulnerability maps with single parameters and performance of parameter removal sensitivity analysis, important considerations were highlighted to develop an adaptive methodology to match regional hydrogeological characteristics. The four methods showed agreement classifying areas as VHV and HV in correspondence to high dolines density areas and strong fissuring. However, coastal areas with shallow water tables (less than 1 m) and sands as soil cover were classified as MV. This classification is considered inaccurate for the region. Similarly, the hill area in the south was also considered as MV despite having an unsaturated zone reaching more than 140 m depth. Regional classification and re-assignment of weights for parameters like depth to ground water, slope, precipitation, and soils are suggested to improve the vulnerability mapping of the region. However, the process to reach this goal requires the development of an adaptive methodology to asses a general scheme based on data availability, regional features classification, and parameters weight adjustment.