First application of the Integrated Karst Aquifer Vulnerability (IKAV) method – potential and actual vulnerability in Yucatán, Mexico
Groundwater vulnerability maps are important decision support tools for water resource protection against pollution and helpful in minimizing environmental damage. However, these tools carry a high subjectivity along the multiple steps taken for the development of such maps. Additionally, the theoretical model on which they are based does not consider other important parameters, such as pollutant concentration or pollutant residence time in a given section of the aquifer, solely focusing on the theoretical travel time of a pollutant particle from a release point towards a target. In this work, an integrated methodology for the evaluation of potential (intrinsic) and actual vulnerability is presented. This integrated method, named Integrated Karst Aquifer Vulnerability (IKAV), was developed after the analysis of several study cases around the world and the application of multiple intrinsic groundwater vulnerability methods in a selected study area. Also, a solute transport model served as the basis to define additional parameters for vulnerability analysis for areas severely affected by anthropogenic practices. However, the focus of the transport model must not be mistaken to be hazards and risk mapping. A defined workflow and several criteria for parameters and attributes selection, rating and weighting, and vulnerability classification are presented here. The first application of the IKAV method was carried out in the Yucatán karst, demonstrating to be a reliable method for vulnerability estimation. Results demonstrated the scope of the IKAV method in highlighting important regional conditions, minimizing the subjectivity, and expanding the analysis of vulnerability.
Since the introduction of the groundwater vulnerability concept by Albinet (1970), several redefinitions and subclassifications for this concept have been continuously proposed. The current groundwater vulnerability conceptual model evaluates the geological, hydrological, and hydrogeological characteristics of a given area by their ability to transport a pollutant particle from a release point towards a specific location in the aquifer (Zwahlen, 2003). The general process to estimate this sectional aquifer vulnerability then follows a release–pathway–target model (Goldscheider, 2005). With focus on the natural characteristics along the pathway, such as soils, vegetation, lithology, and slope, among others, the travel time of a theoretical and immutable pollutant particle is estimated. This “intrinsic” (natural) vulnerability does not consider the possible changes that a pollutant particle can experience along the pathway; when such changes are also evaluated, the vulnerability analysis is then referred to as “specific” (Vrba and Zaporožec, 1994). However, the extensive cost to obtain data and carry out fieldwork to evaluate the possible changes of a given pollutant along the pathway make specific vulnerability analysis less applicable in comparison with intrinsic-vulnerability maps.
Given the high heterogeneity and anisotropy of karst, groundwater vulnerability methods for porous aquifers are not efficient in estimating karst aquifer vulnerability. According to the literature, EPIK is the first proposed method to evaluate the vulnerability of karst aquifers with the aim to define protection areas and to fulfil Swiss water regulations (Dörfliger et al., 1999). The EPIK, the DRASTIC (Aller et al., 1987), and the European framework (Zwahlen, 2003) are the basis for multiple vulnerability methods for karst that have arisen during the last 2 decades (Iván and Mádl-Szőnyi, 2017; Parise et al., 2018). The recurrent appearance of new methods and/or adaptations to estimate groundwater vulnerability in karst exhibits the complexity of such systems. The existence of multiple methodologies following the same purpose complicates the selection of an appropriate method to be applied for a given karst area. Application of several methods over the same region usually displays a considerable mismatch among final vulnerability classifications (Gogu et al., 2003; Ravbar and Goldscheider, 2009; Marín et al., 2012; Stevanović, 2015; Kavousi et al., 2018). However, recent studies have demonstrated that, under some regional conditions, outcomes from multiple vulnerability methods can display a remarkable agreement in vulnerability classification; nevertheless, a high correlation among methods must never be taken as a confirmation of reliability (Moreno-Gómez, 2021).
In general, two main problems are directly related to current methodologies in estimating intrinsic vulnerability: the high subjectivity and the uncertainty of the current conceptual model. For the former, several factors influence the high subjectivity along the steps to estimate vulnerability such as personal interpretations, the use of dissimilar standards, inconsistent parameters and attributes, and discretional rates and weights (Fig. 1). Regarding the latter, the sole evaluation of the theoretical travel time of an immutable pollutant particle can severely mislead the vulnerability analysis and decisions based on it; this problematic has been pointed out by the COST Action 620 (European Cooperation in Science and Technology) final report (Zwahlen, 2003). Therefore, the reliability of the intrinsic-vulnerability analysis is arguable for karst aquifers that are already affected by anthropogenic activities.
Given the high heterogeneity of karst, either locally, regionally, or continentally, the standardization of parameters (including attributes, rates, and weights) utilized to estimate intrinsic groundwater vulnerability is a very complex task (Ford and Williams, 2007; Parise et al., 2015). However, a standardized process can help to minimize the subjectivity of the current conceptual model. From previous studies with the aim to develop an Integrated Karst Aquifer Vulnerability approach, important considerations have been highlighted for map selection, parameter filtering, attribute discretization, and rating schemes (Moreno-Gómez et al., 2018, 2019; Martínez-Salvador et al., 2019; Moreno-Gómez, 2021).
Nowadays, the anthropogenic stress that karst aquifers around the world are experiencing is undeniable. Pollution generated by agricultural practices, cattle raising, wastewater disposal, and dumping sites, among others, have already affected karst groundwater quality in some regions (Parise et al., 2004; Pronk et al., 2007; Ravbar and Kovacic, 2015). On the other hand, large well fields, with considerable water extraction volumes, can diverge groundwater flow from its natural course. This increases the uncertainty of current intrinsic-vulnerability methods (either for source or resource estimations) with regard to their application on already anthropogenically affected karst areas. In order to evaluate current vulnerability scenarios, solute transport models can be beneficial additional criteria to enhance the role of vulnerability maps as decision support tools (Martínez-Salvador, 2018; Martínez-Salvador et al., 2019).
With the aim of solving the aforementioned problems and enhancing the scope of groundwater vulnerability as a decision support tool, the Integrated Karst Aquifer Vulnerability (IKAV) method is presented here as the groundwork for further vulnerability studies. This integrated strategy was developed based on three research steps: (1) the detailed and critical analysis of numerous applications of current intrinsic-vulnerability methods around the world, (2) the application and comparison of eight selected methods on a karst study area, (3) the application and evaluation of a transport model in a severely polluted karst region to determine additional criteria for vulnerability estimation. The state of Yucatán, an important karst region located in Mexico, was selected to carry out steps 2 and 3. The IKAV method demonstrates the necessity to expand the vulnerability analysis beyond the two-dimensional intrinsic scheme and the importance of considering the current vulnerability scenario of the evaluated aquifer. This work presents the first application of IKAV; the study area is the Mexican state of Yucatán.
The state of Yucatán (39 500 km2) is located in the northern part of the Yucatán Peninsula, a transboundary limestone platform (approximately 160 000 km2) covering parts of Mexico, Guatemala, and Belize (Fig. 2a). According to Weidie (1985), the peninsula is formed by limestone, dolomites, and evaporites with thicknesses reaching up to 1500 m. Rocks at the surface correspond to northward sequences from the Upper Cretaceous (Paleogene period) to Holocene (Quaternary period) epochs (Butterlin, 1958; Bonet and Butterlin, 1962; López-Ramos, 1975). The Yucatán Peninsula is classified as a well-developed karst given the existence of systems of solutional conduits of considerable diameter in the subsurface, extending in the range of kilometres (Bauer-Gottwein et al., 2011). The state of Yucatán (hereinafter Yucatán) has four hydrogeological regions: Coastal Area, Inner Cenote Ring, Central Plain, and Valleys and Hills. Yucatán presents interesting characteristics such as areas of high doline (cenote) density, regional faults, and a nearly flat topography for most of the state (Fig. 2b); other important characteristics of Yucatán, such as the shallow water table, soil distribution, spatial precipitation pattern, and lithology are visually displayed in Appendix A.
The nearly flat topography and the considerable secondary porosity do not allow surface streams to generate, making diffuse infiltration dominant in Yucatán (Fig. 2c). For more detailed information regarding the characteristics of Yucatán, the works of Bonet and Butterlin (1962), Doehring and Butler (1974), Lesser (1976), Pope et al. (1991, 1993), Hildebrand et al. (1995), and Lugo-Hubp and García (1999) are highly recommended.
Yucatán is divided into 106 municipalities of variable extension and population density (Fig. 3a). An important area in Yucatán is the Mérida metropolitan area (MMA), located in the Inner Cenote Ring region and composed of six municipalities (Fig. 3b). This highly urbanized area with around 1.1 million inhabitants (approximately 52 % of Yucatán's population) presents several environmental problems related to fast urbanization and the practices derived from it (Pacheco and Cabrera, 1997; Marín et al., 2000; Pacheco et al., 2002; González-Herrera et al., 2014; Rojas-Fabro et al., 2015). It is estimated that more than 55 % of the water extraction for human consumption in Yucatán takes place in the MMA (Fig. 3c). Wastewater volumes are estimated from 75 % to 80 % of the extracted water for human consumption. With approximately 80 % of the population of the MMA utilizing artisanal septic tanks, which are permeable, a pollution plume has been generated from the continuous leaking of wastewater towards the aquifer (Graniel et al., 1999; Marín et al., 2000). The pollution plume has made the upper 20 m of the aquifer below the city of Mérida unsuitable for human consumption (Luis Marín, personal communication, July 2017).
The continuous disposal of untreated wastewater into the aquifer, the use of nitrogen-rich fertilizers used in agriculture, and large-scale pig farming has put Yucatán under a severe nitrate () pollution scenario (Pacheco and Cabrera, 1997; Pacheco et al., 2001, 2002; Drucker et al., 2003; Pacheco-Ávila et al., 2004b; Delgado et al., 2010). Given that intrinsic vulnerability solely evaluates the probable advection of a theoretical pollutant, no other characteristics reflecting the current vulnerable state of the aquifer (or water extraction wells) are analysed. In the Yucatán karst, as well as many other already affected karst areas, the inclusion of additional parameters for vulnerability estimation is critical.
The proposed IKAV method aims to minimize the subjectivity of the current process in order to estimate intrinsic resource vulnerability and to provide additional parameters, such as pollutant concentration from solute transport, in order to enhance the vulnerability analysis. Given that intrinsic vulnerability evaluates an “IF” condition (vulnerability under an incidental pollution scenario), while solute transport aims to simulate a current vulnerable scenario, IKAV is presented as a complementary analysis of these two conditions: potential-vulnerability (IKAV-P) and actual-vulnerability (IKAV-A) maps are proposed. IKAV is based on three interconnected principles to be fulfilled: infiltration distinction, regionalization, and representative activity/pollutant (Fig. 4).
For the infiltration distinction principle, the objective is to define the goal of the analysis, settling either a point infiltration or diffuse infiltration condition; this selection dictates how some parameters will be evaluated. For example, if a point infiltration scenario is going to be investigated, high slopes and fine-textured soils will represent a more vulnerable condition due to their runoff generation capacity; however, these intrinsic characteristics will have an opposite role in vulnerability for diffuse infiltration conditions.
Regionalization settles the individual ranges for the selected parameters; here it is critical how the rates will be assigned to attributes in congruence with the criteria from the infiltration distinction procedure. This step aims to evaluate vulnerability according to existing conditions in a given study area, avoiding the indirect evaluation of external characteristics proposed by current methodologies but not present in the area of interest. The link of this principle to the infiltration distinction principle relies on the rating scheme of some parameters and their further influence to be represented by weights.
The representative-pollution principle defines the contaminant to be evaluated in correspondence with anthropogenic activities carried out in a given study area. The representative pollutant or the activity from which it is derived is directly linked to the infiltration condition. For example, the use of fertilizer in hillside farming could represent a more vulnerable condition for a swallow hole catchment than sewage leakage occurring in the same area.
In addition to these principles, a group of multiple criteria is proposed in order to minimize the high subjectivity and the discretional definition of vulnerability classes for IKAV-P. This additional criterion is also beneficial for the establishment of vulnerability classes for IKAV-A, according to the permissible maximum of the studied pollutant or the purpose of the extracted water at a given point of the aquifer (e.g. human consumption or irrigation).
IKAV-P is proposed as a multi-criteria decision analysis (MCDA) following a rating–weighting system. As proposed by several authors, this type of analysis must evaluate attributes following certain characteristics for their classification (Carver, 1991; Heywood et al., 1995; Malczewski, 1999, 2006). Therefore, with a focus on karst vulnerability, the attributes must be the following.
Measurable. Given the fact that it is challenging to measure some attributes in karst areas, such as the degree of karstification, statistical tools can help to better define notable differences in the study area.
Non-redundant. According to current intrinsic-vulnerability methods, some parameters can be defined as utilizing the same basemap (e.g. karstification and epikarst development from karst surface features according to EPIK). The inclusion of multiple parameters derived from the same evaluation process must be avoided to prevent over-/underestimations of vulnerability.
Minimal. Considering that some degree of subjectivity is inherently attached to each evaluated characteristic, the use of an extensive number of parameters can complicate the process. Therefore, the number of evaluated parameters must be kept to a minimum, only considering the relevant ones to characterize and evaluate specific karst features.
Available. The objective must be achieved with the available data from the region. Under conditions of low-resolution data, such data must be excluded, even if this has a significant influence on the objective.
From the analysis of multiple study cases around the world and the application of eight intrinsic groundwater vulnerability methodologies in the Yucatán karst (see Moreno-Gómez et al., 2018, 2019), it was possible to analyse and redefine additional criteria to be added to the previous aspects. In order to improve the intrinsic-vulnerability analysis with IKAV-P, the following attributes' criteria must also be considered.
Variable. Given that the goal of potential intrinsic vulnerability (or any other MCDA) is to provide a range of conditions to evaluate their relevance for a given objective, a map layer displaying homogeneity will have no purpose for the research. Although a homogeneous map layer could be useful when comparing vulnerability maps of different areas, homogeneous layers must be excluded from the analysis to achieve a regionalization principle.
Unambiguous. The objective must be clear and well defined. Either for point or for diffuse infiltration, the characteristics influencing such processes must be considered separately, i.e. assessing the performance of an individual analysis for each scenario and avoiding the use of “rest-of-the-area” conditions.
Distinctive. The influence of some parameters is different for point and diffuse infiltration; therefore, rates and weights assigned to individual parameters must be evaluated according to the objective. This is clear for slope and soil texture, the influence of which must be evaluated as oppositional for each scenario.
Territorial. The number of attributes, their evaluation, and their rating must be solely dependent upon the local or regional characteristics but follow a rating pattern linked to the objective. Attributes proposed by other methods but not present in the study area must be avoided as vulnerability indicators.
IKAV-A is proposed based on solute transport modelling to represent an approximation of a current pollution scenario. By taking the most influential anthropogenic activities, a representative pollutant can be utilized for a better interpretation of resource (or source) vulnerability by pollution levels (permissible maximum). Activities affecting groundwater quality in karst areas are variable, as well as the pollutants derived from such activities. Therefore, several IKAV-A maps can be produced according to the regional anthropogenic activities. The definition of actual-vulnerability rates will depend on the type of pollutant, regional regulations, and the purpose of a given extraction well. Having defined the three principles of IKAV and a general criterion for map layer selection (including rating and discretization), a workflow is presented as the groundwork to evaluate an integrated groundwater vulnerability (Fig. 5).
The process is considered an improved guideline to obtain a potential (intrinsic) groundwater vulnerability map with the addition of an actual-vulnerability map from solute transport based on the regional anthropogenic activities. For IKAV-P, the data utilized to generate the necessary map layers were obtained from public sources (see Appendix B). For IKAV-A, a solute transport model of the MMA, presented by Martínez-Salvador et al. (2019), was utilized with slight modifications.
Following the previously presented criteria in order to fulfil the IKAV principles, IKAV-P was tested in Yucatán. Given the Yucatán hydrogeological characteristics, the infiltration scenario was regionally defined as diffuse due to the low relief of the area and the high fissuring conditions not allowing for surface streams to generate. Therefore, the vulnerability analysis and the rating of attributes strictly followed a vertical-advection point of view. A filtering process for the available data, according to the attribute's selection criteria, helped to exclude parameters not contributing to the analysis. Yucatán's lithology is considered to be homogeneous (limestone) at a regional scale; since no significant differences can be evaluated for our purposes, the lithology map was excluded from the process. Having defined a regionally diffuse infiltration condition, high slopes represent the most protective attribute from the topographic parameter. Nevertheless, given the low variability displayed by the slope map (see Appendix A) the slope parameter was also excluded. Vegetation is commonly associated with infiltration conditions and its influence on runoff generation (see EPIK, COP, PI, and the Slovene approach). Having defined a diffuse infiltration condition for our study area, vegetation maps were excluded from the analysis. Similarly, due to the lack of data regarding the spatial distribution of recharge, this parameter was also excluded from IKAV-P but included as a boundary condition in IKAV-A.
Given the lack of indirect data in Yucatán, such as spring flow, doline density maps can be used as representative of karstification, epikarst, or aquifer development. However, given the fact that the use of three map layers derived from doline density could lead to over-/underestimations, doline density was selected to generally represent a karstification map (Fig. 6a); fissuring density maps were excluded from the analysis due to their similar discretization with those from doline density (see Appendix A). Given the contrasting depths at which groundwater can be found in Yucatán, the thickness of the unsaturated zone was selected as a vulnerability parameter. The attribute discretization of this parameter was carried out to highlight the shallowness of the water table in the flat plain and its deeper location in the Valleys and Hills region (Fig. 6b). Soils were evaluated based on clay content percentage for the relationship between hydraulic conductivity (K) and fine particles content in soils (Fig. 6c). This evaluation of soils as a vulnerability parameter is more sensible and serves to avoid the subjectivity derived from the multiple standards commonly applied to classify vulnerability from soil texture. For this work, it was decided not to include soil thickness as part of the analysis given the extension of the study area and the low spatial resolution of data from boreholes. Arbitrarily, an inverse relationship was defined for rating purposes: the lower the rating, the higher the vulnerability.
In total, three parameters, representing a notable variability, were selected from the data-filtering step to represent karstification, vadose zone, and soils. In order to avoid misrepresentations from the usually applied discretional approach and to fulfil the regionalization rule, the attribute classification was performed statistically. The Jenks classification method (natural breaks) was utilized for this purpose; this method minimizes the variance within a given attribute and maximizes the variance between them (Jenks, 1967). Parameters' attributes were rated solely in reference to vertical advection (infiltration distinction rule), allowing pollutants to migrate from the surface to the water table. The ranges assigned by the Jenks classification are intended to fulfil the regionalization rule; hence, these ranges are only applicable for our study area (Table 1).
Regionally, the most vulnerable attributes are areas of high doline density, a shallow water table, and soils with low clay content. Previous tests of IKAV-P indicate that a correlative rating system is beneficial for a statistical classification of the vulnerability index and the determination of vulnerability classes. After the classification of attributes and the assignment of rates, a procedure to define the importance of each map layer was carried out. The analytical hierarchy process presented by Saaty (2008), one of the several methods utilized for MCDA, was selected to determine the ranking or importance of each parameter for IKAV-P. After the standardization of the pairwise comparison matrix, the eigenvalues (priority) were obtained (Table 2).
With the consistency index (CI) being lower than the random index (RI), the weights were accepted. The final vulnerability index for the Yucatán case was then defined by Eq. (1):
where K, VZ, and S correspond to the karstification, vadose zone, and soil map layers, respectively; numbers inside brackets are the weights calculated from the AHP.
For the development of IKAV-A, a previously presented solute transport model was slightly modified. The conceptual model's setup defines the main urban area from each municipality as an infiltration basin for pollution (Martínez-Salvador et al., 2019). Four large well fields (JAPAY I to IV; Junta de Agua Potable y Alcantarillado de Yucatán) supplying water for the city of Mérida were defined as stressors of the aquifer. Piezometric data were obtained from six monitoring wells from the metropolitan monitoring network to serve for the calibration process. Precipitation values for monthly and yearly averages were obtained from historic datasets publicly available from the CLICOM (CLImate COMputing) CICESE (Centro de Investigación Científica y de Educación Superior de Ensenada) website (SMN, 2017). As boundary conditions, recharge was estimated from precipitation by the application of the APLIS (altitud, pendiente, litología, infiltración y suelo; altitude, slope, lithology, infiltration, and soils) methodology (see Martínez-Salvador, 2018); recharge was subdivided as coastal, metropolitan, and rest of the area. Neumann conditions were settled for eastern and western boundaries, and Dirichlet conditions were settled for the southern and northern (sea level) boundaries (Fig. 7).
The subsurface was defined with four layers, one of them representing a preferential flow layer; this was defined according to drilling works for monitoring wells indicating cavities (Fig. 8). The model was run in MODFLOW-2005 (Harbaugh, 2005) utilizing the graphical user interface (GUI) ModelMuse version 188.8.131.52 (Winston, 2009).
In order to simulate the turbulent conditions of the preferential layer, the Conduit Flow Process (CFP) package for MODFLOW-2005 was utilized in mode 2 (Shoemaker et al., 2008). After calibration, horizontal hydraulic conductivities (Kx and Ky) were defined as 1, 0.5, 1.115, and 1 m s−1 for the up to bottom layers; vertical hydraulic conductivity (Kz) was defined homogeneously as 1.115 m s−1. The Modular 3-Dimensional Transport model with multi-Species structure (MT3DMS) was utilized to simulate advection and dispersion (Zheng and Wang, 1999); given the anthropogenic impact in the MMA, was selected as the species. The initial concentration was settled as 80 mg L−1 in order to approximate previously reported concentrations in the Mérida subsurface (Pacheco-Ávila et al., 2004a; Rojas-Fabro et al., 2015). The contaminant load was settled as continuous for the major urban settlements of the MMA in an attempt to study the effect of the constant leaking of wastewater from septic tanks that actually occurs in the MMA.
Additionally, 130 extraction wells were included solely as representative given the lack of data regarding their extraction volumes and operational times. For more detailed information regarding the model setup, considerations, and calibration, the work of Martínez-Salvador et al. (2019) is recommended.
In order to define degrees of vulnerability, the permissible maximum of for drinking water was used as a base; according to Mexican standards, this value is settled as 45 mg L−1 (Diario Oficial de la Federación, 2000). Taking this as a base to determine highly vulnerable water sources, the subsequent vulnerability classes were defined according to previous water quality monitoring campaigns performed in Yucatán (Pacheco and Cabrera, 1997; Pacheco, 2003; Pacheco-Ávila et al., 2004b; Pérez-Ceballos, 2004; Rojas-Fabro et al., 2015). Given that concentrations can occur naturally, concentrations below 9 mg L−1 were defined as a low vulnerable condition (Table 3). Although the natural occurrence of is dependent on characteristics such as geology and soils, the proposed value for low vulnerability approximates the estimates for non-inhabited areas in Yucatán as presented by Pacheco and Cabrera (1997).
∗ Indicates no influence of the pollution source.
The IKAV-P map indicates that areas of high doline density, suggesting advanced karstification, are the most vulnerable under an incidental pollution scenario in agreement with European methods previously tested in Yucatán (see Moreno-Gómez et al., 2018, 2019). However, IKAV-P also classifies the near coastal areas for very high vulnerability (VHV) where very shallow water tables, fissuring, and sands are present. Despite karstification representing the most influential characteristic, given its assigned weight, a distinctive arrangement of some characteristics (coarse soils and a shallow water table) also indicates a VHV regional condition (Fig. 9). A vulnerability reduction pattern, from high vulnerability (HV) to moderate vulnerability (MV), is displayed in relation to an increment of the vadose zone in the Inner Cenote Ring and the Central Plain areas; in the Valleys and Hills region, the considerable depth of the unsaturated zone promoted a class of low vulnerability (LV). This manifests a consistent interpretation of this coastal aquifer in relation to the contrasting conditions with respect to the flat-plain and the southern-hill area; both are important regional conditions that were not highlighted by the previously applied intrinsic groundwater vulnerability methods in Yucatán (see Appendix C).
Despite the weight assigned to the soil map being the lowest, it exhibits an important role under different conditions. Soils are depicted as protective when their clay content is above 30 %. In the southern area of the Inner Cenote Ring, a clay-rich soil promotes said area for moderate vulnerability, despite the shallow water table and fissuring. Similarly, in the south-eastern Central Plain, clay-rich soils provide some protection in areas were karstification is at its maximum. Soils are the main promotor for very low vulnerability (VLV) in the Valleys and Hills area where karstification is low and the water table is found at its deepest. In general, results from IKAV-P display percentual vulnerability classes as 30.5, 32.1, 27.2, 7, and 3.3 for VHV, HV, MV, LV, and VLV, respectively.
Results from the solute transport model indicate that pollution moves northward, following the regional groundwater flow. Due to the fact that the city of Mérida is the most populated urban area, representing an elevated number of septic tanks leaking wastewater into the aquifer, the effect of the plume in the northern area of the city is more pronounced in comparison with other urban areas. Model layers 1 and 2 (upper and lower epikarst) present similarities regarding the spatial distribution of concentration; however, model layer 2 is more representative of the IKAV-A source vulnerability given the large number of extraction wells located at this depth (Fig. 10a). The constant leaking of wastewater from septic tanks has already affected the underlying aquifer in the MMA; however, results indicate that concentrations decrease with depth, not representing an immediate threat for extraction wells fields located at these depths (Fig. 10b and c).
Comparing results from IKAV-P for the MMA, it is can be clearly seen that no further discretization of vulnerability is possible from the intrinsic characteristics of Yucatán and the data spatial resolution (Fig. 10d). Nevertheless, IKAV-A provides important insights and an additional parameter (pollutant concentration) to expand the analysis in order to improve decisions regarding protection and corrective strategies (Table 4).
Although approximately 85 % of the extraction wells are classified as low vulnerable, 11 % could be severely affected by concentrations higher than the permissible maximum for drinking water. According to the estimated number of inhabitants who receive water supply from these highly vulnerable wells, approximately 10 % of the MMA population (∼100 000 inhabitants) can be considered under risk given the simulated concentrations at these extraction points.
Due to the lack of tracer test data, IKAV-P was compared by spatial correlation with the regional intrinsic-vulnerability method IVAKY (índice de la vulnerabilidad del acuífero kárstico yucateco; Aguilar-Duarte et al., 2016). Given that IVAKY categorizes vulnerability into six classes, including the additional class “extreme vulnerability”, this class was merged with the VHV category for comparative purposes (Fig. 11a). Despite the considerable differences in the number of used parameters, attributes, rating schemes, and assigned weights, IKAV-P and IVAKY show a very similar percentual tendency to classify vulnerability in the Yucatán karst (Fig. 11b). The percentual similarity is not a definitive indication of spatial relationship; therefore, to investigate the spatial correlation in vulnerability classes between IVAKY and IKAV-P, an overlapping process was carried out utilizing ArcGIS tools. A remarkable total correlation, above 50 %, is displayed by IKAV-P and IVAKY; this correspondence in the assignation of vulnerabilities is outstanding given the fact that the best correlated European methods displayed less than 30 % of agreement with the IVAKY method (see maps in Appendix C). With these results, the plausibility of the IKAV-P method to display potential vulnerability in accordance with the principles of regionalization and infiltration scenario was demonstrated.
Unfortunately, no water quality data were available to validate outcomes from IKAV-A; however, results are consistent with previous studies and water sampling campaigns that took place in Mérida and peripheral shallow wells. Outcomes from IKAV-A are, to some degree, consistent with the previous studies (see among others Pacheco et al., 2001; Pacheco-Ávila et al., 2004a; Pacheco, 2004; and Rojas-Fabro et al., 2015). In general, the Mérida subsurface experiences a continuous pollution condition, in which the upper aquifer layers seem to have a permanently high concentration. The pollution plume generated in Mérida moves northward according to the natural groundwater flow, temporarily increasing levels in northern areas during the high-rainfall season. In the case of other cities, the pollution seems to be locally generated, with highly transient concentration levels; the temporal variability, reported by preceding water sampling campaigns, is also demonstrated by the model. The increment of levels can be associated with a flush effect, incrementing the pollution in the north of Mérida; therefore, the seasonal variations of the representative pollutant should also be considered in further protection strategies.
The empirical interpretation of IKAV-P showed a better approximation of the natural conditions in the Yucatán karst when compared with results from well-established European methodologies previously applied in the same area. By following the steps proposed by the IKAV method, it was possible to obtain a representative potential-vulnerability map, highlighting the regional characteristics and their influence on the migration of accidental pollution from the surface to the water table. Regionally, a combination of sands, shallow water table, and fracturing also represent a VHV condition, despite the infiltration scenario in Yucatán settled as diffuse. On the other hand, the most protective natural characteristics in the region are those found in the Valleys and Hills area, where soils with a high clay content, the deepest groundwater table, and a low fracturing, are the dominant natural characteristics. It is important to note that these contrasting conditions were not highlighted by any of the previous European methods applied in Yucatán.
IKAV provides a general guideline to estimate groundwater vulnerability from a regional point of view but allows for necessary flexibility to fulfil the principles and the criteria presented previously. Flexibility is always necessary, since interpretations can vary according to data, infiltration scenario, pollutant type, and objective. The combined evaluation of potential and actual vulnerability is demonstrated to be necessary for vulnerability as a decision support tool for karst areas already affected by anthropogenic practices. Therefore, IKAV can provide an enhanced analysis, beneficial for decision support and the development of strategies for protective, preventative, and corrective measures.
Given the necessity to minimize the subjectivity of current intrinsic karst groundwater vulnerability methods and to expand the analysis of vulnerability including the anthropogenic influence, the IKAV method is proposed. Following a defined criterion for map layer selection and filtering, a correlative rating system, and a set of principles to be fulfilled, a potential (intrinsic) vulnerability map can be obtained; this procedure aims to minimize the subjectivity and to present a vulnerability map in accordance with the characteristics of the region being evaluated. By quantifying solute transport emerging from anthropogenic activities, the actual vulnerability of the aquifer (or water sources) can be estimated, expanding the decision scope of vulnerability maps as a decision support tool. Some important contributions of IKAV for vulnerability studies are understood to be as follows.
The over-/underestimation of the vulnerability outcomes is reduced by performing individual analyses for point and diffuse infiltration conditions (infiltration distinction principle).
Independent of the number of parameters or attributes, a correlative rating system is favourable for further vulnerability classification.
The vulnerability index partition is not discretional but dependent on statistical distribution, allowing for a better representation of the study area's characteristics (regionalization principle).
Vulnerability results are more consistent with regional characteristics when the attributes display a pronounced variability (regionalization principle).
The evaluation of the anthropogenic influence via solute transport enhances the vulnerability analysis, depicting conditions not displayed by intrinsic methods (representative-pollutant principle).
Actual-vulnerability maps can represent individual aquifer layers, providing additional criteria for cost–benefit judgements.
The IKAV method is not only a vulnerability indicator, but it is also capable of revealing possible solutions for endangered water sources.
In conclusion, the IKAV method is an improved scheme to estimate potential groundwater vulnerability, integrating solute transport to evaluate an aquifer's current state of vulnerability. The IKAV method reduces the inevitable subjectivity of other vulnerability methods by proposing a workflow with well-defined principles, rules, and systemic attribute evaluation criteria. The IKAV-P map provides decisive insights for protective–preventive procedures, and the IKAV-A map focuses on presently vulnerable sections of the aquifer in order to implement corrective actions and maintain optimal groundwater quality. IKAV was developed considering how parameters, attributes, and values can be assigned to highlight regional intrinsic differences according to infiltration scenarios (regionalization and infiltration distinction rules). IKAV could be applied for different karst areas inasmuch as the proposed steps and considerations to develop vulnerability maps are followed. This integrated methodology can be taken as the groundwork to expand further vulnerability studies and their role as decision support tools.
a Data publicly available at https://www.inegi.org.mx/app/mapa/espacioydatos/default.aspx (last access: 6 September 2018) (vector maps at 1 : 50 000 scale). b Digital elevation model (ASTER GDEM, 30 m resolution) from https://earthexplorer.usgs.gov/ (last access: 6 September 2018). c Data publicly available at http://clicom-mex.cicese.mx/ (last access: 6 September 2018).
Datasets and maps are available upon request to the authors.
MM-G conceptualized the project, designed the methodology, validated and curated the data, conducted the investigation, prepared the original draft of the paper, and acquired funding. MM-G and CM-S conducted the formal analysis and visualized the data. MM-G, CM-S, RL, CS, and JP reviewed and edited the paper. RL, CS, and JP supervised the project.
The contact author has declared that neither they nor their co-authors have any competing interests. The funders had no role in the design of the study; in the collection, analysis, or interpretation of data; in the writing of the manuscript; or in the decision to publish the results.
Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
This research has been supported by the Consejo Nacional de Ciencia y Tecnología (grant no. CVU: 466945).
This open-access publication was funded by the Technische Universität Dresden (TUD).
This paper was edited by Mario Parise and reviewed by two anonymous referees.
Aguilar-Duarte, Y., Bautista, F., Mendoza, M. E., Frausto, O., Ihl, T., and Delgado, C.: IVAKY: índice de la vulnerabilidad del acuífero kárstico yucateco a la contaminación [IVAKY: vulnerability index to pollution of the Yucatecan karst aquifer], Rev. Mex. Ing. Quim., 15, 913–933, 2016.
Albinet, M.: Les cartes de vulnérabilité des nappes d'eau souterraine á la pollution [Maps of groundwater vulnerability to pollution], Bulletin du Bureau de Recherches Géologiques et Minières (BRGM), Orléans, France, http://infoterre.brgm.fr/rapports/70-SGN-325-HYD.pdf (last access: 19 February 2019), 1970.
Aller, L., Bennett, T., Lehr, J. H., Petty, R. J., and Hackett, G.: DRASTIC: A standardized system for evaluating ground water pollution potential using hydrogeologic settings, Report EPA/600/2-87/035, United States Environmental Protection Agency, Dublin, USA, https://nepis.epa.gov/Exe/ZyPDF.cgi/20007KU4.PDF?Dockey=20007KU4.PDF (last access: 12 February 2018), 1987.
Bauer-Gottwein, P., Gondwe, B. R., Charvet, G., Marín, L. E., Rebolledo-Vieyra, M., and Merediz-Alonso, G.: Review: the Yucatan Peninsula karst aquifer, Mexico, Hydrogeol. J., 19, 507–524, https://doi.org/10.1007/s10040-010-0699-5, 2011.
Bonet, F. and Butterlin, J.: Stratigraphy of the northern part of the Yucatan Peninsula, in: Guide book. Field trip to Peninsula of Yucatan, New Orleans Geological Society, New Orleans, USA, 52–57, ISBN MSU:31293013771849, 1962.
Butterlin, J.: Reconocimiento geologico preliminar del territorio de Quintana Roo [Preliminary geological reconnaissance of the Quintana Roo territory], Asociación Mexicana de Geólogos Petroleros (AMGP), DF, Mexico, https://archives.datapages.com/data/amgp/pdf-content/1958/1958_Sep_Oct_04_X.htm (last access: 9 March 2018), 1958.
Carver, S. J.: Integrating multi-criteria evaluation with geographical information systems, Int. J. Geogr. Inf. Syst., 5, 321–339, https://doi.org/10.1080/02693799108927858, 1991.
Delgado, C., Pacheco, J., Cabrera, A., Batllori, E., Orellana, R., and Bautista, F.: Quality of groundwater for irrigation in tropical karst environment: The case of Yucatan, Mexico, Agr. Water Manage., 97, 1423–1433, https://doi.org/10.1016/j.agwat.2010.04.006, 2010.
Diario Oficial de la Federación: Modificacion de la Norma NOM-127-SSA1-1994: Límites permisibles de calidad y tratamientos a que debe de someterse el agua para su potabilización [Modification of the NOM-127-SSA1-1994 Standard: Permissible quality limits and treatments to which the water must be subjected for its purification], NOM, 127-SSA1, 73, https://www.dof.gob.mx/nota_detalle.php?codigo=2063863&fecha=22/11/2000 (last access: 9 July 2018), 2000.
Doehring, D. O. and Butler, J. H.: Hydrogeologic constraints on Yucatán's development, Science, 186, 591–595, https://doi.org/10.1126/science.186.4164.591, 1974.
Dörfliger, N., Jeannin, P.-Y., and Zwahlen, F.: Water vulnerability assessment in karst environments: a new method of defining protection areas using a multi-attribute approach and GIS tools (EPIK method), Environ. Geol., 39, 165–176, https://doi.org/10.1007/s002540050446, 1999.
Drucker, A. G., Escalante-Semerena, R., Gómez-González, V., and Magaña-Rueda, S.: La industria porcina en Yucatán: un análisis de la generación de aguas residuales [Pig farming industry in Yucatan: an analysis of wastewater generation], Probl. Desarrollo, 34, 105–124, https://doi.org/10.22201/iiec.20078951e.2003.135.7505, 2003.
Ford, D. C. and Williams, P. D.: Karst hydrogeology and geomorphology, Wiley, Chichester, UK, 562 pp., ISBN 978-0-470-84996-5, 2007.
Gogu, R. C., Hallet, V., and Dassargues, A.: Comparison of aquifer vulnerability assessment techniques. Application to the Néblon river basin (Belgium), Environ. Geol., 44, 881–892, https://doi.org/10.1007/s00254-003-0842-x, 2003.
Goldscheider, N.: Karst groundwater vulnerability mapping: application of a new method in the Swabian Alb, Germany, Hydrogeol. J., 13, 555–564, https://doi.org/10.1007/s10040-003-0291-3, 2005.
González-Herrera, R., Martínez-Santibañez, E., Pacheco-Avila, J., and Cabrera-Sansores, A.: Leaching and dilution of fertilizers in the Yucatan karstic aquifer, Environ. Earth Sci., 72, 2879–2886, https://doi.org/10.1007/s12665-014-3192-y, 2014.
Graniel, C. E., Morris, L. B., and Carrillo-Rivera, J. J.: Effects of urbanization on groundwater resources of Merida, Yucatan, Mexico, Environ. Geol., 37, 303–312, https://doi.org/10.1007/s002540050388, 1999.
Harbaugh, A. W.: MODFLOW-2005, the US Geological Survey modular ground-water model: the ground-water flow process, US Geological Survey, Reston, USA, https://doi.org/10.3133/tm6A16, 2005.
Heywood, I., Oliver, J., and Tomlinson, S.: Building an exploratory multi-criteria modelling environment for spatial decision support, in: Innovations in GIS. Selected papers from the Second National Conference on GIS Research UK, vol. 2, edited by: Fischer, P., Taylor & Francis, London, UK, 127–136, ISBN 0748402683, 1995.
Hildebrand, A. R., Pilkington, M., Connors, M., Ortíz-Alemán, C., and Chávez, R. E.: Size and structure of the Chicxulub crater revealed by horizontal gravity gradients and cenotes, Nature, 376, 415, https://doi.org/10.1038/376415a0, 1995.
INEGI: Conjunto de datos vectoriales geológicos serie I [Geology vector data set series I], Instituto Nacional de Estadistica y Geografia (web portal), Mexico, https://www.datos.gob.mx/busca/dataset/conjunto-de-datos-geologicos-vectoriales-escala-1-250-000 (last access: 14 August 2018), 1984.
Iván, V. and Mádl-Szőnyi, J.: State of the art of karst vulnerability assessment: overview, evaluation and outlook, Environ. Earth Sci., 76, 112, https://doi.org/10.1007/s12665-017-6422-2, 2017.
Jenks, G. F.: The Data Model Concept in Statistical Mapping, International Yearbook of Cartography, 7, 186–190, 1967.
Kavousi, A., Javadi, S., Kardan Moghadam, H., and Mirarabi, A.: Evaluation of Water Resources Exploitation in a Karst Region Using Intrinsic Vulnerability Assessment, Water Harvesting Research, 3, 53–68, https://doi.org/10.22077/JWHR.2019.1055, 2018.
Lesser, J. M.: Resumen del estudio hidrogeológico e hidrogeoquímico de la Peninsula de Yucatán [Summary of the hydrogeological and hydrochemical study of the Yucatan Peninsula], Secretaria de Recursos Hidraulicos (SARH), DF, Mexico, http://www.lesser.com.mx/files/76-4-Resumen-Estudio-Geohidrologico-Yucatan.pdf (last access: 28 July 2018), 1976.
López-Ramos, E.: Geological summary of the Yucatan Peninsula, in: The Gulf of Mexico and the Caribbean, vol. III, edited by: Nair, A. and Stehli, F. G., Springer, Boston, Massachusetts, USA, 257–282, https://doi.org/0.1007/978-1-4684-8535-6_7, 1975.
Lugo-Hubp, J. and García, M.: El relieve de la península de Yucatán [The Yucatan Peninsula relief], in: Atlas de procesos territoriales de Yucatán, Facultad de Arquitectura, Universidad Autonoma de Yucatan (UADY), Merida, Mexico, 159–162, ISBN 9687556927, 1999.
Malczewski, J.: GIS and multicriteria decision analysis, Wiley, New York, USA, 392 pp., ISBN 978-0-471-32944-2, 1999.
Malczewski, J.: GIS-based multicriteria decision analysis: a survey of the literature, Int. J. Geogr. Inf. Sci., 20, 703–726, https://doi.org/10.1080/13658810600661508, 2006.
Marín, A., Dörfliger, N., and Andreo, B.: Comparative application of two methods (COP and PaPRIKa) for groundwater vulnerability mapping in Mediterranean karst aquifers (France and Spain), Environ. Earth Sci., 65, 2407–2421, https://doi.org/10.1007/s12665-011-1056-2, 2012.
Marín, L. E., Steinich, B., Pacheco, J., and Escolero, O. A.: Hydrogeology of a contaminated sole-source karst aquifer, Mérida, Yucatán, Mexico, Geofis. Int., 39, 359–365, 2000.
Martínez-Salvador, C.: Estimation of pollutant residence time and its inclusion in vulnerability assessment approaches for the Yucatan karst system in Mexico, Master thesis, Joint Master Programme in Groundwater and Global Change – Impacts and Adaptation (GroundwatCH), Dresden, Germany, 84 pp., https://fenix.tecnico.ulisboa.pt/downloadFile/844820067126346/MSc 2018 TUD Carolina Martinez Thesis_V2.pdf (last access: 6 February 2020), 2018.
Martínez-Salvador, C., Moreno-Gómez, M., and Liedl, R.: Estimating pollutant residence time and NO3 concentrations in the Yucatan karst aquifer; considerations for an integrated karst aquifer vulnerability methodology, Water, 11, 1431, https://doi.org/10.3390/w11071431, 2019.
Moreno-Gómez, M.: Development of an Integrated Methodology to Estimate Groundwater Vulnerability to Pollution in Karst Areas, 1st edn., Eigenverlag des Forums für Abfallwirtschaft und Altlasten e.V., Dresden, Germany, 181 pp., https://nbn-resolving.org/urn:nbn:de:bsz:14-qucosa2-778534, last access: 24 November 2021.
Moreno-Gómez, M., Pacheco, J., Liedl, R., and Stefan, C.: Evaluating the applicability of European karst vulnerability assessment methods to the Yucatan karst, Mexico, Environ. Earth Sci., 77, 682, https://doi.org/10.1007/s12665-018-7869-5, 2018.
Moreno-Gómez, M., Martínez-Salvador, C., Moulahoum, A.-W., Liedl, R., Stefan, C., and Pacheco, J.: First Steps into an Integrated Karst Aquifer Vulnerability Approach (IKAV). Intrinsic Groundwater Vulnerability Analysis of the Yucatan Karst, Mexico, Water, 11, 1610, https://doi.org/10.3390/w11081610, 2019.
NASA, METI, AIST, Japan Spacesystems and US/Japan ASTER Science Team: ASTER Global Digital Elevation Model V003, NASA EOSDIS Land Processes DAAC [data set], https://doi.org/10.5067/ASTER/ASTGTM.003, 2019.
Pacheco, J.: Determinación y prueba de un índice de contaminación por nitratos en el acuífero cárstico de Yucatán, México [Determination and testing of a nitrate contamination index in the Yucatan krast aquifer, Mexico], Merida, Yucatan, Mexico, 2003.
Pacheco, J.: Delimitación de una zona de reserva hidrogeológica para el abastecimiento de agua potable a la ciudad de Mérida [Delimitation of a hydrogeological reserve area for the supply of drinking water to the city of Mérida], Sistema de Investigación Justo Sierra del CONACYT, Merida, Yucatan, Mexico, 2004.
Pacheco, J. and Cabrera, A.: Groundwater contamination by nitrates in the Yucatan Peninsula, Mexico, Hydrogeol. J., 5, 47–53, https://doi.org/10.1007/s100400050113, 1997.
Pacheco, J., Marín, L., Cabrera, A., Steinich, B., and Escolero, O.: Nitrate temporal and spatial patterns in 12 water-supply wells, Yucatan, Mexico, Environ. Geol., 40, 708–715, https://doi.org/10.1007/s002540000180, 2001.
Pacheco, J., Cabrera, A., Steinich, B., Frías, J., Coronado, V., and Vázquez, J.: Efecto de la aplicación agrícola de la excreta porcina en la calidad del agua subterránea [Effect of agricultural application of porcine excreta on groundwater quality], Ingeniería, 6, 7–17, 2002.
Pacheco-Ávila, J., Calderón-Rocher, L., and Cabrera-Sansores, A.: Delineación de la zona de protección hidrogeológica para el campo de pozos de la planta Mérida I, en la ciudad de Mérida, Yucatán, México [Delineation of the hydrogeological protection zone for the well field “Mérida I”, in the city of Mérida, Yucatán, México], Ingeniería, 8, 7–16, 2004a.
Pacheco-Ávila, J., Cabrera-Sansores, A., and Pérez-Ceballos, R.: Diagnóstico de la calidad del agua subterránea en los sistemas municipales de abastecimiento en el Estado de Yucatán, México [Diagnosis of groundwater quality from municipal supply systems in the State of Yucatan, Mexico], Ingeniería, 8, 165–179, 2004b.
Parise, M., Qiriazi, P., and Sala, S.: Natural and anthropogenic hazards in karst areas of Albania, Nat. Hazards Earth Syst. Sci., 4, 569–581, https://doi.org/10.5194/nhess-4-569-2004, 2004.
Parise, M., Ravbar, N., Živanović, V., Mikszewski, A., Kresic, N., Mádl-Szőnyi, J., and Kukurić, N.: Hazards in karst and managing water resources quality, in: Karst aquifers – Characterization and engineering, Springer, 601–687, https://doi.org/10.1007/978-3-319-12850-4_17, 2015.
Parise, M., Gabrovsek, F., Kaufmann, G., and Ravbar, N.: Recent advances in karst research: from theory to fieldwork and applications, Geological Society, London, Special Publications, 466, 1–24, 2018.
Pérez-Ceballos, R.: Vulnerabilidad del agua subterránea a la contaminación de nitratos en el estado de Yucatán [Groundwater vulnerability to Nitrate pollution in the Yucatan State], Master thesis, Universidad Autónoma de Yucatán, Merida, Yucatan, Mexico, 94 pp., 2004.
Pope, K. O., Ocampo, A. C., and Duller, C. E.: Mexican site for K/T impact crater?, Nature, 351, 105–105, https://doi.org/10.1038/351105a0, 1991.
Pope, K. O., Ocampo, A. C., and Duller, C. E.: Surficial geology of the Chicxulub impact crater, Yucatan, Mexico, Earth Moon Planets, 63, 93–104, https://doi.org/10.1007/BF00575099, 1993.
Pronk, M., Goldscheider, N., and Zopfi, J.: Particle-size distribution as indicator for fecal bacteria contamination of drinking water from karst springs, Environ. Sci. Technol., 41, 8400–8405, https://doi.org/10.1021/es071976f, 2007.
Ravbar, N. and Goldscheider, N.: Comparative application of four methods of groundwater vulnerability mapping in a Slovene karst catchment, Hydrogeol. J., 17, 725–733, https://doi.org/10.1007/s10040-008-0368-0, 2009.
Ravbar, N. and Kovacic, G.: Vulnerability and protection aspects of some Dinaric karst aquifers: a synthesis, Environ. Earth Sci., 74, 129–141, https://doi.org/10.1007/s12665-014-3945-7, 2015.
Rojas-Fabro, A. Y., Pacheco-Ávila, J. G., Esteller-Alberich, M. V., Cabrera-Sansores, S. A., and Camargo-Valero, M. A.: Spatial distribution of nitrate health risk associated with groundwater use as drinking water in Merida, Mexico, Appl. Geogr., 65, 49–57, https://doi.org/10.1016/j.apgeog.2015.10.004, 2015.
Saaty, T. L.: Decision making with the analytic hierarchy process, International Journal of Services Sciences, 1, 83–98, 2008.
SARH: Sinópsis Geohidrológica del Estado de Yucatán [Geohydrologic synopsis of the State of Yucatan], Secretaría de Agricultura y Recursos Hidráulicos, Subsecretaría de Infraestructura Hidráulica, Dirección General de Administración y Control de Sistemas Hidrológicos, Mexico, D.F., 1989.
SMN: Datos climáticos diarios del CLICOM del SMN a través de su plataforma web del CICESE [Daily climate data from the CLICOM of the SMN through its CICESE web platform], http://clicom-mex.cicese.mx, last access: 6 October 2017.
Shoemaker, W. B., Kuniansky, E. L., Birk, S., Bauer, S., and Swain, E. D.: Documentation of a conduit flow process (CFP) for MODFLOW-2005, US Geological Survey, Reston, USA, https://doi.org/10.3133/tm6A24, 2008.
Stevanović, Z.: Characterization of Karst Aquifer, in: Karst aquifers-characterization and engineering, edited by: Stevanović, Z., Springer International Publishing, 692, https://doi.org/10.1007/978-3-319-12850-4_3, 2015.
Vrba, J. and Zaporožec, A.: Guidebook on Mapping Groundwater Vulnerability, Heise, Hannover, Germany, 131 pp., ISBN 978-3922705970, 1994.
Weidie, A. E.: Geology of Yucatan Platform, in: Geology and hydrogeology of the Yucatan and Quaternary geology of northeastern Yucatan Peninsula, edited by: Ward, W. C., Weidie A. E., and Back, W., New Orleans Geological Society, New Orleans, USA, 1–19, 1985.
Winston, R. B.: ModelMuse: a graphical user interface for MODFLOW-2005 and PHAST, US Geological Survey, Reston, USA, https://pubs.usgs.gov/tm/tm6A29 (last access: 12 September 2018), 2009.
Zheng, C. and Wang, P. P.: MT3DMS: a modular three-dimensional multispecies transport model for simulation of advection, dispersion, and chemical reactions of contaminants in groundwater systems; documentation and user's guide, University of Alabama, Department of Geological Sciences, Strategic Environmental Research and Development Program (USA), Alabama, USA, https://hydro.geo.ua.edu/mt3d/index.htm (last access: 8 August 2018), 1999.
Zwahlen, F. (Ed).: EUR 20912 – COST 620 – Vulnerability and risk mapping for the protection of carbonate (karst) aquifers, European Commission, Directorate-General for Research, Brussels, Belgium, https://op.europa.eu/en/publication-detail/-/publication/be3c99bf-1a0a-4213-b35d-c3faffcd355b (last access: 12 September 2018), 2003.