1. Introducción

Según la ^{Organización Mundial de la Salud (OMS)(2002)} sólo un 2.5% aproximadamente del agua total de nuestro planeta corresponde al agua dulce. Esta se encuentra en ríos, lagos y debajo de la tierra y está amenazada por el impacto negativo de las actividades humanas, el cambio climático y el crecimiento poblacional, los cuales ejercen presión con respecto al aumento del consumo y calidad de agua dulce y pone en peligro el suministro para la población.

Por su parte, la ^{Comisión Nacional del Agua (2015)}, señala que las estimaciones más aceptadas indican que el 97% del volumen de agua es salada y apenas el 3% es agua dulce o de baja salinidad, de los cuales el 75% se concentra en casquetes polares, nieves y glaciares; el 21% está almacenado en el subsuelo y el 4% corresponde a los cuerpos de agua superficial como lagos y ríos. Con respecto al agua subterránea, el 38.7% del agua utilizada en el país proviene de sus 653 acuíferos distribuidos en el territorio (^{Instituto Mexicano de Tecnología del Agua, 2019}). En este sentido, la importancia del agua subterránea es incuestionable a nivel mundial, ya que proporciona más de la mitad del agua para el abastecimiento humano, de ahí la necesidad de protegerlo de la contaminación, pues su recuperación suele ser, en muchos casos, un proceso muy largo, costoso y a veces irreversible o irrealizable (^{Valcarce et al., 2001}). Como una medida para la conservación de la calidad de las aguas subterráneas a nivel mundial, se sugiere utilizar herramientas cartográficas sobre la vulnerabilidad del acuífero a la contaminación, como una propiedad intrínseca de los sistemas de agua subterránea, independientemente del contaminante que se trate (^{Foster e Hirata, 1991}). De acuerdo con ^{Bessouat et al., (2001)}, la vulnerabilidad depende de: 1). la función de factores hidrológicos que determinan la poca accesibilidad de la zona saturada a los contaminantes, 2). la capacidad de atenuación de ésta y 3). los estratos por encima de ella. Las propiedades del medio varían de un punto a otro, lo que hace variable el potencial de un acuífero para protegerse, razón por la que algunas áreas son más vulnerables que otras.

La vulnerabilidad se define como el grado de susceptibilidad a que las aguas subterráneas se contaminen con algún elemento, en concentraciones por encima de los valores recomendados por la ^{Organización Mundial de la Salud (OMS)(2002)}, para la calidad del agua de consumo humano (^{Foster e Hirata, 1991}).

El grado de vulnerabilidad intrínseca se puede expresar mediante un índice, y para determinarlo existen diversas metodologías, las cuales tienen por lo general un enfoque hidrogeológico, y consideran principalmente parámetros como la profundidad de la zona saturada y no saturada, la litología, la recarga neta, la infiltración y la pendiente del terreno.

Entre las metodologías más comunes destacan DRASTIC, AVI, GOD, SINTACS y Pi, que no son específicos para acuíferos kársticos (^{Báez, 2001}; ^{Espinoza y Ramírez, 2002}; ^{Vrba y Zoporezec, 1994}). ^{Bolio Barrios et al., (2011)}, utiliza el método EPIK, desarrollado para acuíferos kársticos, particularmente europeos. Toma en cuenta variables como la alta permeabilidad en roca caliza, la profundidad del suelo, la infiltración y el grado de desarrollo del karst (^{Aguilar et al., 2016a}, ^2016b; ^{Delgado-Carranza et al., 2011}). Sin embargo, existe una gran subjetividad en las diversas metodologías que tienden a sesgar los resultados esperados cuando diferentes modelos se aplican en una misma área (^{Moreno-Gómez, et al., 2018}). La interpretación de los datos, el peso de los factores y las tasas de vulnerabilidad, varían entre los diferentes métodos, desplegando resultados espaciales diferentes cuando se aplican para un mismo sitio (^{Andreo et al., 2006}).

Una adecuación realizada para el Estado de Yucatán es el Índice de la Vulnerabilidad del Acuífero Kárstico Yucateco (IVAKY), considerando las formas del relieve, suelos y la precipitación total para identificar los diferentes niveles de la vulnerabilidad del agua subterránea a la contaminación, como una derivación del método EPIK (^{Aguilar et al., 2016a}). Se cuenta también con el Atlas de Peligros por Fenómenos Naturales del Estado de Yucatán, del ^{Servicio Geológico Mexicano (2013)} que ofrece información sobre densidad de dolinas, espesor de suelos, fallas, fracturas, susceptibilidad de colapsos y hundimientos, entre otros temas.

Para el caso de Yucatán, los resultados demuestran que las diversas metodologías utilizadas, sólo son aplicables en casos concretos y deben ser adaptadas para ser congruentes con las características kársticas regionales de Yucatán (^{Moreno-Gómez et al., 2018}).

Por lo tanto, el objetivo de este trabajo es generar un mapa de vulnerabilidad intrínseca a la contaminación del acuífero para el estado de Yucatán mediante la selección y el análisis minucioso de los atributos presentes en el Estado, con respecto a las formaciones geológicas, variables climáticas, usos y tipos de la vegetación, entre otros, incorporando el mapa de las Anomalías de Bouguer y considerando los métodos DRASTIC y EPIK, para identificar las zonas de mayor vulnerabilidad.

La hipótesis es demostrar la reconciliación de los métodos DRASTIC de ^{Arcega (2018)} y EPIK, en su versión IVAKY de ^{Aguilar et al., (2016)}, vinculados con la información de las anomalías gravitacionales de Bouguer, evaluadas por ^{Arellano Catalán (2017)}. En este sentido, el presente trabajo ofrece una mayor congruencia con las características kársticas regionales, así como una alternativa para evaluar y analizar las zonas con mayor vulnerabilidad en el Estado. De igual forma, el mapa resultante podrá ser utilizado para investigaciones, trabajos y proyectos posteriores, ya que es una herramienta muy útil para incorporar en las evaluaciones, diagnósticos y resoluciones en materia de impacto social y ambiental.

2. Descripción del área de estudio

La Figura 1, muestra las mediciones del campo gravitacional de la tierra en varios puntos de la superficie del Estado de Yucatán que genera el mapa de contornos de las fluctuaciones en las fuerzas de gravedad en un Datum dado. Estas fluctuaciones son causadas por las variaciones en la densidad de los diferentes tipos de roca, y permite identificar la localización y configuración de ciertos depósitos, en donde la estructura de la roca pueda presentar relación con el movimiento y almacenamiento del agua, como zonas puntuales de infiltración o áreas de infiltración difusa, u otros que afecten la ocurrencia de las aguas subterráneas, como conductos preferenciales de agua subterránea (^{Arellano Catalán, 2017}).

Figura 1 Mapa del gradiente horizontal de las anomalías gravimétricas de Bouguer (modificado de ^{Urrutia-Fucugauchi et al., 2011}). 

Con el fin de profundizar sobre las características regionales del karst de Yucatán, tomaremos el ejemplo del Anillo de Cenotes que es el borde del conocido Cráter de Chicxulub. En este punto se destaca el trabajo de ^{Urrutia-Fucugauchi et al., (2011)} sobre cráter de impacto multi-anillo de Chicxulub, en la plataforma carbonatada.

Los valores gravitacionales de los estudios geofísicos muestran un patrón de anomalías concéntricas semicirculares amplias bien definidas por arriba de los 30 mGals con máximo en un diámetro central de 40 km y dos mínimos a casi 70 km y 120 km de diámetro. Las anomalías semicirculares concéntricas fueron inferidas a 52.5 ± 5 km, 77.1 ± 6.3 km, 99.6 ±km y 139 ± 11 km de distancia radial, las cuales definen una estructura multianillos con la cumbre del anillo central limitado por el anillo 1, los anillos 2 y 3 se interpretan como los límites interior y exterior de la pared transicional del cráter, con su depresión externa al borde del cráter dado por el anillo 4.

La Península de Yucatán carece de afloramientos superficiales de litologías de impacto, sólo la presencia de la estructura enterrada que se indica por una depresión topográfica semicircular de baja amplitud (aproximadamente de 5 m) que puede observarse mediante imágenes de satélite con radar interferométrico de alta resolución (^{Hernández et al., 2011}). Esta ligera característica topográfica coincide con uno de los más conspicuos aspectos kársticos reconocidos en la Península de Yucatán, el arreglo semicircular de cenotes referidos como el anillo de cenotes, o en términos geomorfológicos como dolinas agrupadas, término utilizado por ^{Aguilar et al., (2016)}. El borde estructural del cráter enterrado se marca en superficie como resultado de la compactación diferencial de casi 5 m de las brechas y fracturamiento dentro del cráter y favorecido por los procesos de disolución a través de las fracturas y el terreno carbonatado (^{Hernández et al., 2011}).

El modelo conceptual de cómo este mecanismo de compactación de la secuencia sedimentaria del Cenozoico, al interior del cráter, dio origen al anillo de cenotes indica que el fallamiento secundario, asociado a la compactación, genera porosidad secundaria, dando origen al sistema kárstico actual. En una siguiente fase, debido a la acción disolutiva de la circulación de agua acidificada con ácido carbónico, la disolución de la roca caliza se acelera dando origen a grietas y cuevas kársticas de mayor tamaño (^{Hernández et al., 2011}), que le confieren importantes atributos hidrogeológicos como una elevada transmisividad y conductividad hidráulica. La interpretación estructural de las secciones sísmicas de la porción norte del cráter de Chicxulub muestra que el intenso fallamiento normal, representa una zona de alta permeabilidad que favorece la circulación de agua subterránea (^{Urrutia-Fucugauchi, 2011}).

Eventualmente, la disolución de la caliza se acelera por la presencia de la intrusión salina en las zonas donde la interfase agua dulce y agua salada llega. Cabe destacar que la acción de la intrusión salina juega un papel primordial, ya que, la disolución de la caliza se presenta en un plano horizontal, y no vertical como se hubiese esperado. Se ha asociado este fenómeno a la acción horizontal de la interfase salina (agua dulce - agua salada), producto de la acción de mareas astronómicas (^{Perry et al., 2002}), de las oscilaciones marinas y de bombeos intensos cerca de la costa. Durante este proceso, la Península de Yucatán experimentó periodos de ascenso y descenso del nivel del mar, y es evidente como demuestra la presencia de espeleotemas subaéreos en algunas cavernas inundadas en la actualidad. Estos períodos de ascenso del nivel del mar, probablemente, también jugaron un papel importante en la formación del anillo de cenotes (^{Hernández et al., 2011}). Los diversos programas de perforaciones, como los realizados por PEMEX, como parte de sus programas de exploración, así como los del Programa Chicxulub de la UNAM, el Proyecto Científico de Perforación “Chicxulub” y el Programa de exploración de la Comisión Federal de Electricidad, han permitido conocer a profundidad la estratigrafía de las secuencias del Mesozoico y Cenozoico dentro y fuera del Cráter de Chicxulub (Figura 2).

Figura 2 Modelo estructural del Cráter de Chicxulub, exagerando la línea vertical y mostrando los pozos S-I, Sacapuc; C-I, Chicxulub-1; Y-6, Yucatán-6; U-2, UNAM; -2; T-I, Ticul-I; U-5, UNAM-5; U-7, UNAM-7; Y-2, Yucatán-2; U-6, UNAM-6; Y-1, Yucatán-1. (modificada de ^{Urrutia-Fucugauchi et al., 2011}). 

De acuerdo con ^{Pérez-Ceballos et al., (2012)}, el Anillo de Cenotes concentra y conduce el agua subterránea hacia la costa de Celestún, al poniente, y Dzilam, por el oriente, por lo que entender y conocer la influencia del agua subterránea en la zona costera es de gran importancia. Hacia el sur de la entidad, el drenaje es mínimo y desaparece en resumideros o en cuerpos de agua superficial llamados aguadas. En la parte norte, la descarga de esta agua se realiza por medio de manantiales y en forma difusa hacia el mar, alimentando las ciénagas y lagunas costeras. Los autores realizaron un estudio de modelación inversa que permitió cuantificar los minerales y su comportamiento en temporada de secas, lluvias y nortes. En la Figura 3 se muestra la ubicación de los diversos cenotes estudiados.

Figura 3 Grupos de cenotes basado en las tendencias espaciales y temporales de las variables físicas y químicas del agua de los cenotes y modelo hidrogeoquímico de los procesos de disolución y precipitación de los minerales del acuífero para la temporada de secas (modificada de ^{Pérez-Ceballos et al., 2011}, ²⁰¹²). 

Con base en las características hidroquímicas como el pH, la conductividad eléctrica y los iones mayores, aplicando técnicas de estadística multivariada, los autores identificaron tres regiones en el Anillo de Cenotes, determinando que el comportamiento hidroquímico entre cada región fue independiente a la temporada climática. Los resultados sugieren que las Regiones 1 (zona poniente) y 2 (zona centro) están influenciadas por un mismo flujo regional y con penetración de aguas azufradas del sur de la entidad, ambas regiones se comportan de manera similar en cuanto a la química el agua y a los índices de saturación, sin embargo, la Región 3 (zona oriente) parece estar influenciada por otra fuente de agua subterránea que se incorpora al Anillo de Cenotes por la parte oriental de Dzilam de Bravo y no del sur del estado de Yucatán. Es decir, los resultados de la modelación hidrogeoquímica inversa, corroboraron la entrada de agua subterránea cercanas a los cenotes 7 y 4 (al oeste), proveniente del sur de Estado de Yucatán; y hacia el lado Este (Dzilam) se identificó otra entrada de agua subterránea cercana al cenote 18 la cual es de un origen diferente a la que ingresa en los cenotes 7 y 4 (^{Pérez-Ceballos, 2012}). Otro anillo interior de cenotes se aprecia también desde el área de Chuburná Puerto - Sisal, por el occidente, bordeando la Ciudad de Mérida y sale por el área de Telchac Puerto - San Crisanto, al oriente (Figura 3). En estos puntos, el agua subterránea aflora a manera de manantiales y fluye hacia la laguna de Chuburná Puerto y la laguna Rosada de Telchac Puerto. Al centro los escurrimientos subterráneos son principalmente por infiltraciones en manto poroso, por lo que su descarga es más estable con aperturas temporales al mar. En la localidad se les conoce con el nombre de ciénagas, como en Progreso y San Bruno (Figura 4).

Figura 4 Modelo Conceptual del agua subterránea (representación hidrogeológica básica del subsuelo de la ciudad de Mérida) (modificada de ^{Martínez et al, 2019}). 

Según los Consultores en Agua Potable, Alcantarillado, Geohidrología e Hidráulica Costera (CAPAGHC) todas las aguas subterráneas del planeta se dirigen hacia las zonas costeras y es el nivel medio del mar (nmm) el que sirve como base para la existencia de los flujos hacia las zonas de descarga, debido a esto se van originando los gradientes hidráulicos, obteniendo mediante las alturas piezométricas la red de flujo del territorio.

3. Metodología

Para la elaboración de mapa de vulnerabilidad intrínseca se consideró la metodología DRASTIC desarrollada por ^{Aller et al., (1987)} y aplicada para Yucatán por ^{Arcega (2018)}, la cual permite el análisis de la vulnerabilidad de un sistema por medio de la superposición de índices mediante el análisis de sus características geográficas. Sus siglas en inglés corresponden a D = depth (profundidad a la lámina de agua freática), R = recharge (recarga neta), A = aquifer (litología del acuífero), S = soil (tipo de suelo), T = topography (topografía), I = impact (naturaleza de la zona no saturada) y C = hydraulic conductivity (conductividad hidráulica del acuífero). El método asume que el contaminante es introducido sobre la superficie de la tierra, que es trasladado al agua subterránea por precipitación, filtración y percolación y que el contaminante es móvil en el agua.

Este modelo es el más utilizado y suele ser muy útil cuando se quiere construir mapas de vulnerabilidad con diversos datos, basándose en la asignación de índices que van de 1 (vulnerabilidad mínima) a 10 (vulnerabilidad máxima), de acuerdo a las características y al comportamiento de las variables consideradas en el acrónimo DRASTIC, además de lo expresado, se pondera su influencia dentro de la evaluación de la vulnerabilidad mediante la asignación de unos pesos de 1 a 5. Ambos índices se multiplican y luego se suman los diferentes resultados, para obtener un valor final o índice de vulnerabilidad, cuyos extremos resultantes en este caso son 354 (mínimo) y 496 (máximo), según se indica en la siguiente fórmula:

El método EPIK, es un método paramétrico desarrollado por Doerfliger y Zwahlen (1998) para evaluar acuíferos kársticos y aplicado a Yucatán en su versión IVAKY por ^{Aguilar, Y. (2016)} y se presenta en la Tabla 3. Sus siglas en inglés equivalen a E = Epikarst, P =Protective cover (cobertura), I = Infiltration conditions (condiciones de infiltración) y K = Karst network development (grado de desarrollo kárstico).

En este estudio se tomaron en cuenta las variables de los modelos DRASTIC utilizadas por ^{Arcega (2018)}, excepto las variables de dureza del agua y topografía, pues resultan muy homogéneas, así como aquellas utilizadas en el método EPIK, en cuanto a las características del relieve kárstico, como la densidad de dolinas, presencia de depresiones húmedas (cenotes).

En este caso en particular, el estudio incorpora el mapa de anomalías gravimétricas de Bouguer, en donde las diferencias de densidad de las rocas subyacentes conllevan a generar incrementos en la conductividad hidráulica de algunas áreas del relieve kárstico, siendo ésta una característica particular regional del cráter multianillos de Chicxulub, sometido a largo plazo a cierto tipo de clima, seco o subhúmedo, con fuertes o ligeras oscilaciones de temperatura. El humedecimiento favorece los procesos de disolución de la roca calcárea, las diferencias fuertes de temperatura pueden intemperizar y resquebrajar la capa de caliche superficial, favoreciendo la infiltración.

En la Tabla 1 y Figuras 7 y 8 se muestran las variables utilizadas y los mapas resultantes, así como la de los otros autores. Para elaborar el mapa de vulnerabilidad del presente estudio se utilizaron un total de 15 variables.

Figura 5 Anomalías gravimétricas de Bouguer de la península de Yucatán, que muestra el Cráter de Chicxulub, el resultado del impacto de un asteroide hace 65 años, que llevó a una extinción masiva. (modificada de Hernández (s/f). 

Figura 6 Mapa del Gradiente Horizontal Total de la anomalía de Bouguer (modificado de ^{Arellano-Catalán, 2017} y ^{Urrutia-Fucugauchi et al., 2011}). 

Figura 7 Conjunto de Mapas para la obtención del cálculo del grado de vulnerabilidad. 

Figura 8 Mapa de vulnerabilidad intrínseca del estado de Yucatán utilizando las anomalías de Bouguer, el tipo de clima y la variación de temperatura. 

Después de haber analizado y seleccionado las variables a utilizar con base a sus características y en relación con la vulnerabilidad intrínseca se procedió a la clasificación y valoración de vulnerabilidad, de acuerdo a los valores de cada variable. Para establecer los valores de vulnerabilidad del 1 al 10, como se indica en la Tabla 2, fue necesario realizar un análisis por cada capa de información y llevar a cabo la ponderación correspondiente. Posteriormente, los valores obtenidos se distribuyeron en rangos y grados de vulnerabilidad para representar el grado de vulnerabilidad en los mapas de manera individual. Finalmente, se obtuvo la sumatoria de los valores de todas las capas generadas obteniendo así, el mapa de vulnerabilidad intrínseca del Estado. Los valores que adoptan los factores de clasificación y de ponderación para cada parámetro se resumen en la Tabla 2 de Matriz de vulnerabilidad y ponderación que se presenta más adelante.

La integración de los pesos de cada factor y obtención del mapa de vulnerabilidad con la variable de la anomalía de Bouguer se realizó mediante una suma ponderada con la herramienta “Spatial analyst” de ArcGis 9.1. Este método es uno de los más fáciles de implementar, obteniéndose el valor de la vulnerabilidad como resultado de la sumatoria ponderada para cada píxel de cada factor (^{Molero et al., 2007}). Por lo tanto, en la Tabla 3 se muestra la suma de los valores distribuidos en “RANGOS” y su grado de vulnerabilidad para generar la clasificación correspondiente de cada una de las variables y su representación en un mapa (Figura 7).

A continuación, se describen las variables adicionales utilizadas para el presente estudio, la información sobre las otras variables se encuentra en ^{Arcega (2018)} y ^{Aguilar et al., (2016a)}.

4. Discusión y resultados

5. Conclusiones

El mapa muestra que más del 50% del Estado se encuentra en una vulnerabilidad entre alta y extrema, distribuyéndose en la región centro-oriental del estado, donde existe la mayor zona de infiltración. La vulnerabilidad muy alta se asocia al sistema multianillos derivado del Cráter de Chicxulub y a las diferencias de densidades de la roca y de mayor conductividad hidráulica. Por su parte, la vulnerabilidad extrema se presenta en el oriente del Estado y está vinculada con la permeabilidad del territorio, derivado de la mayor presencia de cenotes, densidad de dolinas, fracturas y fallas, y la mayor precipitación, generando procesos kársticos muy dinámicos.

En este trabajo se resalta la incorporación de las anomalías de Bouguer, en donde las diferencias de densidad de las rocas subyacentes conllevan a generar incrementos en la conductividad hidráulica de algunas áreas del relieve kárstico, lo que incrementa el nivel de vulnerabilidad en el estado. En el mapa resultante de este estudio se pueden observar las variaciones de la vulnerabilidad intrínseca a la contaminación de acuerdo con las metodologías consideradas. La vulnerabilidad de extrema a media en el territorio se desarrolla proveniente del sureste hacia el noroeste, y esto es congruente con el método DRASTIC, mientras que se presentan oleadas concéntricas de muy alta y extrema vulnerabilidad que se extienden sobre este gradiente, cubriendo muy bien los factores geomorfológicos y edáficos, humedad de suelo y climáticos principalmente, descritos en el método IVAKY.

El presente trabajo ofrece una mayor congruencia con las características kársticas regionales, así como una alternativa para evaluar y analizar las zonas con mayor vulnerabilidad en el estado en la perspectiva de elaborar políticas públicas para su conservación.

Contribuciones de los autores

En este trabajo la contribución es: (1) Conceptualización: Batllori, S. (2) Análisis o adquisición de datos: Batllori, S y Canto, M. (3) Desarrollo metodológico/técnico: Batllori, S y Canto, M. (4) Redacción del manuscrito original: Batllori, S y Canto, M. (5) Redacción del manuscrito corregido y editado: Canto, M. (6) Diseño gráfico: Canto, M.

Financiamiento

Este trabajo se realizó sin financiamiento.

Conflictos de interes

Este escrito no presenta algún conflicto de intereses.